Model of Reaction progress with Sulphur concentration

[jt2020]

Hello,

I am looking to design a reaction path model with reaction progress (mole of reaction progress) on the x axis and sulphur concentration (mol) on the y axis.

I have previously attempted to model the reaction pathway from mackinawite, to greigite, to the end-member, pyrite. I managed to produce a phase diagram following Example 6 in the manual but when trying to design a pathway it seems that pyrite forms instantly as oppose to following the pathway through the phases

I thought it might show the path better by using reaction progress and sulphur content. I know the latter can be done in Phreeqc but can reaction progress?

Thank you,

Jack

[dlparkhurst]

Maybe this will help. Here is a rewrite of example 6 that has the moles of K-feldspar reacted as the X axis. It uses INCREMENTAL_REACTIONS, so the reaction increment at each point (Basic function RXN) needs to be accumulated through the GET/PUT functions to get the total amount of feldspar reacted at each step. If you do not use INCREMENTAL_REACTIONS, then RXN is the total amount added at each step.

Code: [Select]

TITLE React K-feldspar
SOLUTION 1  PURE WATER
            pH    7.0 charge
          temp   25.0
PHASES
        Gibbsite
                Al(OH)3 + 3 H+ = Al+3 + 3 H2O
                log_k          8.049
                delta_h        -22.792 kcal
        Kaolinite
                Al2Si2O5(OH)4 + 6 H+ = H2O + 2 H4SiO4 + 2 Al+3
                log_k          5.708
                delta_h        -35.306 kcal
        K-mica
                KAl3Si3O10(OH)2 + 10 H+ = 3 Al+3 + 3 H4SiO4 + K+
                log_k          12.970
                delta_h        -59.377 kcal
        K-feldspar
                KAlSi3O8 + 4 H2O + 4 H+ = Al+3 + 3 H4SiO4 + K+
                log_k          0.875
                delta_h        -12.467 kcal
END
INCREMENTAL_REACTIONS
REACTION
K-feldspar 1
20*1e-8 20*1e-7  20*1e-6  20*1e-5  20*1e-4  20*1e-3  20*1e-2 
USE solution 1
EQUILIBRIUM_PHASES 1
        Gibbsite        0.0     0.0
        Kaolinite       0.0     0.0
        K-mica          0.0     0.0
        K-feldspar      0.0     0.0
USER_GRAPH 1 Simulation 6
    -headings               Log(K/H) Gibbsite Kaolinite K-mica K-feldspar
    -axis_titles            "Feldspar reacted, moles" "Log([K+] / [H+])" "Mineral, moles"
    -chart_title            "K-Feldspar Reaction Path"
    -axis_scale sy_axis     auto auto auto auto log
    -initial_solutions      false
    -connect_simulations    false
    -plot_concentration_vs  x
  -start
 5 PUT(GET(1) + RXN, 1)
  10 PLOT_XY log10(GET(1)),(LA("K+")-LA("H+")), color = Red, line_w = 0, \
  symbol = Circle, symbol_size = 10
  20 PLOT_XY log10(GET(1)), EQUI("Gibbsite") color = Green, line_w = 0, \
  symbol = Plus, symbol_size = 5, y_axis = 2
  30 PLOT_XY log10(GET(1)), EQUI("Kaolinite") color = Blue, line_w = 0, \
  symbol = Plus, symbol_size = 5, y_axis = 2
  40 PLOT_XY log10(GET(1)), EQUI("K-mica") color = Cyan, line_w = 0, \
  symbol = Plus, symbol_size = 5, y_axis = 2
  50 PLOT_XY log10(GET(1)), EQUI("K-feldspar") color = Magenta, line_w = 0, \
  symbol = Plus, symbol_size = 5, y_axis = 2
  -end
    -active                 true
END


[jt2020]

Ah great thank you,

I will give it a go.

Thanks again,

Jack

[dlparkhurst]

Here is a fancy version that has evenly spaced log-scale reaction increments.

Code: [Select]

PHASES
Gibbsite
    Al(OH)3 + 3H+ = Al+3 + 3H2O
    log_k     8.049
    delta_h   -22.792 kcal
Kaolinite
    Al2Si2O5(OH)4 + 6H+ = 2Al+3 + H2O + 2H4SiO4
    log_k     5.708
    delta_h   -35.306 kcal
K-mica
    KAl3Si3O10(OH)2 + 10H+ = 3Al+3 + 3H4SiO4 + K+
    log_k     12.97
    delta_h   -59.377 kcal
K-feldspar
    KAlSi3O8 + 4H+ + 4H2O = Al+3 + 3H4SiO4 + K+
    log_k     0.875
    delta_h   -12.467 kcal
END
TITLE Simulation 6A.--React to phase boundaries.
SOLUTION 1 PURE WATER
    temp      25
    pH        7 charge
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    -water    1 # kg
SELECTED_OUTPUT 100
-file increments.pqi
USER_PUNCH 100
10 low = -8
20 up = -3
30 steps = 200
35 size = (up -low)/steps
40 FOR i = low TO up STEP size
50 PUNCH 10^i, EOL$
60 NEXT i
END
SELECTED_OUTPUT 100
-active false
END
REACTION 1
    K-feldspar 1
INCLUDE$ increments.pqi
USE solution 1
EQUILIBRIUM_PHASES 1
    Gibbsite  0 0
    K-feldspar 0 0
    K-mica    0 0
    Kaolinite 0 0
USER_GRAPH 1 Simulation 6
    -headings               Log(K+/H+) Gibbsite Kaolinite K-mica K-feldspar
    -axis_titles            "Feldspar reacted, moles" "Log([K+] / [H+])" "Mineral, moles"
    -chart_title            "K-Feldspar Reaction Path"
    -axis_scale x_axis      -8 -3 auto auto
    -axis_scale sy_axis     1e-9 0.001 auto auto log
    -initial_solutions      false
    -connect_simulations    false
    -plot_concentration_vs  x
  -start
#5 PUT(GET(1) + RXN, 1)
 10 PUT(RXN, 1)
 20 PLOT_XY LOG10(GET(1)), (LA("K+")-LA("H+")), color = red, line_w = 2, symbol = circle, symbol_size = 0
 30 IF EQUI("Gibbsite") > 0 THEN y = EQUI("Gibbsite") ELSE y = 1e-12
 40 PLOT_XY LOG10(GET(1)), y, color = green, line_w = 2, symbol = plus, symbol_size = 0, y_axis = 2
 50 IF EQUI("Kaolinite") > 0 THEN y = EQUI("Kaolinite") ELSE y = 1e-12
 60 PLOT_XY LOG10(GET(1)), y, color = blue, line_w = 2, symbol = plus, symbol_size = 0, y_axis = 2
 70 IF EQUI("K-mica") > 0 THEN y = EQUI("K-mica") ELSE y = 1e-12
 80 PLOT_XY LOG10(GET(1)), y, color = cyan, line_w = 2, symbol = plus, symbol_size = 0, y_axis = 2
 90 IF EQUI("K-feldspar") > 0 THEN y = EQUI("K-feldspar") ELSE y = 1e-12
100 PLOT_XY LOG10(GET(1)), y, color = magenta, line_w = 2, symbol = plus, symbol_size = 0, y_axis = 2
  -end
    -active                 true
END


[jt2020]

Hello again,

Thanks for the log scaled diagram. It looked great and I have got it to work for the mackinawite-pyrite pathway.

I've tried, more recently, to input reactions using elemental S but when trying to run the model there are two errors.

ERROR: Elements in species have not been tabulated, Fe3S4.
ERROR: Reaction for species has not been defined, Fe3S4.
ERROR: Elements in species have not been tabulated, FeS.
ERROR: Reaction for species has not been defined, FeS.
ERROR: Elements in species have not been tabulated, S.
ERROR: Reaction for species has not been defined, S.
ERROR: Calculations terminating due to input errors.

I think that there is something wrong with the database I've used. Im not sure if these reactions are suitable as they show a transformation from the old phase, i.e. mackinawite, to the new phase, i.e. pyrite. Plus, I'm not sure if elemental S is included correctly (I think it needs a '0' next to it).

I have attached the database and input file.

Thanks

Jack

[dlparkhurst]

I would not change the database file (minteq.dat). Instead, simply put any new definitions or redefinitions in the input file. I think that will be more clear as to what reactions and log Ks you are using in your calculations.

Your reaction for pyrite does not balance.

Typically, the reaction for a phase is simply the dissolution reaction for the phase in terms of aqueous species. You can use other PHASE definitions in your reactions, but you must identify them with (s). In addition, the coefficient of the phase--the first entity on the left-hand side of the equation--must be 1.0. The definition for Sulfur is from minteq.dat. I am not sure this is what you intended, but the following definitions will work:

Code: [Select]

PHASES
Sulfur
    S + H+ + 2e- = HS-
    log_k     -2.11
    delta_h   -4.2 kcal
Mackinawite
    FeS + H+ = Fe+2 + HS-
    log_k     -2.137
Greigite
    Fe3S4 = 3Mackinawite(s) + Sulfur(s)
    log_k     -6.475
