Guidance of evaporation simulation results

[2012HHU]

hello everyone
Please take a look at my evaporation simulation. The following is the ion concentration:
SOLUTION 1
    temp      25
    pH        7
    pe        4
    redox     pe
    units     mg/l
    density   1
    Na        162.22
    K         8.74
    Ca        42.49
    Mg        45.13
    Cl        253.32
    S(6)      123.5
    C(4)      17.81
    Alkalinity 195.22 as HCO3-
    -water    1 # kg
EQUILIBRIUM_PHASES 1
    Anhydrite 0 0
    Aragonite 0 0
    Bischofite 0 0
    Borax     0 0
    Boric_acid,s 0 0
    Calcite   0 0
    Carnallite 0 0
    Dolomite  0 0
    Gaylussite 0 0
    Glaserite 0 0
    Glauberite 0 0
    Gypsum    0 0
    Halite    0 0
    Huntite   0 0
    K2B4O7:4H2O 0 0
    KB5O8:4H2O 0 0
    Kalicinite 0 0
    Magnesite 0 0
    Mirabilite 0 0
    NaB5O8:5H2O 0 0
    NaBO2:4H2O 0 0
    Pirssonite 0 0
    Polyhalite 0 0
    Sylvite   0 0
    Syngenite 0 0
    Trona     0 0
    Bloedite  0 0
    Leonite   0 0
    Kainite   0 0
REACTION 1
    H2O(g)     -1
    53.334008  55.2424418  55.4073065  55.4333962  55.4356166  55.48457642  55.49551189 
    55.49567842  55.49717719  55.5022286  55.502627268  55.50517063 moles

The correct result is：

Evaporate   53.334008   mol of water   Gypsum+Dolomite                     precipitated out
Evaporate   55.2424418 mol of water   Magnesite+Gypsum+Dolomite    precipitated out
Evaporate   55.4073065 mol of water   Glauberite+Magnesite+Gypsum  precipitated out
Evaporate   55.4333962 mol of water   Bloedite+Glauberite                   precipitated out
Evaporate   55.4356166 mol of water   Halite +Bloedite+Glauberite       precipitated out
..........
                   
The attachment shows all the correct results, but my simulation results are not in the picture. I hope you can help me see if there is any mistake in my steps

Thank you very much for your guidance！！！

[2012HHU]

       
The picture cannot be uploaded. The correct result can only be displayed through this reply method

Moles of evaporated water（The number represents the number of moles of water evaporated，The latter minerals evaporate and precipitate）：
             
Evaporate  53.334008     mol of water      Dolomite+Gypsum      precipitated out
                 55.2424418            Magnesite+Dolomite+Gypsum      precipitated out
                 55.4073065            Glauberite+Magnesite+Gypsum     precipitated out
                 55.4333962            Bloedite+Glauberite
                 55.4356166            Halite+Bloedite+Glauberite
                 55.48457642          Polyhalite+Halite+Bloedite+Glauberite
                 55.49551189          Leonite+Polyhalite+Halite+Bloedite
                 55.49567842          Epsomite+Leonite+Polyhalite+Halite
                 55.49717719          Kainite+Epsomite+Leonite+Polyhalite+Halite
                 55.5022286            Kieserite+Kainite+Halite
                 55.502627268        Carnallite+Kieserite+Kainite+Halite
                 55.50517063          Bischofite+Carnallite+Kieserite+Kainite+Halite

[dlparkhurst]

I'm not sure what the "correct" answer is; there are a lot of possibilities.

First, I think you probably want to change your SOLUTION definition. By including both C(4) and Alkalinity, PHREEQC adjusts the pH to charge balance. If you remove C(4) and retain Alkalinity, the charge balance of the solution is close. However, if your are going to simulate extreme evaporation, I would also charge balance your solution (add "charge" to one of the analytes or use C(4) and charge).

Then there is the question of CO2(aq). If you are evaporating into air, it would make sense to include CO2(g) with a log10 partial pressure of -3.4. If you do not include this in EQUILIBRIUM_PHASES, you could generate large CO2 partial pressures in the solution (hundreds of  atmospheres).

