Simulating methane and CO2 gas in deep porefluids

[GeeqC]

Degassing of methane and CO2 in deep pore fluid could be implemented by a GAS_PHASE in the PhreeqcRM input. This gas phase would then form and equilibrate according to pressure, temperature and aqeous fluid composition.

However, because methane shows a peculiar behaviour at high pressure it may be more accurately calculated by the empirical formula of Tishchenko et al. (2005).

So far, I have simulated the gas phase separately in the C++ code, by equilibrating fluid and gas phase in the transport step, via a source/sink term.

Would it be possible to implement this long formula somewhere in PhreeqcRM, in order to run the equilibration synchronously during the speciation step?

Reference:
Tishchenko, P., Hensen, C.,Wallmann, K., and Wong, C. S. (2005). Calculation of
the Stability and Solubility of Methane Hydrate in Seawater. Chem. Geology.
219, 37?52. doi:10.1016/j.chemgeo.2005.02.008

[dlparkhurst]

I would first check how different the methane Peng-Robinson formulations in PHREEQC are from the equation of state that you want to implement. I think PHREEQC's EOS should be adequate to some T and P.

[GeeqC]

The two models can well be compared.

But unfortunately, I encoutered a further problem with setting the boundary conditions in PhreeqcRM:

std::vector<double> bc_conc;
std::vector<int> bc1(2,0), bc2(2,0);
bc1[0] = 1; //Upper boundary condition, fixed to SOLUTION 1
bc1[1] = 2; //Lower boundary condition, fixed to SOLUTION 2
status = rm.InitialPhreeqc2SpeciesConcentrations(bc_conc, bc1); //no mixing

The boundary conditions are correctly set in bc_conc, but they do not appear in the workers (c_species vector). From the example in advect_cpp.cpp it is not clear to me how the boundary conditions are set initially and during each time step.

[dlparkhurst]

The boundary conditions are up to the transport modeler.

You can have a cell in the transport model and  a corresponding chemistry cell that represents the boundary. You would set that boundary condition as part of SpeciesConcentrations2Module or SetIthSpeciesConcentration. You can use  InitialPhreeqc2SpeciesConcentrations to get the concentrations to set boundary cells in the concentration arrays used to set PhreeqcRM concentrations.

Often, the boundary is not treated as a separate cell. In PHAST the solution compositions for fluxes or constant concentrations are applied strictly by the transport model, and the values to use for the boundary calculations are taken from InitialPhreeqc2Concentrations (InitialPhreeqc2SpeciesConcentrations if you are working with species).

[GeeqC]

This is what I did, which seems correct.
But I became aware of another more fundamental problem, which is better discussed in another thread.

I will come back to the gas model later, providing a comparison between the Tishchenko and Peng-Robinson equations. Thanks so far!

[GeeqC]

Having clarified the redox-problem, I come now back to the calculation of the methane gas equilibrium. It seems, phreeqc correctly calculates the methane saturation concentration. However, the values (CH4 conc. [mM]) and the overall trend differ between the two models:

200 m waterdepth:
Phreeqc:  43.3 mM       Tishchenko: 36.7 mM

2000 m waterdepth:
Phreeqc: 227.0 mM       Tishchenko: 362.6 mM

I suppose the deviation is owed to the anomalous behaviour of methane at high pressure, but I am not able further comment on this.

There is a further issue: degassing and re-dissolution is often kinetically limited. Work by Boudreau, Dale, Mogollon and others have treated this topic (which involves complex models of bubble size) in detail. The question is, how it is possible to treat gas equilibration kinetics in Phreeqc?

This is essentially the analogue question as for redox equilibra. Can the equilibration between CH4(aq) and CH4(g) be switched off?

[dlparkhurst]

I believe the model in pitzer.dat and phreeqc.dat are fit to the data of Qin et al., 2008, J. Chem. Eng. Data 53, 1246-1249. You will have to decide if the methane model is adequate.

Both databases have a definition for Mtg, which is an unreactive (redox-decoupled) definition for methane. It looks like your reference at least starts with a Pitzer model, so you could adjust the Pitzer parameters, and log K for solubility.

As for partial equilibrium, the following will kinetically transform Mtg(g) to CH4(aq)/CH4(g) equilibrium. The rate is arbitrary.

Code: [Select]

