Author Topic: How to add Fertilizer formular in the Phreeqc.dat Database for Inverse modelling (Read 1339 times)

evansmanu · « **on:** March 28, 2022, 11:53:34 AM »

Dear Phreeqc Users,

I come your way again with this task. I am running an inverse modelling and would want to account for the influence of fertilizer on the chemistry of the water under investigation. Here I provide you with the composition of the type of the fertilizer.

The formula to produce the ASAASE WURA Fertilizer is given as bellow

P = 22% (0.22), K = 18% (0.18), CaO = 9% (0.09), MgO = 6% (0.06), S = 7% (0.07)

All the above percentages sum up to 62% (0.62). The other nutrients which I have been informed so far is chlorine which comes in the form of KCl. There are also other micro nutrients which am not preview to. In the below equation, I attempted to find a balance for the main constituents of the fertilizer plus dissolving with H2O.

##Fertilizer
## K0.18O0.8P0.2Ca0.09Mg0.06(SO4)0.07 + 0.26 e - = 0.18 K+ + 0.2 PO4-3 + 0.09 Ca+2 + 0.06 Mg+2 + 0.07 SO4-2
## -log_k 0.0 for inverse modelling only

I included this equation in the data base and used fertilizer as part of the minerals by constraining it to dissolve only. I had only one model solution at the end of the inverse modelling which did not load Fertilizer as part of the modeled solution. I am not sure if the approach I used to form the equation above is right or wrong.

Thanks and hope to get input from you all

Evans

dlparkhurst · « **Reply #1 on:** March 28, 2022, 03:59:00 PM »

It is difficult to know how to formulate the reaction for your fertilizer. I would probably try to use some well-defined formulas for compounds. I'm not sure what would be right, but formulas like Ca3(PO4)2, KCl, CaCO3, MgCO3, Na2SO4 would make it easy to write a balanced reaction. I expect the percentages are weight percentages, so the stoichiometric coefficients of the compounds need to be adjusted so that you get the correct mass percentages of the elements.

I don't know why your fertilizer did not appear in a model. You must have had other sources for all of the elements so that fertilizer was not needed. The e- in the formula may be part of the problem because it would need an electron donor to dissolve the fertilizer. If you use the balanced compounds as suggested, you would not need e- and it may make a difference.

evansmanu · « **Reply #2 on:** March 29, 2022, 10:18:18 PM »

Thanks very much. I will work from that angle and update the forum

evansmanu · « **Reply #3 on:** May 25, 2022, 09:59:47 AM »

Dear Phreeqc users,

In my previous post, I tried to formulate a balanced stoichiometric equation for fertilizer which I used in my inverse model. I have tried to update my equation as shown below:

Fertilizer
   K0.18O0.72P0.22Ca0.09Mg0.06Na0.18(SO4)0.07Cl0.18 + 0.16H2O = 0.18K+ + 0.09Ca+2 + 0.06Mg+2 +
0.18Na+ + 0.32H+ + 0.22PO4-3 + 0.18Cl- + 0.07SO4-2
   -log_k 0.0

The inverse model produced only one model solution without fertilizer. First of all, I would like a comment on my fertilizer formula if it is correct according to the information in my earlier post. I have tried to come out with this balanced equation for the fertilizer. Can I say that in the given model solution fertilizer is not an important component in the development of groundwater? or that my stoichiometric equation is wrong? This is an area where fertilizers are widely used.

Any comments from you will be appreciated

Find below Phreeqc script

SOLUTION 1 Rain_Water Chemistry - here we evaporated the rainwater until the Cl concentration of the Northern zone is matched.
temp 25
pH 3.749
pe 3.017
redox pe
units mol/kgw
density 1
C(4) 1.147e-03 as HCO3
Ca 1.504e-03 charge
Cl 1.321e-03
K 1.559e-04
Mg 1.081e-04
Na 1.692e-04
S(6) 1.216e-03 as SO4
-water 1 # kg
SOLUTION 2 Northern Zone CaHCO3 - Water type
temp 25
pH 6.4635
pe 4
redox pe
units mg/l
density 1
Al 4.419314e-03
C(4) 1.085457e+02 as HCO3
Na     1.190286e+01
Ca     2.635481e+01
Cl 8.449450e+00 charge
K 7.085174e-01
Mg 6.433950e+00
S(6) 1.133550e+00 as SO4
Si 2.394551e+01 as SiO2
Fe     1.918470e-01
-water 1 # kg
PHASES
Phlogopite
KMg3(AlSi3)O10(OH)2 + 10 H+ = 1 Al+3 + 1 K+ + 3 Mg+2 + 3 H4SiO4
log_k 41.082
-delta_H -360.123 kJ/mol # References :92cir/nav
-analytic -1.7201279e+03 -2.6579576e-01 1.0718208e+05 6.1999929e+02 -4.7275095e+06
Plagioclase
    Na0.62Ca0.38Al1.38Si2.62O8 + 8H2O = 0.62Na+ + 0.38Ca+2 + 1.38Al(OH)4- + 2.62H4SiO4
    -log_k -18.65256
    -delta_h 20.45592 kcal
    -Vm 103.6288
Formalin
   CH2O + H2O = 4 e- + 4 H + + CO2
   -log_K = 4.8
Fertilizer
   K0.18O0.72P0.22Ca0.09Mg0.06Na0.18(SO4)0.07Cl0.18 + 0.16H2O = 0.18K+ + 0.09Ca+2 + 0.06Mg+2 + 0.18Na+ + 0.32H+ + 0.22PO4-3 + 0.18Cl- + 0.07SO4-2
   -log_k 0.0
Formalin
   CH2O + H2O = 4 e - + 4 H + + CO2
   -log_K = 4.8

INVERSE_MODELING 1
-solutions 1 2
-uncertainty 0.05 0.05
-phases
Albite dis
Phlogopite dis
Pyrite
K-mica
Fertilizer dis
Calcite
CO2(g)
Fe(OH)3(a)
K-feldspar #pre
Formalin
-balances
Al 0.1 0.1
Na 0.1 0.1
-range 1000
-tolerance 1e-10
-mineral_water true
END

