Processes > Reactive transport modelling

Questions about KINETICS in PHREEQCRM with MODFLOW coupling

(1/1)

ZYWTDaEO6:
Hello dlparkhurst,

I am currently working with PHREEQCRM coupled with MODFLOW, and I have two questions regarding the KINETICS module and the use of the .pqi file.

1. Regarding -steps in KINETICS

In PHREEQC, when defining KINETICS, I normally need to specify time steps using the -steps keyword (e.g., -steps 100 in 100 steps). However, when I write a .pqi file for PHREEQCRM and couple it with MODFLOW, I am not sure whether -steps is still required.
From my tests, the reactive transport model seems to run both with and without -steps. Does PHREEQCRM automatically ignore the -steps definition when it is coupled with a transport model (MODFLOW)? Should I explicitly include it, or is it unnecessary in this context?

2. About my RATES and KINETICS definition

Here is a simplified example of what I wrote:


--- Code: ---KINETICS 1
Halite
    -m0  18.380
    -parms  1.0e-3
INCREMENTAL_REACTIONS true
END

RATES 1
Halite
    -start
    10 k0 = PARM(1)
    20 IF M0 <= 0 THEN M0 = 1E-20
    30 IF M > 0 THEN k_t = k0 * (M/M0)^(2/3) ELSE k_t = 0
    80 moles = k_t * (1 - SR("Halite")) * TIME
    100 IF moles > 0 AND moles > M THEN moles = M
    110 SAVE moles
    -end
END
--- End code ---

I noticed that when -parms is set too small, the mineral hardly reacts. But when I increase it slightly, the simulation often crashes with errors. Is there something wrong in the way I defined the RATES expression, or is this behavior expected due to the sensitivity of the parameterization?

Thank you very much for your time and support!

dlparkhurst:
The steps definition is used only when the input file is initially processed. Normally, you probably won't perform a KINETICS calculation in the input file calculations, but if you do react KINETICS with a solution, then the steps definition is used. After the input file is processed, the steps definition is not used.

Parm(1) is the rate constant that has units of 1/second. So the rate of approach to halite equilibrium is inversely proportional to the rate constant. The smaller the constant, the longer it takes to reach equilibrium. A rate constant of 1 will reach equilibrium in a few seconds.

In general, I try to avoid "if" statements in RATES definitions. Especially if the implicit integrator CVODE is used, it is best to have a differentiable rate definition. "if" statements tend to produce singularities and non-differentiable rates.  The reaction will stop automatically when the moles of reactant goes to zero. I would use the following, which runs smoothly for me:


--- Code: ---RATES 1
Halite
    -start
    10 k0 = PARM(1)
    #20 IF M0 <= 0 THEN M0 = 1E-20
    30 IF M > 0 THEN k_t = k0 * (M/M0)^(2/3) ELSE k_t = 0
    80 moles = k_t * (1 - SR("Halite")) * TIME
    #100 IF moles > 0 AND moles > M THEN moles = M
    110 SAVE moles
    -end
END
--- End code ---

Navigation

[0] Message Index