Issue with Saturation field update (using PhreeqcRM).

[SaiP]

Hello everyone,

I could establish communication between PHREEQC and my CFD package using PhreeqcRM. I am interested in looking at two phase flows.

My test case is a 1D column with no inflow or outflow (can be seen as a batch reactor) at 25C temperature. I initially have the column filled completely with water having Na. I expect the final solution to be identical to my initial solution as there are no chemical reactions involved and theres no phase change due to low temperature.
For the case, porosity is 1 and representative volume (rv) is 1L.

I use the SetSaturation(sat) and GetSaturation(sat) functions for updating the saturation field every time step (in future I would like to use it when there is phase change). In my CFD package I restrict saturation to be bounded between 0.001 - 0.999 (as 0 - 1 can be complicated for some numerical calculations). When I set saturation within the column as 0.999 (entirely filled with aqueous solution), and run the case for several time steps, I notice that the saturation value goes on decreasing every time step (and for some cases it exceeds 1 but because I restrict the value to 0.999 the simulation doesnt crash). That is, if I run the case for very long, the entire aqueous solution would probably vanish.

I verified that the saturation value is sent correctly from my CFD solver to PHREEQC. It is only after the reaction step, the solution volume and consequently the saturation decreases (as my porosity = 1 and rv = 1L, the values of solution volume and saturation match).       

Any clue why is this happening and advice how to cure this? I notice that if I use much smaller time step, the saturation plumments faster. It looks like the error is aggrevating by performing more time steps. 

For reference, I am using Phreeqc.dat and my Phreeqc input file is attached. 

Many thanks.

[dlparkhurst]

My guess is that you are not handling the concentration of H2O correctly.

Your choice of units (SetUnitsSolution) makes a difference. Assuming you are using mol/L, then the moles of elements and water used in the Phreeqc calculation are calculated by c*sat*RV. My guess is that the c(H2O) you are using results in a slightly smaller solution volume when the Phreeqc calculation is performed. Saturation is calculated as soln_vol*RV.

Now, if you use the results from the calculation (GetConcentrations) for the next calculation, you should not get a continuous decrease in saturation. Are you continually using a fixed mol/L for H2O, rather than the result from GetConcentrations?

The intended sequence is to use InitialPhreeqc2Module to set the initial concentrations, RunCells, and then use GetConcentrations to get the concentrations for transport. If you skip the transport step and simply reintroduce the concentrations (from GetConcentrations),  after another RunCells, GetConcentrations should return the same concentrations and the same saturation (assuming no KINETIC reactions).

If you try this approach and still have a problem, then make another post.

[SaiP]

Thanks David for your reply.

My "SetUnitsSolution" is 'mol/l'. I use "GetConcentration()" to update component concentration at run time.

As a check, I output concentrations of Hydrogen and Oxygen during run time. The Hydrogen concentration all through the run remains at 110.684 mol/l and Oxygen concentration at all times is 55.3421 mol/l. So, the water concentration is 55.34 mol/l. To make things easier to debug I didn't take into account any minerals (so no kinetics). As you mentioned - concentration of the components stay the same but not saturation (saturation keeps reducing by 4e-6-5e-6 every time step, that is why if I use smaller time steps I noticed more aqueous solution vanish).

I do follow the same procedure you mentioned - (1) initialize components; (2) SetSaturation(sat) & SetConcentration(c); (3) RunCells and (4) GetSaturation(sat) & GetConcentration(c).     

Many thanks!

[dlparkhurst]

Try setting this parameter to true:

Code: [Select]

phreeqc_rm.UseSolutionDensityVolume(true);


[SaiP]

Thanks David!!

"UseSolutionDensityVolume(true)" works for this case but there are few "interesting" issues that got me confused.

1. I output the solution density using "GetDensity". Output is 0.997048. I use this as initial value for "SetDensity()".
I initialise representative volume (RV) as "SetRV(nxyz,1.003974)".

I got the RV value as follows:
S = (sol_vol)/(RV*por);
0.999 = (1.00297)/(RV*1); //-> 1.00297 collected from PHREEQC batch reaction run
RV = 1.003974;
 