Best regards
Evans

dlparkhurst · « **Reply #4 on:** May 25, 2022, 09:54:54 PM »

The fertilizer formula has both sulfur and phosphorus. However, sulfur actually decreases from initial to final solution, and phosphorus is zero in both solutions. Therefore, there is not really a need for fertilizer in an inverse model to account for the change in water chemistry. Sulfur is removed in pyrite, and you cannot add fertilizer (without an additional P sink) and end up with zero phosphorus.

So, consider whether your analyses are representative of the overall evolution of the water chemistry, and whether you should see a fertilizer signal.

evansmanu · « **Reply #5 on:** May 25, 2022, 11:36:36 PM »

Thank you for your thoughtful comments.

Yes, phosphorus concentration in groundwater was below the detection limit. This means that its concentration in the groundwater system is relatively low to be detected, and assuming that fertilizer is a source of phosphorus may not be a step in the right direction.

In general, sulphate is higher in the source solution (rainwater) affected by marine aerosols. The concentration in groundwater is far lower than that measured in rainwater. It may be that the groundwater system is a reduced environment causing sulphate to be reduced by other processes. A preliminary estimate of the pe was found to be -5, which is favourable for sulphate reduction to occur.

From here it is only fair to suppress the influence of fertilizers in my evolutionary analysis as field data may not support such a conceptual model.

evansmanu · « **Reply #6 on:** May 26, 2022, 08:49:52 PM »

Dear Phreeqc users

I've been using Phreeqc for a while now and just realized I still have more to learn. I am currently running a forward reaction pathway model, and the final step is to include a cation exchange reaction to assess the effect of cation exchange on the Na, K, Mg and Ca concentrations. The modeling process consists of the following.

1. The inverse modeling was performed using a RedModRphree package, which leverages the r programming language to perform combinatorial inverse modeling. It uses the PHREEQC as a simulator to perform combinatorial inverse modelling. The combinatorial inverse model assumes an equilibration of primary minerals and precipitation of secondary minerals. Based on the rainwater chemistry, I was able to determine plausible mineral compositions that are responsible for the chemical development of the groundwater.

2. I then insert the identified mineral compositions into a forward model to reproduce the end member solution. After going through a series of reaction steps, I got very close to the observed concentrations for most of the major ions. However, there is still a significant discrepancy between the observed and modeled concentrations for Ca, Mg, Na and K, which I think could be due to cation exchange.

Incorporating the cation exchange reaction has become a problem, and for now, I've used a few examples that I don't have much confidence in their implementation. After performing the cation exchange, I found that all other concentrations, including Al, Fe, HCO3 and others, were significantly reduced and further deviated from the observed concentrations.

My original model shows that Na and K are underestimated relative to the observed concentrations, while Ca and Mg are overestimated.

How can I best implement cation exchange in my forward model to ensure that Na and K increase while Ca or Mg decrease the cation deviations?

Please find the input script below.

TITLE Northern evolution with new conceptual model
SOLUTION   0
temp 25
pH 3.749
# pe 14 #-4.944438
redox O(0)/O(-2)
units mol/kgw
density 1
O(0) 1.0 O2(g) -0.7
C(4) 1.147e-03
Ca 1.504e-03 charge
Cl 1.321e-03 #charge
K 1.559e-04
Mg 1.081e-04
Na 1.692e-04
S(6) 1.216e-03
EQUILIBRIUM_PHASES 1
#   O2(g)      -0.7
   CO2(g)   -3.5
   Kaolinite   0.0
#   Hematite   0.0

SAVE Solution 1
END
USE Solution 1
REACTION 1
   Formalin 1
   0.001 moles
SAVE Solution 2
END
USE Solution 2
EQUILIBRIUM PHASES 2
   Albite 0 8.017e-04
   Phlogopite 0 1.669e-04
   Pyrite 0 2.166e-09
   K-feldspar 0 0
   Kaolinite 0 0
   Fe(OH)3(a) 0
   pe_Fix -4.944438 O2(g) 0.5
SAVE Solution 3
END
USE Solution 3
EXCHANGE_SPECIES
X- + K+ = KX; log_k 0.7
2X- + Ca+2 = CaX2; log_k 0.8
Na+ + X- = NaX; log_k 0
EXCHANGE 1
NaX 0.417 # exchangeable Na and K in mol
KX 0.1871
MIX 3; 1 1e-9 #
SAVE Solution 4
END

dlparkhurst · « **Reply #7 on:** May 27, 2022, 02:51:01 AM »

In your final calculation, you have USE solution 3 and MIX 3; 1 1e-9. These are conflicting definitions, and PHREEQC will use MIX 3 in preference to USE solution 3. The solution volume will be 1e-9 liters and the system will have only a tiny amount of Ca and Mg. I don't think that is what you intended.

I think you should reconsider the initial condition of your exchanger. If you start with fairly large amounts of just KX and NaX you are trying too hard to fix the system. In most cases, you would expect the exchanger to be in equilibrium with the pre-existing groundwater composition (say, EXCHANGE 1; -eq 1 ...). Then if the infilling water composition changes or there are reactions, the exchange composition will respond.

The equilibrium constants for CaX2, NaX, and KX will determine the relative affinity for the ions on the exchanger, and in turn the concentrations remaining in solution.

Finally, I do not like fixing the pe. You can force the system out of the stability field of water if you are not careful. Better to fix the O2(g) partial pressure.

evansmanu · « **Reply #8 on:** May 27, 2022, 12:22:36 PM »

