phreeqcRM: initialization of boundary for non-diluted solution

[cyprien]

Dear PHREEQC users and developpers,


I coupled phreeqcRM with a transport library. I am facing some problems with the initialization of the boundary condition at the inlet in the case of non-diluted solutions.


My domain is simple: a one-dimensional column of 10 cells, with an inlet and an outlet. Initially, the column contains pure water, and I inject a concentration of 3 mol / L of NaCl:


SOLUTION 0
   units mol/L
   Na   3
   Cl   3
END

SOLUTION 1

END


The code returns (Na and Cl correspond to labels 3 and 4):

bc_conc[0] = 103.07456
bc_conc[1] = 51.537281
bc_conc[2] = -3.5986734e-08
bc_conc[3] = 3.3776924
bc_conc[4] = 3.3776924
bc_conc[5] = 8.7449619e-322


On the other hand, if I change the units in SOLUTION 0 by mol/kgw I got the following results:

bc_conc[0] = 104.42797
bc_conc[1] = 52.213986
bc_conc[2] = -3.3800986e-08
bc_conc[3] = 2.8220615
bc_conc[4] = 2.8220615
bc_conc[5] = 8.7449619e-322


In both cases, I do not get the expected value!

In the code, I specify that I am using mol / L for transporting the solution :

status = phreeqc_.SetUnitsSolution(2); // 1, mg/L; 2, mol/L; 3, kg/kgs

I also define:
status = phreeqc_.SetComponentH2O(false);
phreeqc_.UseSolutionDensityVolume(false);

The reference volume is defined as :

Info << "Set volume, initial porosity and saturation... ";
std::vector<double> rv,por,sat;
rv.resize(nxyz_, 1.0); // 1.0=1 L
por.resize(nxyz_, 1.0);
sat.resize(nxyz_, 1.0);

// Set representative volume
status = phreeqc_.SetRepresentativeVolume(rv);
// Set initial porosity
status = phreeqc_.SetPorosity(por);
// Set initial saturation
status = phreeqc_.SetSaturation(sat);


And to get the values in my inlet boundary:

Info << "Get boundary conditions... ";
//Get boundary conditions
std::vector<double> bc_conc;
std::vector<int> bc1;
int nbound = 1;
// bc1 is solution 0
bc1.resize(nbound,0);
phreeqc_.InitialPhreeqc2Concentrations(bc_conc,bc1);
Info << "OK"<<nl<<endl;



I am not sure what I am doing wrong during the initialization process. Is anybody have a clue about the problem?

Best,
Cyprien

[dlparkhurst]

Yeah, the conversions are a little tricky.

In GetConcentrations, you can chose whether you want to use the PHREEQC calculated density and solution volume, but in InitialPhreeqc2Concentrations, you do not get a choice. The mol/L results from InitialPhreeqc2Concentrations always use the PHREEQC calculated volume (moles divided by volume).

Short answer is that you need to use the correct density in your SOLUTION definition. You should use -density 1 calc, with will iterate to find the density for SOLUTION that results in the same calculated density in the output.

Here is a script that reproduces your results, and adds a calculation that calculates density.

Code: [Select]

USER_PRINT
10 PRINT "Moles:   ", TOTMOL("Na")
20 PRINT "Volume:  ", SOLN_VOL
30 PRINT "Density: ", RHO
40 PRINT "Mol/L:   ", TOTMOL("Na") / SOLN_VOL

SOLUTION 0
   -density 1
   units mol/L
   Na   3
   Cl   3
END
SOLUTION 0
   -density 1 calc
   units mol/kgw
   Na   3
   Cl   3
END
SOLUTION 0
   -density 1 calc
   units mol/L
   Na   3
   Cl   3
END


[cyprien]

Hi David,

Thanks a lot for the help! It's working well! That was a tough one. We spent a lot of time trying to figure out what was wrong...

Cheers,
Cyprien

[cyprien]

Hi,

I have a follow-up question. If I consider a more complex solution (SOLUTION 0) that is injected in a column full of pure water (SOLUTION 1) at pH 4.576 as described below, I always recover a pH in my first cell that is 4.69.

Code: [Select]


SOLUTION 0
temp                     119.18
    pressure              263.055
    pH                          4.576     
    pe                          -1.715
    -density         1.0 calc 
water           1.00
    -units mol/kgw
    Ba                1.065e-04
    C                 1.090e-02
    Ca                2.370e-01
    Cl                4.513e+00
    Fe                2.220e-03
    K                 2.072e-02
    Li                4.482e-03
    Mg                4.936e-02
    Mn                1.875e-04
    Mtg               9.594e-03
    Na                3.912e+00
    Ntg               1.732e-03
    S                 4.038e-03
    Si                1.559e-04
    Sr                5.022e-03
END

SOLUTION 1

END



Do you have any idea about the origin of this difference?


Best,
Cyprien

[dlparkhurst]

I assume you do not have any other reactants in cell 1, like equilibrium phases, that could change the pH.

A temperature or pressure change could change the pH.

It may be that you do not have pure advection in your model. There could be a little dispersion--specified, or numerical. If you are starting with pure water in the column, a little mixing between pure water and your boundary solution will tend to raise the pH. It could also depend on the way your boundary condition is implemented. There would be a difference between a flux boundary condition and a specified concentration boundary condition.

If you run for multiple pore volumes of cell 1, concentrations should approach the boundary solution by the second pore volume.

[cyprien]

Thank you David for the advice.

I found the issue: the PHREEQC module was not reading the pressure from the flow solver properly. This is fixed now.


Best,
Cyprien