Next question is whether minerals are allowed to redissolve. By default, the EQUILIBRIUM_PHASES definition will find the stable phase assemblage at each calculation. In other words, minerals may precipitate but then possibly dissolve when a more stable phase assemblage can be found. If minerals are conceptually removed from the system when they precipitate, you should add "pre", for precipitate only, at the end of each mineral in EQUILIBRIUM_PHASES.

PHREEQC will ultimately have trouble finding a numerical solution as more and more water is removed; with hydrated minerals, I think the numerical method has trouble distributing the remaining water between solution and minerals. I attach a file that uses a slightly different approach to removing the water. Water is removed by a first-order KINETICS reaction. The reaction takes out a fraction of the free water, which hopefully is a bit more stable than removing a greater and greater fraction of the initial water (because some water is tied up in minerals or removed by hydrolysis reactions).

Finally, the calculations depend strongly on the log Ks of all the phases. If different log Ks are selected, the set and sequence of minerals could change.

The attached file assumes equilibrium with atmospheric CO2 and uses all of the appropriate minerals from pitzer.dat for the elements in the solution, except Goergeyite (because you did not list it in any of your assemblages). I have adjusted C(4) to produce charge balance, which changes the alkalinity of the solution from your analyzed value. I did not allow minerals to redissolve. The set of phases seems to match the ones you listed, but maybe not the sequence or the amount of water removed to reach each assemblage.

Code: [Select]

SOLUTION 1
    temp      25
    pH        7
    pe        4
    redox     pe
    units     mg/l
    density   1
    C(4)      17.81 charge
    Ca        42.49
    Cl        253.32 #charge
    K         8.74
    Mg        45.13
    Na        162.22
    S(6)      123.5
    #Alkalinity 195.22 as HCO3-
    -water    1 # kg
END
USE solution 1
EQUILIBRIUM_PHASES 1
Anhydrite   0 0   pre
Aragonite   0 0   pre       
Arcanite    0 0   pre       
Artinite    0 0   pre     
Bischofite  0 0   pre     
Bloedite    0 0   pre   
Brucite     0 0   pre   
Burkeite    0 0   pre   
Calcite     0 0   pre   
Carnallite  0 0   pre   
CO2(g)      -3.4 100   
Dolomite    0 0   pre 
Epsomite    0 0   pre   
Gaylussite  0 0   pre   
Glaserite   0 0   pre   
Glauberite  0 0   pre   
#Goergeyite  0 0   pre   
Gypsum      0 0   pre   
H2O(g)      0 0   pre   
Halite      0 0   pre   
Hexahydrite 0 0   pre   
Huntite     0 0   pre   
Kainite     0 0   pre   
Kalicinite  0 0   pre   
Kieserite   0 0   pre   
Labile_S    0 0   pre   
Leonhardite 0 0   pre   
Leonite     0 0   pre   
Magnesite   0 0   pre   
MgCl2_2H2O  0 0   pre   
MgCl2_4H2O  0 0   pre   
Mirabilite  0 0   pre   
Misenite    0 0   pre   
Nahcolite   0 0   pre   
Natron      0 0   pre   
Nesquehonite  0 0   pre
Pentahydrite  0 0   pre
Pirssonite    0 0   pre
Polyhalite    0 0   pre
Portlandite   0 0   pre
Schoenite     0 0   pre
Sylvite       0 0   pre
Syngenite     0 0   pre
Thenardite    0 0   pre
Trona         0 0   pre
RATES
Evap
10 k = 0.5/86400 # 0.5 per day
20 rate = k * (TOT("water") * 1000) / GFW("H2O")
30 moles = rate * TIME
40 SAVE -moles
INCREMENTAL_REACTIONS
KINETICS
#-cvode
Evap
-formula H2O 1
-m 0
-tol 1e-10
-steps 100*0.25  days
USER_PRINT
10 t = SYS("equi", count , name$ , type$ , moles)
20 FOR i = 1 TO count
30   IF EQUI(name$(i)) > 0 THEN PRINT "xxx ", name$(i)
40 NEXT i
50 PRINT "H2O removed, moles: ", STR_F$(KIN("evap"), 15, 8)
60 END
USER_GRAPH 1
    -headings               frac Bischofite Carnallite Dolomite Epsomite Goergeyite Gypsum \
                            Halite Magnesite Bloedite Glauberite Polyhalite Leonite Kainite \
                            Kieserite H2O_evap
    -axis_titles            "Concentration factor" "Precipitate, log10(moles)" "H2O removed, moles"
    -axis_scale x_axis      auto auto auto auto log
    -axis_scale sy_axis      53 56 auto auto

    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
 10 GRAPH_X 1 / TOT("water")
 20 GRAPH_Y LOG10( EQUI("Bischofite") )
 30 GRAPH_Y LOG10( EQUI("Carnallite") )
 40 GRAPH_Y LOG10( EQUI("Dolomite") )
 50 GRAPH_Y LOG10( EQUI("Epsomite") )
 60 GRAPH_Y LOG10( EQUI("Goergeyite") )
 70 GRAPH_Y LOG10( EQUI("Gypsum") )
 80 GRAPH_Y LOG10( EQUI("Halite") )
 90 GRAPH_Y LOG10( EQUI("Magnesite") )
