CO2(g) partial pressure / molar concentration in gas phase

[bashuisman]

Hi all,

As a context: I'm looking at waters from carbonate reservoirs at relatively high P&T (around 500 ATM and 75C). Since there can sometimes be 3 different C(4) concentrations for which a water is in equilibrium with Calcite (unphysically low C(4) and high pH, mid range C(4) and mid range pH, unphysically high C(4) and low pH) I add CO2 in many steps to a water using the REACTION keyword, and take the most physical crossing of SI_CO2(g) with SI_Calcite=0. I then use that C(4) to generate a water based on my (room temperature) samples that is in equilibrium with Calcite. I can automate this also for larger number of water samples. However, for each of the waters I would like to know the CO2(g) concentration that would generate this type SI_CO2(g) at this P&T. However, I'm failing to be able to calculate the partial pressure of CO2(g).

If I simply speciate a water at a high P&T, including C(4) as a species, I always get as output an SI_CO2(g) with a PR_PHI("CO2(g)") of 1.

As soon as I USE that solution and equilibrate it with CO2(g) using EQUILIBRIUM_PHASES with that output SI_CO2(g) as new input, I get a much lower SI_CO2(g), and a completely different carbonate/pH/C concentration. Now it does calculate a reasonable PR_PHI.

There's 3 problems with this for me:
1) If CO2(g) is supplied as an input phase (either via GAS_PHASES or EQUILIBRIUM_PHASES) it calculates a PR_PHI, if it is not, it doesn't calculate PR_PHI and assumes it to be 1.
2) It seems that output SI_CO2(g) is the log10( CO2(g) fugacity), and the input SI_CO2(g) is the log10( partial pressure CO2(g)). This inconsistency can cause mistakes very easily.
3) These 2 facts combined make that it is near impossible for me to calculate the partial pressure of CO2(g) for a specific solution. The only thing I can do is USE my speciated solution manually and iteratively change the input SI_CO2(g) until the output SI_CO2(g) matches the output SI_CO2(g) of my speciated solution. This is fine for 1 solution, but terrible work if you want to do this for lots of water samples.

Perhaps it's not well defined what the PR_PHI of CO2(g) is if you don't know the concentration of the other elements in the gas phase, but then it may be better to assume that CO2(g) is the only species, so as to at least know if CO2(g) is under- of oversaturated at the current P&T, if water was the only phase present.

I guess that two changes in phreeqc would allow me to do this more easily:
1) Phreeqc could calculate PR properties for CO2 even if it is not an "input", then I could just read out the partial pressure of CO2(g) using USER_PUNCH.
2) As a work around if I could choose to input real fugacities instead of partial pressures for CO2(g), then I could just take all my solutions, use the output SI as input fugacity, and read off the PR_P("CO2(g)") again.

Do any of you maybe know a workaround that would not require a code change to Phreeqc? I don't really want to iteratively change SI_CO2(g), although I guess I can also script around that somehow.

Regards,
Bastiaan Huisman

[dlparkhurst]

Peng-Robinson applies considers a non-ideal GAS_PHASE, so, I think you need to use a GAS_PHASE to be able to consider the effects of temperature and pressure on the fugacities of GAS_PHASE components.

I did not entirely follow all of your post, it seems like you are trying to infer gas-phase compositions from solution compositions, but the phi of the gas phase is not well defined (because the composition of the gas phase is not known?). Perhaps adding a fixed-volume GAS_PHASE with a small volume (<= 1e-3 L?) would be a way to determine gas compositions that would be in equilibrium with your solution, but have no effect on the aqueous calculation.

If you know the GAS_PHASE composition that you want, then define it and give it a large volume, so that it is not affected by a small transfer to the aqueous solutions.

[bashuisman]

Hi David,

many thanks for your quick reply. Sorry that I wasn't more clear.

I had tried the GAS_PHASES, but found some behavior I didn't understand (all of a sudden the water was also calculated at the same pressure as the partial pressure of the CO2, which e.g. was 10% of the original pressure), so I abandoned that. I'll re-research this option, including Methane in the mix, and using a much smaller gas volume.

Basically, I indeed want to infer a gas composition from a solution composition, to see what kind of CO2 gas composition I'd need for this water to be in equilibrium with calcite at a certain P&T. It's not that I'm going to "estimate" the CO2 composition from a water, more that I want to know up-front the region of stability of calcite w.r.t. CO2 gas composition. In this case it is a hydrocarbon reservoir, so the rest of the gas would be hydrocarbons, but I guess using Methane and CO2 together should give me a close enough approximation. At a minimum I would like to know if the water is over-saturated with CO2 (it's my understanding that you can't have PR_P(CO2) larger than log10(pressure), because CO2 would spontaneously form a gas phase then).