I anticipated results to match independent of "UseSolutionDensityVolume(bool)" choice. But, the saturation still keeps falling when the 'bool is false'. Why? Ideally, I am using same data for the solution density and volume in both cases right?

2. From my understanding, UseSolutionDensityVolume(true) will not work when porosity lies between 0 - 1 as PHREEQC calculations per se dont take into account porosity. Is my guess true? So, when porosity not equal to 1, we need to always use "UseSolutionDensityVolume(false)"?

3. I did a run considering phase change. Water converts to steam at 125C and 1.1 atm pressure (input file attached, I use phreeqc.dat as my database). When I use "UseSolutionDensityVolume(true)" the saturation field is not updated correctly. That is, I get saturation (liquid phase) of 0.999.
However, when I use "UseSolutionDensityVolume(false)" - I get correct results i.e. saturation of liquid is 0.001.
I think that, when we use "UseSolutionDensityVolume(true)" the initial solutions density and volume (before the chemical reaction) are taken into account but not the final solutions data. When there is complete phase change, the volume of liquid (and saturation) beomes 0. In this case how can we use "UseSolutionDensityVolume(true)"?

Many thanks for your time David.

[dlparkhurst]

1. I think the key is to set the saturation before you get the concentrations. When bool = false, v = saturation * porosity * rv is used to calculate mol/L; the number of moles in the solution is divided by v. So, unless you reset saturation before you get the concentrations, the saturations from the beginning of the calculation are used. If you do the following, the volumes resulting from RunCells are used to calculate mol/L.

Code: [Select]

status = phreeqc_rm.RunCells();
status = phreeqc_rm.GetSaturation(sat);
status = phreeqc_rm.SetSaturation(sat);
status = phreeqc_rm.GetConcentrations(c);

2. No, PhreeqcRM should work with non-unit porosity. PHREEQC itself does not use the porosity, but these conversions from mol/L to moles in the SOLUTION and from moles in solution to mol/L do use the porosity, saturation, and RV.

3. The volume of the solution should be just the volume of solution, and the saturation will be based on that. I don't know how you are handling the phase change. If you have a GAS_PHASE, where water partitions into the GAS_PHASE, then the volume of solution should decrease after the calculation, and the saturation should decrease; some of the water will end up in the GAS_PHASE. If you include only the liquid phase and no gas phases (no GAS_PHASE and no gases in EQUILIBRIUM_PHASES) when you SetConcentrations, then you will get the volume of the liquid phase after interactions with other reactants (EXCHANGE, SURFACE, non-gas EQUILIBRIUM_PHASES, etc).

[SaiP]

Hi David.

I have a two-phase isothermal system where I consider a column of 1 meter (1D, 400 cells). I consider “inlet @ top”, “outlet @ bottom”.  The temperature at the top half of the column is 25C and for the bottom half its 125C (Fig. 1). I don't solve the energy equation. I set the pressure for chemistry to 1.1 atm. The temperature, pressure remains constant all through my run while solving chemistry.

What I was expecting?
The top half of the column will always contain water (25C & 1.1 atm). After solving for flow, water enters the “temperature transition cell” (see fig1) which is at 125C & 1.1 atm. The water that enters this transition cell evaporates completely.

NOTE: When water completely evaporates in a cell (Saturation is around 1e-8), I find that PHREEQC is not happy and crashes. Hence, I restrict the saturation to be bounded between 0.001 – 0.999.

My results:
Though I notice evaporation taking place as water enters the transition cell (start of 125C), I observe the concentration of vapour linearly increasing with time in this specific cell (fig2, the colour shows saturation of water, the values show vapour concentration in the gas phase). I don't see any physical reason why this shall happen.
Trying to backtrack things: I notice the solution volume in this transition cell gradually increasing (a consequence of more H, O in this cell after transport). But I don't understand why H, O is accumulating over time in this cell after the chemistry time step? Ideally, all the water that has entered the transition cell after the flow has to evaporate.

I follow the procedure below:
Flow → Transport of components → SetSaturaton() → SetConcentration(c) → RunCells() → GetSaturation(sat) → bound Saturation b/w 1e-3 to 0.999 → SetSaturation(sat) → GetConcentration(c)