Pyrite
    FeS2 = 0.3333333333Greigite(s) + 0.6666666666Sulfur(s)
    log_k     -31.1184
END


[jt2020]

Hello,

Thanks so much for getting back to me. I just need a model showing a pathway of phases by adding elemental S so I think the reaction here should be OK, considering the coefficient on the left-hand side has to be 1.

Im having to change the saturation index to get a pathway from greigite to pyrite but an SI of 12 seems very high. How realistic is this in nature and how does a delay in precipitation happen like that? Id think that as soon as a solution becomes supersaturated with respect to a mineral, that mineral would then precipitate out straight away to reach equilibrium.

Thanks again

Jack

[jt2020]

In addition, I realized after running the model that logIAP for greigite and pyrite remains at 0.0 (probably because the reactions only contain solids and not ions such as Fe2+ and HS-). However, pyrite, and on removing pyrite, greigite, start to precipitate.

It is my understanding that the mineral should precipitate when logIAP = logK therefore SI = 0.0. If logIAP is at 0.0 and cant be calculated, what is causing the mineral to precipitate in the model?

Thanks

Jack

[dlparkhurst]

Although the reactions are written using solids [(s)], that is just a shorthand to simplify the definitions. For the reaction

Code: [Select]

Fe3S4 = 3Mackinawite(s) + Sulfur(s)

PHREEQC rewrites the equation by adding the following reactions for Mackinawite(s) and Sulfur(s):

Code: [Select]

    Mackinawite(s) + H+ = Fe+2 + HS-
    Sulfur(s) + H+ + 2e- = HS-

using the appropriate stoichiometry and log Ks. So, the ultimate equation is a dissolution reaction for Greigite to aqueous species. The SI is calculated by using  this rewritten equation; so, the SI can always be calculated.

Here's my take on your system. Pyrite is always more stable than greigite or mackinawite, so any calculations that produce greigite or mackinawite represent a metastable situation. I suspect mackinawite forms in rapidly sulfate reducing environments, but would alter to pyrite with time. I'm not sure when you would expect to find Greigite, but maybe under a different redox regime. You are setting up the metastable situation in your simulation by including only greigite and mackinawite in EQUILIBRIUM_PHASES which may be reasonable under some conditions, but the SI for pyrite is always supersaturated.

[jt2020]

Hi,

Ah I see, so I have to substitute those reactions into the greigite reaction, then use the concentrations in the model to give me logIAP?

Your take is in keeping with the literature. Experimentally, it has been shown that by adding elemental S, there is a pathway from mackinawite to greigite, which then ultimately forms pyrite. Greigite is seen in nature so it would be nice to see what conditions help it be preserved, whether that's a delay to the formation of pyrite (where inhibitors may prevent the transformation until a certain level of supersaturation with respect to pyrite- in this case 12), or greigite may form as a result of a different oxidant being used (as it seems using HS- in a reaction gives a much lower logK for greigite, then the reaction including S0).

In terms of length of time, it would be nice to show the moles of mack/greig/pyrite in a system with time- this may perhaps be a better model. I presume PHREEQC has these capabilities?

Thanks again for all your help,

Jack

[dlparkhurst]

You can model a kinetically controlled reaction with KINETICS and RATES, where RATES is used to define the rate expression for the reaction. Without some experimental data or some rate estimates from a chemical system, the rate expression is completely unconstrained.

[jt2020]

Hello Dr Parkhurst


Thank you so much for your help with my models.


I am producing some nice models of the pyrite reaction pathway using different sulphur species as the oxidant. I have come across some articles about pyrite being formed via a polysulfide pathway which I hope to look at. I wonder if you could guide me in the correct direction.


I have found some reactions and balanced them but unfortunately, the model won't work. It says the reactions aren't balanced but I'm pretty sure they are. I wonder if you had any suggestions? Do I need to add some phases? I have attached my input file and below are my reactions.





SOLUTION 1  PURE WATER
            pH    7.0 charge
          temp   25
PHASES
Mackinawite
    FeS + H+ = Fe+2 + HS-
    log_k     -2.137
    delta_h   -3.578 kcal
Greigite
    Fe3S4 + S5-2 = 3Mackinawite(s) + S6-2   
    log_k     -6.721
    delta_h   38.521 kcal
Pyrite
    FeS2 + 0.33333333333333S3-2 = 0.3333333333333Greigite(s) + 0.3333333333333S5-2
    log_k     -10.075
    delta_h   7.23 kcal


Thanks again,


Jack

[dlparkhurst]

I don't do polysulfides.

The problem with the reaction goes back to the definitions of S5-2 and S6-2. In wateq4f.dat, the reactions are

Code: [Select]

SOLUTION_SPECIES
#S5-2                505
        HS- = S5-2 + H+                 # (lhs) +4S
        log_k           -9.595
        delta_h 9.3 kcal
        -no_check
        -mass_balance   S(-2)5
        -gamma  12.0    0.0
#S6-2                506
        HS- = S6-2 + H+                 # (lhs) +5S
        log_k           -9.881
        -no_check
        -mass_balance   S(-2)6
        -gamma  14.0    0.0

The reactions are based on the assumption that the solution is in equilibrium with native sulfur, and native sulfur has been left out of the equation. When these equations are substituted into your PHASES definitions for greigite and pyrite, the equations for the phases are not mole balanced on sulfur.

You can get the script to run by using -no_check in the definitions for greigite and pyrite.

llnl.dat has different definitions for the polysulfide species that are mole balanced, but assume the polysulfides are in equilibrium based on the solution pe. However, S6-2 is not defined for llnl.dat, so you would have to add it. Note that pe or O2 does not appear in the wateq4f.dat definitions. I think you can redefine the polysulfides for wateq4f.dat to be like llnl.dat (equilibrium based on solution pe) if you rewrite the equations, for example,

Code: [Select]

S-2 + 5Sulfur = S6-2
        log_k           -9.881

and then rewrite using the reaction for native sulfur:

Code: [Select]

        S + 2e- = S-2
        log_k           -15.026
        delta_h 7.9 kcal

You would have to rewrite all of the polysulfide SOLUTION_SPECIES definitions.

Whether -no_check, llnl.dat, or rewriting wateq4f.dat makes any sense you will have to decide for yourself.

[jt2020]

Hi Dr Parkhurst,

Thanks for all your help over the past few months. I have produced a series of really good models showing the pyrite reaction pathway.

I have come across another pathway that may lead to the preservation of metastable phases. It involves nitrate but Im not sure how to incorporate it into the models. I am using the MINTEQ database and have added some N2 phases from PHREEQ but whenever I run the input file is says "ERROR: Elements in species have not been tabulated, NO2". Its the same of N2 as well.

I have rearanged the reactions so that the phase of interest is on the left. I've also made sure everything in the database matches everything in the input file. I wonder if you could help me to understand the issue? Below is the input file.

DATABASE C:\Users\jt2020\Dropbox\PC\Desktop\llnlpyrrminteq.dat
DATABASE C:\Users\jt2020\Dropbox\PC\Desktop\llnlpyrrminteq.dat
DATABASE C:\Users\jt2020\Dropbox\PC\Desktop\llnlpyrrminteq.dat
DATABASE C:\Users\jt2020\OneDrive - Imperial College London\Data\Thermodynamics\Thermodynamics\Phreeq\Reaction Progress\Nitrate and Fe loss\llnlpyrrminteq.dat
DATABASE C:\Program Files (x86)\USGS\Phreeqc Interactive 3.7.1-15876\database\phreeqc.dat
DATABASE C:\Users\jt2020\OneDrive - Imperial College London\Data\Thermodynamics\Thermodynamics\Phreeq\Reaction Progress\Nitrate and Fe loss\llnlpyrrminteq.dat
DATABASE C:\Users\jt2020\OneDrive - Imperial College London\Data\Thermodynamics\Thermodynamics\Phreeq\Reaction Progress\Nitrate and Fe loss\llnlpyrrminteq.dat
TITLE React K-feldspar
SOLUTION 1  PURE WATER
            pH    7.0 charge
          temp   25.0
             
PHASES
       pH_fix
                H + = H +
                log_k = 0.0
       N2
             N2 = N2
             log_k       -3.1864
       NO2
                NO2 = 0.5N2 + O2
                log_k   8.9857
                delta_h -7.91421  kcal
       Mackinawite   
                FeS + H+ = Fe+2 + HS-
                log_k          -2.1371
                delta_h        -3.578 kcal
        Greigite
                Fe3S4 + NO2 + Fe+2 + H2O = 4Mackinawite(s) + NO3- + 2H+
                log_k   -19.62791199
                delta_h  48.2874405    kcal
        Pyrite
                FeS2 + 0.5Fe+2 + NO2 + H2O = 0.5Greigite(s) + NO3- + 2H+
                log_k          -22.44229642
                delta_h        11.59814325    kcal
END
INCREMENTAL_REACTIONS
REACTION
Fe4S1 1
20*1e-12 20*1e-10  20*1e-8  20*1e-6  20*1e-4  20*1e-3  20*1e-2 mol
USE solution 1
EQUILIBRIUM_PHASES 1
        Mackinawite     0.0     0.0
        Greigite        0.0     0.0
        Pyrite          0.0     0.0
        pH_fix         -6.999     HCl      10.0
         
