On the right usage of KNOBS

[Flopi]

Dear Forum,
I am having some convergence issues (I am pushing the poor PHREEEQC to its limits I'm afraid) and I know that I should be using KNOBS to help converging. However, it is not clear to me how the values shoud be modified and on which direction. How to know if I should be increasing of decreasing the values or -iter -convergence_tolerance -tolerance -step_size -pe_step_size and so on ?

I have been playing with the numbers but haven't found a proper method to be successful.

Thanks for your help

(This post is related to a post I wrote in the 'Precipitatation and dissolution' section called "multiple questions on hydrothermal modeling"

[dlparkhurst]

I use KNOBS only occasionally. The most common use is to set -step 5, -pe 2. These settings allow the numerical method to take smaller steps to avoid possible overshoot and undershoot in finding the solution. If you set values smaller than these, you may need to increase the number of iterations.

Diagonal_scale may be useful if you are dealing with very small concentrations (< 1e-12).

Often the problem is not with convergence parameters, but a lack of redox buffering. In this case, adding this aqueous species can help:

SOLUTION_SPECIES
H2O + 0.01e- = H2O-0.01
log_k -9

It adds a redox species with a concentration of about 1e-9 regardless of redox conditions, which makes the equations easier to solve.

If you provide a simple example of your convergence problem, I will try to give you an analysis of the numerical issues.

[Flopi]

The problem I'm currently on is the following:

Code: [Select]

PHASES

Olivine-Fo90
Mg1.8Fe0.2SiO4 + 4 H+ = 1.0H4SiO4 + 1.8Mg+2 + 0.2Fe+2
log_k 27.646
-analytic -7.3366336e2 -0.1113871 4.90932603e4 2.63819198e2 -2.0395157e5

Cronstedtite-7A #From Carbfix
Fe2Fe2SiO5(OH)4 + 10 H+ = H4SiO4 + 2 Fe+2 + 2 Fe+3 + 5 H2O
log_k 16.2603
-delta_H -244.266 kJ/mol
-analytic -2.3783e2 -7.1026e-2 1.7752e4 8.7147e1 2.7707e2
-Vm 110.9 # HDN+78

Opx-En90
Mg0.9Fe0.1SiO3 + 2 H+ + 1.0H2O = 0.9Mg+2 + 0.1Fe+2 + 1.0H4SiO4
log_k 11.46
-analytic -3.7593221e2 -0.0573594 2.41700954e4 1.35369045e2 -1.0253772e5

Cpx-Di90
Mg0.90Fe0.1CaSi2O6 + 4H+ + 2H2O = 0.9Mg+2 + 0.1Fe+2 + 1Ca+2 + 2H4SiO4
log_k 21.555
-analytic -6.75640e2 -0.10647 4.251246e4 2.451167e2 -1.79344068e6

#SOLUTION_SPECIES
#H2O + 0.01e- = H2O-0.01
#log_k -9

SOLID_SOLUTIONS 1
Olivine_II
-comp Forsterite 0
-comp Fayalite 0

Opx_II
-comp Enstatite(alpha) 0
-comp Ferrosilite(alpha) 0

Cpx_II
-comp Diopside 0
-comp Hedenbergite 0

Serpentine-II
-comp Lizardite 0
-comp Greenalite 0
-comp Cronstedtite-7A 0
-comp Kaolinite 0

Brucite_Fe
-comp Brucite 0
-comp Fe(OH)2 0

Talc_Fe
-comp Talc 0
-comp Minnesotaite 0

Chlorite
-comp Clinochlore 0
-comp Chamosite(Daphnite) 0

Tremolite_SS
-comp Tremolite 0
-comp Ferrotremolite 0

END

SOLUTION 1
    temp      25
pressure 500
    pH        7.8
    pe        0
    redox     pe
    units     mmol/kgw
    density   1
    Al        0.000037
    C(4)   2.34
    Ca        10.2
    Cl        546.0
    Fe(2)     0.0000015
    K         9.8
    Mg        53.0
    Na        464.0
    O(0)      0.25
    Si        0.11
    -water   1 #kg
SAVE SOLUTION 1
END

USE SOLUTION 1
USE SOLID_SOLUTION 1

REACTION_TEMPERATURE 1
    25 300 in 12 steps

EQUILIBRIUM_PHASES 
Olivine-Fo90 0 6.12 dissolve_only
Opx-En90 0 0 dissolve_only
Cpx-Di90 0 0.455 dissolve_only
Magnetite 0 0
Hematite 0 0
Chalcedony 0 0
H2(g) 2.7 0

KNOBS
-step 5
-pe 2
-iter 500
-logfile
-tolerance 1e-6

I'm working with the thermoddem database. With a bit of tweaking of the tolerance and the -step and -pe you gave, it works most of the time, but for this one, I can't manage to get a convergence at 300°C.

[dlparkhurst]

PHREEQC is not so good at determining when a solid solution should form. It seems to do better with finding equilibrium with a solid solution that is already present. I added a small amount of each solid solution component and decreased the step sizes and it struggled through.

Simplifying the problem by removing solid solutions that do not form would help convergence; same with equilibrium phases. We must be on some other planet if you specified a maximum P(H2) of nearly 500 atm.

BTW, the thermoddemv database that I have does not have any molar volumes defined, so the pressure specification does not have much, if any, effect.


Also, I do not like to increase the -tolerance and -convergence_tolerance because it may produce a numerical solution that is not a good representation of the original problem definition.

Code: [Select]

PHASES

Olivine-Fo90
Mg1.8Fe0.2SiO4 + 4 H+ = 1.0H4SiO4 + 1.8Mg+2 + 0.2Fe+2
log_k 27.646
-analytic -7.3366336e2 -0.1113871 4.90932603e4 2.63819198e2 -2.0395157e5

Cronstedtite-7A #From Carbfix
Fe2Fe2SiO5(OH)4 + 10 H+ = H4SiO4 + 2 Fe+2 + 2 Fe+3 + 5 H2O
log_k 16.2603
-delta_H -244.266 kJ/mol
-analytic -2.3783e2 -7.1026e-2 1.7752e4 8.7147e1 2.7707e2
-Vm 110.9 # HDN+78

Opx-En90
Mg0.9Fe0.1SiO3 + 2 H+ + 1.0H2O = 0.9Mg+2 + 0.1Fe+2 + 1.0H4SiO4
log_k 11.46
-analytic -3.7593221e2 -0.0573594 2.41700954e4 1.35369045e2 -1.0253772e5

Cpx-Di90
Mg0.90Fe0.1CaSi2O6 + 4H+ + 2H2O = 0.9Mg+2 + 0.1Fe+2 + 1Ca+2 + 2H4SiO4
log_k 21.555
-analytic -6.75640e2 -0.10647 4.251246e4 2.451167e2 -1.79344068e6

SOLUTION_SPECIES
H2O + 0.01e- = H2O-0.01
log_k -9

SOLID_SOLUTIONS 1
Olivine_II
-comp Forsterite 1e-10
-comp Fayalite 1e-10

Opx_II
-comp Enstatite(alpha) 1e-10
-comp Ferrosilite(alpha) 1e-10

Cpx_II
-comp Diopside 1e-10
-comp Hedenbergite 1e-10

Serpentine-II
-comp Lizardite 1e-10
-comp Greenalite 1e-10
-comp Cronstedtite-7A 1e-10
-comp Kaolinite 1e-10

Brucite_Fe
-comp Brucite 1e-10
-comp Fe(OH)2 1e-10

Talc_Fe
-comp Talc 1e-10
-comp Minnesotaite 1e-10

Chlorite
-comp Clinochlore 1e-10
-comp Chamosite(Daphnite) 1e-10

Tremolite_SS
-comp Tremolite 1e-10
-comp Ferrotremolite 1e-10

END

SOLUTION 1
    temp      25
pressure 500
    pH        7.8
    pe        0
    redox     pe
    units     mmol/kgw
    density   1
    Al        0.000037
    C(4)   2.34
    Ca        10.2
    Cl        546.0
    Fe(2)     0.0000015
    K         9.8
    Mg        53.0
    Na        464.0
    O(0)      0.25
    Si        0.11
    -water   1 #kg
SAVE SOLUTION 1
END

USE SOLUTION 1
USE SOLID_SOLUTION 1

REACTION_TEMPERATURE 1
    25 300 in 12 steps

EQUILIBRIUM_PHASES
Olivine-Fo90 0 6.12 dissolve_only
Opx-En90 0 0 dissolve_only
Cpx-Di90 0 0.455 dissolve_only
Magnetite 0 0
Hematite 0 0
Chalcedony 0 0
H2(g) 2.7 0

KNOBS
-step 1.2
-pe 1.2
-iter 1000


[Flopi]

We must be on some other planet if you specified a maximum P(H2) of nearly 500 atm

I actually tried to verify if the H2 would exsolve from the solution at 500 bar (usual pressure for hydrothermal reactions). Was that the correct way ?

[dlparkhurst]

Yes, you set it to exsolve at 500 atm. Still seems unlikely; you would be outside the stability field of water and the limits of PHREEQC.

[Micheal]

Could you please share details of how to plot the results from this code. Especially plots of hydrogen generation and others. I am new to PHREEQC and I am trying to figure out how to generate hydrogen from a reaction between olivine (forsterite and fayalite). Other reaction products are Magnetite, Brucite, and silica. My model has 80% olivine, 15% orthopyroxene, and 5% clinopyroxene. It's somehow similar to this code only that I don't know how to proceed with the plotting. Thanks

[dlparkhurst]

You will have to provide a simple script and describe a specific plot that you want to make.