running suddenly stopped

[nurlan]

Hello,
I simulated RTM with PHREEQC, running initiated but suddenly stopped at ~3500 years. Command line showed some warnings regarding bad time max. But increasing it can't solve the problem. I am not sure if I can rely on simulation results got for 3500 years.
Please, can you help with it? I can't find the problem though.
Thank you a lot
The code is

Code: [Select]

DATABASE c:\phreeqc\database\llnl.dat
SOLUTION 0 #solution injected into the 1st cell 

temp 40
pH 3.37
pe -3.1
redox pe
units mol/kgw
density 1.1596
pressure 88.8

Al 1.15E-08
C 9.88E-01
Ca 3.86E-02
Cl 5.24E+00
Fe 6.84E-04
K 4.02E-02
Mg 6.32E-02
Na 4.90E+00
S 1.17E-06
Si 1.84E-04
Sr 7.82E-04
   
water 1 

#All rate data from Palandri and Kharaka, except for dawsonite...
 Equilibrium_phases 1
 CO2(g) 1.75 1000
 calcite 0 21.9
  Equilibrium_phases 2-30
calcite 0 21.9
SOLUTION 1-30 #solution into which CO2 is dissolved

temp 40
pH 6.55
pe -3.1
redox pe
units mol/kgw
density 1.1596
pressure 88.8

Al          1.153e-08 
C                 9.941e-03 
Ca                3.857e-02 
Cl                5.237e+00 
Fe                6.842e-04 
K                 4.024e-02 
Mg                6.319e-02 
Na                4.903e+00 
S                 1.170e-06 
Si                1.843e-04 
Sr                7.824e-04 
water 1 
KINETICS 1-30
-cvode

Quartz
-m0    158.2
-parms  60.08  0.111 #https://escholarship.org/content/qt3h77d5bj/qt3h77d5bj_noSplash_925369b5813359fe46019e017e739492.pdf
Kaolinite
-m0    3.1
-parms  258.16   13.2 #https://escholarship.org/content/qt3h77d5bj/qt3h77d5bj_noSplash_925369b5813359fe46019e017e739492.pdf

Albite
-m0 9
-parms 263.02  0.164 #https://escholarship.org/content/qt3h77d5bj/qt3h77d5bj_noSplash_925369b5813359fe46019e017e739492.pdf
Dawsonite
-m0     0.0
        -parms 144.0  9.8 #hallevang
Magnesite
          -m0     0.177
          -parms   84.31  0.8 # https://www.sciencedirect.com/science/article/pii/S0883292714002315#b0230
 
Muscovite
-m0    0
-parms   398.71 0.68
INCREMENTAL_REACTIONS true

RATES
Quartz
-start
10 reacS = parm(1)*parm(2)*m
13 n = 0.0
15 logkH = -13.99
    16 Ea = 87700
    17 R = 8.3145
    20 logkHT = logkH - (Ea/(2.303*R))*(1/Tk - 1/298.15)
    21 kHT = 10^logkH
22 k_p = (act("H+")^n)*kHT/1
23 if SR("Quartz") > 1 then expTerm = 1/((Tk^3)*(LOG(SR("Quartz")))^2)
    24 k_n = 1
25 gamma_n = 5e10
26  lamda = 1
    30  if SR("Quartz") <= 1 then moles = reacS*kHT*(act("H+")^n)*(1 - SR("Quartz"))*time
    40  if SR("Quartz") > 1 then moles = (-reacS*lamda*k_p*(SR("Quartz") - 1)^2 - k_n*exp(-gamma_n*expTerm))*time
    100 save moles
-end
Kaolinite
-start
10 reacS = parm(1)*parm(2)*m
13 n = 0.777
15 logkH = -11.31
    16 Ea = 65900
    17 R = 8.3145
    20 logkHT = logkH - (Ea/(2.303*R))*(1/Tk - 1/298.15)
    21 kHT = 10^logkH
22 k_p = (act("H+")^n)*kHT/10
23 if SR("Kaolinite") > 1 then expTerm = 1/((Tk^3)*(LOG(SR("Kaolinite")))^2)
    24 k_n = 1
25 gamma_n = 2e10
26  lamda = 1
    30  if SR("Kaolinite") <= 1 then moles = reacS*kHT*(act("H+")^n)*(1 - SR("Kaolinite"))*time
    40  if SR("Kaolinite") > 1 then moles = (-reacS*lamda*k_p*(SR("Kaolinite") - 1)^2 - k_n*exp(-gamma_n*expTerm))*time
    100 save moles
-end


Albite
-start
10 reacS = parm(1)*parm(2)*m
13 n = 0.457
15 logkH = -10.16
    16 Ea = 65000
    17 R = 8.3145
    20 logkHT = logkH - (Ea/(2.303*R))*(1/Tk - 1/298.15)
    21 kHT = 10^logkH
22 k_p = (act("H+")^n)*kHT/1
23 if SR("Albite") > 1 then expTerm = 1/((Tk^3)*(LOG(SR("Albite")))^2)
    24 k_n = 1
25 gamma_n = 2e10
26  lamda = 1
    30  if SR("Albite") <= 1 then moles = reacS*kHT*(act("H+")^n)*(1 - SR("Albite"))*time
    40  if SR("Albite") > 1 then moles = (-reacS*lamda*k_p*(SR("Albite") - 1)^2 - k_n*exp(-gamma_n*expTerm))*time
    100 save moles
-end

Dawsonite
-start
10 reacS = parm(1)*parm(2)*m
13 n = 0.98
15 logkH = -4.5
    16 Ea = 63820
    17 R = 8.3145
    20 logkHT = logkH - (Ea/(2.303*R))*(1/Tk - 1/298.15)
    21 kHT = 10^logkH