USER_GRAPH 1 Simulation 6
    -headings               Pyrite Greigite Mackinawite
    -axis_titles            "Fe and S added, Log moles" "Mineral, moles"
    -chart_title            "Pyrite pathway in an iron dominant system with HS- as the oxidant"
    -initial_solutions      false
    -connect_simulations    false
    -axis_scale y_axis      auto auto auto auto Log
    -axis_scale sy_axis       auto auto auto auto Log
    -plot_concentration_vs  x
  -start
 5 PUT(GET(1) + RXN, 1)
  20 PLOT_XY log10(GET(1)), EQUI("Pyrite") color = Green, line_w = 0, \
     symbol = Plus, symbol_size = 5, y_axis = 2
  30 PLOT_XY log10(GET(1)), EQUI("Greigite") color = Blue, line_w = 0, \
     symbol = Plus, symbol_size = 5, y_axis = 2
  40 PLOT_XY log10(GET(1)), EQUI("Mackinawite") color = Cyan, line_w = 0, \
     symbol = Plus, symbol_size = 5, y_axis = 2
  -end
    -active                 true
END

Thanks and Merry Christmas!

Jack

[dlparkhurst]

In the reaction equation in PHASES, the default is to assume all of the entities are aqueous species. If one of the entities is defined in PHASES, it is necessary to append either (g) or (s) so that the correct list (PHASES instead of SOLUTION_SPECIES definitions) is searched when rewriting the equation.

Still, I'm not sure whether you intended NO2- or NO2(g) in your definitions of Greigite and Pyrite. Charge is balanced with NO2(g) with Pyrite and NO2- for Greigite.

 
Although not required, the convention in the main databases for PHASES  is to identify gas components with (g). Whether you use the convention or not, the (g) or (s) is required as described above.

Code: [Select]

       N2(g)
             N2 = N2
             log_k       -3.1864
       NO2(g)
                NO2 = 0.5N2 + O2
                log_k   8.9857
                delta_h -7.91421  kcal


[jt2020]

Hi,

Thanks so much for the quick reply. It looks like its working. I had to change the N2 in the database to N2(g) and had to rebalance the equation.

Thanks so much!

Jack

[jt2020]

Hello Dr Parkhurst,

I have been developing a Fe sulphide model for the past few months as I want to test the stability of different minerals under specific conditions, including Fe vs S concentration, pH, and pe. I am using buffers so that pH and pe remain fixed and adding Fe and S to the solution at different ratios to see what precipitates out.

I have a question regarding the fixing of pH and pe. I used O2(g) to fix pe but with pH it appears a little more complicated. I have found suitable buffers using pKa values but some of the buffers are not in the database. In order to get around this, I have inserted the formula of the correct buffer into the input file. I understand this means that the buffer won't complex like it should (as its not specified in the database), however, it seems to work in fixing the pH.

In the example, I've used glyclyglycine as the buffer to fix a pH of 6 and kept pe at -4.23 using O2. As I cant see glyclyglycine in the minteq database I typed in its formula C4H8N2O3. Up to a certain point, the pH remains fixed but after that, it starts to change (see step 90 in output file). Up to this point (step 89), C and N remain at 2.48 M and 1.24 M respectively, in solution, however, at step 90 C and N disappear from solution or become very low, and the pH decreases. I cant see where C and N is being used up or where they go? Why does this sudden change occur?

I have used -force to avoid this problem but I'm not sure how accurate this method of buffering is and whether I should actually be using buffers seen in the database (as opposed to writing formulas). I'd like your opinion on how legitimate this method of buffering is? I assume I probably need to use buffers, such as Citrate, that are in the database in order for them to complex with the other ions present as this would be more applicable to natural systems.

Thanks for your time

Jack

Here is the input file. For some reason the output file was not 'writable' even when I tried to attach a CSV file. I've had to paste in step 89 and 90 of the output file to show the data.

TITLE React K-feldspar
SOLUTION 1  PURE WATER
            pH    6 charge
          temp   25.0
           pe     -4.23       #-250 mV
PHASES
       pH_fix
                H + = H +
                log_k = 0.0
       pe_fix
                e- = e-
                log_k 0.0
Mackinawite
    FeS + 2H+ = Fe+2 + H2S
    log_k     4.850
    delta_h        -8.901 kcal
Greigite
    Fe3S4 + H2 = 3Mackinawite(s) + H2S
    log_k     -1.633
    delta_h        28.408 kcal
Pyrite
    FeS2 + 0.6666666666H2 = 0.3333333333Greigite(s) + 0.6666666666H2S
    log_k     -7.144
    delta_h        0.907783 kcal
END
INCREMENTAL_REACTIONS
REACTION
Fe1S3 1
10*1e-12 10*1e-11 10*1e-10 10*1e-9 10*1e-8 10*1e-7 10*1e-6 10*1e-5 10*1e-4  10*1e-3  10*1e-2 mol
USE solution 1
EQUILIBRIUM_PHASES 1
        Mackinawite     0.0     0.0
        Greigite        0.0     0.0
        Pyrite          0.0     0.0   
        pH_fix          6       C4H8N2O3  0.6    #formula for Glyclyglycine
        pe_fix          4.23    O2(g)     10.0

USER_GRAPH 1 Simulation 6
    -headings               Pyrite Greigite Mackinawite
    -axis_titles            "Fe and S added, Log moles" "Mineral, moles"
    -chart_title            "Sulfur dominant system with H2S as the oxidant. Fix both -250mV at pH 6"
    -axis_scale sy_axis     auto auto auto auto log
    -initial_solutions      false
    -connect_simulations    false
    -axis_scale y_axis      auto auto auto auto Log
    -axis_scale sy_axis    auto auto auto auto Log
    -plot_concentration_vs  x
  -start
 5 PUT(GET(1) + RXN, 1)
  20 PLOT_XY log10(GET(1)), EQUI("Pyrite") color = Green, line_w = 0, \
  symbol = Plus, symbol_size = 5, y_axis = 2
  30 PLOT_XY log10(GET(1)), EQUI("Greigite") color = Blue, line_w = 0, \
  symbol = Plus, symbol_size = 5, y_axis = 2
  40 PLOT_XY log10(GET(1)), EQUI("Mackinawite") color = Cyan, line_w = 0, \
  symbol = Plus, symbol_size = 5, y_axis = 2
  -end
    -active                 true
SELECTED_OUTPUT
    -file                 250 S both fix.csv
    -saturation_indices   Mackinawite  Greigite  Pyrite
    -reaction           
    -totals               Fe S C N
    -equilibrium_phases   Mackinawite  Greigite  Pyrite
END


Reaction step 89.

Using solution 1. PURE WATER
Using pure phase assemblage 1.
Using reaction 1.

Reaction 1.

 1.000e-04 moles of the following reaction have been added:

                Relative
Reactant            moles

Fe1S3                1.00000

                Relative
Element             moles
Fe                   1.00000
S                    3.00000

-------------------------------Phase assemblage--------------------------------

                                                      Moles in assemblage
Phase               SI  log IAP  log K(T, P)   Initial       Final       Delta

Greigite          5.43     3.79     -1.63    0.000e+00           0   0.000e+00
Mackinawite     -11.48    -6.63      4.85    0.000e+00           0   0.000e+00
Pyrite            7.14     0.00     -7.14    9.111e-04   1.011e-03   1.000e-04
pH_fix           -6.04    -6.04      0.00
C4H8N2O3           is reactant       0.000e+00           0   0.000e+00
pe_fix            4.23     4.23      0.00
O2                 is reactant       8.200e+00   8.201e+00   4.997e-05

-----------------------------Solution composition------------------------------

Elements           Molality       Moles

C                 2.480e+00   2.400e+00
Fe                3.201e-15   3.097e-15
N                 1.240e+00   1.200e+00
S                 1.045e-03   1.011e-03

----------------------------Description of solution----------------------------

                                       pH  =   6.036      Charge balance
                                       pe  =  -4.230      Adjusted to redox equilibrium
                        Activity of water  =   0.937
                 Ionic strength (mol/kgw)  =   1.240e+00
                       Mass of water (kg)  =   9.676e-01
                 Total alkalinity (eq/kg)  =   1.240e+00
                       Total CO2 (mol/kg)  =   2.480e+00
                         Temperature (°C)  =  25.00
                  Electrical balance (eq)  =  -1.302e-11
 Percent error, 100*(Cat-|An|)/(Cat+|An|)  =  -0.00
                               Iterations  =  18
                                  Total H  = 1.158124e+02
                                  Total O  = 6.090521e+01

----------------------------Distribution of species----------------------------

                                               Log       Log       Log    mole V
   Species          Molality    Activity  Molality  Activity     Gamma    cm³/mol

   H+              1.249e-06   9.208e-07    -5.903    -6.036    -0.132      0.00
   OH-             1.813e-08   1.022e-08    -7.741    -7.991    -0.249     (0) 
   H2O             5.551e+01   9.367e-01     1.744    -0.028     0.000     18.07
