EWASG

The message tells you there is a phase transition from L to L+G at the very first iteration, which is not allowed. The appearance of gas phase occurs because of a wrong initialization of the third primary variable X3=10.0, which makes no sense, as it is interpreted as the mass fraction of H2 in the liquid phase.

$$$$$$$$$$$ GAS PHASE EVOLVES AT ELEMENT *000 0*  $$$$$    XH2  = 0.100000E+02   PX = 0.106144E+08 Pa   PBOIL = 0.817076E+10 Pa
 }}}}}}}}}}} ERRONEOUS DATA INITIALIZATION }}}}}}}}}}           STOP EXECUTION----------

I guess you want to initialize a single-liquid system  for which X3 should be the mass fraction of hydrogen dissolved in the aqueous phase. This value must be lower than H2 solubility at local P&T conditions.

Have a look to the user's guide and the informative printout in the output file.

 ***************************************************************************************

 THE PRIMARY VARIABLES ARE
  P - PRESSURE          XSM - SALT MASS FRACTION IN LIQUID (XS)  OR  SOLID SATURATION  (SS+10.)
  X3 - NCG MASS FRACTION  OR  SG+10. - GAS PHASE SATURATION + 10.         T - TEMPERATURE

 ***************************************************************************************
 *         COMPONENTS         *                                *   FLUID PHASE CONDITION        PRIMARY VARIABLES   * ************************************************************************************** *
  *       # 1  -  WATER     *                                        *   SINGLE-PHASE GAS             P,     XSM,          X3,   T *
  *       # 2  -  NaCl         *                                        *   SINGLE-PHASE LIQUID      P,     XSM,          X3,   T *
  *       # 3  -  NCG         *                                        *   TWO-PHASE                           P,     XSM,  SG+10.,   T *
  *       # 4  -  HEAT        *                                        

 ***************************************************************************************

Be aware that the thermodynamic approach in EWASG is not accurate at very high partial pressures of the NCG, such as in your simulation of H2 storage in an aquifer with pressures around 100 bar.

Regards,

Alfredo

Dear Gao,

H2 was included as an optional NCG inside EWASG with a treatment of water-H2 mixture very similar to what is done in EOS5.

The density of gaseous H2 is evaluated with the ideal gas law in both EOS.

The viscosity of gaseous H2 is evaluated with the same routine using data from Vargaftik (1975).

The enthalpy of gaseous H2 is computed using a constant specific heat in EOS5, while is computed according to Irvine and Liley (1984) correlation in EWASG.

One of the major differences is in the H2 solubility in liquid water. In EOS5 there is a linear relationship  for Henry's constant  between 0 and 25°C, then the 25°C value is kept constant at T>25°C.

In EWASG H2 Henry's constant is a function of T from 0 to 370°C. The two Henry's constants are shown in the figure below.

In addition, in EOS5 the contribution of dissolved H2 on the enthalpy of aqueous phase is neglected, while in EWASG is accounted for. But H2 solubility is very low, so related effects should not be very important.

Another difference is related to the effects of dissolved salts, in EWASG simulated using NaCl as a proxy fro TDS. If the salinity in your aquifer is low (less than a couple tens of thousands ppm), the pure water properties of EOS5 are almost equivalent to the brine properties computed by EWASG. At higher salinities, you may start to see the differences between water and brine properties.