22 k_p = (act("H+")^n)*kHT
23 if SR("Dawsonite") > 1 then expTerm = 1/((Tk^3)*(LOG(SR("Dawsonite")))^2)
    24 k_n = 1
25 gamma_n = 2e10
26  lamda = 1
    30  if SR("Dawsonite") <= 1 then moles = reacS*kHT*(act("H+")^n)*(1 - SR("Dawsonite"))*time
    40  if SR("Dawsonite") > 1 then moles = (-reacS*lamda*k_p*(SR("Dawsonite") - 1)^2 - k_n*exp(-gamma_n*expTerm))*time
    100 save moles
-end
Magnesite
-start
10 reacS = parm(1)*parm(2)*m
13 n = 1.0
15 logkH = -6.38
    16 Ea = 14400
    17 R = 8.3145
    20 logkHT = logkH - (Ea/(2.303*R))*(1/Tk - 1/298.15)
    21 kHT = 10^logkH
22 k_p = 7.36E-14
23 if SR("Magnesite") > 1 then expTerm = 1/((Tk^3)*(LOG(SR("Magnesite")))^2)
    24 k_n = 1
25 gamma_n = 4e10
26  lamda = 1
    30  if SR("Magnesite") <= 1 then moles = reacS*kHT*(act("H+")^n)*(1 - SR("Magnesite"))*time
    40  if SR("Magnesite") > 1 then moles = (-reacS*lamda*k_p*(SR("Magnesite") - 1)^2 - k_n*exp(-gamma_n*expTerm))*time
    100 save moles
-end
Muscovite  #Rate parameters from muscovite
-start
10 reacS = parm(1)*parm(2)*m
13 n = 0.37
15 logkH = -11.85
    16 Ea = 22000
    17 R = 8.3145
    20 logkHT = logkH - (Ea/(2.303*R))*(1/Tk - 1/298.15)
    21 kHT = 10^logkH
22 k_p = (act("H+")^n)*kHT/1
23 if SR("Muscovite") > 1 then expTerm = 1/((Tk^3)*(LOG(SR("Muscovite")))^2)
    24 k_n = 1
25 gamma_n = 2e10
26  lamda = 1
    30  if SR("Muscovite") <= 1 then moles = reacS*kHT*(act("H+")^n)*(1 - SR("Muscovite"))*time
    40  if SR("Muscovite") > 1 then moles = (-reacS*lamda*k_p*(SR("Muscovite") - 1)^2 - k_n*exp(-gamma_n*expTerm))*time
    100 save moles
-end
KNOBS
       -iterations             150
       -convergence_tolerance  1e-8
       -tolerance              1e-14
      -step_size              10
     -pe_step_size           5
    -diagonal_scale         TRUE
    -debug_diffuse_layer    TRUE
      -debug_inverse          TRUE
      -debug_model            TRUE
      -debug_prep             TRUE
      -debug_set              TRUE
      -logfile                TRUE
TRANSPORT
 -cells 30 # total length = cells X length
 -length 1.0 # total depth = 100x1.0=100m
 -shifts 1000 # =
 -time_step 3.1536e+8 # in seconds
 -flow_direction forward
 -boundary_cond flux flux
 -diffc 3e-9
 -dispersivity 0.1 # (m)
 -correct_disp true
 -punch_cells 1 15 30
 -punch_frequency 1
 -print_cells 1-30
 -print_frequency 1
SELECTED_OUTPUT
 -file rtm nurlan ebeity.sel
 -pe true
 -totals Br Ca Mg Na K Si C(4) Al Fe O(0) Fe(3) Fe(2) S(6)
  -si O2(g) CO2(g)
 -kinetic_reactants Quartz Kaolinite  Albite  Dawsonite Magnesite  muscovite
-equilibrium_phases CO2(g) calcite
END

[dlparkhurst]

You have more patience than I. The calculation runs very slowly. I think the main problem is that the implicit CVode method relies on derivatives of the rate to extrapolate. Your rates are non-differentiable at SR=1, in particular Dawsonite, which is nearly always near equilibrium. Consequently, the rate oscillates from rapid precipitation to slow dissolution as equilibrium is approached.

The simplest fix is to switch Dawsonite from KINETICS to EQUILIBRIUM_PHASES.

If you don't trust that Dawsonite should be near equilibrium, you could move the switch in rate ror Dawsonite from SR=1 to, say, SR=1.1. That will make an oscillation near equilibrium less likely.

More generally, the calculation will near steady state after a few pore volumes (30 shifts per pore volume). At least until you are satisfied with the chemistry, most of the information--concentrations vs distance and rates of dissolution and precipitation--will be available after 10s of shifts, rather than 1000 shifts. Ultimately, you may run out of one of your reactants, but you can estimate that from the rates of dissolution after a couple of pore volume.

[nurlan]

Dear Mr. Parkhurst,
Thank you a lot I ran it and it is much faster now,
Wish you a great weekend.

[nurlan]

Dear Mr Parkhurst,
I got the results for 20000 years for this and the problem is that pH shows stepwise decrease affecting the mineral rates, interestingly after each 2400 years. I can't find the exact reason for it. I read some literature but there is no mention of it.
Can it be because of numerical fail? I tried to increase tolerance but the problem is still there.
Please can you suggest the reason?
Thanks

[dlparkhurst]

Always blaming the model.

No, I think the step decreases in pH are the result of exhausting kinetic or equilibrium reactants in a cell, most likely calcite.

Diffusion is probably more important than advection in your calculation, so exhausting a mineral in a cell can affect an upstream cell.

[nurlan]

Thanks,
This means unless there is calcite in the system, pH would drop as it is,
Seems like there's not much I can do about that.
Got it,
Thanks again