C(4)          2.480e+00
   H2CO3           1.241e+00   1.651e+00     0.094     0.218     0.124     (0) 
   HCO3-           1.239e+00   7.985e-01     0.093    -0.098    -0.191     (0) 
   CO3-2           2.355e-04   4.059e-05    -3.628    -4.392    -0.764     (0) 
Fe(2)         3.201e-15
   Fe(HS)2         2.201e-15   2.929e-15   -14.657   -14.533     0.124     (0) 
   Fe+2            9.489e-16   1.846e-16   -15.023   -15.734    -0.711     (0) 
   Fe(HS)3-        5.104e-17   4.255e-17   -16.292   -16.371    -0.079     (0) 
   FeOH+           9.426e-20   5.937e-20   -19.026   -19.226    -0.201     (0) 
   Fe(OH)2         3.863e-25   5.140e-25   -24.413   -24.289     0.124     (0) 
   Fe(OH)3-        3.085e-29   1.943e-29   -28.511   -28.712    -0.201     (0) 
   FeSO4           1.684e-31   2.240e-31   -30.774   -30.650     0.124     (0) 
Fe(3)         3.529e-27
   Fe(OH)2+        3.466e-27   2.234e-27   -26.460   -26.651    -0.191     (0) 
   FeOH+2          4.213e-29   6.631e-30   -28.375   -29.178    -0.803     (0) 
   Fe(OH)3         2.006e-29   2.669e-29   -28.698   -28.574     0.124     (0) 
   Fe(OH)4-        4.215e-31   2.716e-31   -30.375   -30.566    -0.191     (0) 
   Fe+3            1.567e-32   1.010e-33   -31.805   -32.996    -1.191     (0) 
   FeSO4+          0.000e+00   0.000e+00   -46.041   -46.242    -0.201     (0) 
   Fe2(OH)2+4      0.000e+00   0.000e+00   -55.663   -56.927    -1.263     (0) 
   Fe(SO4)2-       0.000e+00   0.000e+00   -61.829   -61.908    -0.079     (0) 
   Fe3(OH)4+5      0.000e+00   0.000e+00   -79.284   -81.258    -1.974     (0) 
H(0)          2.602e-07
   H2              1.301e-07   1.731e-07    -6.886    -6.762     0.124     (0) 
N(-3)         1.240e+00
   NH4+            1.240e+00   6.262e-01     0.093    -0.203    -0.297     (0) 
   NH3             2.913e-04   3.876e-04    -3.536    -3.412     0.124     (0) 
   NH4SO4-         8.742e-17   5.506e-17   -16.058   -16.259    -0.201     (0) 
N(0)          4.009e-09
   N2              2.005e-09   2.667e-09    -8.698    -8.574     0.124     (0) 
N(3)          0.000e+00
   NO2-            0.000e+00   0.000e+00   -67.782   -67.860    -0.079     (0) 
N(5)          0.000e+00
   NO3-            0.000e+00   0.000e+00   -92.576   -92.847    -0.271     (0) 
O(0)          0.000e+00
   O2              0.000e+00   0.000e+00   -79.038   -78.913     0.124     (0) 
S(-2)         1.045e-03
   H2S             8.074e-04   1.074e-03    -3.093    -2.969     0.124     (0) 
   HS-             2.368e-04   1.334e-04    -3.626    -3.875    -0.249     (0) 
   S5-2            7.620e-08   3.682e-08    -7.118    -7.434    -0.316     (0) 
   S4-2            4.446e-08   2.148e-08    -7.352    -7.668    -0.316     (0) 
   S6-2            3.944e-08   1.906e-08    -7.404    -7.720    -0.316     (0) 
   S-2             1.112e-10   1.750e-11    -9.954   -10.757    -0.803     (0) 
   S3-2            1.567e-11   7.570e-12   -10.805   -11.121    -0.316     (0) 
   S2-2            8.891e-13   4.296e-13   -12.051   -12.367    -0.316     (0) 
   Fe(HS)2         2.201e-15   2.929e-15   -14.657   -14.533     0.124     (0) 
   Fe(HS)3-        5.104e-17   4.255e-17   -16.292   -16.371    -0.079     (0) 
S(6)          1.514e-16
   NH4SO4-         8.742e-17   5.506e-17   -16.058   -16.259    -0.201     (0) 
   SO4-2           6.398e-17   6.825e-18   -16.194   -17.166    -0.972     (0) 
   HSO4-           9.997e-22   6.099e-22   -21.000   -21.215    -0.215     (0) 
   FeSO4           1.684e-31   2.240e-31   -30.774   -30.650     0.124     (0) 
   FeSO4+          0.000e+00   0.000e+00   -46.041   -46.242    -0.201     (0) 
   Fe(SO4)2-       0.000e+00   0.000e+00   -61.829   -61.908    -0.079     (0) 

------------------------------Saturation indices-------------------------------

  Phase               SI** log IAP   log K(298 K,   1 atm)

  CH4(g)            9.28    -30.82  -40.10  CH4
  CO2(g)            1.73    -16.43  -18.16  CO2
  Fe2(SO4)3      -121.07   -117.49    3.58  Fe2(SO4)3
  Fe3(OH)8        -53.89    -33.67   20.22  Fe3(OH)8
  Ferrihydrite    -19.86    -14.97    4.89  Fe(OH)3
  FeS(ppt)         -9.66    -13.57   -3.92  FeS
  Goethite        -15.45    -14.95    0.50  FeOOH
  Greigite        -29.02    -43.87  -14.85  Fe3S4
  Hematite        -25.85    -29.86   -4.01  Fe2O3
  Jarosite-H      -91.24   -103.34  -12.10  (H3O)Fe3(SO4)2(OH)6
  Lepidocrocite   -16.32    -14.95    1.37  FeOOH
  Mackinawite     -11.48     -6.63    4.85  FeS
  Maghemite       -36.25    -29.86    6.39  Fe2O3
  Magnetite       -37.29    -33.55    3.74  Fe3O4
  Melanterite     -30.63    -33.10   -2.47  FeSO4:7H2O
  N2(g)            -5.40     -8.57   -3.18  N2
  NO2(g)          -92.19    -83.20    8.99  NO2
  O2(g)           -75.95      7.17   83.12  O2
  pe_fix            4.23      4.23    0.00  e-
  pH_fix           -6.04     -6.04    0.00  H+
  Pyrite            0.00    -16.72  -16.72  FeS2
  Pyrrhotite       -9.85    -13.57   -3.72  FeS
  Siderite         -9.58    -20.13  -10.55  FeCO3
  Sulfur           -4.19     -6.30   -2.11  S

**For a gas, SI = log10(fugacity). Fugacity = pressure * phi / 1 atm.
  For ideal gases, phi = 1.



Reaction step 90.

WARNING: pH_fix, Pure phase with add formula has not converged.
SI may be a local minimum. Residual: 2.521359e+01

Using solution 1. PURE WATER
Using pure phase assemblage 1.
Using reaction 1.

Reaction 1.

 1.000e-04 moles of the following reaction have been added:

                Relative
Reactant            moles

Fe1S3                1.00000

                Relative
Element             moles
Fe                   1.00000
S                    3.00000

-------------------------------Phase assemblage--------------------------------

                                                      Moles in assemblage
Phase               SI  log IAP  log K(T, P)   Initial       Final       Delta

Greigite          3.26     1.63     -1.63    0.000e+00           0   0.000e+00
Mackinawite      -9.32    -4.47      4.85    0.000e+00           0   0.000e+00
Pyrite            7.14     0.00     -7.14    1.011e-03   1.111e-03   1.000e-04
pH_fix           -4.95    -4.95      0.00
C4H8N2O3           is reactant       0.000e+00   6.000e-01   6.000e-01
pe_fix            4.23     4.23      0.00
O2                 is reactant       8.201e+00   1.000e+01   1.800e+00

-----------------------------Solution composition------------------------------

Elements           Molality       Moles

C                 0.000e+00   0.000e+00
Fe                4.373e-12   4.373e-12
N                 7.243e-14   7.243e-14
S                 1.111e-03   1.111e-03

----------------------------Description of solution----------------------------

                                       pH  =   4.950      Charge balance
                                       pe  =  -4.230      Adjusted to redox equilibrium
                        Activity of water  =   1.000
                 Ionic strength (mol/kgw)  =   1.126e-05
                       Mass of water (kg)  =   1.000e+00
                 Total alkalinity (eq/kg)  =   2.176e-11
                         Temperature (°C)  =  25.00
                  Electrical balance (eq)  =  -1.294e-11
 Percent error, 100*(Cat-|An|)/(Cat+|An|)  =  -0.00
                               Iterations  =  93
                                  Total H  = 1.110124e+02
                                  Total O  = 5.550508e+01

----------------------------Distribution of species----------------------------

                                               Log       Log       Log    mole V
   Species          Molality    Activity  Molality  Activity     Gamma    cm³/mol

   H+              1.126e-05   1.122e-05    -4.948    -4.950    -0.002      0.00
   OH-             8.991e-10   8.956e-10    -9.046    -9.048    -0.002     (0) 
   H2O             5.551e+01   1.000e+00     1.744    -0.000     0.000     18.07