As your manual says, PR_PHI depends on the composition of entire gas phase. Yet if I do an experiment with one EQUILIBRIUM_PHASES with SI_CO2=log10(60) and SI_Mtg(g)=log10(540), and another with only EQUILIBRIUM_PHASES : SI_CO2=log10(60) for a solution that is at pressure 600 ATM, I find an identical PR_PHI(CO2(g)) in both cases of 7.820724789488e-01 (P=600, T=60, with high_precision true). So at least for the EQUILIBRIUM_PHASES it doesn't seem to use the gas composition, but still reports a PR_PHI for each gas. I was assuming that Phreeqc would at least calculate the same PR_PHI it calculates for EQUILIBRIUM_PHASES (possibly a pure gas phase fugacity coefficient) if it didn't know the gas composition.

I will research your GAS_PHASE suggestion, and let you know my results. I think it might still be valuable to either not calculate PR_PHI for equilibrium_phases if it doesn't actually calculate the influence of other gas components (or warn the user that it's not accurate), or if that's considered a minor approximation, to also calculate this PR_PHI for gas phases for which you state an SI.

[bashuisman]

OK, I had a little go with the GAS_PHASE.
I start with a tiny fixed volume (1e-8L seems to not change the volumes of C in the solution enough to change the solution). I first use only CO2(g), left empty, and same -pressure/-temperature as I used for the speciation.
I again see the same problem I had before, the pressure in the entire solution drops to the (partial) pressure of CO2(g), even though the gas volume is so tiny that it shouldn’t change the pressure in the solution (increasing the bubble 10 fold has no influence on the pressure drop). More worryingly, the SI_CO2(g) drops by 0.2.
I then used EQUILIBRIUM_PHASES and hand-tuned the input SI_CO2(g) to get the same C in solution. I then get the same SI_CO2(g), a much higher PR_P(CO2(g)) and a much lower PR_PHI(CO2(g)). The difference is large. The GAS_PHASE solution results in a partial pressure of 245 ATM, the EQUILIBRIUM_PHASES solution results in 648 ATM, higher than the solution pressure, which is 600 ATM.
Is there anything I’m doing wrong with the GAS_PHASE command?
Here’s my input file:

Code: [Select]

DATABASE  ../databases/pitzer.dat
SELECTED_OUTPUT
   -reset false
   -totals C
   -saturation_indices CO2(g)
   -file tTP.txt
   USER_PUNCH
   -headings PRESSURE PHICO2 lPCO2
   10 PUNCH  PRESSURE
   20 PUNCH  PR_PHI("CO2(g)")
   30 PUNCH  log10(PR_P("CO2(g)"))
END

SOLUTION_SPREAD
-units mol/kgw
NUMBER     pH    temp   Alkalinity  C        Ca    Cl      Na  pressure 
                                                    charge                  
   1       3.95   60   1.75e-02    0.772    0.38  3.14    2.4 600       
END

USE SOLUTION 1
   GAS_PHASE 1
   -fixed_volume
   -pressure 600
   -temperature 60
   -volume 0.00000001
   CO2(g)
END

USE SOLUTION 1
   EQUILIBRIUM_PHASES 1
   CO2(g) 2.695425805
END(you can see I had fun iterating to find the correct SI_CO2(g) ;-)). The output was:

Code: [Select]

         C    si_CO2(g)     PRESSURE       PHICO2        lPCO2
7.7200e-01       2.1922   6.0000e+02   1.0000e+00         -inf
7.7200e-01       1.9939   2.4509e+02   4.0232e-01   2.3893e+00
7.7200e-01       2.1922   6.0000e+02   3.1389e-01   2.6954e+00

Any comments are greatly appreciated!


PS. I also tried with Mtg. However, without Mtg in the initial solution, it obviously will not calculate a concentration for Mtg(g) in a small gas bubble, since the only source of Mtg comes from the solution. However, adding Mtg to the solution doesn't seem right either, since I don't actually know the Mtg concentration, and in my simple problem it is 1-M_CO2(g), which is the one thing I'm trying to determine. I could saturate the solution with Mtg at speciation, with
Mtg 0.1 Mtg(g) 2.415 #log10(600)
But that doesn't converge, I guess because Mtg is already at 100% saturation.

[dlparkhurst]