I made some changes to the exchange equation and this time put it in equilibrium as shown in the script below. My intention is to equilibrate the final solution 3 with the exchanger and observes the changes in the cation concentrations (Ca2+, Na+, Mg2+ and K+. I do not see many changes in the ion concentrations after bringing the final solution 3 in equilibrium with the exchanger.

I may be wrong, however, if am right then the model may suggest that cation exchange may not have much influence on the evolution of groundwater. I also do not understand the interpretation of the exchange composition below. Other exchangers like FeX2 and AlOHX2 are also included and would like to understand their role.

The final pH was also reduced and this suggests the influence of H+ ion, how can this be controlled

-----------------------------Exchange composition------------------------------

X 1.000e-03 mol

     Equiv- Equivalent Log
   Species Moles alents Fraction Gamma

   CaX2 4.647e-04 9.293e-04 9.293e-01 -0.000
   FeX2 2.026e-05 4.053e-05 4.053e-02 -0.138
   NaX 1.095e-05 1.095e-05 1.095e-02 -0.000
   MgX2 8.856e-06 1.771e-05 1.771e-02 -0.139
   KX 1.503e-06 1.503e-06 1.503e-03 -0.000
   AlOHX2 2.118e-14 4.235e-14 4.235e-11 0.036

The objective is to use ascertain the influence of cation exchange in the system.

TITLE Northern evolution with new conceptual model
SOLUTION   0
temp 25
pH 3.749
# pe 14 #-4.944438
redox O(0)/O(-2)
units mol/kgw
density 1
O(0) 1.0 O2(g) -0.7
C(4) 1.147e-03
Ca 1.504e-03 charge
Cl 1.321e-03 #charge
K 1.559e-04
Mg 1.081e-04
Na 1.692e-04
S(6) 1.216e-03
EQUILIBRIUM_PHASES 1
   O2(g)      -0.7
   CO2(g)   -3.5
   Kaolinite   0.0
#   Hematite   0.0

SAVE Solution 1
END
USE Solution 1
REACTION 1
   Formalin 1
   0.001 moles
SAVE Solution 2
END
USE Solution 2
EQUILIBRIUM PHASES 2
   Albite 0 8.017e-04
   Phlogopite 0 1.669e-04
   Pyrite 0 2.166e-09
   K-feldspar 0 0
   Kaolinite 0 0
   Fe(OH)3(a) 0
#   pe_Fix -4.944438 O2(g) 0.5
SAVE Solution 3
END
USE Solution 3
EXCHANGE_SPECIES
X- + K+ = KX; log_k 0.7
2X- + Ca+2 = CaX2; log_k 0.8
Na+ + X- = NaX; log_k 0
EXCHANGE 1
   -equilibrate 1
   X   1.e-3
SAVE Solution 4
END

dlparkhurst · « **Reply #9 on:** May 28, 2022, 12:15:49 AM »

A few comments.

You have no Ca source or sink in your set of minerals other than calcite. At one point you had a definition of plagioclase, which seems a likely reactant if you have K and Na feldspar. The reacted water of your last post is supersaturated with calcite, which would mean calcite is likely to precipitate.

The "Formalin" definition is not wrong but seems misleading. I prefer CH2O or Organic_matter, where CH2O is a simplification of something like cellulose.

If this is a soil-zone or saturated zone weathering, soil-zone pCO2 is usually greater than atmospheric--in the range 10-3 to 10-1 depending on productivity. If you are modeling soil zone, you should probably include CO2(g) in the equilibrium phases and perhaps O2(g).

I think it is difficult to assess the effect of exchange. I think a TRANSPORT calculation may be more appropriate, where the exchange composition evolves as water percolates through the system. Still it is hard to know how to define the system. If you start with rain and a pure primary mineral assemblage and weather to secondary clays and minerals, then I think you might have a well-defined system that would exchange ions as they were created by weathering. But if you start at some intermediate weathering point, I am not sure what the correct exchange composition would be.

I might go back to the inverse modeling approach to try to get a net reaction. Before, we got hung up on the fertilizer composition. Now, you might reconsider the minerals that are weathered and the secondary minerals that form. You could start with pure water (rainwater) and run to one of your water compositions. It still may be difficult to arrive at anything definitive, but I think it would help to formalize what you think the important reactions are.

evansmanu · « **Reply #10 on:** May 30, 2022, 12:02:42 PM »

Dear David,

Thanks very much for your support. First of all, I conducted petrographic studies on samples of rock outcrops from the study area and established primary and secondary minerals pertaining to the geology of the area. With the issue of Plagioclase, the mineralogical studies revealed the presence of Na-rich (albite) and no Ca-rich (Anorthite). The reason I constrained my model to only include albite as a primary mineral which weathers into the clay minerals such as kaolinite.

Calcite, from the mineral studies, was found to be a secondary mineral as it is always found in the features and not as a bulk mineral. In the inverse modelling, it was constrained to precipitate only. These are the breakdown of the minerals from the petrographic analysis, Primary includes Albite, Phlogopite, K-mica, K-feldspar, Quartz, Pyrite and secondary, Kaolinite, Chlorite and Goethite.

You suggested a Transport model to help understand the cation exchange process. I will appreciate it if you can give me a lead or a simple script to start with.

Here is the updated script

TITLE Northern evolution with new conceptual model
SOLUTION   0
temp 25
pH 3.749
redox O(0)/O(-2)
units mol/kgw
density 1
O(0) 1.0 O2(g) -0.7
C(4) 1.147e-03
Ca 1.504e-03 charge
Cl 1.321e-03
K 1.559e-04
Mg 1.081e-04
Na 1.692e-04
S(6) 1.216e-03

PHASES
Phlogopite
KMg3(AlSi3)O10(OH)2 + 10 H+ = 1 Al+3 + 1 K+ + 3 Mg+2 + 3 H4SiO4
log_k 41.082
-delta_H -360.123 kJ/mol # References :92cir/nav
-analytic -1.7201279e+03 -2.6579576e-01 1.0718208e+05 6.1999929e+02 -4.7275095e+06
Plagioclase
Na0.62Ca0.38Al1.38Si2.62O8 + 8H2O = 0.62Na+ + 0.38Ca+2 + 1.38Al(OH)4- + 2.62H4SiO4
-log_k -18.65256
-delta_h 20.45592 kcal
-Vm 103.6288
Organic_matter
CH2O + H2O = 4 e- + 4 H + + CO2
-log_K = 4.8
Fertilizer
K0.18O0.72P0.22Ca0.09Mg0.06Na0.18(SO4)0.07Cl0.18 + 0.16H2O = 0.18K+ + 0.09Ca+2 + 0.06Mg+2 + 0.18Na+ + 0.32H+ + 0.22PO4-3 + 0.18Cl- + 0.07SO4-2
-log_k 0.0

EQUILIBRIUM_PHASES 1
   O2(g)      -0.7
   CO2(g)   -3.5
   Kaolinite   0.0
SAVE Solution 1
END

USE Solution 1
REACTION 1
   Organic_matter 1
   0.001 moles in 20 steps
SAVE Solution 2
END

USE Solution 2
EQUILIBRIUM PHASES 2
   CO2(g)   -3.0
   O2(g)      -0.7
   Albite 0 8.017e-04
   Phlogopite 0 1.669e-04
   Pyrite 0 2.166e-09
   K-feldspar 0 0
   Kaolinite 0 0
   Fe(OH)3(a) 0
   Calcite 0 0
SAVE Solution 3
END

USE Solution 3
EXCHANGE_SPECIES
X- + K+ = KX; log_k 0.7
2X- + Ca+2 = CaX2; log_k 0.8
Na+ + X- = NaX; log_k 0
EXCHANGE 1
   -equilibrate 1
   X   1.e-3
SAVE Solution 4
END

I am not able to include any attachments. Is there anything for me to do?

dlparkhurst · « **Reply #11 on:** May 30, 2022, 08:33:40 PM »

If there are no Ca sources, it is hard to make Ca+2 the dominant cation and to make secondary calcite. Right now, all of the Ca comes from your initial solution, which begs the question, Where did it come from?

So, I think there must be a source of Ca. There can be calcite in the feldspar matrix or even if it is low-Ca feldspar, it may still be the source of the Ca.

Here is a simple model that generates a CaX2 clay related to Ca-Montmorillonite, and then allows ion exchange equilibrium with major cations (I removed FeX2). It is a little forced in that I did not allow the most stable phase, kaolinite to form. Pyrite is also removed because all of it will dissolve completely in the first step given the partial pressure of O2(g).

I don't think the calculation can be made very reliable. The equilibrium constants of the aluminosilicate minerals are uncertain; moreover, equilibrium is probably not reasonable. Uncertainty of stoichiometry of the reactive minerals plus reaction kinetics, and there are just too many degrees of freedom.

I hope you have gone back and tried some inverse modeling again. With your mineralogy and analyses of water samples, seems like you should be able to convince yourself of some reasonable reactions. One of the major uses of inverse modeling is helping to establish that your preconceived notions are plausible.

Code: [Select]

PHASES
Phlogopite
    KMg3(AlSi3)O10(OH)2 + 10 H+ = 1 Al+3 + 1 K+ + 3 Mg+2 + 3 H4SiO4
    log_k        41.082
    -delta_H        -360.123    kJ/mol    # References :92cir/nav
    -analytic    -1.7201279e+03    -2.6579576e-01    1.0718208e+05    6.1999929e+02    -4.7275095e+06
Plagioclase
          Na0.62Ca0.38Al1.38Si2.62O8 + 8H2O = 0.62Na+ + 0.38Ca+2 + 1.38Al(OH)4- + 2.62H4SiO4
          -log_k  -18.65256
          -delta_h 20.45592 kcal
          -Vm 103.6288
Organic_matter
   CH2O + H2O = 4 e- + 4 H + + CO2
   -log_K = 4.8
END
EXCHANGE_SPECIES
 X- + K+ = KX; log_k 0.7
 2X- + Ca+2 = CaX2; log_k 0.8
 Na+ + X- = NaX; log_k 0
 Fe+2 + 2X- = FeX2; log_k -20
TITLE Northern evolution with new conceptual model
SOLUTION   0
    temp      25
    pH        3.749
    redox     O(0)/O(-2)
    units     mol/kgw
    density   1
    O(0)      1.0     O2(g) -0.7
    C(4)      1.147e-03
    Ca        1.504e-03 charge
    Cl        1.321e-03
    K         1.559e-04
    Mg        1.081e-04
    Na        1.692e-04
    S(6)      1.216e-03
END
SOLUTION 0-3
END
EQUILIBRIUM_PHASES 1
   #Pyrite 0 #2.166e-09
   #Fe(OH)3(a) 0 0
   K-feldspar 0 1 dis
   #Albite 0 8.017e-04
   Plagioclase  0 10 
   Phlogopite 0 #1.669e-04
   CO2(g)   -2.0
   O2(g)    -0.7
   Calcite    0 0
   Ca-Montmorillonite 0 1e-8
   #Kaolinite 0 0
   Quartz    0 0
END
USE solution 1
USE equilibrium_phases 1
SAVE solution 1
END
EXCHANGE 1
   CaX2  Ca-Montmorillonite equilibrium_phase 0.165
END 
TRANSPORT
-cells 1
-shifts 10
USER_GRAPH 1
    -headings               step CaX2 MgX2 NaX KX
    -axis_scale y_axis      auto auto auto auto log
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X STEP_NO
20 GRAPH_Y MOL("CaX2"), MOL("MgX2"), MOL("NaX"), MOL("KX"), 
  -end
    -active                 true
USER_GRAPH 2
    -headings               step Ca Mg Na K
    -axis_scale y_axis      auto auto auto auto log
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X STEP_NO
20 GRAPH_Y MOL("Ca+2"), MOL("Mg+2"), MOL("Na+"), MOL("K+"), 
  -end
    -active                 true

evansmanu · « **Reply #12 on:** May 31, 2022, 12:16:42 PM »

Hello David,

I performed a new inverse modelling using pure water instead of rainwater. The final inverse modelling resulted in a major mineral assemblages as shown bellow.

INVERSE MODELLING USING RAINWATER (Old phases used): Other minor mineral assemblages also includes: Calcite, Ca-Mont, Goethite and K-feldspar

delta
Albite -8.017e-04
K-mica 3.228e-04
Phlogopite -1.669e-04
Pyrite -2.166e-09

INVERSE MODELLING USING PURE WATER (New phases): Other minor mineral phases also includes: Calcite, Fe(OH)3, and Kaolinite

delta
Chalcedony -3.168e-04
Hematite 2.321e-08
K-mica 3.568e-06
Phlogopite -1.079e-05
Pyrite -4.652e-08

From these two results, dissolution of phlogopite, pyrite and the precipitation of k-mica, calcite and some iron oxides are common in both models. In-congruent dissolution of Phlogopite is responsible for the precipitation of K-mica whiles oxidation of pyrite also causes the formation of iron oxides such as Fe(OH)3 and Hematite.

in the coming days i will test the models with these two mineral assemblages from 1) using a pure water and 2) using rainwater.