C(4)          0.000e+00
   H2CO3           0.000e+00   0.000e+00   -44.729   -44.729     0.000     (0) 
   HCO3-           0.000e+00   0.000e+00   -46.128   -46.130    -0.002     (0) 
   CO3-2           0.000e+00   0.000e+00   -51.503   -51.509    -0.007     (0) 
Fe(2)         4.373e-12
   Fe+2            3.938e-12   3.877e-12   -11.405   -11.412    -0.007     (0) 
   Fe(HS)2         4.346e-13   4.346e-13   -12.362   -12.362     0.000     (0) 
   Fe(HS)3-        5.328e-16   5.307e-16   -15.273   -15.275    -0.002     (0) 
   FeOH+           1.097e-16   1.093e-16   -15.960   -15.961    -0.002     (0) 
   Fe(OH)2         8.293e-23   8.293e-23   -22.081   -22.081     0.000     (0) 
   Fe(OH)3-        2.758e-28   2.747e-28   -27.559   -27.561    -0.002     (0) 
   FeSO4           8.697e-38   8.697e-38   -37.061   -37.061     0.000     (0) 
Fe(3)         3.746e-25
   Fe(OH)2+        3.618e-25   3.603e-25   -24.442   -24.443    -0.002     (0) 
   FeOH+2          1.240e-26   1.221e-26   -25.907   -25.913    -0.007     (0) 
   Fe(OH)3         3.774e-28   3.774e-28   -27.423   -27.423     0.000     (0) 
   Fe+3            2.197e-29   2.121e-29   -28.658   -28.674    -0.015     (0) 
   Fe(OH)4-        3.378e-31   3.365e-31   -30.471   -30.473    -0.002     (0) 
   Fe2(OH)2+4      0.000e+00   0.000e+00   -50.370   -50.397    -0.027     (0) 
   FeSO4+          0.000e+00   0.000e+00   -52.651   -52.653    -0.002     (0) 
   Fe3(OH)4+5      0.000e+00   0.000e+00   -72.477   -72.520    -0.043     (0) 
   Fe(SO4)2-       0.000e+00   0.000e+00   -79.050   -79.052    -0.002     (0) 
H(0)          5.138e-05
   H2              2.569e-05   2.569e-05    -4.590    -4.590     0.000     (0) 
N(-3)         7.243e-14
   NH4+            7.242e-14   7.214e-14   -13.140   -13.142    -0.002     (0) 
   NH3             3.665e-18   3.665e-18   -17.436   -17.436     0.000     (0) 
   NH4SO4-         1.177e-40   1.172e-40   -39.929   -39.931    -0.002     (0) 
N(0)          0.000e+00
   N2              0.000e+00   0.000e+00   -43.137   -43.137     0.000     (0) 
N(3)          0.000e+00
   NO2-            0.000e+00   0.000e+00   -89.426   -89.428    -0.002     (0) 
N(5)          0.000e+00
   NO3-            0.000e+00   0.000e+00  -116.556  -116.558    -0.002     (0) 
O(0)          0.000e+00
   O2              0.000e+00   0.000e+00   -83.200   -83.200     0.000     (0) 
S(-2)         1.111e-03
   H2S             1.100e-03   1.100e-03    -2.959    -2.959     0.000     (0) 
   HS-             1.126e-05   1.121e-05    -4.949    -4.950    -0.002     (0) 
   S5-2            2.581e-10   2.540e-10    -9.588    -9.595    -0.007     (0) 
   S4-2            1.506e-10   1.482e-10    -9.822    -9.829    -0.007     (0) 
   S6-2            1.336e-10   1.315e-10    -9.874    -9.881    -0.007     (0) 
   Fe(HS)2         4.346e-13   4.346e-13   -12.362   -12.362     0.000     (0) 
   S-2             1.227e-13   1.208e-13   -12.911   -12.918    -0.007     (0) 
   S3-2            5.306e-14   5.223e-14   -13.275   -13.282    -0.007     (0) 
   S2-2            3.011e-15   2.964e-15   -14.521   -14.528    -0.007     (0) 
   Fe(HS)3-        5.328e-16   5.307e-16   -15.273   -15.275    -0.002     (0) 
S(6)          1.283e-28
   SO4-2           1.281e-28   1.261e-28   -27.892   -27.899    -0.007     (0) 
   HSO4-           1.378e-31   1.373e-31   -30.861   -30.862    -0.002     (0) 
   FeSO4           8.697e-38   8.697e-38   -37.061   -37.061     0.000     (0) 
   NH4SO4-         1.177e-40   1.172e-40   -39.929   -39.931    -0.002     (0) 
   FeSO4+          0.000e+00   0.000e+00   -52.651   -52.653    -0.002     (0) 
   Fe(SO4)2-       0.000e+00   0.000e+00   -79.050   -79.052    -0.002     (0) 

------------------------------Saturation indices-------------------------------

  Phase               SI** log IAP   log K(298 K,   1 atm)

  CH4(g)          -27.07    -67.17  -40.10  CH4
  CO2(g)          -43.25    -61.41  -18.16  CO2
  Fe2(SO4)3      -144.62   -141.04    3.58  Fe2(SO4)3
  Fe3(OH)8        -49.38    -29.16   20.22  Fe3(OH)8
  Ferrihydrite    -18.71    -13.82    4.89  Fe(OH)3
  FeS(ppt)         -7.50    -11.41   -3.92  FeS
  Goethite        -14.32    -13.82    0.50  FeOOH
  Greigite        -24.70    -39.54  -14.85  Fe3S4
  Hematite        -23.64    -27.65   -4.01  Fe2O3
  Jarosite-H     -104.97   -117.07  -12.10  (H3O)Fe3(SO4)2(OH)6
  Lepidocrocite   -15.19    -13.82    1.37  FeOOH
  Mackinawite      -9.32     -4.47    4.85  FeS
  Maghemite       -34.03    -27.65    6.39  Fe2O3
  Magnetite       -32.89    -29.16    3.74  Fe3O4
  Melanterite     -36.84    -39.31   -2.47  FeSO4:7H2O
  N2(g)           -39.96    -43.14   -3.18  N2
  NO2(g)         -113.75   -104.77    8.99  NO2
  O2(g)           -80.24      2.88   83.12  O2
  pe_fix            4.23      4.23    0.00  e-
  pH_fix           -4.95     -4.95    0.00  H+
  Pyrite           -0.00    -16.72  -16.72  FeS2
  Pyrrhotite       -7.69    -11.41   -3.72  FeS
  Siderite        -52.37    -62.92  -10.55  FeCO3
  Sulfur           -6.35     -8.46   -2.11  S

**For a gas, SI = log10(fugacity). Fugacity = pressure * phi / 1 atm.
  For ideal gases, phi = 1.

[dlparkhurst]

Although you entered the formula of glyclyglycine, without definitions of SOLUTION_MASTER_SPECIES and SOLUTION_SPECIES, molecular glyclyglycine will not form, and it will not buffer the solution the way that you think. Even if you find a log K for C4H8N2O3, it will not be thermodynamically stable relative to  other C,N,H, and O aqueous species. If you wanted a buffer that does not decompose, you will need to define something like the following:

Code: [Select]

SOLUTION_MASTER_SPECIES
Glygly HGlygly 2 1 1
SOLUTION_SPECIES
HGlygly = HGlygly
log_k 0
HGlygly = Glygly- + H+
-log_k -3.133
END

That said, if you are simply doing a numerical experiment, you do not need to use realistic buffers, you can just use HCl or NaOH to fix the pH. Note there is a trick in a footnote of example 8 of the manual that allows acid or base to be added as needed to fix the pH.

Be careful when you fix both pe and pH for your experiments. It is possible to step into unreasonable redox environments, particularly, you can step outside the stability range for water.

[jt2020]

Thanks for getting back to me so quickly.

Yes I was worried about the complexation with other C, N species. In the long term, I hope to test out the models on natural systems so I think it would be best using realistic buffers but will try the trick in example 8 first.

Thanks again,

Jack

[jt2020]

Hello Dr Parkhurst.

I feel I am now very close to getting this to work. Thanks for the previous suggestions. By adding Glygly as a phase it seems to be working well, however, I had one or two problems I hope you can help with.

It seems using Glyclyglycinate as a buffer works well for Fe dominated ratios (Fe3 S1), however, when inverting the ratio to S dominated, the model doesnt complete due to a numerical failure. At Fe1S2 the model runs fully, however, at a slightly more S dominated ratio of Fe1S2.1 the model fails. I've tried changing the buffer logK and amount but this seems to have no effect. I've looked at the output file and I assume the problem arises from the amount of Fe relative to S at the higher increments?  Is this because at the highest amounts of S the solution wants to be more acidic, and there are too many protons so the buffer becomes incompatible?

I also followed your suggestion and used example 8. I got a full model without numerical errors, however, the pH would not fix at 6 and fluctuates from 7 to 6 throughout the model. I wonder if you could help me.