Additional information:
I was doing some number crunching. For the injection velocity (U = 1e-5 m/s) and a time step of 50 seconds, I notice the saturation in the transition cell after the flow is 0.2. So, in 10-time steps i.e. 500 seconds 2 litres of water enters the hot cells (125C). To make sense, I notice the vapour concentration in the transition cell increase by 111 for every 500 seconds. 

Thanks for all the help!!

[dlparkhurst]

I don't think I have completely grasped the details of your calculation.

It is a fundamental limitation of PHREEQC that there must be a solution; it cannot calculate a system without liquid water.

At temperatures above 125, the phreeqc.dat database will normally set the pressure to the pressure where liquid and gas are in equilibrium. At 125 C, that is about 2.29 atm. 

I'm not sure how you have used GAS_PHASE or EQUILIBRIUM_PHASES. Seems like you should use GAS_PHASE. If you use a fixed_volume gas phase, then H2O(g) will accumulate in the gas phase up to a pressure of 2.29, at which point, it will be in equilibrium with any water that is advected into the hot column. Not sure how it ran, but the attached file and graph shows the accumulation of gas in a large gas reservoir with time.

Not sure how this calculation compares to yours, but maybe it is a start for discussion.

[SaiP]

Thanks David.

I use "GAS_PHASE_MODIFY". A snippet from my output during runtime looks like this:

Code: [Select]

GAS_PHASE_MODIFY 198 ; \ #Contain gas
-type     0 ; \
-total_p  1.10492 ;\
 GAS_PHASE_MODIFY 199 ; \ #Contain gas - temperature transition cell
-type     0 ; \
-total_p  1.10494 ;\
 GAS_PHASE_MODIFY 200 ; \ #Contain water
-type     0 ; \
-total_p  1.10497 ;\

For -total_p: I use pressure data from my CFD solver.

I understand what you mean by the equilibrium pressure using -fixed_volume. But in my case, I use "-total_p/ fixed_pressure" rather.

I used "selected_output -gas" to see if I copy the correct pressure and for the above cells, I have the output as below:

Code: [Select]

#Cell 198
so op 0 (pressure) = 1.10492
so op 1 (moles of gas)= 52.3284
so op 2 (gas vol)= 1534.83
#Cell 199
so op 0 (pressure) = 1.10494
so op 1 (moles of gas)= 2385.94
so op 2 (gas vol)= 69980.3
#Cell 200
so op 0 (pressure) = 1.10497
so op 1 (moles of gas)= 0
so op 2 (gas vol)= 0

As you can see only in "cell 199", the moles of gas and gas volume are increasing. From Cells #0 - #198, I have same gas moles (52.3284) but the gas volume slightly varies due to the pressure variation.

Cell #199, is the transition cell where the temperature is 125C. All cells above 199 (200 and above) have a temperature of 25C and contain water only. It looks like the vapour is accumulating only in this cell (#199) as if there is no transport taking place. But I indeed solve for the transport of the gas components as well.

Do you think its physical that the vapour can accumulate in the transition cell even for this set-up? 

Many thanks!!

[dlparkhurst]

In cell 199, the temperature is 125 C. The partial pressure of H2O(g) for the solution is about 2.2 atm. So, basically, all of the water will evaporate. I'm not sure what your conceptual model is for a temperature of 125 and a pressure of 1 atm. You are boiling water into an infinite atmosphere; the volume and moles of gas will increase indefinitely. However, PHREEQC will not be happy if there is no liquid water.

If there are solutes in the water, it is possible that PHREEQC will calculate a very low activity of water at very high concentrations (although the numerical method is unstable at these high concentrations).

If gas-liquid equilibrium occurs, then K = a(H2O(aq))/f(H2O(g)),

If f(H2O(g))~P(H2O(g))~1, K ~0.5,  then if activity of water is about 0.5, equilibrium could obtain. For NaCl, phreeqc.dat calculates a(H2O(aq)) ~ 0.5 at a concentration of about 15 mol/kgw.

So, PHREEQC might calculate gas-liquid equilibrium at 1 atm, but it would be at high solute concentrations that are beyond model reliability, both geochemical and numeric.


