A forward model in WRI

[1106485097@qq.com]

Hello everybody
I'm a fresh PHREEQC user

present,I have studied a paper about WATER-ROCK INTERACTION by forward model of CHILLER
general speaking，there are tow sample data，a water sample and a rock sample ，as shown in the figure below

During the modelling processes, 0.01 or 0.1grams of rocks were titrated into 1 kg groundwater until the limit given, making it
convenient to control the water/rock ratio，modeling under 150℃ and 250℃ .
The result as show below.

My question is 1) how the rock data (in concentration form) are used as  a “titrated rock ", by REACTION or EQUILIBRIUM_PHASES ?
 2) is it possible to make the same model using PHREEQC ? and how to do ? Could you give me any advice?

If somebody can help me with some advice it would be much appreciated.

 

[dlparkhurst]

Attached is a file that does similar calculations.

There are a lot of caveats. I chose similar minerals from the llnl.dat database, but the formulas and log Ks (especially at 120 C) may be very different from Chiller. The aqueous model for species and log Ks may also be very different from Chiller. Not knowing the groundwater composition, I simply generated a groundwater in equilibrium with calcite and dolomite at PCO2 10^2. Your description of adding 0.01 or 0.1 g of rock is inconsistent with a water rock ratio of 10 (which would be 100 g of rock); I added up to 1 g of rock. Finally, it is probably not reasonable to dissolve a rock at constant composition (maybe a glass?).

The file first calculates the number of moles per gram of rock for each element, then it calculates the carbonate groundwater, and finally the titration of the rock while generating three different graphs with the rock/water ratio to keep the axis similar to the plots that you showed.

[1106485097@qq.com]

Thanks Ms Parkhurst very much,you're really great and devoted （sincerely）
I have modeled after adding the groundwater data(lack for my negligence）with your remarkable model，the result is in
agreement with the previous study.

however，Some problems around me still exist.

1）You've guessed right, the rock in the paper research is basalt glass ,so it may be  reasonable to dissolve a rock at constant composition ,if not a glass,such as granite，is it meaninful to do same work？and how should i deal with？

2) what is the reasonable rock/water ratio interval ，or just according to guess?

I look forward to your answer ，Many thanks in advance.

ps：I really hope the forum can run all the time，i have learned so much ， but I haven't been able to open the donation website
past few days,may be because of my location(CNA,you kown),so if any other ways for donation？thanks again.

[dlparkhurst]

For a granite, I don't think it makes sense to use constant composition. In concept, you need a rate expression for each primary mineral in the granite, and either rates or equilibrium for appropriate secondary minerals. Example 6 in the manual is a simple example of incongruent dissolution of K-spar.

I don't think a water/rock ratio is meaningful; it is the minerals present and their rates of reaction. A start into the literature is the work of Art White, who has looked at weathering sequences and detailed weathering processes for granites.

Another approach is inverse modeling. Example 16 is an inverse model of granite weathering based on the classic work of Garrels and Christ.

[1106485097@qq.com]

 Thanks dlparkhurst a lot ,you are not only a programming expert, but also an excellent geologist.

In fact ,I have some geothermal water(also groudwater) data and rock data(as shown in the attachment） ,I want to reveal the characteristics of water rock interaction and interprete of the chemical evolution of hot spring waterin this field ,I have finished a inverse modeling，and want to set up a forward model for mutual demonstration,is it meaningful ?

My programming ability is very poor，there are other questions looking forward to your answer，thanks again.
1）If  RATE is used ,how to control the Rock/water rate ?
2) Should I consider ADVECTION or TRANSPORT ？I'm really confused.

I would appreciate it if you could give me some more advice.

[dlparkhurst]

If you have an inverse model, you should be able to start with the initial solution and add the reactants from the inverse model in REACTION. That said, provided your inverse model does indeed generate your final water, you have not really learned anything.

First you should check that the minerals that dissolve are undersaturated or saturated, and that minerals that precipitate are supersaturated or saturated. It would also be important to have more waters that corroborate the set of reactions.