PHASES
Mtg(g)
Mtg = Mtg
#-analytic 10.44 -7.65e-3 -6669 0 1.014e6 # CH4 solubilities 25 - 100?C
log_k -100
T_c 190.6; -P_c 45.4; -Omega 0.008
END
RATES
Mtg2CH4
10 k = 1/86400 # 1/day
20 rate = k * GAS("Mtg(g)")
30 moles = rate * TIME
40 SAVE moles
SOLUTION 1
END
GAS_PHASE 1
Mtg(g) 100
CH4(g) 0
END
KINETICS
Mtg2CH4
-formula CH4 +1 Mtg -1
-m0 100
-time  30*8640
END
INCREMENTAL_REACTIONS
END
USER_GRAPH 1
    -headings               time CH4(aq) CH4(g) Mtg(g)
    -axis_titles            "Days" "Molality" "Moles in gas phase"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME / 86400
20 GRAPH_Y MOL("CH4")
30 GRAPH_SY GAS("CH4(g)"), GAS("Mtg(g)")
  -end
    -active                 true
USE solution 1
USE gas_phase 1
USE kinetics 1
END

[GeeqC]

As you say, the rate is arbitrary.
Nevertheless, in order to implement this it would be important to calculate the rate as a function of departure from equilibrium:

rate = k * f[Mtg/Mtg(eq)]

Is it possible to retrieve the calculated equilibrium gas content [Mtg(eq)]?

If this procedure is used for an array of workers, the parameters would have to be calculated according to pressure and temperature for each cell. Is there a way the parameters are updated for each cell in each time step?

[dlparkhurst]

Code: [Select]

PHASES
Mtg(g) # No aqueous
Mtg = Mtg
#-analytic 10.44 -7.65e-3 -6669 0 1.014e6 # CH4 solubilities 25 - 100?C
log_k -100
T_c 190.6; -P_c 45.4; -Omega 0.008
CH4(g) # No gas
CH4 = CH4
log_k 100
#-log_k -2.8
#-analytic 10.44 -7.65e-3 -6669 0 1.014e6 # CH4 solubilities 25 - 100?C
-T_c 190.6; -P_c 45.4; -Omega 0.008
CH4(g)_logk # for logk calculation
CH4 = CH4
-log_k -2.8
-analytic 10.44 -7.65e-3 -6669 0 1.014e6 # CH4 solubilities 25 - 100?C
-T_c 190.6; -P_c 45.4; -Omega 0.008
END
RATES
Mtg2CH4
10 k = 1e-7 # 1/d
20 logk = LK_PHASE("CH4(g)_logk")
30 sr_calc = ACT("CH4") / (PR_PHI("Mtg(g)") * PR_P("Mtg(g)") * 10^logk)
40 rate = k * (1 - sr_calc)
50 moles = rate * TIME
60 SAVE moles
END
SOLUTION 1
END
GAS_PHASE 1
Mtg(g) 10
END
KINETICS
Mtg2CH4
-formula CH4 +1 Mtg -1
-m0 100
-time  20*8640
END
INCREMENTAL_REACTIONS
END
USER_GRAPH 1
    -headings               time CH4(aq) CH4(g) Mtg(g)
    -axis_titles            "Days" "Molality" "Moles in gas phase"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME / 86400
20 GRAPH_Y MOL("CH4")
30 GRAPH_SY GAS("CH4(g)"), GAS("Mtg(g)")
  -end
    -active                 true
USE solution 1
USE gas_phase 1
USE kinetics 1
END


[GeeqC]

Sorry for my slow reponse. I am trying to complete this thread now.

I do not understand in your response how the rate is calculated:
sr_calc = ACT("CH4") / (PR_PHI("Mtg(g)") * PR_P("Mtg(g)") * 10^logk)

What are these parameters? Are they arbitrary or real Phreeqc syntax?

Would it be possible to calculate a rate using parameters defined elsewhere in the C++ code?
For example:

sr_calc = k * ACT("CH4") / CH4_sat(i)

where CH4_sat(i) is the saturation concentration (or activity) of methane at each depth. I presume, the difficulty will be to transfer the information from the vector variable CH4_sat(i) to the Phreeqc array.

[dlparkhurst]

They are real Basic functions. Look at the section on The Basic Interpreter in the manual for PHREEQC version 3.

You can extract the information--ACT("CH4"), PR_PHI("Mtg(g)"), PR_P("Mtg(g)")--by using SELECTED_OUTPUT and USER_PUNCH (defined in an input file) and processing these selected output data with the methods available in PhreeqcRM.

Of course, you should start by considering the equilibrium case, thereby avoiding the need for Mtg and KINETICS.

Code: [Select]

GetCurrentSelectedOutputUserNumber
GetNthSelectedOutputUserNumber
GetSelectedOutputColumnCount
GetSelectedOutputCount
GetSelectedOutputHeading
GetSelectedOutputOn
GetSelectedOutputRowCount
SetCurrentSelectedOutputUserNumber
SetNthSelectedOutput
SetSelectedOutputOn