FEW COMMENTS ON YOUR CODE
1) what is the meaning SOLUTION 0-3
2) I do not understand this "CaX2 Ca-Montmorillonite equilibrium_phase 0.165" in the EXCHANGE data block. Is the number 0.165 arbitrary?

Best

dlparkhurst · « **Reply #13 on:** May 31, 2022, 03:45:21 PM »

For ADVECTION and TRANSPORT modeling you need solutions defined for solution 0 (for forward flow direction) and each cell of the column. SOLUTION 0-3 defines solutions numbered 0, 1, 2, and 3. I ended up using only one cell, so solutions 2 and 3 were not used.

0.165 is the stoichiometry of Ca in Ca-Montmorillontite. This definition makes all of the Ca that precipitates in Ca-Montmorillonite exchangeable.

evansmanu · « **Reply #14 on:** July 27, 2022, 04:05:23 PM »

Dear all,

I am uncertain about the ion exchange component of my script. What I want to do is to equilibrate the preceding solution with an ion exchange given an exchanger amount.

For instance, after running a code that saves the resulting reaction into Solution 3 and then applies the exchange with solution 3. Must I change the number after the equilibrate to 3 as well? Find the entire script below.

.
.
.
.
USE Solution 3
EXCHANGE 1
equilibrate 1
X 0.015
SAVE Solution 4