Minerals near equilibrium are reacting quickly. For the others, if you know the transport time, or you know groundwater ages, you have a start on rates, but you need a lot more data to start to pin down rates.

[1106485097@qq.com]

thank dlparkhurst very much ！
In fact，I want to set up a forward model using the rock and water data by  KINETICS and RATES ,and research how rock-water carry on , finally , I want to explore  the variety under  different temperature.
 
I get some RATES data from research results of others ，but the result of my model  Seems unreasonable ，is it right ？thank for your reply .

Code: [Select]

TITLE forward model
SOLUTION 1
    temp      25
    pH        7.2
    pe        4
    redox     pe
    units     mg/l
    density   1
    Cl        11
    Alkalinity 36.94 as HCO3
    S(6)      16.6
    Ca        12.4
    K         4.5
    Na        16.06
    Mg        2.46
    F         0.86
    N         4.75
    Ba        0.06
    Si        11.32
    Sr        0.1
    Al        0.058
    -water    1 # kg
REACTION_TEMPERATURE 1
   100
EQUILIBRIUM_PHASES 1
    Albite    0 0
    Anorthite 0 0
    Ca-Montmorillonite 0 0
    Calcite   0 0
    Chlorite(14A) 0 0
    Dolomite  0 0
    Illite    0 0
    K-feldspar 0 0
    Quartz    0 0
 RATES
Albite
-start
10  DATA 11.5, 0.5, 4e-6, 0.4, 500e-6, 0.2, 13.7, 0.14, 0.15, 11.8, 0.3
20  RESTORE 10
30  READ pK_H, n_H, lim_Al, x_Al, lim_BC, x_BC, pK_H2O, z_Al, z_BC, pK_OH, o_OH
40  DATA 3500, 2000, 2500, 2000
50  RESTORE 40
60  READ e_H, e_H2O, e_OH, e_CO2
70  pk_CO2 = 13
80  n_CO2 = 0.6
100 REM Generic rate follows
110 dif_temp = 1/TK - 1/281
120 BC       = ACT("Na+") + ACT("K+") + ACT("Mg+2") + ACT("Ca+2")
130 REM rate by H+
140 pk_H     = pk_H + e_H * dif_temp
150 rate_H   = 10^-pk_H * ACT("H+")^n_H / ((1 + ACT("Al+3") / lim_Al)^x_Al * (1 + BC / lim_BC)^x_BC)
160 REM rate by hydrolysis
170 pk_H2O   = pk_H2O + e_H2O * dif_temp
180 rate_H2O = 10^-pk_H2O / ((1 + ACT("Al+3") / lim_Al)^z_Al * (1 + BC / lim_BC)^z_BC)
190 REM rate by OH-
200 pk_OH    = pk_OH + e_OH * dif_temp
210 rate_OH  = 10^-pk_OH * ACT("OH-")^o_OH
220 REM rate by CO2
230 pk_CO2   = pk_CO2 + e_CO2 * dif_temp
240 rate_CO2 = 10^-pk_CO2 * (SR("CO2(g)"))^n_CO2
250 rate     = rate_H + rate_H2O + rate_OH + rate_CO2
260 area     = PARM(1) * M0 *(M/M0)^0.67
270 rate     = PARM(2) * area * rate * (1-SR("Albite"))
280 moles    = rate * TIME
290 SAVE moles
 -end
    Anorthite
-start
 10  DATA 6.9, 1.0, 4e-6, 0.4, 500e-6, 0.25, 13.2, 0.14, 0.25, 12.0, 0.25
 20  RESTORE 10
 30  READ pK_H, n_H, lim_Al, x_Al, lim_BC, x_BC, pK_H2O, z_Al, z_BC, pK_OH, o_OH
 40  DATA 3500, 2000, 2500, 2000
 50  RESTORE 40
 60  READ e_H, e_H2O, e_OH, e_CO2
 70  pk_CO2 = 13
 80  n_CO2 = 0.6
