Issue during initializing component concentrations.

[SaiP]

Hi David,

I have a question regarding initialising of component concentrations before commencing the run. That is, I try to explain what is happening before and during the 1st time step. I attach the corresponding PHREEQC files (water.ph, air.ph).

What I have done:
I use the initial conditions keyword 'ic1' to activate the concentration of 'SOLUTION' and also read data of the 'GAS_PHASE'. I provide information of the SOLUTION and GAS_PHASE in the PHREQC file and run 'RunCells()' for the 1st time here to initialize the components concentration.
 
NOTE:
1. there are no mineral interactions within the system and I only model fluid flows,
2. column is completely water saturated (S_w = 1)
3. the representative volume, porosity are 1 (always).
4. Units: mol/l of cell.

Corresponding PHREEQC file attached: 'water.ph'
I initially, run a case with Temperature, T = 25C and Pressure, p(inlet) = 25 atm and p(outlet) = 15 atm with spatial discritisation of 20 cells. At these conditions, the liquid shall continue to be a liquid (i.e. no phase change to gas). I vary the bool option of "UseSolutionDensityVolume(true/false)" and my results are as follows:

When UseSolutionDensityVolume(true):
The saturation during the initialization has values as follows:

Code: [Select]

saturation = 0.998929
saturation = 0.998952
saturation = 0.998974
saturation = 0.998997
saturation = 0.999019
saturation = 0.999042
saturation = 0.999064
saturation = 0.999087
saturation = 0.999109
saturation = 0.999132
.
.
.

and the Hydrogen concentration is:

Code: [Select]

110.803
110.8
110.798
110.795
110.793
110.79
110.788
110.785
110.783
110.78
.
.
.

When I commence the run with saturation (of water) 1, I wonder why the saturation reduces below 1 considering there is no phase change? This is what happens during the initialization stage.

After running the first 'RunCells' at the 1st time step, the saturation goes back to 1 but the Hydrogen concentration remains the same as above. I wonder if by using UseSolutionDensityVolume(true), we consider initial solution volume as a reference volume and calculate the saturation for the rest of the simulation with respect to this reference saturation. Is my hypothesis true? 

On the contrary, when I run the same case by using UseSolutionDensityVolume(false),
The saturation and hydrogen concentrations during initialization are as follows:
Saturation:

Code: [Select]

saturation = 0.998929
saturation = 0.998952
saturation = 0.998974
saturation = 0.998997
saturation = 0.999019
saturation = 0.999042
saturation = 0.999064
saturation = 0.999087
saturation = 0.999109
saturation = 0.999132
.
.
.
and Hydrogen concentration is:

Code: [Select]

110.684
110.684
110.684
110.684
110.684
110.684
110.684
110.684
110.684
110.684
.
.
.
After entering the first time step and using 'RunCells', the saturation doesnt go to 1 but continues to reduce at the inlet from 0.998929 slowly. And the Hydrogen concentration stays at:

Code: [Select]

110.803
110.8
110.798
110.795
110.793
110.79
110.788
110.785
110.783
110.78
.
.
Here, as I set the bool of UseSolutionDensityVolume to 'false', perhaps, the initial saturation is not being taken as a reference?
--------------------------------------------
Corresponding PHREEQC file: air.ph
I noticed another interesting issue about using UseSolutionDensityVolume(true). When I consider the column to be at T=150C, p=1 atm(approx), the fluid configuration shall be gaseous. During the initialization stage, I have the saturation as follows:

Code: [Select]

saturation = 1.74697e-06
saturation = 1.74697e-06
saturation = 1.74697e-06
saturation = 1.74697e-06
saturation = 1.74697e-06
saturation = 1.74697e-06
saturation = 1.74697e-06
saturation = 1.74697e-06
saturation = 1.74697e-06
saturation = 1.74697e-06
.
.
.
After entering the 1st time step and using 'RunCells', the saturation becomes 1 everywhere in the column. The reference saturation is then set to 1.75e-6 and hence perhaps the saturation stays at 1 according to the observation above? Ideally, the saturation here shall be 1e-6 at all times. But when I set the UseSolutionDensityVolume(false) the saturation is as expected, at 1e-6.

I dont know if my hypothesis is true but is there a solution that can be framed either we have liquids or gases in the system?

Thanks!!

[dlparkhurst]

I'm not sure what SOLUTION and GAS_PHASE definitions you intended.

For water.ph, PHREEQC reads the file until it hits an END. It then does all of the initial SOLUTION calculations and calculates the number of moles of gas components in each GAS_PHASE definition. At this point there are definitions for SOLUTIONs 0-50 and GAS_PHASE 0 and 1. All the solutions are the same and both the GAS_PHASES have 0 moles of components.

Here is where it gets tricky because of the lack of END statements. The last "USE solution 0" will be reacted with the first GAS_PHASE definition--GAS_PHASE 0 and the resulting solution will be saved as solutions 1-50. There are no SAVE gas_phase statements, so the GAS_PHASES will still have 0 moles of components.

The remainder of the file has no SAVE statements, so effectively you are defining pressures for the GAS_PHASES with no components present.

At 25 C, a gas phase will not be present, so the definitions are probably what you want, although the logic of the USE and SAVES is not correct.

I think this is what you should define

Code: [Select]

SOLUTION 0
temp 25
pH 7
Na 10 charge
END

GAS_PHASE 0
 -fixed_pressure
 -pressure 24.75
 H2O(g) 0
END

You would use these definitions (SOLUTION 0 and GAS_PHASE 0) to populate your cells and then, before RUN_CELLS, set the pressures with with the PhreeqcRM SetPressure method. Alternatively, you can define the pressures with the series of GAS_PHASE definitions in the input file and not call SetPressure before RunCells.

It is complicated, so with the uncertainty in your initial definitions I might not get it right the first time, but with a value of true, you have defined solutions at 1 atm pressure and 25 C. When you use RunCells, the pressure from the GAS_PHASE definitions is used (15-25 atm) in RunCells and the solution is compressed to a smaller volume. When that smaller volume is sent back to your transport code, you get saturations less than 1. The concentration of H is thus greater, the moles of H is the same but the volume is smaller, so mol/L is greater.

With a value of false, the results of RunCells are scaled to 1 L of pore space based on saturation (1), RV (1), and porosity (1), but also a density of 1. I think the difference between the calculated density and the assumed density may be the reason the saturation is not exact, but it would help if the initial condition was well defined (i.e. pressure for the solutions).