TITLE Simulation of groundwater chemical evolution along the flow path considering mineral reactions (evaporated rainwater as the initial solution)
SOLUTION   0 Precipitation from the literature (evaporation until groundwater Cl concentration is matched)
temp 25
pH 3.749
pe 4
redox O(0)/O(-2)
redox pe
units mol/kgw
density 1
O(0) 0.01 O2(g) -0.7
C(4) 1.147e-03 as HCO3 #CO2(g) -3.5
Ca 1.504e-03 charge
Cl 1.321e-03
K 1.559e-04
Mg 1.081e-04
Na 1.692e-04
S(6) 1.216e-03
EQUILIBRIUM_PHASES 1 ##----------- equilibrated with the potential phases in the unsaturated zone (Saprolite)
   Kaolinite   0.0
   Hematite 0.0
SAVE Solution 1
END

SELECTED_OUTPUT
-file organic_matter.dat
-reaction true
-totals C(4) S(6) Fe
-molalities SO4-2 Ca+2 Na+ Cl- Mg+2 K+

USE Solution 1
REACTION 1 ## --------- decomposition of organic matter
   Organic_matter 1
    2.680e-03 moles ##0.004 moles in 100 steps #2.800e-03 moles (1.56E-05) ## --------- manual input through try and error
SAVE Solution 2
END

USE Solution 2 ## --------- Equilibration with mineral phases identified from the combinatorial inverse model in the Northern Zone.
EQUILIBRIUM PHASES 2
   CO2(g)   -1.7 ## -------- considered for soil zone CO2
   Albite 0 8.017e-02
#   Anorthite 0
#   Phlogopite 0 1.669e-02
#   Pyrite 0 #2.166e-04
   K-feldspar 0
# K-mica 0 0
#   Kaolinite 0 0
   Fe(OH)3(a) 0
Chalcedony 0
#   Calcite 0 0
   pe_Fix -4.672 O2(g) 0.5 ## ---- for redox optimization derived from the combinatorial inverse calculation
SAVE Solution 3
END
USE Solution 3
EXCHANGE 1 ## ----- Ion exchange reaction
equilibrate 1
X 0.015 ## Exchange capacity estimated (currently, there is no value from literature)
SAVE Solution 4
END

TITLE Simulation of evolution from Northern zone through Central zone ---------------------------------------------------------

USE Solution 4
EQUILIBRIUM_PHASES   3 ## Potential phase assemblages from the Combinatorial inverse calculation in the Central zone
CO2(g)   -1.75 ## considered for soil zone CO2 (estimated)
   Chalcedony 0 #2.806e-04
#   Anorthite 0
#   Pyrite 0 #1.38e-08
   K-mica 0 2.044e-04
#   K-feldspar 0 0
#   Fe(OH)3(a) 0 0
   Kaolinite 0 0
   pe_Fix -1.900623 O2(g) 0.5
SAVE Solution 5
END
USE Solution 5
EXCHANGE 2
equilibrate 2
X 0.015 #working 0.015
SAVE Solution 6
END

TITLE Simulation of evolution from Central zone to Southern zone -- here we equilibrate with the aquifer mineralogy --------------------

USE Solution 6
EQUILIBRIUM_PHASES   3 ##----- Mineral assemblages from the combinatorial inverse model in the Southern zone
#   CO2(g)   -2.0
   Calcite 0 0
#   K-feldspar 0
K-mica 0
   Phlogopite 0 10
   Chalcedony 0 2.82e6-04
   Fe(OH)3(a) 0
   pe_Fix -4.35 O2(g) 0.5
SAVE Solution 7
END
USE Solution 7
EXCHANGE 3
equilibrate 3
X 0.00009
SAVE Solution 8
END