100 dif_temp = 1/TK - 1/281
110 BC       = ACT("Na+") + ACT("K+") + ACT("Mg+2") + ACT("Ca+2")
120 pk_H     = pk_H + e_H * dif_temp
130 rate_H   = 10^-pk_H * ACT("H+")^n_H / ((1 + ACT("Al+3") / lim_Al)^x_Al * (1 + BC / lim_BC)^x_BC)
140 pk_H2O   = pk_H2O + e_H2O * dif_temp
150 rate_H2O = 10^-pk_H2O / ((1 + ACT("Al+3") / lim_Al)^z_Al * (1 + BC / lim_BC)^z_BC)
160 REM rate by OH-
170 pk_OH    = pk_OH + e_OH * dif_temp
180 rate_OH  = 10^-pk_OH * ACT("OH-")^o_OH
190 pk_CO2   = pk_CO2 + e_CO2 * dif_temp
200 rate_CO2 = 10^-pk_CO2 * (SR("CO2(g)"))^n_CO2
210 rate     = rate_H + rate_H2O + rate_OH + rate_CO2
220 rate = rate * Parm(1) * m0 * (m/m0)^0.67 * (1 - SR("Anorthite"))
240 moles    = rate * TIME
250 SAVE moles
260 PUT(rate,1)
-end
    Quartz
-start
 1  REM  Specific rate k from Rimstidt and Barnes, 1980, GCA 44,1683
 2  REM  k = 10^-13.7 mol/m2/s (25 C), Ea = 90 kJ/mol
 3  REM  sp. rate * parm(2) due to salts (Dove and Rimstidt, MSA Rev. 29, 259)
 4  REM  PARM(1) = Specific area of Quartz, m^2/mol Quartz
 5  REM  PARM(2) = salt correction: (1 + 1.5 * c_Na (mM)), < 35
10 dif_temp = 1/TK - 1/298
20 pk_w = 13.7 + 4700.4 * dif_temp
40 rate = PARM(1) * M0 *  (M/M0)^0.67 * 10^-pk_w * (1 -  SR("Quartz"))
50 moles = rate * time
60 save moles
-end
K-feldspar
 -start
1   REM Sverdrup and Warfvinge, 1995, mol m^-2 s^-1
2   REM PARM(1) = Specific area of Kspar m^2/mol Kspar
3   REM PARM(2) = Adjusts lab rate to field rate
4   REM temp corr: from A&P, p. 162. E (kJ/mol) / R / 2.303 = H in H*(1/T-1/281)
5   REM K-Feldspar parameters
10  DATA 11.7, 0.5, 4e-6, 0.4, 500e-6, 0.15, 14.5, 0.14, 0.15, 13.1, 0.3
20  RESTORE 10
30  READ pK_H, n_H, lim_Al, x_Al, lim_BC, x_BC, pK_H2O, z_Al, z_BC, pK_OH, o_OH
40  DATA 3500, 2000, 2500, 2000
50  RESTORE 40
60  READ e_H, e_H2O, e_OH, e_CO2
70  pk_CO2 = 13
80  n_CO2 = 0.6
100 REM Generic rate follows
110 dif_temp = 1/TK - 1/281
120 BC       = ACT("Na+") + ACT("K+") + ACT("Mg+2") + ACT("Ca+2")
130 REM rate by H+
140 pk_H     = pk_H + e_H * dif_temp
150 rate_H   = 10^-pk_H * ACT("H+")^n_H / ((1 + ACT("Al+3") / lim_Al)^x_Al * (1 + BC / lim_BC)^x_BC)
160 REM rate by hydrolysis
170 pk_H2O   = pk_H2O + e_H2O * dif_temp
180 rate_H2O = 10^-pk_H2O / ((1 + ACT("Al+3") / lim_Al)^z_Al * (1 + BC / lim_BC)^z_BC)
190 REM rate by OH-
200 pk_OH    = pk_OH + e_OH * dif_temp
210 rate_OH  = 10^-pk_OH * ACT("OH-")^o_OH
220 REM rate by CO2
230 pk_CO2   = pk_CO2 + e_CO2 * dif_temp
240 rate_CO2 = 10^-pk_CO2 * (SR("CO2(g)"))^n_CO2
250 rate     = rate_H + rate_H2O + rate_OH + rate_CO2
260 area     = PARM(1) * M0 *(M/M0)^0.67
270 rate     = PARM(2) * area * rate * (1-SR("K-feldspar"))
280 moles    = rate * TIME
290 SAVE moles
-end
 KINETICS 1