100 GRAPH_Y LOG10( EQUI("Bloedite") )
110 GRAPH_Y LOG10( EQUI("Glauberite") )
120 GRAPH_Y LOG10( EQUI("Polyhalite") )
130 GRAPH_Y LOG10( EQUI("Leonite") )
140 GRAPH_Y LOG10( EQUI("Kainite") )
150 GRAPH_Y LOG10( EQUI("Kieserite") )

210 GRAPH_SY KIN("Evap")
  -end
END


[2012HHU]

Thank you very much for your reply. Your steps are very reasonable.
I am not familiar with the KINETICS module and RATES module in PHREEQC. What I want to ask is how to define logks and how to choose logks for salt lake water with very high TDS?
Thank you again for your reply!!

[dlparkhurst]

Now you have an example of KINETICS and RATES that removes water from the solution over time. Rates are usually a function of solution composition (in this case the mass of water present) and are integrated over time; in this case, the length and total number of time steps are defined in the KINETICS data block.

Log Ks for aqueous species are defined with the SOLUTION_SPECIES data block. Temperature dependence of log Ks are defined either with a delta H (enthalpy of reaction) or an analytical expression which is a function of temperature. Databases have this data block, but species can be added or modified by including a data block in the input file. The pitzer.dat database has few aqueous species, and most aqueous species are simply bare ions that have log K of zero.

Log Ks for minerals and gases are defined in the PHASES data block. Again, the database contains log K; temperature dependence is delta H or analytical expression. A PHASES data block may be included in the input file to replace or add a new mineral definition.

The log Ks (both aqueous and phases) are independent of ionic strength. All of the nonideality correction is in the activity coefficients for the aqueous species.

[2012HHU]

Thank you very much for your patience!

 I would like to know which is better for simulating water evaporation in water body and studying salt precipitation sequence through RECTION module or KINETICS  and RATESmodule?

I only understand the methods described in the example of PHREEQC instruction.

In addition, for the method you provide, if I only change the ion concentration in the solution (without adding new ion types), equilibrium_ Mineral types in phases module remain unchanged.  Can the code be referenced？

[dlparkhurst]

You can use either the REACTION approach or the KINETICS approach. Both are removing water from the solution.

I tried the KINETICS approach before I realized the problems with Alkalinity and charge balance. If they converge, both approaches should give the same results. I think the KINETICS approach may work better at very large concentration factors.

You can use the approach that you feel comfortable with.

I think I included all of the minerals in pitzer.dat that are in the system defined by your SOLUTION definition (except Goergeyite) . If you add other elements, there will be more minerals that should be considered, or possibly there are other minerals that are not included in pitzer.dat.

I do not expect to be cited for the script that was posted. All of the calculations would be included in a citation of PHREEQC.

[2012HHU]

Thank you again for your patience! Thank you.