Time unit difference between HP1 and RATE expressions

[mpwickham]

I am using rate expressions to manage dissolution kinetics and two potential issues.

First, I'm using print functions to generate verbose output from the rate expressions (for debugging) and have noticed that it always reports "time" as 1. "Time" is used to compute the number of moles reacted for each time step called (see the example below), so it appears to always be returning the moles for a time step size of one second (the rate convention in PHREEQC is seconds) without any adjustment for the time step size. How is the time step size "Time" passed from HP1 to PHREEQC for RATE calculations? Shouldn't this be variable?

Second, I would like to confirm that I need to directly account in RATE expressions for different time units between HP1 (if I specify hours, for example, in HYDRUS) and RATES (which are in seconds) in PHREEQC, by multiplying "time" by 86400. Or is that time conversion automatically managed by the HP1 code?

For example, Calcite dissolution where Parm(3) is 86400.

Calcite
-start
1 rem M = current number of moles of calcite
2 rem M0 = number of moles of calcite initially present
3 rem PARM(1) = A/V, cm^2/L
4 rem PARM(2) = exponent for M/M0
5 rem PARM(3) = time conversion for HP1
10 si_cc = SI("Calcite")
20 if (M <= 0 and si_cc < 0) then goto 200
30 k1 = 10^(0.198 - 444.0 / TK )
40 k2 = 10^(2.84 - 2177.0 / TK)
50 if TC <= 25 then k3 = 10^(-5.86 - 317.0 / TK )
60 if TC > 25 then k3 = 10^(-1.1 - 1737.0 / TK )
70 t = 1
80 if M0 > 0 then t = M/M0
90 if t = 0 then t = 1
100 area = PARM(1) * (t)^PARM(2)
110 rf = k1*ACT("H+")+k2*ACT("CO2")+k3*ACT("H2O")
120 rem 1e-3 converts mmol to mol
130 rate = area * 1e-3 * rf * (1 - 10^(2/3*si_cc))
140 moles = rate * TIME * PARM(3)
200 save moles
-end

Thanks,
M

[diederik]

Dear,

the basic variable TIME in the rate equation does have the same units as defined in the GUI of HYDRUS. It should have the value of the current time step in hydrus (e.g. 0.02 units). I do not understand why you have always 1 (e.g. use 10 punch time in user_punch)?

Because TIME is in the same units as defined in the GUI, all parameters with unit [T] should be adapted accordingly. Thus, if your rate constant is in s-1, you should multiply it with 3600 to have the value in the same unit as the GUI of HYDRUS.

Diederik

[mpwickham]

D,

Thank you for confirming the requirement for a time unit conversion in RATE from the units of the rate expression to what units were selected in HP1. The second matter raised I've isolated to the use of -cvode. Printing TIME from within the RATE expression returns a "1" when using -cvode for integration - computations using TIME with -cvode appear to be using the correct time step size. Printing TIME from within the RATE expression returns the time step size when using the default Runge-Kutta.

M

[diederik]

Hi,

So the rate equation is solved correctly, but the TIME basic function prints something unexpected. Correct?
I have looked in the source code, and I do not see immediately why it returns 1 with -cvode true. Note that this is also the case when running PHREEQC only (e.g. PHREEQCI). In following program, I also get 1 as an output:

SOLUTION_MASTER_SPECIES
Pola Pola 0.0 Pola 1.0

SOLUTION_SPECIES
Pola = Pola; log_k 0.0


RATES
degradation
-start
10 rem parm(1) first-order degradation coefficient (day-1)
30 moles=-parm(1)*TIME*TOT("water")*MOL("Pola")
35 print total_time TIME
40 SAVE moles
-end

solution 3001 upper boundary solution
-units mol/kgw
Pola 1

KINETICS 3001
degradation
-formula Pola 1
-m 1
-m0 1
-parms 0.02
-tol 1e-08
-steps 1 2 3 4 5
-step_divide 1
-runge_kutta 3
-bad_step_max 500
-cvode true
-cvode_steps 100
-cvode_order 5