Conceptual Models > Design of conceptual models

Granite water-rock reaction fails to converge

(1/1)

duanjiabin:
  I am trying to model a water-rock reaction in granite. It works when the water/rock ratio is large, and fails to converge when the water/rock ratio is small

Reaction step 1.

WARNING: Maximum iterations exceeded, 100

WARNING: Numerical method failed with this set of convergence parameters.

WARNING: Trying smaller step size, pe step size 10, 5 ...

WARNING: Maximum iterations exceeded, 200

WARNING: Numerical method failed with this set of convergence parameters.

WARNING: Trying reduced tolerance 1e-16 ...

WARNING: Maximum iterations exceeded, 200

WARNING: Numerical method failed with this set of convergence parameters.

WARNING: Trying increased tolerance 1e-14 ...

WARNING: Maximum iterations exceeded, 200

WARNING: Numerical method failed with this set of convergence parameters.

WARNING: Trying diagonal scaling ...

WARNING: Maximum iterations exceeded, 200

WARNING: Numerical method failed with this set of convergence parameters.

WARNING: Trying diagonal scaling and reduced tolerance 1e-16 ...

WARNING: Maximum iterations exceeded, 200

WARNING: Numerical method failed with this set of convergence parameters.

WARNING: Trying scaling pure_phase columns 1e-10 ...

WARNING: Maximum iterations exceeded, 200

WARNING: Numerical method failed with this set of convergence parameters.

WARNING: Trying scaling pure_phase columns and diagonal scale 1e-10 ...

WARNING: Maximum iterations exceeded, 200

WARNING: Numerical method failed with this set of convergence parameters.

WARNING: Trying increased scaling 1e-09 ...

WARNING: Maximum iterations exceeded, 200

WARNING: Numerical method failed with this set of convergence parameters.

WARNING: Skipping optimize equations for first 5 iterations ...

WARNING: Maximum iterations exceeded, 100

WARNING: Numerical method failed with this set of convergence parameters.

WARNING: Adding inequality to make concentrations greater than zero.

WARNING: Maximum iterations exceeded, 100

WARNING: Numerical method failed with this set of convergence parameters.

WARNING: Trying reduced tolerance 1e-17 ...

WARNING: Maximum iterations exceeded, 200

WARNING: Numerical method failed with this set of convergence parameters.

WARNING: Trying reduced tolerance 1e-18 ...

WARNING: Maximum iterations exceeded, 200

WARNING: Numerical method failed with this set of convergence parameters.

WARNING: The program has failed to converge to a numerical solution.

The following equations were not satisfied:
ERROR:                   Al has not converged. Total: 3.308424e-01   Calculated: 2.982420e-02   Residual: 3.010182e-01

ERROR:                   Ca has not converged. Total: 6.536400e-02   Calculated: 7.318785e-03   Residual: 5.804522e-02

ERROR:                   Cs has not converged. Total: 1.005600e-04   Calculated: 2.783960e-05   Residual: 7.272040e-05

ERROR:                   Fe has not converged. Total: 4.994480e-02   Calculated: 6.429787e-03   Residual: 4.351501e-02

ERROR:                    K has not converged. Total: 1.012304e-01   Calculated: 9.190688e-03   Residual: 9.203971e-02

ERROR:                   Li has not converged. Total: 8.044800e-04   Calculated: 8.661147e-04   Residual: -6.163469e-05

ERROR:                   Mg has not converged. Total: 4.039160e-02   Calculated: 6.778091e-03   Residual: 3.361351e-02

ERROR:                   Mn has not converged. Total: 8.380000e-04   Calculated: 1.059996e-03   Residual: -2.219956e-04

ERROR:                   Na has not converged. Total: 1.025712e-01   Calculated: 1.120810e-02   Residual: 9.136310e-02

ERROR:                   Rb has not converged. Total: 6.704000e-05   Calculated: 1.930427e-05   Residual: 4.773573e-05

ERROR:                   Si has not converged. Total: 1.276442e+00   Calculated: 1.521808e-01   Residual: 1.124261e+00

ERROR:                   Mu Ionic strength has not converged.    Residual: 1.041370e-01

ERROR:               A(H2O) Activity of water has not converged.    Residual: 4.275897e-04

ERROR:                   pH Charge balance has not converged.    Residual: 2.066848e-02

ERROR:             Hydrogen Mass of hydrogen has not converged.    Residual: -5.875737e+00

ERROR:               Oxygen Mass of oxygen has not converged.    Residual: 2.938330e+00