I have copied the main areas of interest from the two models. I wanted to attach is but for some reason it says the attachments "are not writable".

Thanks again and sorry for the long post.

Jack

Here is the input file for the Fe1S2.1 model

TITLE React K-feldspar
SOLUTION 1  PURE WATER
            pH    6 charge
          temp   25.0
           pe     -4.23       #-250 mV
PHASES
       pH_fix
                H + = H +
                log_k = 0.0
       pe_fix
                e- = e-
                log_k 0.0
Mackinawite
    FeS + 2H+ = Fe+2 + H2S
    log_k     4.850
    delta_h        -8.901 kcal
Pyrrhotite
    FeS + 2H+ = Fe+2 + H2S
    log_k     3.74177
    delta_h        -7.977 kcal
Greigite
    Fe3S4 + H2 = 3Mackinawite(s) + H2S
    log_k     -1.633
    delta_h        28.408 kcal
Pyrite
    FeS2 + 0.6666666666H2 = 0.3333333333Greigite(s) + 0.6666666666H2S
    log_k     -7.144
    delta_h        0.907783 kcal
HGlygly
     HGlygly = Glygly- + H+
     -log_k -3.133
END
INCREMENTAL_REACTIONS
REACTION
Fe1S2.1  1
10*1e-12 10*1e-11 10*1e-10 10*1e-9 10*1e-8 10*1e-7 10*1e-6 10*1e-5 10*1e-4  10*1e-3  10*1e-2 mol
USE solution 1
EQUILIBRIUM_PHASES 1
        Mackinawite     0.0     0.0
        Pyrrhotite      0.0     0.0
        Greigite        0.0     0.0
        Pyrite          0.0     0.0   
        pH_fix          -6      Glygly   10.0
        pe_fix          4.23   O2(g)     10.0


USER_GRAPH 1 Simulation 6
    -headings               Pyrite Greigite Mackinawite
    -axis_titles            "Fe and S added, Log moles" "Mineral, moles"
    -chart_title            "Iron dominant system with H2S as the oxidant. Fix both -250mV at pH 6"
    -axis_scale sy_axis     auto auto auto auto log
    -initial_solutions      false
    -connect_simulations    false
    -axis_scale y_axis      auto auto auto auto Log
    -axis_scale sy_axis       auto auto auto auto Log
    -plot_concentration_vs  x
  -start
 5 PUT(GET(1) + RXN, 1)
  20 PLOT_XY log10(GET(1)), EQUI("Pyrite") color = Green, line_w = 0, \
     symbol = Plus, symbol_size = 5, y_axis = 2
  30 PLOT_XY log10(GET(1)), EQUI("Greigite") color = Blue, line_w = 0, \
     symbol = Plus, symbol_size = 5, y_axis = 2
  40 PLOT_XY log10(GET(1)), EQUI("Mackinawite") color = Cyan, line_w = 0, \
     symbol = Plus, symbol_size = 5, y_axis = 2
  60 PLOT_XY log10(GET(1)), EQUI("Pyrrhotite") color = Red, line_w = 0, \
     symbol = Plus, symbol_size = 5, y_axis = 2
  -end
    -active                 true
SELECTED_OUTPUT
    -file                 250 Fe both fix.csv
    -saturation_indices   Mackinawite  Greigite  Pyrite Pyrrhotite
    -reaction             
    -totals               Fe S
    -equilibrium_phases   Mackinawite  Greigite  Pyrite Pyrrhotite
END

OUTPUT FILE- I've selected the point where it fails.

Reaction step 78.

Using solution 1.   PURE WATER
Using pure phase assemblage 1.   
Using reaction 1.   

Reaction 1.   

     1.000e-05 moles of the following reaction have been added:

                    Relative
   Reactant            moles

   Fe1S2.1              1.00000

                    Relative
   Element             moles
   Fe                   1.00000
   S                    2.10000

-------------------------------Phase assemblage--------------------------------

                                                      Moles in assemblage
Phase               SI  log IAP  log K(T, P)   Initial       Final       Delta

Greigite          3.24     1.60     -1.63    0.000e+00           0   0.000e+00
Mackinawite      -9.29    -4.44      4.85    0.000e+00           0   0.000e+00
Pyrite            7.14     0.00     -7.14    8.111e-05   9.111e-05   1.000e-05
Pyrrhotite       -8.18    -4.44      3.74    0.000e+00           0   0.000e+00
pH_fix           -6.00    -6.00      0.00
    Glygly             is reactant       1.000e+01   1.000e+01   1.030e-07
pe_fix            4.23     4.23      0.00
    O2                 is reactant       1.000e+01   1.000e+01   4.741e-07

-----------------------------Solution composition------------------------------

   Elements           Molality       Moles

   Fe                4.461e-12   4.461e-12
   Glygly            5.435e-08   5.435e-08
   S                 9.111e-06   9.111e-06

----------------------------Description of solution----------------------------

                                       pH  =   6.000      Charge balance
                                       pe  =  -4.230      Adjusted to redox equilibrium
                        Activity of water  =   1.000
                 Ionic strength (mol/kgw)  =   1.002e-06
                       Mass of water (kg)  =   1.000e+00
                 Total alkalinity (eq/kg)  =   1.087e-07
                         Temperature (°C)  =  25.00
                  Electrical balance (eq)  =  -6.019e-16
 Percent error, 100*(Cat-|An|)/(Cat+|An|)  =  -0.00
                               Iterations  =   5
                                  Total H  = 1.110124e+02
                                  Total O  = 5.550621e+01

----------------------------Distribution of species----------------------------

                                               Log       Log       Log    mole V
   Species          Molality    Activity  Molality  Activity     Gamma    cm³/mol

   H+              1.001e-06   1.000e-06    -5.999    -6.000    -0.001      0.00
   OH-             1.006e-08   1.005e-08    -7.997    -7.998    -0.001     (0) 
   H2O             5.551e+01   1.000e+00     1.744    -0.000     0.000     18.07
Fe(2)         4.461e-12
   Fe+2            4.456e-12   4.435e-12   -11.351   -11.353    -0.002     (0) 
   Fe(HS)2         3.454e-15   3.454e-15   -14.462   -14.462     0.000     (0) 
   FeOH+           1.404e-15   1.403e-15   -14.853   -14.853    -0.001     (0) 
   Fe(HS)3-        3.520e-19   3.516e-19   -18.453   -18.454    -0.001     (0) 
   Fe(OH)2         1.194e-20   1.194e-20   -19.923   -19.923     0.000     (0) 
   Fe(OH)3-        4.441e-25   4.435e-25   -24.353   -24.353    -0.001     (0) 
   FeSO4           2.331e-29   2.331e-29   -28.632   -28.632     0.000     (0) 
Fe(3)         5.270e-23
   Fe(OH)2+        5.193e-23   5.187e-23   -22.285   -22.285    -0.001     (0) 
   Fe(OH)3         6.094e-25   6.094e-25   -24.215   -24.215     0.000     (0) 
   FeOH+2          1.574e-25   1.566e-25   -24.803   -24.805    -0.002     (0) 
   Fe(OH)4-        6.101e-27   6.094e-27   -26.215   -26.215    -0.001     (0) 
   Fe+3            2.452e-29   2.426e-29   -28.610   -28.615    -0.005     (0) 
   FeSO4+          0.000e+00   0.000e+00   -44.224   -44.224    -0.001     (0) 
   Fe2(OH)2+4      0.000e+00   0.000e+00   -48.172   -48.180    -0.008     (0) 
   Fe(SO4)2-       0.000e+00   0.000e+00   -62.253   -62.254    -0.001     (0) 
   Fe3(OH)4+5      0.000e+00   0.000e+00   -68.132   -68.145    -0.013     (0) 
Glygly        5.435e-08
   Glygly-         5.428e-08   5.422e-08    -7.265    -7.266    -0.001     (0) 
   HGlygly         7.364e-11   7.364e-11   -10.133   -10.133     0.000     (0) 
H(0)          4.083e-07
   H2              2.042e-07   2.042e-07    -6.690    -6.690     0.000     (0) 
O(0)          0.000e+00
   O2              0.000e+00   0.000e+00   -79.000   -79.000     0.000     (0) 
S(-2)         9.111e-06
   H2S             8.173e-06   8.173e-06    -5.088    -5.088     0.000     (0) 
   HS-             9.358e-07   9.347e-07    -6.029    -6.029    -0.001     (0) 
   S5-2            2.386e-10   2.375e-10    -9.622    -9.624    -0.002     (0) 
   S4-2            1.392e-10   1.386e-10    -9.856    -9.858    -0.002     (0) 
   S6-2            1.235e-10   1.229e-10    -9.908    -9.910    -0.002     (0) 
   S-2             1.134e-13   1.129e-13   -12.945   -12.947    -0.002     (0) 
   S3-2            4.906e-14   4.883e-14   -13.309   -13.311    -0.002     (0) 
   Fe(HS)2         3.454e-15   3.454e-15   -14.462   -14.462     0.000     (0) 
   S2-2            2.784e-15   2.771e-15   -14.555   -14.557    -0.002     (0) 
   Fe(HS)3-        3.520e-19   3.516e-19   -18.453   -18.454    -0.001     (0) 
