PhreeqcUsers Discussion Forum
Conceptual Models => Design of conceptual models => Topic started by: nickvz on August 28, 2014, 07:34:23 AM

Hello
I have been struggling to understand how PHREEQC solves differential equations in the RATES block. I am trying to replicate a code that was published in the article by TirutaBarna ( Using PHREEQC for modelling and simulation of dynamic leaching tests and scenarions). The author models a stirred tank reactor by solving the following equation in the RATES block:
dC'/dt = (Q/V_leachate*(c'_in c'))
This equation describes the change in concentration of over time for a CSTR. Initially the reactor contains 0.001 moles of reactant, in this case NaCl. Pure water enters the reactor at flow rate Q, and the solution containing NaCl exist the reactor at flow rate Q. The analytical solution to the above mentioned equation is the following:
c'= c'_0 exp(tQ/V_leachate)
The PHREEQC code is as follows:
RATES
FLOW_Na
start
10 total_element = TOT("Na")
20 IF (total_element <= 0) THEN GOTO 60
30 rates = parm(1) /parm(2) * (parm(3) total_element)
40 moles = rate*TIME
50 IF (total_element+moles < 0) THEN moles = total_element
60 SAVE moles
end
SOLUTION 1
KINETICS 1
FLOW_Na
formula Na+ 1
m 1
m0 1
parms 0.001 1 0 # parm(1) = Q, parm(2) = V_leachate, Parm(3) = c"_in
tol 1e008
steps 3600 in 10 steps
step_divide 1
runge_kutta 3
USER_GRAPH 1 Compare PHREEQC and analytical solution
axis_titles "time" " c' mol/l"
axis_scale x_axis 0 4000
axis_scale y_axis 0 2e03
initial_solutions true
start
10 plot_xy total_time tot("Na"), symbol_size = 0, color = Red
end
END
The amount of Na+ that is in solution stay unchanged throughout the simulation. Am I missing something here? Please advise.
Kind regards
Nicolaus

Try the following. The changes I made were (1) changed "rates" to "rate" in line 30, (2) removed negative in line 40 (because parm(3)  total_element) was already negative, (3) removed the IF statements (optional), the IF statements can add discontinuities in the rate expression, (4) added 1 mmol NaCl to solution, (5) made the formula NaCl in kinetics, otherwise it is equivalent to removing sodium metal, which causes redox reactions, the plus sign is ignored, and (6) added a comma in line 10 of USER_GRAPH.
RATES
FLOW_Na
start
10 total_element = TOT("Na")
#20 IF (total_element <= 0) THEN GOTO 60
30 rate = parm(1) /parm(2) * (parm(3) total_element)
40 moles = rate*TIME
#50 IF (total_element+moles < 0) THEN moles = total_element
60 SAVE moles
end
SOLUTION 1
Na 1
Cl 1
KINETICS 1
FLOW_Na
formula NaCl 1
m 1
m0 1
parms 0.001 1 0 # parm(1) = Q, parm(2) = V_leachate, Parm(3) = c"_in
tol 1e008
steps 3600 in 10 steps
step_divide 1
runge_kutta 3
USER_GRAPH 1 Compare PHREEQC and analytical solution
axis_titles "time" " c' mol/l"
axis_scale x_axis 0 4000
axis_scale y_axis 0 2e03
initial_solutions true
start
10 plot_xy total_time, tot("Na"), symbol_size = 0, color = Red
end
END

Hello
Thank you, I did not notice all the simple mistakes. It must be the heat. I am in Mozambique at the moment and is hot as hell. Thank you very much for the help. I appreciate the corrections.
Kind regards
Nicolaus

Hello
I just want to make sure about the addition of NaCl to the SOLUTION block. I was under the impression if you define initial moles in the kinetic block it PHREEQC reads it as the amount of moles present in the system. Is this not the case?
Kind regards
Nicolaus

Yes, you do define the number of moles of a kinetic reaction in the KINETICS data block, particularly if it is a certain amount of a solid that could react completely. However, your rate expression requires Na in solution to have a nonzero rate, and, for your case, the amount of the kinetic reactant is irrelevant. It is the amount of Na in solution that drives the reaction.

Hello David
Can you please advise on the following:
The link below is to an example from Tony's website. It is to solve a general ODE. I am just
having some trouble to plot the variable.
http://www.hydrochemistry.eu/exmpls/RK.html
Please correct me if I am not understanding this correctly. To solve any ODE using the RATES block in PHREEQC, the following applies:
dy/dt = ((1C)/y)*D, or in the case of the example dy/dx = x * exp(x)
The constants are defined in either the KINETICS block or the RATES block. The expression is integrated by multiplying the expression by TIME.
How do you plot this variable as a function of time etc. to calculate for instance the conversion as a function of time for a plug flow reactor.
Kind regards
Nicolaus

