Question about gas solubility calculation in Phreeqc

[Christinali91]

Hi David,

I'm looking into a problem with CO2 in water equilibrium at relatively high P&T and have some questions about gas solubility calculation using Peng-Robinson EOS in the Phreeqc source code.

First, does Phreeqc use the equation given in the Phreeqc manual version 3.0 page 445 to calculate gas solubility?

The solubility of gas i is given by m_i = K_H * (φ_i*P_i)/γ_i,where m is the molality, γ is the activity coefficient in water, K_H is the equilibrium constant, P is the partial pressure, and φ is the fugacity coefficient (the activity coefficient in the gas phase).

From my understanding, the corresponding code for this part should be in function calc_fixed_volume_gas_pressures, gases.cpp, line 647. The code is attached below. However, I didn't find equilibrium constant K_H in the equation. If Phreeqc use equilibrium constant K_H, could you help me point out where K_H is defined?

Code: [Select]

if (pr_done)
{
lp = phase_ptr->p_soln_x / gas_phase_ptr->Get_total_p() *
gas_phase_ptr->Get_volume() / gas_phase_ptr->Get_v_m();
phase_ptr->moles_x = lp;Second, in gases.cpp line 342, it mentioned a correction of the solubility of gas. Is there a reference on how we get this correlation? Is this correlation specific for water?

Code: [Select]

  3) correct the solubility of gas i with:
  pr_si_f = log10(phi_i) -  Delta_V_i * (P - 1) / (2.303 * R * TK);
Thanks so much for your help!

Best,
Christina

[dlparkhurst]

Yes, PHREEQC uses that equation. When the numerical method has converged, that equation will obtain--the ratio of the aqueous activity to the gas phase fugacity will equal the Henry's law constant.

I think you need to look at the code just above the section you cited. The first line below sets lp equal to the log10(KH), and the rest of that block before if(pr_done) calculates the partial pressure of the gas from the aqueous activities. The code that you cite calculates the moles of gas in the gas phase from the total pressure, volume, and partial pressure of the gas (calculated from activities). The moles of gas in the aqueous phase is calculated elsewhere, probably in "molalities". As the numerical iterations proceed, the values of the aqueous activities, moles of gas, and total gas pressure are revised until convergence is satisfied.

Code: [Select]

        lp = -phase_ptr->lk;
for (rxn_ptr = phase_ptr->rxn_x->token + 1; rxn_ptr->s != NULL;
rxn_ptr++)
{
lp += rxn_ptr->s->la * rxn_ptr->coef;
}
phase_ptr->p_soln_x = exp(LOG_10 * (lp - phase_ptr->pr_si_f));

if (pr_done)
{
lp = phase_ptr->p_soln_x / gas_phase_ptr->Get_total_p() *
gas_phase_ptr->Get_volume() / gas_phase_ptr->Get_v_m();
phase_ptr->moles_x = lp;
}

The equation for the fugacity coefficient is given in "Notes 1" of the PHASES keyword in the version 3 manual. You can to Appelo, Parkhurst, and Post, 2014 (http://www.lmcg.ufpe.br/~leo/Modelagem_Hidrogeoquimica/docs/AppeloHydrogeochem.pdf) for more discussion, or to the original Peng and Robinson 1976.