S(6)          2.970e-20
   SO4-2           2.970e-20   2.956e-20   -19.527   -19.529    -0.002     (0) 
   HSO4-           2.872e-24   2.868e-24   -23.542   -23.542    -0.001     (0) 
   FeSO4           2.331e-29   2.331e-29   -28.632   -28.632     0.000     (0) 
   FeSO4+          0.000e+00   0.000e+00   -44.224   -44.224    -0.001     (0) 
   Fe(SO4)2-       0.000e+00   0.000e+00   -62.253   -62.254    -0.001     (0) 

------------------------------Saturation indices-------------------------------

  Phase               SI** log IAP   log K(298 K,   1 atm)

  Fe2(SO4)3      -119.40   -115.82    3.58  Fe2(SO4)3
  Fe3(OH)8        -40.81    -20.58   20.22  Fe3(OH)8
  Ferrihydrite    -15.51    -10.62    4.89  Fe(OH)3
  FeS(ppt)         -7.47    -11.38   -3.92  FeS
  Goethite        -11.12    -10.62    0.50  FeOOH
  Greigite        -24.64    -39.49  -14.85  Fe3S4
  Hematite        -17.22    -21.23   -4.01  Fe2O3
  HGlygly         -10.13    -13.27   -3.13  HGlygly
  Jarosite-H      -82.80    -94.90  -12.10  (H3O)Fe3(SO4)2(OH)6
  Lepidocrocite   -11.99    -10.62    1.37  FeOOH
  Mackinawite      -9.29     -4.44    4.85  FeS
  Maghemite       -27.62    -21.23    6.39  Fe2O3
  Magnetite       -24.32    -20.58    3.74  Fe3O4
  Melanterite     -28.41    -30.88   -2.47  FeSO4:7H2O
  O2(g)           -76.04      7.08   83.12  O2
  pe_fix            4.23      4.23    0.00  e-
  pH_fix           -6.00     -6.00    0.00  H+
  Pyrite            0.00    -16.72  -16.72  FeS2
  Pyrrhotite       -8.18     -4.44    3.74  FeS
  Sulfur           -6.38     -8.49   -2.11  S

**For a gas, SI = log10(fugacity). Fugacity = pressure * phi / 1 atm.
  For ideal gases, phi = 1.

Reaction step 79.

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 delay removal of equilibrium phases 1 ...

WARNING: Maximum iterations exceeded, 100

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:
WARNING: pH_fix, Pure phase with add formula has not converged.
    SI may be a local minimum.   Residual: -2.492053e-02

Using solution 1.   PURE WATER
Using pure phase assemblage 1.   
Using reaction 1.   

Reaction 1.   

     1.000e-05 moles of the following reaction have been added:

                    Relative
   Reactant            moles

   Fe1S2.1              1.00000

                    Relative
   Element             moles
   Fe                   1.00000
   S                    2.10000

-------------------------------Phase assemblage--------------------------------

                                                      Moles in assemblage
Phase               SI  log IAP  log K(T, P)   Initial       Final       Delta

Greigite          3.26     1.63     -1.63    0.000e+00           0   0.000e+00
Mackinawite      -9.32    -4.47      4.85    0.000e+00           0   0.000e+00
Pyrite            7.14     0.00     -7.14    9.111e-05   1.011e-04   1.000e-05
Pyrrhotite       -8.21    -4.47      3.74    0.000e+00           0   0.000e+00
pH_fix           -5.99    -5.99      0.00
    Glygly             is reactant       1.000e+01   1.000e+01   5.435e-08
pe_fix            4.23     4.23      0.00
    O2                 is reactant       1.000e+01   1.000e+01   4.916e-07

-----------------------------Solution composition------------------------------

   Elements           Molality       Moles

   Fe                3.982e-12   3.982e-12
   Glygly            0.000e+00   0.000e+00
   S                 1.011e-05   1.011e-05

----------------------------Description of solution----------------------------

                                       pH  =   5.989      Charge balance
                                       pe  =  -4.230      Adjusted to redox equilibrium
                        Activity of water  =   1.000
                 Ionic strength (mol/kgw)  =   1.027e-06
                       Mass of water (kg)  =   1.000e+00
                 Total alkalinity (eq/kg)  =   7.965e-12
                         Temperature (°C)  =  25.00
                  Electrical balance (eq)  =  -6.019e-16
 Percent error, 100*(Cat-|An|)/(Cat+|An|)  =  -0.00
                               Iterations  = 201 (2414 overall)
                                  Total H  = 1.110124e+02
                                  Total O  = 5.550621e+01

----------------------------Distribution of species----------------------------

                                               Log       Log       Log    mole V
   Species          Molality    Activity  Molality  Activity     Gamma    cm³/mol

   H+              1.026e-06   1.025e-06    -5.989    -5.989    -0.001      0.00
   OH-             9.811e-09   9.799e-09    -8.008    -8.009    -0.001     (0) 
   H2O             5.551e+01   1.000e+00     1.744    -0.000     0.000     18.07
Fe(2)         3.982e-12
   Fe+2            3.978e-12   3.959e-12   -11.400   -11.402    -0.002     (0) 
   Fe(HS)2         3.630e-15   3.630e-15   -14.440   -14.440     0.000     (0) 
   FeOH+           1.222e-15   1.221e-15   -14.913   -14.913    -0.001     (0) 
   Fe(HS)3-        4.015e-19   4.010e-19   -18.396   -18.397    -0.001     (0) 
   Fe(OH)2         1.014e-20   1.014e-20   -19.994   -19.994     0.000     (0) 
   Fe(OH)3-        3.678e-25   3.674e-25   -24.434   -24.435    -0.001     (0) 
   FeSO4           1.804e-29   1.804e-29   -28.744   -28.744     0.000     (0) 
Fe(3)         4.475e-23
   Fe(OH)2+        4.410e-23   4.405e-23   -22.356   -22.356    -0.001     (0) 
   Fe(OH)3         5.048e-25   5.048e-25   -24.297   -24.297     0.000     (0) 
   FeOH+2          1.370e-25   1.364e-25   -24.863   -24.865    -0.002     (0) 
   Fe(OH)4-        4.929e-27   4.923e-27   -26.307   -26.308    -0.001     (0) 
   Fe+3            2.189e-29   2.165e-29   -28.660   -28.664    -0.005     (0) 
   FeSO4+          0.000e+00   0.000e+00   -44.335   -44.336    -0.001     (0) 
   Fe2(OH)2+4      0.000e+00   0.000e+00   -48.292   -48.301    -0.008     (0) 
   Fe(SO4)2-       0.000e+00   0.000e+00   -62.426   -62.427    -0.001     (0) 
   Fe3(OH)4+5      0.000e+00   0.000e+00   -68.324   -68.337    -0.013     (0) 
Glygly        0.000e+00
   Glygly-         0.000e+00   0.000e+00   -41.917   -41.918    -0.001     (0) 
   HGlygly         0.000e+00   0.000e+00   -44.774   -44.774     0.000     (0) 
H(0)          4.292e-07
   H2              2.146e-07   2.146e-07    -6.668    -6.668     0.000     (0) 
O(0)          0.000e+00
   O2              0.000e+00   0.000e+00   -79.043   -79.043     0.000     (0) 
S(-2)         1.011e-05
   H2S             9.093e-06   9.093e-06    -5.041    -5.041     0.000     (0) 
   HS-             1.016e-06   1.014e-06    -5.993    -5.994    -0.001     (0) 
   S5-2            2.526e-10   2.514e-10    -9.598    -9.600    -0.002     (0) 
   S4-2            1.474e-10   1.467e-10    -9.832    -9.834    -0.002     (0) 
   S6-2            1.308e-10   1.301e-10    -9.884    -9.886    -0.002     (0) 
   S-2             1.201e-13   1.195e-13   -12.921   -12.923    -0.002     (0) 
   S3-2            5.193e-14   5.169e-14   -13.285   -13.287    -0.002     (0) 
   Fe(HS)2         3.630e-15   3.630e-15   -14.440   -14.440     0.000     (0) 
   S2-2            2.947e-15   2.933e-15   -14.531   -14.533    -0.002     (0) 
   Fe(HS)3-        4.015e-19   4.010e-19   -18.396   -18.397    -0.001     (0) 
S(6)          2.576e-20
   SO4-2           2.575e-20   2.563e-20   -19.589   -19.591    -0.002     (0) 
   HSO4-           2.553e-24   2.550e-24   -23.593   -23.593    -0.001     (0) 
   FeSO4           1.804e-29   1.804e-29   -28.744   -28.744     0.000     (0) 
   FeSO4+          0.000e+00   0.000e+00   -44.335   -44.336    -0.001     (0) 
   Fe(SO4)2-       0.000e+00   0.000e+00   -62.426   -62.427    -0.001     (0) 