When you do a SOLUTION calculation, you calculate the activity of CO2(aq). Using Henry's law, you can then calculate the fugacity of CO2(g) that is in equilibrium with that solution.

However, the gas phases that have that fugacity is not unique. There may be multiple gas-phase compositions that produce the same fugacity of CO2(g).

If you assume that you want a gas phase that is pure CO2(g), then the calculation with a tiny GAS_PHASE that contains only CO2 gives you the pressure and phi for that gas phase.

If you want to define an EQUILIBRIUM_PHASE that requires no transfer of CO2 to the aqueous phase, you need to use the log pressure determined from the GAS_PHASE calculation defined in the last paragraph. You noted the "feature" that the EQUILIBRIUM_PHASES target SI is the log10 pressure, whereas the Saturation Index list in the output is log10 fugacity.

The following defines such an EQUILIBRIUM_PHASE where 2.39 is the log pressure in the output from the GAS_PHASE description.

Code: [Select]

SELECTED_OUTPUT
   -reset false
   -totals C
   -saturation_indices CO2(g)
   -file tTP.txt
   USER_PUNCH
   -headings PRESSURE PHICO2 lPCO2
   10 PUNCH  PRESSURE
   20 PUNCH  PR_PHI("CO2(g)")
   30 PUNCH  log10(PR_P("CO2(g)"))
END

SOLUTION_SPREAD
-units mol/kgw
NUMBER     pH    temp   Alkalinity  C        Ca    Cl      Na  pressure 
                                                    charge                  
   1       3.95   60   1.75e-02    0.772    0.38  3.14    2.4 600       
END

USE SOLUTION 1
GAS_PHASE 1
   -fixed_volume
   -pressure 600
   -temperature 60
   -volume 0.00000001
   CO2(g) 0
SAVE solution 1
END

USE SOLUTION 1
   EQUILIBRIUM_PHASES 1
   CO2(g) 2.39   
END


[bashuisman]

Your suggestion not only uses 245ATM for the gas phase, but also for the water phase, and that initially confused me quite a bit, as I always kept my solution at 600ATM. But now I thought about it longer, your suggestion is actually very useful to me, so thank you very much! It tells me easily, and without the need for scripting, how much the fluid is under-/oversaturated w.r.t pure CO2(g). I would have to drop the pressure of the entire vessel to 245ATM before pure CO2(g) comes bubbling out. This is a useful number and something easy to communicate to "peers".

However, if I want to know what gas composition at 600 ATM is in equilibrium with my solution, let's say in a methane gas reservoir, then I guess X_CO2=245/600=41% is probably not the right answer, because as you say, the fugacity coefficient of CO2 depends on the presence of the other constituents of the gas phase.

In my little experiments I kept the solution pressure at 600 ATM, and then, using EQUILIBRIUM_PHASES, I needed to use a much higher partial pressure of CO2 (83%), to equilibrate with my solution. Adding methane to EQUILIBRIUM_PHASES doesn't really help because EQUILIBRIUM_PHASES doesn't calculate a different fugacity coefficient for CO2 in the presence of methane. I'm wondering if the difference between 83% CO2 and 41% CO2 is entirely due to a difference in calculating the fugacity in the gas phase though, because at higher pressures also e.g. the pH is different, and I can imagine that the different activity coefficients in the water phase at higher pressures allow for different amounts of CO2 to dissolve into that solution.

I tried using a fixed pressure GAS_PHASE with CO2 and Mtg, but I haven't been able to calculate a CO2/Mtg composition in the gas phase that at 600 ATM is in equilibrium with this solution without it either changing the solution or failing to converge. E.g. if I use 83% CO2 and 17% Mtg, a lot of CO2 (and Mtg, but I expected that, since there's 0 Mtg in the initial solution) goes from the gas phase into the water. If I increase the gas phase volume it fails to converge, and if I increase the Mtg/CO2 ratio in the gas phase it also fails to converge.

Any idea how I could calculate which gas phase (Mtg+CO2) composition my solution is in equilibrium with? Using your fixed_volume with a very small volume trick doesn't work, because there's no Mtg initially in the solution, and so it ends up giving me the 245ATM answer again with 0 Mtg in the gas phase.

I tried using the result of an EQUILIBRIUM_PHASES with 83% CO2 and 17% Mtg. But then the result in the selected output looks very strange to me. si_CO2 is again 2.19, but with a partial pressure of only 2.6ATM, a PR_PHI close to 1, and 63% CO2 in the gas phase. So the equation SI=log10(PR_P*PR_PHI) doesn't hold anymore here, and that makes me a bit scared to use that 63% number, although it is nicely in between 41 and 83% :-).