dlparkhurst:
Please include a simple script (using the # button) that demonstrates the problem.

duanjiabin:

--- Quote from: dlparkhurst on February 22, 2024, 04:32:57 PM ---Please include a simple script (using the # button) that demonstrates the problem.

--- End quote ---

SOLUTION_MASTER_SPECIES
    Rb            Rb+              0     85.47           85.47
    Cs            Cs+              0     132.9054        132.9054
    Al            Al+3             0     Al2O3           27
    F             F-               0     18.9984         18.9984
    H+            H+               0     1               1
    Si            H4SiO4           0     SiO2            28.1


SOLUTION_SPECIES
Rb+ = Rb+
    log_k     0
Cs+ = Cs+
    log_k     0
Cl- + Li+ = LiCl
    log_k     0.2
    delta_h   -700 kJ
F- = F-
    log_k     0
CO3-2 + 2Li+ = Li2CO3
    log_k     0.64
    delta_h   5 kJ
H4SiO4 = SiO2 + 2H2O
    log_k     2.71
    delta_h   -15.469 kJ
H+ = H+
    log_k     0
H2O = OH- + H+
    log_k     -14
    delta_h   79.92 kJ
Fe+2 + H2O = FeO + 2H+
    log_k     -13.39
    delta_h   76.45 kJ
2Fe+3 + 3H2O = Fe2O3 + 6H+
    log_k     -2.79
    delta_h   15.963 kJ
Ca+2 + H2O = CaO + 2H+
    log_k     -32.7
    delta_h   186.65 kJ
Na+ = Na+
    log_k     0
    delta_h   3 kJ
H2O + Na+ = NaOH + H+
    log_k     -14.45995
H2O + 2Na+ = Na2O + 2H+
    log_k     -67.46
2Cs+ + H2O = Cs2O + 2H+
    log_k     -89.68
    delta_h   511.91 kJ
2H+ + 2e- = H2
    log_k     -3.15
    delta_h   -1.759 kJ
Cl- + K+ = KCl
    log_k     -0.5
    delta_h   2.854 kJ
Cl- + Rb+ = RbCl
    log_k     -1.01
    delta_h   5.739 kJ
H2O + Mn+2 = MnO + 2H+
    log_k     -17.9
    delta_h   102.173 kJ
2Cl- + Fe+2 = FeCl2
    log_k     -8.93
    delta_h   50.966 kJ
Al+3 + 4H2O = Al(OH)4- + 4H+
    log_k     -22.87
    delta_h   0.62 kJ
H2O + Rb+ = RbOH + H+
    log_k     -14.26
    delta_h   81.411 kJ
2Cs+ + SO4-2 = Cs2SO4
    log_k     -0.58
    delta_h   3.334 kJ
Cl- + Cs+ = CsCl
    log_k     -1.55
    delta_h   8.867 kJ
CO3-2 + 2Cs+ = Cs2CO3
    log_k     -9.9
    delta_h   56.51 kJ
2Al+3 + H2O + 2H4SiO4 = Al2Si2O5(OH)4 + 6H+
    log_k     -6.5
    delta_h   37.086 kJ
H4SiO4 + Mn+2 = Mn(H2SiO4) + 2H+
    log_k     -12.44
    delta_h   71.008 kJ
H2O + K+ = K(OH) + H+
    log_k     -14.46
    delta_h   82.538 kJ
H2O + 2K+ = K2O + 2H+
    log_k     -84.11
    delta_h   480.073 kJ
Ca+2 + H4SiO4 = CaH3SiO4+ + H+
    log_k     -8.83
    delta_h   50.402 kJ
11Al+3 + 49H4SiO4 + 11K+ = K11Si49Al11O120:27H2O + 44H+ + 49H2O
    log_k     1.23
    delta_h   -7.021 kJ
3Ca+2 + 2H2O + 2H4SiO4 = Ca3Si2O4(OH)6 + 6H+
    log_k     -49.42
    delta_h   282.09 kJ

PHASES
Li2O
    Li2O + H2O = 2Li+ + 2OH-
    log_k     -5.31
    delta_h   -484.3 kJ
Rb2O
    Rb2O + H2O = 2OH- + 2Rb+
    log_k     -18.745
    delta_h   -772 kJ
SiO2
    SiO2 + 2H2O = H4SiO4
    log_k     -3.98
    delta_h   5.99 kJ
Al(OH)3
    Al+3 + 3H2O = Al(OH)3 + 3H+
    log_k     -16.42
    delta_h   -700 kJ
Li2CO3
    Li2CO3 = CO3-2 + 2Li+
    log_k     1.5625
    delta_h   -5 kJ
FeO
    FeO + 2H+ = Fe+2 + H2O
    log_k     13.39
    delta_h   -76.45 kJ
Fe2O3
    Fe2O3 + 6H+ = 2Fe+3 + 3H2O
    log_k     2.79
    delta_h   -15.963 kJ
CaO
    CaO + 2H+ = Ca+2 + H2O
    log_k     32.7
    delta_h   -186.65 kJ
Na2O
    Na2O + 2H+ = H2O + 2Na+
    log_k     67.46
    delta_h   -385.046 kJ
Cs2O
    Cs2O + 2H+ = H2O + 2Cs+
    log_k     89.68
    delta_h   -511.91 kJ
MgO
    MgO + 2H+ = H2O + Mg+2
    log_k     21.5841
    delta_h   -151.23 kJ
KCl
    K+ + Cl- = KCl
    log_k     -0.5
    delta_h   2.854 kJ
MnO
    Mn+2 + H2O = 2H+ + MnO
    log_k     -17.9
    delta_h   102.173 kJ
Al2O3
    Al2O3 + 6H+ = 2Al+3 + 3H2O
    log_k     19.6224
    delta_h   -258.59 kJ
CsCl
    Cl- + Cs+ = CsCl
    log_k     -1.55
    delta_h   8.867 kJ
Mn(H2SiO4)
    Mn+2 + H4SiO4 = 2H+ + Mn(H2SiO4)
    log_k     -12.44
    delta_h   71.008 kJ
K2O
    K2O + 2H+ = H2O + 2K+
    log_k     84.11
    delta_h   -480.073 kJ


SOLUTION 2
    temp      25
    pH        6.5
    pe        4
    redox     pe
    units     mg/l
    density   1
    Cs        0
    Li        0
    Rb        0
    S(6)      0
    -water    0.003 # kg
REACTION 1
    Al2O3      9.87
    CaO        3.9
    Cs2O       0.003
    Fe2O3      1.49
    FeO        0
    K2O        3.02
    Li2O       0.024
    MgO        2.41
    MnO        0.05
    Na2O       3.06
    Quartz     76.16
    Rb2O       0.002
    0.0151 moles in 1 steps
END

Here's my input. Thanks.

dlparkhurst:
Here are a few comments about why your script fails.

First, SOLUTION_SPECIES is used to define aqueous species. You have defined some aqueous species and some minerals in this data block. You should either use another database that has the elements that you need (like llnl.dat), or remove all of the minerals from SOLUTION_SPECIES and include all relevant aqueous species.

Second, you do not need to define a PHASE for each item that you have in REACTION. You can use chemical formulas like MgO in REACTION without defining MgO in PHASES.

Third (and most importantly), with REACTION, you are adding moles of elements to the solution, that is, the dissolved concentrations. You are trying to add tens or hundreds of moles of solutes without any sinks (EQUILIBRIUM_PHASES), and PHREEQC cannot deal with concentrations that high.

One calculation that would make sense is to find the stable mineral assemblage that would result from complete equilibration of the rock with water. The calculation would assume geologic time and essentially complete dissolution of the rock and reprecipitation of the stable phases (so even this is not too reasonable).

Here is a script that blindly considers all of the phases in the llnl.dat database as potential stable phases. I have modified the solution volume to be 3 liters (actually kg water) because PHREEQC works best with a solution volume near 1 liter. I have increased the reaction by a factor of 1000 as well. This is a difficult problem for PHREEQC with the many moles and many phases involved. I have incrementally added the entire reaction to help with the numerics of the calculation.

I doubt the calculation makes too much sense because all minerals were considered, many of which are unlikely to form. The final solution has only 1e-3 kg of water remaining after the stable phases have formed. Hopefully, the calculation will give you an idea about how you want to proceed. I don't see what your were trying to do with your original calculation. Perhaps you meant to include the minerals that you had in SOLUTION_SPECIES in EQUILIBRIUM_PHASES.

This calculation runs for me with llnl.dat, but I don't think it is very robust.

--- Code: ---SOLUTION 2
    temp      25
    pH        6.5
    pe        4
    redox     pe
    units     mg/l
    density   1
    Cs        0
    Li        0
    Rb        0
    S(6)      0
    #-water    0.003 # kg
    -water 3
END
USE solution 2
REACTION 1
    Al2O3      1
    CaO        1
    Cs2O       1
    Fe2O3      1
    #FeO        1
    K2O        1
    Li2O       1
    MgO        1
    MnO        1
    Na2O       1
    Quartz     1
    Rb2O       1
1e-9 moles
SELECTED_OUTPUT 2
-file equilibrium_phases.pqi
USER_PUNCH 2
10 p = SYS("phases", count , name$ , type$ , moles )
20 s$ = "EQUILIBRIUM_PHASES 1" + EOL$
30 FOR i = 1 TO count
40   s$ = s$ + name$(i) + " 0 0 " + EOL$
50 NEXT i
60 s$ = s$ + "END" + EOL$
70 PUNCH s$
END
SELECTED_OUTPUT 2
-active false
END
INCREMENTAL_REACTIONS true
USE solution 2
REACTION 1
    Al2O3      9.87
    CaO        3.9
    Cs2O       0.003
    Fe2O3      1.49
    #FeO        0
    K2O        3.02
    Li2O       0.024
    MgO        2.41
    MnO        0.05
    Na2O       3.06
    Quartz     76.16
    Rb2O       0.002
    #0.0151e3
    0.0001e3 0.005e3 0.01e3
INCLUDE$ equilibrium_phases.pqi
END

--- End code ---

duanjiabin:

--- Quote from: dlparkhurst on February 23, 2024, 05:30:54 AM ---Here are a few comments about why your script fails.

First, SOLUTION_SPECIES is used to define aqueous species. You have defined some aqueous species and some minerals in this data block. You should either use another database that has the elements that you need (like llnl.dat), or remove all of the minerals from SOLUTION_SPECIES and include all relevant aqueous species.

Second, you do not need to define a PHASE for each item that you have in REACTION. You can use chemical formulas like MgO in REACTION without defining MgO in PHASES.

Third (and most importantly), with REACTION, you are adding moles of elements to the solution, that is, the dissolved concentrations. You are trying to add tens or hundreds of moles of solutes without any sinks (EQUILIBRIUM_PHASES), and PHREEQC cannot deal with concentrations that high.

One calculation that would make sense is to find the stable mineral assemblage that would result from complete equilibration of the rock with water. The calculation would assume geologic time and essentially complete dissolution of the rock and reprecipitation of the stable phases (so even this is not too reasonable).

Here is a script that blindly considers all of the phases in the llnl.dat database as potential stable phases. I have modified the solution volume to be 3 liters (actually kg water) because PHREEQC works best with a solution volume near 1 liter. I have increased the reaction by a factor of 1000 as well. This is a difficult problem for PHREEQC with the many moles and many phases involved. I have incrementally added the entire reaction to help with the numerics of the calculation.

I doubt the calculation makes too much sense because all minerals were considered, many of which are unlikely to form. The final solution has only 1e-3 kg of water remaining after the stable phases have formed. Hopefully, the calculation will give you an idea about how you want to proceed. I don't see what your were trying to do with your original calculation. Perhaps you meant to include the minerals that you had in SOLUTION_SPECIES in EQUILIBRIUM_PHASES.

This calculation runs for me with llnl.dat, but I don't think it is very robust.

--- Code: ---SOLUTION 2
    temp      25
    pH        6.5
    pe        4
    redox     pe
    units     mg/l
    density   1
    Cs        0
    Li        0
    Rb        0
    S(6)      0
    #-water    0.003 # kg
    -water 3
END
USE solution 2
REACTION 1
    Al2O3      1
    CaO        1
    Cs2O       1
    Fe2O3      1
    #FeO        1
    K2O        1
    Li2O       1
    MgO        1
    MnO        1
    Na2O       1
    Quartz     1
    Rb2O       1
1e-9 moles
SELECTED_OUTPUT 2
-file equilibrium_phases.pqi
USER_PUNCH 2
10 p = SYS("phases", count , name$ , type$ , moles )
20 s$ = "EQUILIBRIUM_PHASES 1" + EOL$
30 FOR i = 1 TO count
40   s$ = s$ + name$(i) + " 0 0 " + EOL$
50 NEXT i
60 s$ = s$ + "END" + EOL$
70 PUNCH s$
END
SELECTED_OUTPUT 2
-active false
END
INCREMENTAL_REACTIONS true
USE solution 2
REACTION 1
    Al2O3      9.87
    CaO        3.9
    Cs2O       0.003
    Fe2O3      1.49
    #FeO        0
    K2O        3.02
    Li2O       0.024
    MgO        2.41
    MnO        0.05
    Na2O       3.06
    Quartz     76.16
    Rb2O       0.002
    #0.0151e3
    0.0001e3 0.005e3 0.01e3
INCLUDE$ equilibrium_phases.pqi
END

--- End code ---

--- End quote ---


Thank you very much for your answer, my aim is to use the software to simulate all the major components of granite dissolution. After reading your reply I will rethink the meaning of this simulation and continue to learn.

Navigation

[0] Message Index