Converting from molality (mol/kgw) to mg/l for high solubility minerals

[Tom]

I read on the PHREEQC Q+A thread that the conversion from molality (mol/kgw) to mg/l for dissolved species requires a correction for the density of the final solution and mass of the solutes.

(here: http://wwwbrr.cr.usgs.gov/projects/GWC_coupled/phreeqc/mail/msg00054.html)

Rewritten here, where X is the molar mass of the phase

mg/L = (mol/kgw) * (X g/mol) * (1000 mg/g) * (1 kgw/L)

For most solutions where the solubility of dissolving phases is low, I imagine the correction is only minor and could even be ignored. For very soluble minerals I would guess the density would change.

Is it possible to output the density of a solution? Or is the mass of water considered to represent the density (ie mass of water (in kg) per litre)?

I have attempted to print density in SELECTED_OUTPUT without any luck. Alternatively could it be output using a BASIC command in USER_PUNCH?

Thanks!

[dlparkhurst]

Version 3 will calculate density when using the Amm.dat, phreeqc.dat, or pitzer.dat databases. The density is printed in the output section "Description of solution". In Basic, the function RHO returns the density (kg/L).

[Tom]

Thanks David, that's a great tip.

I had been using version 2.18 for MacOS so will experiment a little with V3 on a Windows PC.

The Basic command worked fine with the MacOS version.

For anyone else looking to output density, I used the following to output to the SELECTED_OUTPUT file:

USER_PUNCH
-headings density
-start
10 punch RHO
-end

[Tom]

Quick tip for output of mg/l without a density consideration (assumes 1kgw/litre):

Examples below for sulfate (S(6)) with molecular mass 96.06 g/mol.

Calculation in BASIC is:

mol/kgw * g/mol * 1000mg/g

The following writes to the SELECTED_OUTPUT file:

USER_PUNCH
-heading sulfate_mg/l
-start
10 PUNCH TOT("S(6)")*96.06*1000
-end

Alternatively the following writes to the general output file:

USER_PRINT
-start
10 PRINT "mg SO4/l", TOT("S(6)")*96.06*1000
-end

[Tom]

For anyone looking to compare the effect of a density correction here is a partial script that may help.

The USER_PUNCH data block outputs density (BASIC 10), then converts mol SO4/kgw to mgSO4/l without a density correction (BASIC 20), then converts mol SO4/kgw to mgSO4/l with a density correction (BASIC 30), then outputs the difference between the uncorrected and corrected value (in mgSO4/l) (BASIC 40). Finally, the percentage difference is calculated (BASIC 50).
--------------
#llnl.dat

SELECTED_OUTPUT
-file out1.xls

USER_PUNCH
-headings density mgSO4/l_uncorrected mgSO4/l_corrected diff %diff
-start
10 PUNCH RHO
20 PUNCH TOT("S(6)")*96.06*1000
30 PUNCH TOT("S(6)")*96.06*1000*RHO
40 PUNCH (TOT("S(6)")*96.06*1000)-(TOT("S(6)")*96.06*1000*RHO)
50 PUNCH (((TOT("S(6)")*96.06*1000)-(TOT("S(6)")*96.06*1000*RHO))/(TOT("S(6)")*96.06*1000*RHO))*100
-end

-----------------
Please let me know if you spot any errors for correction or identify any simplifications that could be made.

I suppose the BASIC could be simplified by defining each punch as a variable? But it probably helps with explanation to leave it written in full...

[dlparkhurst]

To do this calculation correctly, you will need to use pitzer.dat, phreeqc.dat, or Amm.dat. These databases have the molar volume information to calculate density accurately. The following gives the molality, mol/L, and mol/kgs concentrations for a solution saturated with gypsum at 25 C.

SOLUTION
END
USE solution 1
EQUILIBRIUM_PHASES 1
   Gypsum 0 10
USER_PRINT
10 PRINT "Density:   ", RHO
20 PRINT "SO4, mol/kgw (molality): ", TOT("S(6)")
30 PRINT "SO4, mol/L:              ", TOTMOL("S(6)") / SOLN_VOL
40 PRINT "SO4, mol/kgs:            ", (TOTMOL("S(6)") / SOLN_VOL) * RHO

Results:

Density:      9.9908e-01
SO4, mol/kgw (molality):    1.5085e-02
SO4, mol/L:                 1.5040e-02
SO4, mol/kgs:               1.5026e-02

[Tom]

Great! Thanks David, that is a considerable improvement which will certainly help myself and hopefully others.


For clarification, am I right in thinking this BASIC line would give mg SO4/L ?

50 PRINT "SO4, mg/l:            ", (TOTMOL("S(6)") * 96.06) / SOLN_VOL


and then this line would give mg SO4/kgs ?

60 PRINT "SO4, mg/kgs:        ', ((TOTMOL("S(6)") * 96.06) / SOLN_VOL) * RHO

Thanks

[dlparkhurst]

Yes. Note that TOT("S(6)") is the sum of all SO4-2 species (SO4-2, NaSO4-, CaSO4, etc), not just the bare ion SO4-2 (which would be MOL("SO4-2").

[dlparkhurst]

Think you also need a factor of 1000 mg/g. 96 is g/mol, and you want mg.

Note you can use GFW("SO4") to calculate the gram formula weight of sulfate.

[Tom]

Ahh yes, great spot, thanks.

I had copied it from the file where I was after simply g/l.

[John]

Thanks for making these conversions - they are helpful. It's unclear to me why you are multiplying by density, instead of 1/density. Shouldn't the conversion be:

(mol/kgw) x (1kgw/L) x (1L/kgs) = mol/kgs

What am I missing?

[dlparkhurst]

You're right. Should divide by rho.

60 PRINT "SO4, mg/kgs:        ', ((TOTMOL("S(6)") * 96.06) * 1000 / SOLN_VOL) / RHO

[GeeqC]

Dear All,

Related to the problem of molality, I think I have missed something basic: If I choose for a seawater solution a total C(4) concentration of 2.4 mmol/kg solution the sum of all species of inorganic carbon gives a concentration larger than 2.4 mmol/kg. Why does this happen?

In contrast, if the input C(4) is given in mmol/kg pure water, the result is correct, even though the concentration per kg pure water should be smaller than molality, as molality should be concentration in the total solution.

In the user manual it says that density of the total solution is not known before the calculation. I understand this problem. And yet in the latter case, the  concentration was calculated correctly.

Best wishes,
GeeqC

[dlparkhurst]

If you take a kilogram of solution, the mass of water is 1 kg minus the sum of the solutes. So if you have 2.4 mmol/kgs you need to divide by kgw/kgs, which is a number less than 1.0. mmol/kgw will be a number greater than 2.4.

I'm not sure what you mean by per kilogram pure water. The definition of molality is moles per kilogram water. Note that the following receipe does not make a 1 molal solution of calcium, even though you take 1 kilogram of water and add 1 mole of Ca. The reaction 2CaO + H2O = 2Ca+2 + 2OH- consumes water, so that the molality of Ca is 1.018 in the reacted solution.

SOLUTION
REACTION
CaO 1
1 mole
END

[GeeqC]

Dear Dr. Parkhurst,

Thanks! This explains parts of the problem. Now the question is, how I can proceed with the following:

I would like to use the output solution from a reaction-transport model as input for PhreeqC. The concentrations are given in Molarity (i.e. mmol per liter of solution). Phreeqc should then do the speciation, but leave the concentrations unchanged, as otherwise this would create a runaway effect in the reaction-transport model. Is there a way to solve this problem?

Best wishes,
GeeqC

[dlparkhurst]

No Dr.

PhreeqcRM is a version of PHREEQC intended to run chemistry for reactive transport models. It requires that H, O, charge, and each additional element in the system be transported by the transport code. (Alternatively, H2O, total H - 2*H2O, and total O - H2O can be transported instead of total H and total O.)

Three options are available for concentration units used by the transport code: mg/L, mol/L, or mass fraction, kg/kgs. PHREEQC uses moles of H, O, charge, and each element in its underlying solution definition (see PHREEQC manual SOLUTION_RAW and SOLUTION_MODIFY in the appendix). Here is the documentation from PhreeqcRM.

SetUnitsSolution

Solution concentration units used by the transport model. Options are 1, mg/L; 2 mol/L; or 3, mass fraction, kg/kgs. PHREEQC defines solutions by the number of moles of each element in the solution.

To convert from mg/L to moles of element in the representative volume of a reaction cell, mg/L is converted to mol/L and multiplied by the solution volume, which is the product of porosity (SetPorosity), saturation (SetSaturation), and representative volume (SetRepresentativeVolume). To convert from mol/L to moles of element in the representative volume of a reaction cell, mol/L is multiplied by the solution volume. To convert from mass fraction to moles of element in the representative volume of a reaction cell, kg/kgs is converted to mol/kgs, multiplied by density (SetDensity) and multiplied by the solution volume.

To convert from moles of element in the representative volume of a reaction cell to mg/L, the number of moles of an element is divided by the solution volume resulting in mol/L, and then converted to mg/L. To convert from moles of element in a cell to mol/L, the number of moles of an element is divided by the solution volume resulting in mol/L. To convert from moles of element in a cell to mass fraction, the number of moles of an element is converted to kg and divided by the total mass of the solution. Two options are available for the volume and mass of solution that are used in converting to transport concentrations: (1) the volume and mass of solution are calculated by PHREEQC, or (2) the volume of solution is the product of porosity (SetPorosity), saturation (SetSaturation), and representative volume (SetRepresentativeVolume), and the mass of solution is volume times density as defined by SetDensity. Which option is used is determined by UseSolutionDensityVolume.

[GeeqC]

Thank you for the detailed information.
If I understand you correctly, this is about the reaction module. However, I still use my own reaction-transport code and IphreeqCOM for data transfer.

But let's consider a more simple example: I would like to know the saturation index of halite in a 500 mmol/L NaCl solution. The following information is entered into the input file:

SOLUTION_SPREAD
    -units    mmol/l
  Cl     Na   
 500    500

The following appears in the output file:

Molality: 5.151e-1
Density (g/cm3): 1.01761
Volume (L): 1.01228

I assume the Molarity can be calculated by dividing the Molality by the volume and multiplying it with the density, which yields 0.518 mol/L. In other words, the speciation was calculated for a 518 mmol/L NaCl solution, i.e. the saturation index for a 500 mmol/L NaCl solution would not be correct. Where does the excess NaCl come from? And how can I modify the input so that the concentration in the resulting solution is 500 mmol/L?

[dlparkhurst]

The conversions are the same whether you use PhreeqcRM or not, so you should look at the response more carefully. The point is to use SOLUTION_RAW or SOLUTION_MODIFY in IPhreeqcCOM to define the new solutions so that you avoid the issues you are having with SOLUTION and SOLUTION_SPREAD.

A necessary condition is to have the correct density for the solution to convert from /L to /kgw with SOLUTION or SOLUTION_SPREAD. There will still be approximations involved with calculating the mass of solutes, so that for more complicated solutions, say with carbonates, you will not be able to make an exact conversion. The option "calc" will iterate on the density so that the density used in conversions is the same as the density calculated in the initial solution calculation.
 

Code: [Select]

SOLUTION
-units mmol/L
-density 1 calc
Na 500
Cl 500
USER_PRINT
10 PRINT "Na, mol/L: ", TOTMOL("Na")/soln_vol
20 PRINT "Cl, mol/L: ", TOTMOL("Cl")/soln_vol


[GeeqC]

Thanks for the explanations. I think there was an error in my calculation above. With the density and volume indeed the correct concentration of 500 mmol/l can be calculated. But as you say, it may not work for more complex solutions.

I tried "-density 1 calc", and it works. However, in the output file now the Molality (not the Molarity) is 0.5.

As far as I understand, SOLUTION_RAW and SOLUTION_MODIFY are written to a "dump" file. The question is then, how these data can be retrieved through IPhreeqCOM, as they are probably not accessible through "GetSelectedOutputArray"?

Having now two different options for renewing the parameters in the reaction-transport model, I wonder what the advantage of "SOLUTION_RAW/SOLUTION_MODIFY" is over "-density 1 calc" (The latter slows down the model run by ca. 30%, perhaps due to the additional iterations?).

[dlparkhurst]

The script in the previous post produces molality of 0.5061 mol/kgw corresponding to 0.5 mol/L using phreeqc.dat.

Total moles of each element can be written with USER_PUNCH, for example for hydrogen, with either the TOTMOL("H") function, or TOT("H")*TOT("water") [oops, previous version of this post had SOLN_VOL].

SOLUTION_RAW/SOLUTION_MODIFY does not require an estimate of the mass of solutes. If density and solution volumes are known, then the conversions should be exact. With SOLUTION, you still need to know density, but the conversion also requires an inexact estimate of solute mass. Further, it is not easy to define a SOLUTION with a given volume; the default is 1 kg water or some fraction of a kilogram of water.

[GeeqC]

Thanks for the helpful comments. The user punch was successful. For Molarity I assume it should say "TOTMOL/soln_vol" (?)

Now, there is one issue that I did not mention yet: for the reaction-transport model I would also need the molar concentrations of the carbonate species. However, I cannot find any output option for the molarity of a species, neither in "selected output" nor in "user punch".

[dlparkhurst]

Molarity of HCO3- is MOL("HCO3-") * TOT("water") / SOLN_VOL

Not sure what you are doing with your transport model, but, unless you are doing multicomponent diffusion, it is only necessary to transport H, O, charge, and each element.

[GeeqC]

Yes, this works! The molarities of partially dissociated species (e.g. CO3-2, HCO3-, CO2) are needed, because these species have different diffusion coefficients and are calculated separately in the reaction-transport model.

[dlparkhurst]

Still think you should take a look at PhreeqcRM. It does all of these conversions, and there is an option to get all aqueous species to allow for multicomponent diffusion. In addition, there is built-in parallelization. Much simpler than writing it all with IPhreeqc.

[GeeqC]

Thanks for the advice. I will consider PhreeqcRM for further modelling.

[GeeqC]

Implementing your methods above has led to a follow-up question. I use the input in the following form in PhreeqcRM to calculate total molarities in the selected output:

SOLUTION
-units mmol/L
-density 1 calc
Na 500
Cl 500

SELECTED_OUTPUT
-user_punch true

USER_PUNCH
-headings Na(mM) Cl(mM)
-start
10 PUNCH TOT("Na")/soln_vol
20 PUNCH TOT("Cl")/soln_vol
END

However, the parameter soln_vol seems not to be calculated. Values appear in the selected output as 1*+30.

[dlparkhurst]

First, TOT returns moles per kilogram water, but you want moles. To work generally, you should either multiply by the mass of water (TOT("water")) or use the function TOTMOL, which returns the total number of moles in the solution.

I cannot reproduce your problem. Using this script with phreeqc.dat

Code: [Select]

SOLUTION
-units mmol/L
-density 1 calc
Na 500
Cl 500
USER_PRINT
10 print "Volume, L:   ", soln_vol
20 print "Na, mol/L:   ", TOTMOL("Na") / soln_vol
30 print "Cl, mol/L:   ", TOTMOL("Cl") / soln_vol
END
produces this output:

Code: [Select]

   Input file: C:\temp\Phrqc7.pqi
  Output file: C:\temp\Phrqc7.pqo
Database file: C:\Program Files (x86)\USGS\Phreeqc Interactive 3.8.7-17149\database\phreeqc.dat

------------------
Reading data base.
------------------

SOLUTION_MASTER_SPECIES
SOLUTION_SPECIES
PHASES
GAS_BINARY_PARAMETERS
EXCHANGE_MASTER_SPECIES
EXCHANGE_SPECIES
SURFACE_MASTER_SPECIES
SURFACE_SPECIES
MEAN_GAMMAS
RATES
END
------------------------------------
Reading input data for simulation 1.
------------------------------------

DATABASE C:\Program Files (x86)\USGS\Phreeqc Interactive 3.8.7-17149\database\phreeqc.dat
SOLUTION
units mmol/L
density 1 calc
Na 500
Cl 500
USER_PRINT
10 print "Volume, L:   ", soln_vol
20 print "Na, mol/L:   ", TOTMOL("Na") / soln_vol
30 print "Cl, mol/L:   ", TOTMOL("Cl") / soln_vol
END
-------------------------------------------
Beginning of initial solution calculations.
-------------------------------------------

Initial solution 1.

----------------------------------User print-----------------------------------

Volume, L:      1.0121e+00
Na, mol/L:      5.0000e-01
Cl, mol/L:      5.0000e-01

-----------------------------Solution composition------------------------------

Elements           Molality       Moles

Cl                5.061e-01   5.061e-01
Na                5.061e-01   5.061e-01

----------------------------Description of solution----------------------------

                                       pH  =   7.000   
                                       pe  =   4.000   
      Specific Conductance (uS/cm,  25oC)  = 46780
                          Density (g/cm3)  =   1.01726 (Iterated)
                               Volume (L)  =   1.01211
                        Viscosity (mPa s)  =   0.93045
                        Activity of water  =   0.983
                 Ionic strength (mol/kgw)  =   5.061e-01
                       Mass of water (kg)  =   1.000e+00
                 Total alkalinity (eq/kg)  =   1.972e-08
                         Temperature (oC)  =  25.00
                  Electrical balance (eq)  =  -1.972e-08
 Percent error, 100*(Cat-|An|)/(Cat+|An|)  =  -0.00
                               Iterations  =   4 (20 overall)
                                  Total H  = 1.110124e+02
                                  Total O  = 5.550622e+01

----------------------------Distribution of species----------------------------

                                               Log       Log       Log    mole V
   Species          Molality    Activity  Molality  Activity     Gamma   cm3/mol

   OH-             1.575e-07   9.947e-08    -6.803    -7.002    -0.200     -2.92
   H+              1.309e-07   1.000e-07    -6.883    -7.000    -0.117      0.00
   H2O             5.551e+01   9.828e-01     1.744    -0.008     0.000     18.07
Cl            5.061e-01
   Cl-             5.061e-01   3.285e-01    -0.296    -0.484    -0.188     18.69
   HCl             6.891e-09   1.132e-08    -8.162    -7.946     0.215     (0) 
H(0)          1.260e-25
   H2              6.301e-26   7.079e-26   -25.201   -25.150     0.051     28.61
Na            5.061e-01
   Na+             5.061e-01   3.630e-01    -0.296    -0.440    -0.144     -0.63
O(0)          0.000e+00
   O2              0.000e+00   0.000e+00   -42.146   -42.095     0.051     30.40

------------------------------Saturation indices-------------------------------

  Phase               SI** log IAP   log K(298 K,   1 atm)

  H2(g)           -22.05    -25.15   -3.10  H2
  H2O(g)           -1.51     -0.01    1.50  H2O
  Halite           -2.49     -0.92    1.57  NaCl
  O2(g)           -39.20    -42.10   -2.89  O2

**For a gas, SI = log10(fugacity). Fugacity = pressure * phi / 1 atm.
  For ideal gases, phi = 1.