How to add brine solution in PHREEQC

[maria]

Dear all,

I would like your help on how to add brine composition in Phreeqc. We are doing some lab work and we would like to use geochemical modelling as well. We have used a brine composition of Na2SO4 2000mg/lt, KCl 500mg/lt and CaCO3 5.7mg/lt as initial solution where we have added CO2 gas to get 4% saturation. In this solution we have injected a marble block (which is mainly of calcite, dolomite). My first question is if there is a way to add this brine composition as initial solution in Phreeqc. Namely, I know that something like that

Solution #temp pH pe Na SO4 K Cl Ca CO3 is what we normally do but I was wondering if this is the only way or perhaps get these salts under equilibrium_phases and have something like this
Solution temp pH pe  Equilibrium_phases Na2SO4  KCl  CaCO3?

Can Phreeqc recognize these salts?

Thank you very much in advance.

Kind regards,
Maria

[dlparkhurst]

You could do it with EQUILIBRIUM_PHASES with some PHASES definitions, but I would prefer to use REACTION. In either case, you need to convert from mg to moles.

You could just convert to moles yourself, but here is one way to do it, that does the conversion with PHREEQC and perhaps avoids some possible calculation errors. The USER_PUNCH block writes a REACTION block that is then imported into the input file. The pure water solution (SOLUTION 1) is not exactly 1 L at 25 C, but it is close. If you are really picky, you can use MIX to adjust the volume.

The last step equilibrates with CO2(g). I'm not sure how you define 4 percent saturation. I have set the P(CO2) to 0.04 atm. You can adjust the definition if you meant something else.

Code: [Select]

SOLUTION 1
SELECTED_OUTPUT 100
-file mg2mol.pqi
USER_PUNCH 100
10 s$ = "REACTION 1" + EOL$
20 s$ = s$ + "  Na2SO4 " + STR_F$(2000*1e-3 / GFW("Na2SO4"), 15, 10) + EOL$
30 s$ = s$ + "  KCl    " + STR_F$(500*1e-3 / GFW("KCl"), 15, 10) + EOL$
40 s$ = s$ + "  CaCO3  " + STR_F$(5.7*1e-3 / GFW("CaCO3"), 15, 10) + EOL$
50 s$ = s$ + "1 mole" + EOL$
60 PUNCH s$
END
SELECTED_OUTPUT 100
-active false
END
USE solution 1
INCLUDE$ mg2mol.pqi
SAVE solution 2
END
USE solution 2
EQUILIBRIUM_PHASES
CO2(g) -1.4 10 # -1.4 = log10(0.04 atm)
SAVE solution 3
END


[maria]

Dear David,

Thank you very much for your email. I tried to use the code but I am getting an error. I am using pitzer database and because the experiment has been set up under T 220oC and P 2000psi (I converted this to atm for the phreeqc), I have used this under the solution 1.

Code: [Select]

SOLUTION 1
#temp 220
#pressure 130
pH 6.8

I am a bit puzzled.

I don't know if the next step that I will need to model with phreeqc will help but after setting up the initial brine solution with CO2, I will need to find out how long (in months) will reach equilibrium with the rock, a marble block, that they have injected in that mixture. I have the volume of brine 310 ml and the volume of co2 6286 ml. They have done XRF and they have found that the rock is mainly consists of calcite/ dolomite. But I don't have mineral composition with concentrations. I am not sure if the next step will be with Kinetics or do I need more information for phreeqc?

Once again thank you very much in advance. Your help will be valuable.

Kind regards,
Maria

[dlparkhurst]

I'm not sure why your calculation failed. This script runs for me. I did change solution 1 to have your temperature and pressure, but I let the pH adjust to charge balance. I also, adjusted the mass of water so that the volume after addition of salts is close to 1 L.

I then used MIX to set the solution volume to 0.310 L.

I'm still not sure how you want to handle the CO2. Do you have a head space of 0.62 L of pure CO2 at 1 atm? You might use a fixed volume GAS_PHASE, or you could simply specify the CO2(g) partial pressure with EQUILIBRIUM_PHASES if the partial pressure is known.

I have added Calcite and Dolomite to EQUILIBRIUM_PHASES to estimate the amounts that dissolve. It is not much, on the order of 1e-5 moles in 310 mL. At 220 C I would think the reaction is pretty fast, on the order of hours. But you have experimental results.

If you want to run the reaction kinetically, you will need to define a rate expression in RATES. phreeqc.dat has a calcite rate expression from Plummer et al., 1978, AJS 278, 179; Appelo et al., AG 13, 257, but I'm not sure for what temperature range it is fit. There are rate expressions in Palandri and Kharaka (implemented for PHREEQC by Chen ZHU at the University of Indiana). Again you will have to determine whether the temperature range extends to 220 C. You will also need the reactive surface area for each mineral. If you have experimental dissolution rates, you may want to fit your own rate parameters.

Code: [Select]