If the temperature is 25 C in cell 200, and the pressure is 1.10497, then there is no gas phase. Don't know if you have any other gases (N2(g), O2(g), CH4(g), CO2(g)) defined for your gas phase, but the sum of the partial pressures of gases must be greater than 1.10497 for a constant-pressure gas phase to exist. If you have only H2O(g), then the temperature will need to be greater than 100 C for the partial pressure to exceed 1 atm so that a constant-pressure gas phase could form. If you transport liquid and H2O(g) into cell 200 (at 25 C), then all of the gas will be added to the liquid.

[SaiP]

You are boiling water into an infinite atmosphere; the volume and moles of gas will increase indefinitely.

Why is this so? Why cant the vapour formed in cell 199 be transported to cell 198 every time step rather, similarly the vapour in cell 198 be transported to cell 197 and so on? Something like advection of any component in a liquid phase?

However, PHREEQC will not be happy if there is no liquid water.

I realised this issue. Hence, I restrict the saturation to not drop below 1e-3. And the runs don't crash for now.

For now, I consider water to have only Na. There are no other gas components than H2O(g) in the gas phase. I am using phreeqc.dat as my dataset.

[dlparkhurst]

PhreeqcRM does not transport anything, so you will be transporting the liquid and gas phase. If you do that, and assuming no other reactions, the liquid and gas phase in cell 199 should reach a steady state, where the liquid (no gas) transported in will equal the gas (and liquid, if any) transported out.

In cell 198 (sorry, previously I thought transport was in the other direction), the gas transported in will distribute between a gas and liquid, but at 125 C, 1 atm, it should all be gas.

[SaiP]

Hi David.

Indeed as you say, even I expect the solution to reach steady-state. But perhaps something is not correct regarding the way I am solving a reaction time step, I guess my error comes from the 'GAS_PHASE'. As you mention in the earlier post, infact I do the transport of components using the CFD package :).

In my case I use:
"SetUnitsSolution(2)" with 2 - mol/l, and for the rest like "SetUnitsPPassemblage, SetUnitsExchange, ...", I use (1) with 1 - mol/l water.
However, I am not sure where is "SetUnitsGasPhase(0/1/2)" used. For now I set it to 1 - mol/l water to be consistent with my SetUnitsSolution(2) choice. I notice if I change between 0 (mol/l cell),1,2 (mol/l rock) there is suttle difference in results even when porosity is unity.

I use "Selected_Output: -gas" to write the output of total gas pressure, total gas volume, gas components concentration and then collect+store data from the third row [gas components concentrations]. But this data is in "moles". I use the same data in "moles" for transport of gas components and later to solve for the chemistry as well. Is this the right strategy? Or could there be an issue if I use mol/l for aqueous components (H, O) and moles for H2O(g)?

I use "Selected_Output: -gas" because the Set/Get_Concentration(c) are for aqueous phase only and do not work for the GAS_PHASE.

Thanks and continuing my debugging...!

[dlparkhurst]

The units choice for solutions (SetUnitsSolution) is used every time GetConcentrations and SetConcentrations are called.

However, the units choice for every other reactant (GAS_PHASE, EQUILIBRIUM_PHASES, SURFACE, etc) are only used when InitalPhreeqc2Module is called. These units define how the moles in the Initial Phreeqc instance, which are generally defined with GAS_PHASE, EQUILIBRIUM_PHASES, etc data blocks that are read from a file, are to be converted to moles in a transport model chemical cell. Once the initial conditions for the model are set, I don't think the SetUnitsX settings have an effect (other than SetUnitsSolution).

PhreeqcRM is missing pieces that would make multiphase transport easier, mainly something like GetGasPhase and SetGasPhase. I basically did not want to think about it. Seems like there are a few possibilities once you have transported liquid and gas.

(1) You could combine the moles in liquid and gas when you SetConcentrations, assuming you keep the conversions from transport units to PHREEQC moles correct (see SetUnitsSolution for the conversions). You would have to reinitialize the gas phase in each cell to have zero moles. A large number of moles in the gas phase relative to solution could cause numerical problems because the gas phase will be far from the ultimate equilibrium.