I think you want the functions KIN and KIN_DELTA, which give the total amount of a kinetic reactant and the change in the amount of a kinetic reactant during the last rate integration.

Hello David’
I am trying to model the shrinking core model equation
d(r) / d(t) = (((1por)*D2)/(eps*rho_p*f_p*(r^2)))*((R*r)/(Rr))*(C_O2/H)
I have solved the ODE using the software Polymath. I used the following input for PHREEQC:
RATES
SCM
10 por = 0.3 # Porosity
20 D2 = 2.6e8 # Ash layer diffusion coefficient
30 eps = 1.75 # ratio of O2 consumed to moles pyrite oxidized
40 rho_p = 4.8 # g/cm3
50 f_p = 0.03 # Weight fraction of pyrite in waste rock
60 R = 0.01 # Initial core radius (cm)
70 dt = time
80 r = 0.0099
90 r = (((1por)*D2)/(eps * rho_p * (r^2)))*((R*r)/(Rr))*(TOT("O2"))
91 IF r= 0 then goto 100
100 Save r*dt
end
END
SOLUTION 1
Na 1
Cl 1
EQUILIBRIUM_PHASES 1
O2(g) 0.67
KINETICS
SCM
formula NaCl 1
m 1
m0 1
time_step 15000 in 50
step_divide 1
INCREMENTAL_REACTIONS true
USER_GRAPH 1 Compare PHREEQC and analytical solution
axis_titles "time" " core radius"
axis_scale x_axis 0 15000
axis_scale y_axis 0 0.01
headings PHREEQC analytical
initial_solutions true
start
10 plot_xy total_time , r, symbol_size = 0, color = Red
end
END
The following error message is shown in the output:
WARNING: Zero divide in BASIC line
90 r = (((1por)*D2)/(eps * rho_p * (r^2)))*((R*r)/(Rr))*(TOT("O2")).
I am guessing r is zero in the expression. I tried to correct this by adding line 80 into the RATEs block as an initial condition, but the problem is persisting so I am assuming that PHREEQC does not read this as an initial condition.
Can you please assist me with the input file and comment as to how initial conditions is defined for r(0) where r(0) = R.
Kind regards
Nicolaus

Variables are case independent, so R  r is zero. Also, TOT("O(0)") is the total moles of oxygen valence zero, which is one half MOL("O2"). TOT("O2") is zero.
You can use PRINT Basic statements to debug your RATES equations. The PRINT statements write to the output file. It is not pretty, but you can determine what is going wrong with your RATES equations.

Hello David
I have a question about the KIN function. Does this function recall the amount of reactant (in the solid phase) that is present in the system. I need to define a value for r (uncreated core radius) of pyrite before I can solve the differential equation for the rate of decrease of r, see code below:
RATES
Pyrite
start
10 por = 0.3 # Porosity
20 D2 = 3e8 # Ash layer diffusion coefficient
30 eps = 1.75 # ratio of O2 consumed to moles pyrite oxidized
40 rho_p = 30 # kg/m3
50 f_p = 0.03 # Weight fraction of pyrite in waste rock
60 r_0 = 10 # Initial core radius,m
70 PI = 3.14159
80 M_pyite = 119.96 # g/mol
90 u0 = 0.265
100 H = 26.3
110 r = (((3*KIN("Pyrite")*M_pyrite)/(4*PI*rho_p))^(1/3))
120 r = (((1por)*D2)/(eps * rho_p * (r^2)))*((r_0*r)/(r_0r))*(u0/H)
130 moles = (((4/3)*PI*(r^3)*rho_p)/M_pyrite)
140 save moles
end
END
SOLUTION 1
EQUILIBRIUM_PHASES 1
O2(g) 0.67
KINETICS
Pyrite
formula FeS2 1
m 0.083
m0 0.083
time_step 15000 in 100 steps
step_divide 1
INCREMENTAL_REACTIONS true
USER_GRAPH 1 Compare PHREEQC and analytical solution
axis_titles "time" " core radius"
axis_scale x_axis 0 15000
axis_scale y_axis 0 0.01
headings PHREEQC SCM
initial_solutions true
start
10 plot_xy total_time , (((3*KIN("Pyrite")*M_pyrite)/(4*PI*rho_p))^(1/3)), symbol_size = 0, color = Red
end
END
I have tried to define the current value for r by using the KIN function. This seems not to be working as I still have errors with division by 0.
Kind regards
Nicolaus