SOLUTION 1
-pressure 130
-temp 220
pH 6.8  charge
-water 0.849
SELECTED_OUTPUT 100
-file mg2mol.pqi
USER_PUNCH 100
10 s$ = "REACTION 1" + EOL$
20 s$ = s$ + "  Na2SO4 " + STR_F$(2000*1e-3 / GFW("Na2SO4"), 15, 10) + EOL$
30 s$ = s$ + "  KCl    " + STR_F$(500*1e-3 / GFW("KCl"), 15, 10) + EOL$
40 s$ = s$ + "  CaCO3  " + STR_F$(5.7*1e-3 / GFW("CaCO3"), 15, 10) + EOL$
50 s$ = s$ + "1 mole" + EOL$
60 PUNCH s$
END
SELECTED_OUTPUT 100
-active false
END
USE solution 1
INCLUDE$ mg2mol.pqi
SAVE solution 2
END
MIX
2 0.310
SAVE solution 3
END
USE solution 3
EQUILIBRIUM_PHASES
CO2(g) -1.4 10 # -1.4 = log10(0.04 atm)
Calcite 0 10
Dolomite 0 10
SAVE solution 4
END


[maria]

Dear David,

Thank you very much. But I am still getting an error. I am getting the following

Code: [Select]

ERROR: Unknown input in SOLUTION keyword
ERROR: - pressure 130
ERROR: Can't open file, mg2mol.pqi
ERROR: Calculations terminating due to input errors.

I don't have dissolution rates. They have set up the experiment with the brine solution-CO2 and the marble rock and now they would like me to calculate how long will reach equilibrium the whole system so they can open the cabinet. From reading, I though that either kinetics or reaction could calculate this time.

I hope this helps.

Kind regards,
Maria

[maria]

Dear David,

I was running the code to the Phreeqc for windows and I was getting these errors. But I thought to try it to Phreeqc Interactive and it actually run.

But I am still confused of how I can model the time duration for reaching equilibrium.

Thank you very much in advance.

Kind regards,
Maria

[dlparkhurst]

Phreeqc for Windows uses PHREEQC 2, whereas PhreeqcI uses PHREEQC 3. Some of the calculations need PHREEQC 3 features, including the updated pitzer.dat database.

Here is an example of a kinetic calcite dissolution calculation. The parameters have been set arbitrarily. I expect the rate expression only applies at less than 100 C. You have a lot of work to develop a rate expression that applies at 220 C, one for calcite and one for dolomite.

Code: [Select]

RATES
########
#Calcite
########
# Example of KINETICS data block for calcite rate,
#   in mmol/cm2/s, Plummer et al., 1978, AJS 278, 179; Appelo et al., AG 13, 257.
# KINETICS 1
# Calcite
# -tol   1e-8
# -m0    3.e-3
# -m     3.e-3
# -parms 1.67e5   0.6  # cm^2/mol calcite, exp factor
# -time  1 day

Calcite
   -start
1   REM   PARM(1) = specific surface area of calcite, cm^2/mol calcite
2   REM   PARM(2) = exponent for M/M0

10  si_cc = SI("Calcite")
20  IF (M <= 0  and si_cc < 0) THEN GOTO 200
30  k1 = 10^(0.198 - 444.0 / TK )
40  k2 = 10^(2.84 - 2177.0 /TK )
50  IF TC <= 25 THEN k3 = 10^(-5.86 - 317.0 / TK)
60  IF TC > 25 THEN k3 = 10^(-1.1 - 1737.0 / TK )
80  IF M0 > 0 THEN area = PARM(1)*M0*(M/M0)^PARM(2) ELSE area = PARM(1)*M
110 rate = area * (k1 * ACT("H+") + k2 * ACT("CO2") + k3 * ACT("H2O"))
120 rate = rate * (1 - 10^(2/3*si_cc))
130 moles = rate * 0.001 * TIME # convert from mmol to mol
200 SAVE moles
   -end

SOLUTION 1
C 1 CO2(g) -2
END
INCREMENTAL_REACTIONS
USE solution 1
KINETICS
Calcite
 -m     0.1
 -parms 1e4   0.8  # cm^2/mol calcite, exp factor
 -time  7200 in 10 steps
USER_GRAPH 1
    -headings               Time SI(Calcite)
    -axis_titles            "Hours" "SI(Calcite)" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME / 3600
20 GRAPH_Y SI("Calcite")
  -end
    -active                 true


[maria]

Dear David,

Thank you so much for the help. It is a good start for me. I did some reading and I was wondering if I get the calcite amount and the surface area, so I can use these instead, I believe that will make the rates slightly closer to the reactions that are taking place in the experiment.

I want also to ask about the tolerance and the exponent factor. Are these default or again they have to do with rate expression and perhaps they will need modification? I noticed that under the rates you have used 0.6 and under kinetics 0.8. I also noticed that in time you used 1 day. To be honest, I tried to change this to 1 year, but there wasn't any difference.

Zhang et al. 2019 have published a paper "A library of Basic scripts of reaction rates for geochemical modelling using PHREEQC" where they have produced a geothermal database, which is free to download. Perhaps this could help as well.

Thank you very much.

Kind regards,
Maria

[dlparkhurst]

The tolerance simply indicates how accurately you want to integrate the rate equation. The smaller the number, the more accurate, but also longer calculation time. 1e-8 to 1e-10 is usually adequate.

The exponential factor is used to adjust the surface area as a mineral dissolves. A value of 1 indicates the ratio of surface area (A/A0) varies directly with mineral volume (mineral mass). If the minerals dissolve as perfect cubes or perfect spheres the exponential factor would be 0.66. If M/M0 remains close to 1, it makes little difference. Reactive surface area is usually difficult to evaluate and often is more or less a fitting factor.

The time is simply the length of time that the rates are integrated. I used 2 days in 10 steps, so the program integrated over 10 0.2 hour intervals. Shorter time intervals are usually more numerically accurate, but usually, the integration to the same time will be similar regardless of the time steps used to get there.