(2) You could use GAS_PHASE_MODIFY to update the PHREEQC GAS_PHASE to be consistent with the transported gas phase, provided you know the gas composition after transport. You would again have to be careful of how you convert to moles in the GAS_PHASE, given the conversion of the solution. 

(3) You could take the difference between the last reacted GAS_PHASE and the transported GAS_PHASE, and either add it in the SetConcentrations, or if that could cause negative concentrations, add the difference as a REACTION. You would not modify the GAS_PHASE definition in PHREEQC directly.

Most of these options require some cell-by-cell manipulation of something that PHREEQC keeps track of (GAS_PHASE, REACTION, possibly something else). So, you would need a RunString with the modifications (GAS_PHASE_MODIFY, REACTION, ...) to update the PHREEQC definitions, in addition to SetConcentrations, before you execute RunCells. I won't go into it now, but the parallelization requires care. For debugging, you should be able to define a string with the modifications for all cells and run it through all the workers, but it may be possible to do it more efficiently in other ways.

[SaiP]

Thanks David,

for now I will proceed with 'GAS_PHASE_MODIFY' as I have already coded it and from the output I write during runtime, I can see that I pass on the correct values. I use 'RunString' to update 'GAS_PHASE_MODIFY' each time step - to update the pressure from the CFD package as you mentioned (I use fixedPressure - varyingVolume for gas phase equilibrium).

I attach a figure to remind the case and also the output of moles during the debugging process when I turn off the RunCells (I explain how I am debugging below). I generalised the results instead of pasting all cell data which might confuse you.

I found something interesting in my results:
Out of curiosity to see what happens, I turned off "RunCells" (indeed I would need this but in the process of debugging I think its a valid try). So, I run the following in order:

1. Transport of aqueous and gaseous compoenents // I verified and the equations are good
2. SetSaturation(sat)                                           // Saturation data from flow solver
3. SetConcentration(c)                                        // c data for aqueous phase after (1. Transport equation)
4. Gas_Phase_Modify                                          // pass information about the pressure
5. GetSaturation(sat)                                          // output saturation
6. SetSaturation(sat)                                          // set saturation to get accurate concentration below
7. GetConcentration(c)
8. Selected_Output -gas                                     // output gas component con. only H2O(g)

I observe that I reach steady-state solution as we expected.
But theres a catch. My Hydrogen and Oxygen concentrations keep on piling to the previous time steps value in the transition cell. By 110 for H and 55 for O. Perhaps I believe this could be the potential reason why we observe the H2O(g) accumulating in the transition cell when I use RunCells.

But when I turn off RunCells, the observations are confusing. How can gas form (52.5 moles) but the Hydrogen and Oxygen concentrations dont reduce? They ideally shall all be consumed as these are the ingridients which form water? I can only think of GetConcentration(c) going wrong for aqueous phase but not sure why. I use SetUnitsSolution(2 - mol/l). Do you see some other potential leaks from my observations to focus on?


Thanks and wish you a nice weekend!!

[dlparkhurst]

If you do not use RunCells, the GAS_PHASE definition will have no effect on the solution composition that is returned with GetConcentrations.

Skipping RunCells may help debugging, but it is not completely straightforward.

The number of moles in solution, mol(n,i), when you use SetConcentrations(c) will be mol(n,i) = c(n,i) * RV * Por * Sat, for each component n and cell i.  I don' think the solution volume will be updated if RunCells is not called, so solution volume does not change from the previous value in each cell. So, the saturation will be the old solution volume v(k-1)(i) / (RV * por), where k-1 is the last time plane, and the values returned from GetConcentrations(c) will be mol(n,i)/v(k-1)(i).

I think results of GetConcentrations without RunCells should be the same with or without GAS_PHASE_MODIFY. I also think the saturation returned from GetSaturation will not change if RunCells is not run. If you run GetSaturation at step 1.5, and at step 5, I think you will get the same values. Let me know if I am wrong.

Your have a comment "// pass information about the pressure". Just to be clear, the key values in GAS_PHASE_MODIFY will be the number of moles, something like

Code: [Select]

  -component                 H2O(g)
    # GAS_PHASE_MODIFY candidate identifiers #
    -moles                   xxx