You misspelled M_pyrite in line 80.
Also, all Basic variables are local to a Basic script. Defining M_pyrite in a RATES definition does not carry over to a USER_GRAPH definition. You must either (1) define all of the needed variables in USER_GRAPH (M_pyrite, PI, rho_p), or (2) you can use PUT and GET to save data in RATES and use it in USER_GRAPH. I think case (1) would be better because RATES is evaluated many times per time step, whereas USER_GRAPH is evaluated only at the end of the time step, and values saved in RATES may actually not be the values at the end of the time step (there is an interpolation step after RATES has been evaluated).

Hello
It seems like I have managed to get the model to work. I am just struggling with the diffusion coefficient. When I increase it it slows down the reaction to a near zero rate. A slower rate is expected but not this slow. The initial radius seems right at 0.00143 meter. Why does it take PHREEQC 195 seconds before it reaches this value?
Please see below is the input:
PRINT
reset false
kinetics
totals
RATES
Pyrite
start
#
# Define input parameters
#
10 por = 0.4 # Porosity
20 D2 = 2.6e3 # Ash layer diffusion coefficient
30 eps = 1.75 # ratio of O2 consumed to moles pyrite oxidized
40 rho_p = 4800 # kg/m3
50 rho_bulk = 1500 # kg/m3
51 v_t = 1 # bulk_liter
52 f_a = 2/3
60 PI = 3.14159
70 M_pyrite = (119.96/1000) # kg/mol
#
# Ensure r is not equel to 0 at step 1
#
100 IF (m < 0 ) then goto 250
110 ratio = 1
120 if m0 > 0 then ratio = m/m0 # initial moles of the kinetic reactant
130 if ratio = 0 then ratio = 1
#
# Calculate volume fraction
#
140 mass_pyrite = m*M_pyrite # kg pyrite
150 v_p = (mass_pyrite/rho_p) # V_pyrite (l) = 0.009375 l at step 1
160 vol_frac = v_p/v_t # = 0.009375 l at step 1
#
# Calculate the initial radius of the unreacted core =r_0
# Calculate the initial radius of the particle =r_c
#
170 v_par = v_p/633.145 # Volume of individuel pyrite particles
180 r_0 = (((3*v_par*0.857)/(4*PI))^(1/3))
190 r_c = r_0/(10.05) # r_c is the initial core radius plus 5% ash layer
200 r = r_0*(m/m0)
210 if r <= 0 then goto 250
#
# Calculate the reaction surface in the bulk volume
#
220 react_a = ((3*(10^3))*(vol_frac/r))
#
# Calculate the rate of pyrite oxidation (mol/liter_bulk.s)
#
230 rate = (10^3)*react_a*(D2/3)*(r_c/((r_cr)*r))*(mol("O2")) # mol/l_bulk.s
240 moles = rate*time
250 save moles
end
END
#
# Define leachate volume
#
SOLUTION 1
water 0.3
#
# Define EQUILIBRIUM_PHASES 1
#
EQUILIBRIUM_PHASES 1
CO2(g) 2
O2(g) 0.67
#
# Define KINETICS 1 for Pyrite
#
KINETICS 1
Pyrite
formula FeS2 1
m 0.375
m0 0.375
time_step 19500 in 100 steps
step_divide 1
INCREMENTAL_REACTIONS true
USER_GRAPH 1 Compare PHREEQC and analytical solution
axis_titles "time" " core radius"
axis_scale x_axis 0 19500
axis_scale y_axis 0 0.00146
headings PHREEQC SCM
initial_solutions true
start
#
# Define parameters
#
1 rho_p = 4800
2 M_pyrite = (119.96/1000)
3 PI = 3.14159
10 t_mol = KIN("pyrite")
#
# Calculate the unreacted core radius = r_unreacted
#
20 vol_pyrite = (t_mol*M_pyrite)/rho_p
30 vol_par = vol_pyrite/633.145
40 r_unreacted = (((3*vol_par*0.857)/(4*PI))^(1/3))
#
# Plot unreacted core radius
#
50 plot_xy total_time , r_unreacted, symbol_size = 0, color = Red
# 80 plot_xy total_time , kin("pyrite"), symbol_size = 0, color = Red
end
END

195 s is the first time step. If you increase the number of steps or decrease the length of time, you will see a smoother transition.
I think you will have to investigate the rate equations yourself. You can use PRINT statements to look at variables to try to understand why the rate is decreasing as it is.

Hello
Thank you very much. I appreciate the help.
kind regards

Hello
I found the problem. Was an input error. Caused the rate to be 3 orders smaller than it is suppose to be. Thank you for all the help.
Kind regards
Nicolaus