dlparkhurst · « **Reply #15 on:** July 27, 2022, 08:33:40 PM »

I'm not sure what you are trying to do, but I will explain a bit how exchange works.

The following definition determines the exchange composition that is in equilibrium with SOLUTION 1. Solution 1 does not change, only the exchange composition is calculated.

Code: [Select]

EXCHANGE 1                       
        equilibrate 1
        X               0.015

When you use EXCHANGE 1 in a reaction as follows:

Code: [Select]

USE Solution 3
EXCHANGE 1                       ## ----- Ion exchange reaction
        equilibrate 1
        X               0.015   ## Exchange capacity estimated (currently, there is no value from literature)    
SAVE Solution 4
END

then new equilibrium solution and exchange compositions are calculated and printed to the output file. The new solution composition is saved as solution 4. However, the new exchange composition is not saved. So, EXCHANGE 1 remains the same composition that is in equilibrium with solution 1. Depending on your conceptual model, you could include SAVE exchange 4 and the exchange calculated exchange composition would be saved as EXCHANGE 4, and it would be in equilibrium with solution 4.

When you define the following:

Code: [Select]

USE Solution 5
EXCHANGE 2
        equilibrate 2
        X         0.015    #working 0.015
SAVE Solution 6
END

you define exchange 2 to be in equilibrium with solution 2 and then react it with solution 5. So, during the reaction, both the solution and the exchanger react to equilibrium. You save the solution, but not the exchanger, so EXCHANGE 2 is still in equilibrium with solution 2.

You seem to have used only one reactant (EXCHANGE, REACTION, or EQUILIBRIUM_PHASES) in each reaction definition. It is possible to use multiple reactants. For instance, a reaction with REACTION with saved solution followed by reacting the solution with EQUILIBRIUM_PHASES would produce the same solution composition as reacting with REACTION and EQUILIBRIUM_PHASES in one reaction (before the same END).

evansmanu · « **Reply #16 on:** July 28, 2022, 02:10:39 PM »

Dear David,

Thanks very much for your reply. Your write-up has provided me with much information on ion exchange. In my previous post, I just wanted to implement an exchange where my saved solution is brought into equilibrium with the exchanger and the final solution saved for the subsequent reaction. I do not intend to react with the exchanger again with the preceding solution, which has been saved.

for instance (refer to the script below), after the first Equilibrium Phase 1, I saved the solution under Solution 1. Then used the keyword USE Solution as the new input solution to react with the organic, and the final solution was saved as Solution 2
Solution 2 was then used as the input solution for the Equilibrium Phase 2, and the final reaction was saved as Solution 3.

Then I wanted to equilibrate Solution 3 with the exchanger. which I wrote like this;

USE Solution 3 ## is it correct to have this line if I just want to equilibrate the solution with the exchanger without reacting with the solution again?

EXCHANGE 1 ## ----- Ion exchange reaction
equilibrate 3 # equilibrate solution 3 with the exchanger
X 0.00009 ## Exchange capacity estimated

SAVE Solution 4
END

This is how I implemented the exchanger. But reading your writeup, I may be wrong if I use the keyword USE Solution 3 then followed by the equilibrate 3. It means that Solution 3 will then react again with the final equilibrated Solution 3 with the exchanger.

Below is the updated script considering the exchanger after each reaction

TITLE Simulation of groundwater chemical evolution along the flow path considering mineral reactions (evaporated rainwater as the initial solution)
SOLUTION   0 Precipitation from the literature (evaporation until groundwater Cl concentration is matched)
temp 25
pH 3.749
pe 4
redox O(0)/O(-2)
redox pe
units mol/kgw
density 1
O(0) 0.01 O2(g) -0.7
C(4) 1.147e-03 as HCO3 #CO2(g) -3.5
Ca 1.504e-03 charge
Cl 1.321e-03
K 1.559e-04
Mg 1.081e-04
Na 1.692e-04
S(6) 1.216e-03
EQUILIBRIUM_PHASES 1 ##----------- equilibrated with the potential phases in the unsaturated zone (Saprolite)
   Kaolinite   0.0
   Hematite 0.0
SAVE Solution 1
END

SELECTED_OUTPUT
-file organic_matter.dat
-reaction true
-totals C(4) S(6) Fe
-molalities SO4-2 Ca+2 Na+ Cl- Mg+2 K+

USE Solution 1
REACTION 1 ## --------- decomposition of organic matter
   Organic_matter 1
    2.680e-03 moles ##0.004 moles in 100 steps #2.800e-03 moles (1.56E-05) ## --------- manual input through try and error
SAVE Solution 2
END

USE Solution 2 ## --------- Equilibration with mineral phases identified from the combinatorial inverse model in the Northern Zone.
EQUILIBRIUM PHASES 2
   CO2(g)   -1.7 ## -------- considered for soil zone CO2
   Albite 0 #8.017e-02
#   Anorthite 0 0
#   Phlogopite 0 1.669e-02
#   Pyrite 0 #2.166e-04
#   K-feldspar 0
# K-mica 0 0
#   Kaolinite 0 0
   Fe(OH)3(a) 0
Chalcedony 0
#   Calcite 0 0
   pe_Fix -4.9 O2(g) 0.5 ## ---- for redox optimization derived from the combinatorial inverse calculation
SAVE Solution 3
SAVE Solution 3
END
USE Solution 3
EXCHANGE 1 ## ----- Ion exchange reaction
equilibrate 3
X 0.00009 ## Exchange capacity estimated (currently, there are no value from literature)

SAVE Solution 4
END

TITLE Simulation of evolution from Northern zone through Central zone ---------------------------------------------------------

USE Solution 4
EQUILIBRIUM_PHASES   3 ## Potential phase assemblages from the Combinatorial inverse calculation in the Central zone
CO2(g)   -1.75 ## considered for soil zone CO2 (estimated)
   Chalcedony 0 #2.806e-04
#   Anorthite 0 0
#   Pyrite 0 #1.38e-08
   K-mica 0 2.044e-04