To have an idea of the reliability of EOS5 and EWASG EOS approaches in the calculation of phase properties at high H2 partial pressure, what I can suggest is to use the two EOS to compute H2 properties as function of P&T and then compare that properties to a suitable source, such as the on-line NIST calculator for fluid properties (https://webbook.nist.gov/chemistry/fluid/). In that way you may chose the more reliable EOS and have an idea about the approximations involved in using it.

Regards,

Alfredo

The T is changing at the injection element T1717:

 T1715 34118 0.10028E+08 0.29773E+02 0.49812E+00 0.50188E+00 0.00000E+00 0.19998E-01 0.99632E+00 0.15183E-03 -.10815E+04 0.10000E+01 0.10125E+04
 T1716 34119 0.10039E+08 0.27806E+02 0.54849E+00 0.45151E+00 0.00000E+00 0.20014E-01 0.99672E+00 0.15324E-03 -.11102E+04 0.10000E+01 0.10131E+04
 T1717 34120 0.10059E+08 0.14330 E+02 0.61696E+00 0.38304E+00 0.00000E+00 0.20113E-01 0.99857E+00 0.16649E-03 -.11549E+04 0.10000E+01 0.10159E+04
 T1718 34121 0.10039E+08 0.27808E+02 0.54845E+00 0.45155E+00 0.00000E+00 0.20014E-01 0.99672E+00 0.15324E-03 -.11102E+04 0.10000E+01 0.10131E+04
 T1719 34122 0.10028E+08 0.29773E+02 0.49803E+00 0.50197E+00 0.00000E+00 0.19998E-01 0.99632E+00 0.15183E-03 -.10814E+04 0.10000E+01 0.10125E+04

The T change is low because the mass injection rate is very low, only 0.015 kg/s

GENER----1----*----2----*----3----*----4----*----5----*----6----*----7----*----8
T1717A0001    0    0    0    1     COM3  1.500E-02 0.000E+00 0.000E+00    

By the way, you have not assigned the enthalpy for the injected H2, which is then taken equal to 0. J/kg by TOUGH2.  You should evaluate the H2 enthalpy at injection conditios and use it in GENER.

I would also suggest using MOP(11)=2.

Dear Gao,

it seems you have assigned a max time of 1.5552E7 s under GENER

GENER----1----*----2----*----3----*----4----*----5----*----6----*----7----*----8
P1717A0001    0    0    0    3     COM3  0.000E+00 0.000E+00 0.000E+00                                        
  0.000000E+00  1.296000E+07  1.555200E+07                                                                    
  1.280000E-01  1.280000E-01  0.000000E+00   

Then your simulation simply goes beyond this time and TOUGH2 does not know which rate to use.

O3232(  17,  4) ST = 0.131071E+08 DT = 0.655360E+07 DX1= 0.155844E+02 DX2= -.346122E-08 T =  30.313 P =  9917425. S = 0.941512E-01
 0.262143E+08 IS OUTSIDE THE RANGE  ( 0.0000E+00, 0.1555E+08)
 AT KCYC =   18 EXCEED GENERATION TIME TABLE AT ELEMENT P1717 (SOURCE A00 1) -- WILL REDUCE DELTE

You should simply extend the time-rate data to the time foreseen in your input file (4.908E+14 s).

You are running a non-isothermal simulation (NEQ=4). Thus, you must assign the enthalpy of injected NCG, otherwise its enthalpy is assumed to be 0 J/kg.

You wrote you are first simulating the injection of an N2 cushion gas, to simulate the injection of H2 afterward. This is not possible in EWASG, as the EOS can simulate just one NCG, either N2 or H2, as it is clearly stated in the user's guide.

You need  an EOS able to simulate at least 2 gases, including N2 and H2. In the TOUGH2 release, I'm afraid no EOS is able to model N2 and H2 together. The TMVOC version released in 2002 can simulate mutiple NCGs, but at the time H2 was not included in the internal NCG databank. A TMVOC version with many more NCGs, including H2, has been developed at Saipem SpA (Modeling Multiphase Organic Spills in Coastal Sites with TMVOC V.2.0, VZJ, 2014) but that version has never been distributed. 

TOUGH3 should not have an EOS with these capability, but others may give a more detailed answer about TOUGH3 capabilities.

N2 and H2 can be simulated with EOS7 and EWASG versions supported by TOUGH4, which is under testing and not yet released. You may look at the invitation to test TOUGH4 (https://tough.forumbee.com/t/x2yfrfj#:~:text=Invitation%20for%20TOUGH4%20testing).

Regards,

Alfredo