Saturation issue in boundary cell when only 1 phase shall exist.

[SaiP]

Hello David and others!!

I am trying to model a two-phase reactive transport case and stuck at an interesting issue. I account for compressibility of fluids and include necessary equations in my flow solver. Flow results are as expected without considering chemistry.
Now I considered coupling the flow solver with PhreeqcRM using the GAS_PHASE_MODIFY data block. The coupling works too for the moment.

I setup a case where, I have a column with two different temperatures. Top half of column has a temperature of 80C and bottom half of column is at 30C. The pressure at inlet (top of column) is 20.27 atm and at the outlet (bottom of column) it is 20.22 atm. Na, K ions are present in the system and no chemical interactions are considered (input file at the end of message). At these values of pressure and temperature, fluid exists only in liquid state and no phase change can occur. So, this means that the saturation of solution shall be, Sw = 1. However, when I run my case I notice that the saturation in the boundary cell at the outlet (bottom end of the column) goes on decreasing after using the RunCells() function. Out of curiosity I ran the simulation long to see it the saturation of solution moves towards 0. It stops at 0.573131. Not sure if this decimal signifies physical meaning but the system reaches steady-state with rest of the cells having solution volume of around 1.

Trying to debug, the saturation of water is computed as: Sw = (solution volume) / (Rep. volume * porosity).
In my case, the 'Rep; Vol.' and 'porosity' are both equal to 1 and do not vary. So, this implies that the solution volume is decreasing,. This potentially means that H, O ions are being consumed. When I try to check the concentration of H, O in the column they are between 106 - 110 moles of H and 54 moles of O at all times. I don't see why the Saturation of the boundary cell is decreasing then. I comment the RunCells() function (i.e. not running chemistry), and the saturation of the entire system is 1. So, there is something I am doing wrong in the PHREEQCs chemistry portion.

I use the following procedure to solve the reactive transport system:
1. Solve flow equation for Saturation, velocity data
2. use data from above to solve transport equation of both phases
3. SetSaturation() and SetConcentration()
4. RunCells()
5. GetSaturation() -> SetSaturation() -> GetConcentration()
6. Then few functions related to the gas_phase which are not active for this case.     

Phreeqc input data that I used for the case:

Code: [Select]

SOLUTION 0 #injection sol
temp 80
pH 7
Na 10 charge
END

SOLUTION 1 #bottom of column
temp 30
pH 7
K 10 charge
END

SOLUTION 2 #top of column
temp 80
pH 7
Na 10 charge
END

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

GAS_PHASE 1
-fixed_pressure
-pressure 20
H2O(g) 0
END

GAS_PHASE 2
-fixed_pressure
-pressure 20
H2O(g) 0
END

Thanks for any leads.

[dlparkhurst]

I think we have been through some of this before. The conversions are pretty tricky. I don't know the answer, but I have a few suggestions for debugging.

First, I would either remove sodium from the solutions or use equal concentrations for your initial solutions. In either case, I would charge balance with pH.  I would also set the pressure as part of the SOLUTION definition (although that should not matter if you are using SetPressure). The solutions you are using have small, but different sodium concentrations.

Code: [Select]

SOLUTION 0 #injection sol
temp 80
-pressure 20.27
pH 7 charge
Na 1
END
SOLUTION 0 #injection sol
temp 30
-pressure 20.22
pH 7 charge
Na 1
END

It seems odd to me that the saturation is changing just in the boundary cell. If it is a general problem with the saturation-concentration-transport cycle--because of concentration, pressure, temperature, or density--I think all of the cells would exhibit a systematic change in saturation. Because it is just the boundary cell, I suspect it is related to the PhreeqcRM calculation interaction with the boundary condition applied in transport.

One thing I would try is to cycle through RunCells two or more times without transport to examine whether the saturation decreases even without the transport calculation. That may help to isolate the problem.

One last thought is that the density calculated by PhreeqcRM may differ from the compressibility you are using in your transport. It sounds as if you may be using compressibility of water in your transport calculations. Transporting pure water would be a good check; however, with solutes, the density will be a function of solute concentration (and temperature and pressure), so, you will need to take that into account in your transport equations. 

The next step I would try is to try to simplify the problem--for example, see if the problem persist with sequentially constant temperature, constant pressure,  constant concentration, fixed concentration boundary, etc. If you can isolate to a PHREEQC or advection.cpp (PhreeqcRM) example, it would be easier for me.