#   K-feldspar 0 0
#   Fe(OH)3(a) 0 0
   Goethite   0 0
   Kaolinite 0 0
   pe_Fix -1.900623 O2(g) 0.5
SAVE Solution 5
END
USE Solution 5
EXCHANGE 2
equilibrate 5
X 0.015 #working 0.015
SAVE Solution 6
END

TITLE Simulation of evolution from Central zone to Southern zone -- here we equilibrate with the aquifer mineralogy --------------------

USE Solution 6
EQUILIBRIUM_PHASES   3 ##----- Mineral assemblages from the combinatorial inverse model in the Southern zone
#   CO2(g)   -1.7
   Calcite 0
#   Pyrite 0
#   K-feldspar 0
K-mica 0
   Phlogopite 0 10
   Chalcedony 0 2.82e6-04
   Fe(OH)3(a) 0
#   pe_Fix -1.7 O2(g) 0.5 #-4.35
SAVE Solution 7
END
USE Solution 7
EXCHANGE 3
equilibrate 7
X 0.00009 #0.015
SAVE Solution 8
END

dlparkhurst · « **Reply #17 on:** July 28, 2022, 03:09:26 PM »

I'll only add one comment.

Code: [Select]

USE Solution 3
EXCHANGE 1                       ## ----- Ion exchange reaction
        equilibrate 3    # equilibrate solution 3 with the exchanger
        X               0.00009   ## Exchange capacity estimated
   
SAVE Solution 4
END

This snippet causes two calculations to be made. (1) The composition of the given amount of exchanger in equilibrium with solution 3 is calculated (initial exchange calculation). (2) The exchanger is then brought together with solution 3 and reacted to equilibrium. However, because the exchanger is in equilibrium with solution 3 by the initial exchange calculation, the reaction calculation changes neither the solution nor the exchange composition. So, solution 3 and solution 4 are equivalent.

You will need to look closely at the output to see if the calculations are doing what you expect.

evansmanu · « **Reply #18 on:** July 28, 2022, 03:43:26 PM »

Dear David,

Thanks again,

From your suggestion so far, I want to have the composition of the exchanger equilibrated with Solution 3 without bringing it further into reaction with the same solution. In that regard, my code should rather be like this by removing USE Solution 3. Then Solution 4 will be the composition ofter equilibrating Solution 3 with the exchanger

I then use the Solution 4 as the input for the next reaction

However, including USE Solution 3 gives me what am looking for. The modelled concentrations match well with the observed. I was not sure if it was a good coding practice to implement the ion exchange

#USE Solution 3
EXCHANGE 1 ## ----- Ion exchange reaction
equilibrate 3 # equilibrate solution 3 with the exchanger
X 0.00009 ## Exchange capacity estimated

SAVE Solution 4
END

dlparkhurst · « **Reply #19 on:** July 28, 2022, 05:48:26 PM »

If this is your code, then there is only one calculation--finding the composition of the exchanger in equilibrium with solution 3. There is no reaction calculation, and without a reaction calculation, there is no solution to save; "SAVE solution 4" does not save a solution, so solution 4 will not exist in subsequent simulations (assuming it did not exist before this simulation).

Code: [Select]

[b]END[/b]
#USE Solution 3
EXCHANGE 1                       ## ----- Ion exchange reaction
        equilibrate 3    # equilibrate solution 3 with the exchanger
        X               0.00009   ## Exchange capacity estimated
   
SAVE Solution 4
END

evansmanu · « **Reply #20 on:** July 29, 2022, 10:01:39 AM »

Thank you David,

Finally, then how do I get the composition of the exchanger in equilibrium with the initial solution to use as an input in the subsequent batch reaction? In the below code, solution 4 will be saved after the reaction with the exchanger composition.

My concept is to use the composition of the exchanger, which is equilibrated with the solution as the input solution for the next subsequent batch reaction. So that solution 4 can be used as an input solution

USE Solution 3
EXCHANGE 1 ## ----- Ion exchange reaction
equilibrate 3 # equilibrate solution 3 with the exchanger
X 0.00009 ## Exchange capacity estimated

SAVE Solution 4
END

dlparkhurst · « **Reply #21 on:** July 29, 2022, 03:55:01 PM »

I don't know how else to explain it. The following code does two calculations which are detailed in the output file. First, the composition of the exchanger is calculated to be that which is in equilibrium with solution 3. Second, solution 3 is equilibrated with the new exchange composition, but, because the exchanger is in equilibrium with solution 3, there is no net reaction; the composition of the exchanger and the solution do not need to change to be in equilibrium. Please look at solution 3 and solution 4 in the output.

I think as you have written your script, the exchanger has no effect on the calculation sequence.

Code: [Select]

USE Solution 3
EXCHANGE 1                       ## ----- Ion exchange reaction
        equilibrate 3    # equilibrate solution 3 with the exchanger
        X               0.00009   ## Exchange capacity estimated
   
SAVE Solution 4
END

Now if you had this

Code: [Select]

USE solution 3
EXCHANGE 1
        equilibrate 1    
        X               0.00009   
SAVE solution 4

The exchanger would be defined to be in equilibrium with solution 1, so that when the reaction occurs, the solution 3 would react with the exchanger. Assuming solution 3 and solution 1 are not the same, then there would be transfers between solution 3 and the exchanger to reach a new equilibrium, and solution 4 would differ from solution 3. Again look at the output, I'm about tapped out.

evansmanu · « **Reply #22 on:** July 29, 2022, 04:52:46 PM »

Thanks alot David

evansmanu · « **Reply #23 on:** August 02, 2022, 03:53:56 PM »

Dear all,

I have a small task which I am struggling with. The issue is that I want to output the amount and the fractional amounts of the ions sorbed in my SURFACE complexation model below. The output file did not return the moles sorbed, nor are the fractional amounts.

Is there something wrong with my script? Thanks in advance for your support.