------------------------------Saturation indices-------------------------------

  Phase               SI** log IAP   log K(298 K,   1 atm)

  Fe2(SO4)3      -119.68   -116.10    3.58  Fe2(SO4)3
  Fe3(OH)8        -41.04    -20.82   20.22  Fe3(OH)8
  Ferrihydrite    -15.59    -10.70    4.89  Fe(OH)3
  FeS(ppt)         -7.49    -11.41   -3.92  FeS
  Goethite        -11.20    -10.70    0.50  FeOOH
  Greigite        -24.69    -39.54  -14.85  Fe3S4
  Hematite        -17.39    -21.39   -4.01  Fe2O3
  HGlygly         -44.77    -47.91   -3.13  HGlygly
  Jarosite-H      -83.13    -95.23  -12.10  (H3O)Fe3(SO4)2(OH)6
  Lepidocrocite   -12.07    -10.70    1.37  FeOOH
  Mackinawite      -9.32     -4.47    4.85  FeS
  Maghemite       -27.78    -21.39    6.39  Fe2O3
  Magnetite       -24.55    -20.82    3.74  Fe3O4
  Melanterite     -28.52    -30.99   -2.47  FeSO4:7H2O
  O2(g)           -76.08      7.04   83.12  O2
  pe_fix            4.23      4.23    0.00  e-
  pH_fix           -5.99     -5.99    0.00  H+
  Pyrite            0.00    -16.72  -16.72  FeS2
  Pyrrhotite       -8.21     -4.47    3.74  FeS
  Sulfur           -6.35     -8.46   -2.11  S

**For a gas, SI = log10(fugacity). Fugacity = pressure * phi / 1 atm.
  For ideal gases, phi = 1.

ERROR: Numerical method failed on all combinations of convergence parameters, cell/soln/mix 1
-------------------------------
End of Run after 1.188 Seconds.
-------------------------------



HERE IS THE SECOND MODEL WHERE IM USING THE LOOP FROM EXAMPLE 8. I've copied the input file only as this post is getting very long. I can post the output file as well if helpful.


TITLE React K-feldspar
SOLUTION 1  PURE WATER
            pH    6 charge
          temp   25.0
           pe     -4.23       #-250 mV
PHASES
       pH_fix
                H + = H +
                log_k = 0.0
       pe_fix
                e- = e-
                log_k 0.0

Mackinawite
    FeS + 2H+ = Fe+2 + H2S
    log_k     4.850
    delta_h        -8.901 kcal
Greigite
    Fe3S4 + H2 = 3Mackinawite(s) + H2S
    log_k     -1.633
    delta_h        28.408 kcal
Pyrite
    FeS2 + 0.6666666666H2 = 0.3333333333Greigite(s) + 0.6666666666H2S
    log_k     -7.144
    delta_h        0.907783 kcal
HGlygly
    HGlygly = Glygly- + H+
    -log_k      -3.133
NaNO3
    NaNO3 = Na+ + NO3-
    log_K -20
NaCl
    NaCl = Cl- + Na+
    log_k     -20
END
INCREMENTAL_REACTIONS
REACTION
Fe1S3 1
10*1e-12 10*1e-11 10*1e-10 10*1e-9 10*1e-8 10*1e-7 10*1e-6 10*1e-5 10*1e-4  10*1e-3  10*1e-2 mol
USE solution 1
EQUILIBRIUM_PHASES 1
        Mackinawite     0.0     0.0
        Greigite        0.0     0.0
        Pyrite          0.0     0.0
        pH_fix          -6      NaOH      10.0   
        pe_fix          4.23    O2(g)     10.0

SELECTED_OUTPUT
    -file                 250 S both fix.csv
    -saturation_indices   Mackinawite  Greigite  Pyrite
    -reaction             
    -totals               Fe S C N
    -equilibrium_phases   Mackinawite  Greigite  Pyrite


USER_PUNCH
  10 FOR i = 5.9 to 6.1                               #start FOR loop
  30 a$ = a$ + "EQUILIBRIUM_PHASES 1" + EOL$
  40 a$ = a$ + "   pH_fix " + STR$(-i) + "NaOH 10" + EOL$
  45 a$ = a$ + " -force"  + EOL$
  50 a$ = a$ + "END" + EOL$                           #end loop
  60 PUNCH a$
  70 NEXT i


USER_GRAPH 1 Simulation 6
    -headings               Pyrite Greigite Mackinawite
    -axis_titles            "Fe and S added, Log moles" "Mineral, moles"
    -chart_title            "Sulfur dominant system with H2S as the oxidant. Fix both -250mV at pH 6"
    -axis_scale sy_axis     auto auto auto auto log
    -initial_solutions      false
    -connect_simulations    false
    -axis_scale y_axis      auto auto auto auto Log
    -axis_scale sy_axis       auto auto auto auto Log
    -plot_concentration_vs  x
  -start
 5 PUT(GET(1) + RXN, 1)
  20 PLOT_XY log10(GET(1)), EQUI("Pyrite") color = Green, line_w = 0, \
     symbol = Plus, symbol_size = 5, y_axis = 2
  30 PLOT_XY log10(GET(1)), EQUI("Greigite") color = Blue, line_w = 0, \
     symbol = Plus, symbol_size = 5, y_axis = 2
  40 PLOT_XY log10(GET(1)), EQUI("Mackinawite") color = Cyan, line_w = 0, \
     symbol = Plus, symbol_size = 5, y_axis = 2
  -end
    -active                 true
INCLUDE$ 250 S both fix.csv
END

[dlparkhurst]

First, I would prefer to use HGlygly in your definition of pH_fix in EQUILIBRIUM_PHASES; it is clearer that you are adding an acid.

The calculation fails at the point that you need to add a base instead of an acid. You can use the same approach as example 8 and define a phase NaGlygly.

A buffer fixes the pH by the ratio of acid to base, iny our case, the ratio of HGlygly(aq)/Glygly-. In your simulation, you can adjust this ratio up to a point by simply adding HGlygly, but eventually, you need to add Glygly- to maintain the proper ratio and the fixed pH.

Note that the dependence of Pyrite and Greigite on Mackinawite in the PHASES definition is not handled very well with PHREEQC. The calculations are done correctly and the saturation indices are correct in the SI output, but the SIs in the Phase Assemblage are not correct. Unless you are continually adjusting the Mackinawite log K, I would rewrite the Pyrite and Greigite reactions to remove Mackinawite so the output is clearer.

[jt2020]

Hello Dr Parkhurst

Ah I see. I tried the trick on example 8 and that worked fantastically. I was going to go down the root of creating new phases for more realistic buffers but as I just want to see numerically what happens at specific conditions, I think I can just use NaOH and HCl (it makes things much more simple).

Yes I noticed the SI not matching as well. I am really interested in what conditions mackinawite converts to greigite and greigite to pyrite, hence why I thought it best to include them in the reactions as solid phases. When I changed phase reactions by omitting the solid phases (see below), the phase assemblage SI's and 'Saturation indices' SI's matched.  However, despite the older models not having matching SI's I found that the SI in 'Saturation indices' between the older models and the new ones weren't that different so I'm not sure it matters too much? But if you have any thoughts Id really appreciate hearing them.

Obviously, in the original models Im adding H2S to a solid phase, whereas I assume in the newer reactions, it's more indicative of a dissolution/reprecipitation reaction? Is that correct?

Thank you

Jack


Here are the new reactions

Mackinawite
    FeS + 2H+ = Fe+2 + H2S
    log_k     4.850
    delta_h        -8.901 kcal
Greigite
    Fe3S4 + 8H+ = 2Fe+3 + Fe+2 + 4H2S
    log_k     -13.17851821
    delta_h        20.9195763 kcal
Pyrite
    FeS2 + 4H+ + Fe+2 = 2Fe+3 + 2H2S
    log_k     -28.93469201
    delta_h        20.7091683 kcal

Here are the old model reactions

Mackinawite
    FeS + 2H+ = Fe+2 + H2S
    log_k     4.850
    delta_h        -8.901 kcal
Greigite
    Fe3S4 + H2 = 3Mackinawite(s) + H2S
    log_k     -1.633
    delta_h        28.408 kcal
Pyrite
    FeS2 + 0.6666666666H2 = 0.3333333333Greigite(s) + 0.6666666666H2S
    log_k     -7.144
    delta_h        0.907783 kcal

[dlparkhurst]

I'm not sure the two sets are consistent. Here, I start with the second equation for Greigite and derive the first using reactions from llnl.dat.

Code: [Select]

Fe3S4 H2 = 3Mac H2S -1.633 28.408
3 FeS 2H+ = Fe+2 H2S 4.85 -8.901
2 H+ Fe+2 0.25O2 = Fe+3 0.5H2O 8.4899 -97.209
H2O 0.5O2 H2 -46.1066 275.588
Fe3S4 8H+ = Fe+2 2Fe+3 4H2S -16.2098 82.875

The log K (-16.21) is quite a bit different than -13.18 that you give for the Greigite in the first set. I could easily made a mistake, but you may want to recheck your log Ks because they will be critical for the relative stability of the phases.

[jt2020]

Yes you're quite right. I calculated my logK's in the new reactions using an excel spreadsheet of thermodynamic data I compiled myself, hence, why they are different from yours. Probably best to stick with just the one database!

Thanks!

Jack