Albite
    -formula  NaAlSi3O8 
    -m        11.4
    -m0       11.4
    -parms    0.62  0.1
    -tol      1e-08
Anorthite
    -formula  CaAl2Si2O8 
    -m        10.8
    -m0       10.8
    -parms    0.65  0.1
    -tol      1e-08
K-feldspar
    -formula  KAlSi3O8 
    -m        29.35
    -m0       29.35
    -parms    0.65 0.1
    -tol      1e-08
Quartz
    -formula  SiO2 
    -m        85
    -m0       85
    -parms    0.14 0.1
    -tol      1e-08
-steps    0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 2 3 4 5 6 7 8 9 10  year
-step_divide 1
-runge_kutta 3
-bad_step_max 500
INCREMENTAL_REACTIONS true
USER_GRAPH 1
    -headings               TOTAL_TIME Calcite Dolomite Quartz Illite K-Feldspar Albite  Ca-Montmorillonite  Chlorite(14A)  Anorthite
    -axis_titles            "TOTAL_TIME" "Secondary mineral, Log (moles)" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
 10 GRAPH_X SIM_TIME/31536000
 20 GRAPH_Y LOG10(EQUI("Calcite"))
 30 GRAPH_Y LOG10(EQUI("Dolomite"))
 40 GRAPH_Y LOG10(EQUI("Quartz"))
 50 GRAPH_Y LOG10(EQUI("Illite"))
 60 GRAPH_Y LOG10(EQUI("K-Feldspar"))
 70 GRAPH_Y LOG10(EQUI("Albite"))
 80 GRAPH_Y LOG10(EQUI("Ca-Montmorillonite"))
 90 GRAPH_Y LOG10(EQUI("Chlorite(14A)"))
100 GRAPH_Y LOG10(EQUI("Anorthite"))
  -end
    -active                 true
USER_GRAPH 2
    -headings               SIM_TIME Na Si K Al Ca  Mg H+
    -axis_titles            "SIM_TIME" "Log Molality, Log Activity H+" ""
    -initial_solutions      false
    -connect_simulations    false
    -plot_concentration_vs  x
  -start
10 GRAPH_X SIM_TIME/31536000
20 GRAPH_Y LOG10(TOT("Na"))
30 GRAPH_Y LOG10(TOT("Si"))
40 GRAPH_Y LOG10(TOT("K"))
50 GRAPH_Y LOG10(TOT("Al"))
60 GRAPH_Y LOG10(TOT("Ca"))
70 GRAPH_Y LOG10(TOT("Mg"))
80 GRAPH_Y LA("H+")
  -end
    -active                 true

USER_GRAPH 3
    -headings               TIME Albite Quartz K-Feldspar Anorthite
    -axis_titles            "TIME" "Reactant transfer" ""
    -chart_title            "TIME-Reactant"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X SIM_TIME/31536000
20 GRAPH_Y -KIN_DELTA("Albite")
30 GRAPH_Y -KIN_DELTA("Quartz")
40 GRAPH_Y -KIN_DELTA("K-Feldspar")
50 GRAPH_Y -KIN_DELTA("Anorthite")
  -end
    -active                 true
END

[dlparkhurst]

As long as the program is functioning correctly, it is up to you at this point to determine what is right.