I've tried to illustrate all the things I've tried in the following code

Code: [Select]

DATABASE  ../databases/pitzer.dat

SELECTED_OUTPUT
   -reset false
   -totals C
   -saturation_indices CO2(g)
   -pH true
   -file tTP2.txt
   USER_PUNCH
   -headings PRESSURE PHICO2 lPCO2
   10 PUNCH  PRESSURE
   20 PUNCH  PR_PHI("CO2(g)")
   30 PUNCH  log10(PR_P("CO2(g)"))
END

SOLUTION_SPREAD
-units mol/kgw
NUMBER     pH    temp   Alkalinity  C        Ca    Cl      Na  pressure 
                                                    charge                  
   1       3.95   60   1.75e-02    0.772    0.38  3.14    2.4 600       
END

#l2):
USE SOLUTION 1
   GAS_PHASE 1
      -fixed_volume
      -pressure 600
      -temperature 60
      -volume 0.00000001
      CO2(g) 0
END

#l3):
USE SOLUTION 1
   EQUILIBRIUM_PHASES 1
      CO2(g) 2.389323232305 10
END

#l4):
USE SOLUTION 1
   REACTION_TEMPERATURE
      60
   REACTION_PRESSURE
      2.450886684349675e2
END

#l5):
USE SOLUTION 1
   REACTION_TEMPERATURE
      60
   REACTION_PRESSURE
      600
   EQUILIBRIUM_PHASES 1
      CO2(g)  2.695425805
END

#l6):
USE SOLUTION 1
   REACTION_TEMPERATURE
      60
   REACTION_PRESSURE
      600
   EQUILIBRIUM_PHASES 1
      CO2(g) 2.695425805
      Mtg(g) 2.017299704488154
      SAVE SOLUTION 2
END

#l7):
USE SOLUTION 1
   GAS_PHASE 1
      -fixed_pressure
      -pressure 600
      -temperature 60
      -volume 5
      CO2(g) 0.826560323782619
      Mtg(g) 0.173439676217381
END


#l8):
USE SOLUTION 2
   REACTION_TEMPERATURE
      60
   REACTION_PRESSURE
      600
   GAS_PHASE 1
      -fixed_volume
      -pressure 600
      -temperature 60
      -volume 0.000000001
      CO2(g) 0
      Mtg(g) 0
END

with output:

Code: [Select]

          pH            C    si_CO2(g)     PRESSURE       PHICO2        lPCO2
     4.01083   7.7200e-01       2.1922   6.0000e+02   1.0000e+00         -inf
      4.1205   7.7200e-01       1.9939   2.4509e+02   4.0232e-01   2.3893e+00
     4.20408   4.8262e-01       1.9939   6.0000e+02   4.0232e-01   2.3893e+00
      4.1205   7.7200e-01       1.9939   2.4509e+02   1.0000e+00         -inf
     4.01083   7.7200e-01       2.1922   6.0000e+02   3.1389e-01   2.6954e+00
     4.01114   7.7200e-01       2.1922   6.0000e+02   3.1389e-01   2.6954e+00
     3.97867   8.3748e-01       2.2257   6.0000e+02   3.0760e-01   2.7377e+00
     4.01114   7.7200e-01       2.1922   6.0000e+02   9.8437e-01   4.1150e-01 l1) Speciation. The reference case
l2) GAS_PHASE block: decrease of solution pressure and increase in pH
l3) EQUILIBRIUM_PHASES with input CO2(g) = PR_P(CO2(g)) of GAS_PHASE in l2). Uses initial solution pressure, CO2 comes out of solution.
l4) EQUILIBRIUM_PHASES with nothing more than a decrease in pressure: same as l2).
l5) EQUILIBRIUM_PHASES with input CO2(g) tuned to get output si_CO2(g)=output si_CO2(g) of l1): same as l1).
l6) EQUILIBRIUM_PHASES with input CO2(g) same as l5, and input Mtg(g) log10(600-10^siCO2(g)): almost same as l1).
l7) GAS_PHASE with input CO2(g) and Mtg(g) equal to l6: a lot of CO2 goes from the gas phase to the solution.
l8) uses l6 solution, GAS_PHASE with 0 initial CO2(g) and Mtg(g) and small fixed volume: very low PCO2 and high PR_PHI(CO2)

[dlparkhurst]

