Getting bulk/species properties from Phreeqc

[vineetbansal]

In our wrappers, to get all bulk/species properties, we're formulating a phreeqc string like below (without the -file so.out which is included here to test against phreeqc UI). This is the main mechanism we use to interact with phreeqc:

SELECTED_OUTPUT
    -file so.out
    -reset false

USER_PUNCH
5 PUNCH CELL_NO, TOT['water'], OSMOTIC, EOL_NOTAB$
10 t = SYS("aq", count, name$, type$, moles)
20 FOR i = 1 to count
30 PUNCH name$(i), MOL(name$(i)), ACT(name$(i))
40 NEXT i
50 PUNCH EOL$
60 p = SYS("phases", count, name$, type$, moles)
70 FOR j = 1 TO count
80 PUNCH name$(j), SI(name$(j))
90 NEXT j

SOLUTION 0
  pH 7.0
  pe 8.5
  redox pe
  temp 25.0
  units mol/kgw
  water 0.9970480319717386
SAVE SOLUTION 0
END

The part between SOLUTION 0 and SAVE SOLUTION 0 is dynamic and generated using the Solution class. For example, the above solution would correspond to the call:

Solution({"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386})

The output we get in so.out is:

0.0000e+00	  9.9705e-01	  0.0000e+00	           
OH-	  1.0127e-07	  1.0123e-07	          H+	  1.0005e-07	  1.0002e-07	          O2	  8.3176e-25	  8.3176e-25	          H2	  7.0795e-35	  7.0795e-35	           
H2O(g)	 -1.5028e+00	       O2(g)	 -2.1188e+01	       H2(g)	 -3.1049e+01	

interpretable as:

<solution_num>, <kg_water>, <osmotic_coeff>
<species1_name>, <species1_amt_moles>, <species1_activity>, <species2_name>, ...
<phase1_name>, <phase1_si>, <phase2_name>, ...

This issue is created to ensure that this indeed gives us the correct values (unit testing seems to indicate that it does), and to discuss how to add add more bulk/species properties to be reported in so.out.

@rkingsbury mentioned that we may want:

- Bulk properties
  - conductivity
  - density
  - mass
  - volume
  - pH
  - pE
  - total dissolved solids (if available)

- Individual species properties
  - molalities
  - activities
  - activity coefficients  <-- currently calculated in pyEQL using activity+molality
  - diffusion coefficients
  - mobility (if available in PHREEQC)

- Equilibrium-related
  - Saturation indices for any minerals present

[vineetbansal]

Some specific notes on properties here. We can open separate issues for this if we want. I got these from the docs tab in phreeqc UI.

Conductivity

SC gives us conductance in uS/cm. This seems to match in dimensions with what pyEQL calculates as the conductivity property of a solution (in S/m).

Density

From the docs:

• RHO - Density of the solution (kg/L).
• RHO_0 - Density of pure water at the current temperature (kg/L).

Mass

From the docs:

• TOT("water") gives kg water. This is what we're using above.

Volume

• SOLN_VOL - Volume of the solution, in liters.

pH

• Using LA("H+") ?

pE

• Using LA("e-")?

Osmotic Coefficient

The docs indicate:

• OSMOTIC - Returns (1) the osmotic coefficient if the Pitzer model is used, or (2) 0.0 if the ion-association model is used.

  And this seems to be the case if we use pitzer.dat or the default db phreeqc.dat respectively - we get 0.99988 and 0.0 respectively. Interesting to note that pyEQL currently defaults this property to 1.0 recognizing that it cannot calculate it, while phreeqc uses 0.0.

In general, are specific properties only calculable if we swap databases in phreeqc when retrieving those properties? This is doable, but just making sure...

[rkingsbury]

Thanks @vineetbansal !

In general, are specific properties only calculable if we swap databases in phreeqc when retrieving those properties? This is doable, but just making sure...

In principle the properties we want should always be calculable no matter what database PHREEQC is using. In practice...I suspect this could be a little messy / rife with edge cases as you identified above with water activity.

To add some possibilities, see the following two excerpts from the PHREEQC docs. I'm not sure if all these variables would be available via SELECTED_OUTPUT or not, but could be a place to start.

This page describes arguments you can set on individual species in the input file. I suspect it's not possible to retrieve the species-specific values that are used in the database for -dw, -gamma, etc., but just adding this here for reference.

[rkingsbury]

But to summarize, from my quick reading of the docs and Vineet's comment above, it looks like the following are straightforward to get:

• conductivity (SC)
• solution volume (SOLN_VOL)
• water mass (TOT('water') or possibly -water)
• temperature ('-temperature`)
• alkalinity (-alkalinity)
• ionic strength ('-ionic_strength)
• charge balance (-charge_balance)
• pH (-pH)
• pE (-pe)
• total element concentrations (-totals')

and I'd like to include all of the above. We won't immediately use all of them, but some will be helpful for debugging, etc. later on.

[vineetbansal]

Thanks @rkingsbury - I saw that in the docs, but my initial attempts to generate any of those SELECTED_OUTPUT properties (in the UI or programatically) failed, so I was just using a PUNCH call like in my original comment above.

I tried again and now see that the following does work (specifically with -reset true):

SELECTED_OUTPUT
    -file so.out
    -reset true

and I do see a whole lot of bulk properties by default:

         sim	       state	        soln	      dist_x	        time	        step	          pH	          pe	    reaction	        temp	         Alk	          mu	    mass_H2O	      charge	     pct_err	
           1	      i_soln	           0	         -99	         -99	         -99	           7	         8.5	         -99	      25.000	 1.21651e-09	 1.00661e-07	    0.997048	-1.21292e-09	   -0.604265	           
         OH-	  1.0127e-07	  1.0123e-07	          H+	  1.0005e-07	  1.0002e-07	          O2	  8.3176e-25	  8.3176e-25	          H2	  7.0795e-35	  7.0795e-35	     

Along with the USER_PUNCH stuff below it. This seems to work both in the UI and through their C++ interface.

So I suppose any of these two approaches (-reset true or PUNCH CELL_NO, TOT['water'], OSMOTIC, EOL_NOTAB$ ...) will work. We can discuss how to select between these two approaches and whether they're telling us different things.