SOLUTION   0 Precipitation from the literature
temp 25
pH 3.749
pe 4
redox O(0)/O(-2)
redox pe
units mol/kgw
density 1
O(0) 0.01 O2(g) -0.7
C(4) 1.147e-03 as HCO3
Ca 1.504e-03 charge
Cl 1.321e-03
K 1.559e-04
Mg 1.081e-04
Na 1.692e-04
S(6) 1.216e-03

EQUILIBRIUM_PHASES 1 ##----------- equilibrated with the potential phases in the unsaturated zone (Saprolite)
   Kaolinite   0.0
   Hematite 0.0
SAVE Solution 1
END

USE Solution 1
REACTION 1 ## --------- decomposition of organic matter
   Organic_matter 1
    2.680e-03 moles ##0.004 moles in 100 steps #2.800e-03 moles (1.56E-05) ## --------- manual input through try and error
SAVE Solution 2
END

USE Solution 2
EQUILIBRIUM_PHASES 2
   CO2(g)   -1.7 ## -------- considered for soil zone CO2
   Albite 0 #8.017e-02
   Fe(OH)3(a) 0
   Chalcedony 0
   pe_Fix -4.9 O2(g) 0.5 ## ---- for redox optimization derived from the combinatorial inverse calculation

SURFACE 1
   -equilibrate with Solution 2
Hfo_sOH 6e-1 600. 1
Hfo_wOH 3e-4

SELECTED_OUTPUT
   -file    Northern_surf.dat
   -reset   true
USER_PUNCH
   -headings   Ca   Na   Mg   K   Fe
   -start
   10 sorbed_Ca = mol("Hfo_sOHCa+2") + mol("HfOH_wOCa+")
   20 sorbed_Na = mol("Hfo_sOHNa+") + mol("HfOH_wONa+")
   30 sorbed_Mg = mol("Hfo_sOHMg+2") + mol("HfOH_wOMg+")
   40 sorbed_K = mol("Hfo_sOHK+") + mol("HfOH_wOK+")
   50 sorbed_Fe = mol("Hfo_sOHFe+2") + mol("HfOH_wOFe+")
   60 frac_Ca = sorbed_Ca/ (sorbed_Ca + TOT("Ca"))
   70 frac_Na = sorbed_Na/ (sorbed_Na + TOT("Na"))
   80 frac_Mg = sorbed_Mg/ (sorbed_Mg + TOT("Mg"))
   90 frac_K = sorbed_K/ (sorbed_K + TOT("K"))
   95 frac_Fe = sorbed_Fe/ (sorbed_Fe + TOT("Fe"))
   -end
END

dlparkhurst · « **Reply #24 on:** August 02, 2022, 04:08:58 PM »

You need to use the PUNCH command for the items you want to appear in the selected-output file.

Code: [Select]

100 PUNCH frac_Ca, frac_Na, frac_Mg, frac_K, frac_Fe

evansmanu · « **Reply #25 on:** August 02, 2022, 07:22:47 PM »

Dear David,

Thanks for the help. Now I have the below as the sorbed ions, and I want to know how it is interpreted. This is my first time running surface complexation. Does it mean that given the initial conditions of Hfo_sOH 6e-1 600 1 and Hfo_wOH 3e-4, an amount of 3.3181e-03 moles of Ca were sorbed onto the surface? I worry about the amount of Na, Mg, K, and Fe, which looks unrealistic. Is there any better interpretation, or did I do something wrong?

My objective is to reduce the Ca concentration in my solution to mimic the reference Ca, which is in the order of e-04

Ca Na Mg K Fe
3.3181e-03 2.0000e-99 2.0000e-99 2.0000e-99 2.0000e-99

Regards,

dlparkhurst · « **Reply #26 on:** August 02, 2022, 10:47:01 PM »

You do not have the correct formulas for the sorbed species.

Here is your surface in equilibrium with an arbitrary solution containing Ca, Mg, Na, K, and Fe. Note that there are no sorption reactions (no surface species) for Na and K. Look at the exact formulas of the surface species.

Code: [Select]

Surface 1.	

Diffuse Double Layer Surface-Complexation Model

Hfo           
	  5.469e-04  Surface charge, eq
	  8.796e-02  sigma, C/m²
	  1.385e-01  psi, V
	 -5.389e+00  -F*psi/RT
	  4.566e-03  exp(-F*psi/RT)
	  6.000e+02  specific area, m²/g
	  6.000e+02  m² for   1.000e+00 g


Hfo_s         
	  6.000e-01  moles
	                                   Mole                     Log
	Species               Moles    Fraction    Molality    Molality

	Hfo_sO-           2.509e-01       0.418   2.509e-01      -0.600
	Hfo_sOFe+         2.506e-01       0.418   2.506e-01      -0.601
	Hfo_sOH           9.752e-02       0.163   9.752e-02      -1.011
	Hfo_sOH2+         8.682e-04       0.001   8.682e-04      -3.061
	Hfo_sOHCa+2       1.265e-04       0.000   1.265e-04      -3.898

Hfo_w         
	  3.000e-04  moles
	                                   Mole                     Log
	Species               Moles    Fraction    Molality    Molality

	Hfo_wO-           2.060e-04       0.687   2.060e-04      -3.686
	Hfo_wOH           8.004e-05       0.267   8.004e-05      -4.097
	Hfo_wOFeOH        1.131e-05       0.038   1.131e-05      -4.947
	Hfo_wOFe+         1.919e-06       0.006   1.919e-06      -5.717
	Hfo_wOH2+         7.126e-07       0.002   7.126e-07      -6.147
	Hfo_wOMg+         6.159e-08       0.000   6.159e-08      -7.210
	Hfo_wOCa+         3.441e-09       0.000   3.441e-09      -8.463

Another way you can calculate the total moles of an element sorbed on the Hfo surface is with the Basic function SURF.

PUNCH SURF("Ca", "Hfo") is equal to the sum of Hfo_wOCa+ and Hfo_sOHCa+2 in moles (not moles/kg water). PUNCH SURF("Ca", "Hfo")/TOT("water") would be moles per kg water.