Here is a file that calculates a methane-CO2 GAS_PHASE that is in equilibrium with the solution. The GAS_PHASE is fixed volume, and Mtg is added incrementally while maintaining equilibrium with calcite. The pressure of the GAS_PHASE increases and at some point is equal to 600 atm, which is shown on the graph. Is that the point that you are looking for?

The CO2(g) fugacity changes a bit because of the reaction with calcite, but is still similar to the original fugacity calculated by the SOLUTION calculation.

Code: [Select]

SELECTED_OUTPUT
   -reset false
   -totals C
   -saturation_indices CO2(g)
   -file tTP.txt
   USER_PUNCH
   -headings PRESSURE PHICO2 lPCO2
   10 PUNCH  PRESSURE
   20 PUNCH  PR_PHI("CO2(g)")
   30 PUNCH  log10(PR_P("CO2(g)"))
END

SOLUTION_SPREAD
-units mol/kgw
NUMBER     pH    temp   Alkalinity  C        Ca    Cl      Na  pressure 
                                                    charge                  
   1       3.95   60   1.75e-02    0.772    0.38  3.14    2.4 600       
END

USE SOLUTION 1
GAS_PHASE 1
   -fixed_volume
   -pressure 600
   -temperature 60
   -volume 0.00000001
   CO2(g) 0
   Mtg(g)    0
KNOBS
-iter 400
USE SOLUTION 1
EQUILIBRIUM_PHASES
Calcite 0 10
REACTION
Mtg(g)  1
0.05 in 10 steps
USER_GRAPH 1
    -headings               Mtg_added Pressure_atm
    -axis_titles            "Mtg added, moles" "Pressure, atm" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 graph_x rxn
20 graph_y pressure
  -end
    -active                 true
END


[bashuisman]

Hi David,

thanks, that's fantastic! I think you've solved it for me!

It's surprising to me that a gas at 600C with 83.9% CO2 and the rest Methane dissolves the same CO2 in water as a water has that needs to drop to 41% in pressure to bubble out pure CO2. I would have expected these percentages to be much closer.

What's also interesting is that the impact of the presence of Mtg is very small on the PR_PHI of CO2(g) (ok, in this case, it's only 16.1% Mtg). This presence of Mtg leads to a difference in gas concentration of only 1.02% w.r.t. just calculating it iteratively with only CO2(g) using EQUILIBRIUM_PHASES. Although that difference could even be smaller if the initial solution would have been perfectly in equilibrium with Calcite.

This does strengthen my opinion that it could be very worthwhile to, for a next version of Phreeqc, add the calculation of PR_PHI and PR_P for gasses also at speciation. This means, like in EQUILIBRIUM_PHASES, assuming that each gas component exists in its own 100% pure bubble. I also still think that it is confusing to have the same variable name to mean two different things: SI for gas phases. SI input being a log partial pressure, and SI output being a log fugacity. I also wonder whether the fact that a GAS_PHASE with 83% CO2 and 17% Mtg gave me partial pressures of CO2(g) of 2.6ATM and a PR_PHI close to ideal is a bug (l7 in my previous code snippet)?

Thanks again, you've solved all my problems, and showed me a new way of calculating gas phase equilibria!

Cheers,
Bastiaan

[dlparkhurst]

I think the formation of bubbles depends on the sum of the partial pressures, so a solution would not bubble out pure CO2, but a mixture of CO2 and methane.

I do not think there are is any interaction between CO2 and Methane in our Peng-Robinson calculations. There are some interaction coefficients between H2O(g) and other gases (in the method calc_PR), but not other interactions. So, I think the phi(CO2) is just a function of pressure and temperature (~0.3 at 60 C and 600 atm, not sure how you got ~1).

I understand your suggestion about calculating phi, but it probably is not going to happen, and there are ways to calculate phi with GAS_PHASE, which is when Peng-Robinson actually applies. I also agree that the use of saturation index has always been confusing for gases, even before the addition of Peng-Robinson. SI is appropriate for minerals used in EQUILIBRIUM_PHASEs, but gases are treated in the same way as minerals internally, which give rise to the confusion. We tried to be clear in that the output specifically states that the SI is log fugacity for gases and the documentation in the 2013 manual says:

saturation index—Target saturation index for the pure phase in the aqueous phase (Line 1a); for gases, this number is the log of the partial pressure (Line 1b).

And the updated version in PhreeqcI states further:

saturation index --Target saturation index for the pure phase in the aqueous phase (Line 1a); for gases, this number is the log of the partial pressure [which is not equal to fugacity for Peng-Robinson gases, P = fugacity / (phi * 1 atm)] (Line 1b).