Isotopes output in PhreeqcRM

[GeeqC]

For the output of isotope ratios the following keywords seem to work:

SELECTED_OUTPUT
-isotopes R(13C)_CO2(aq) R(13C)_HCO3-

Of course, the delta value can be calculated as: Delta = (R - Rstd)/Rstd * 1000

But I wonder if this is necessary, since the database already contains Rstd and delta values may already have been calculated internally. Is there a keyword to directly select delta values as selected output?

[dlparkhurst]

You should look at the results, and read the manual for the -isotopes option of SELECTED_OUTPUT.

You could also look at the Basic functions ISO and ISO_UNITS.

[GeeqC]

The User manual says:
"The -isotopes identifier allows isotopic values in the units of the isotope standard to be written to the selected output file. The isotopic values are identified by the names of CALCULATE_VALUES Basic functions that have been identified in the ISOTOPE_RATIOS data block."

If these "isotopic values" are 'in the units of the standard', I suppose they are deltas?
But what confuses me is that they are labelled as 'I_R...' in the output.

[dlparkhurst]

They are permit for 13C. You will have to deal with it.

[GeeqC]

You mean:
delta13C = (R - Rstd)/Rstd * 1000  in permille? So, this is calculated internally in Phreeqc?

This is what I need.

[dlparkhurst]

Sorry for the misspelling, permit (not permit). You can do the he calculation to check for yourself.

[GeeqC]

In the meantime, I could investigate the matter but then encountered several further issues:

1.) I recalculated the delta values (from the concentrations and the Rstd given in the database). The I_R-values in the selected output indeed seem to be delta values. But there is some unexplained difference.

2.) The isotope values are correctly set in the initial solution, also used as upper boundary condition. However, all isotopologue concentrations are zero in the bc_conc vector. As a result they also appear as zero in the output.

3.) Under CALCULATE_VALUES, isotope ratios are calculated as: ratio = total_13C/total_C
This may be fine as an approximation. But would it not be better to change the calculations to: 
ratio = total_13C/(total_C - total_13C)

[dlparkhurst]

2. I think you probably need to do a calculation in the initial phreeqc file to get the complete isotopic composition of a solution for the boundary condition. The initial solution calculation calculates the moles of 13C, but not the actual isotopic species distribution. If you add something like

Code: [Select]

END
MIX 1
2 1.0  # 2 is boundary solution
SAVE solution 2
END

then the MIX calculation will populate the entire distribution of isotopic species, and you should be able to use solution 2 for boundary conditions.

3. You are pretty much missing the entire point of the approach to isotopes. [13C] is an entirely separate element once the initial solution calculation is completed. It was convenient to use "C" to represent 12C. So, [13C] is not part of the mass balance for C; C is only 12C and [13C] is only 13C.

[GeeqC]

The MIX calculation worked. Thank you!

My question concerning the isotope ratio calculation was not accurate. Please let me re-phrase. The iso.dat database calculates the R's as in the following example:

R(13C)_CO2(aq)
     -start
10 ratio = -9999.999
20 if (TOT("[13C]") <= 0) THEN GOTO 100
30 total_13C = sum_species("[13C]O2","[13C]")
40 total_12C = sum_species("CO2","C")
50 if (total_12C <= 0) THEN GOTO 100
60 ratio = total_13C/total_12C
100 save ratio
     -end

Line 40 defines the concentration of 12C as concentration of the total species. However, this is only an approximation. The isotope ratio should be R = 13C/12C (not 13C/C). While this small difference may not matter in most cases, it may possibly amplify as part of a numerical solution.

Therefore my suggestion to calculate:
total_12C = sum_species("CO2","C") - sum_species("[13C]O2","[13C]")

[dlparkhurst]

No changes are needed. The PHREEQC calculations are correct and exact.

When you define SOLUTION such as the following:

Code: [Select]

SOLUTION 1
pH 4 charge
C 1
[13C] 10 # permil
END
MIX
1 1.0
SAVE solution 2
END

The concentration given for C is the sum of 12C, 13C, (and 14C). The amounts of each isotope are calculated (based on the standard definition in ISOTOPES) and are printed. In this case, C is 1.00000e-3 mol/kgw, 13C is 9.88834e-04 mol/kgw, and 13C is 1.11659e-05 mol/kgw. After this calculation is done, C now represents only 12C. You will note in the subsequent MIX calculation that the concentration of C is 9.88834e-04 mol/kgw, not 1.0e-3.

Sorry if it is confusing but once the initial solution calculation is done, C is exactly equal to 12C, and does not include any contribution from 13C or 14C. The programming was complicated enough without adding a whole new set of species named [12C]. The manual has the following statement in the ISOTOPES description: "Reaction calculations with isotopes are performed by assuming each isotope is a separate thermodynamic component. Thus, in addition to the principle isotope of an element, which typically is named by the standard element nomenclature (for example, C for carbon), each isotope also is defined as an element in a SOLUTION_MASTER_SPECIES data block."

[GeeqC]

In this case, it would be correct.

But I still do not understand the syntax ("CO2","C"). CO2 and C are also the whole species/element. How can phreeqc know when the principle isotope is meant and when the whole species?

And there is still the discrepancy between I_R-values in the selected output and d13C values calculated from concentrations in c_species:

R = {H[13]CO3-} / {HCO3-} - {H[13]CO3-}
d13C = (R/Rstd - 1) * 1000

Could it be a problem that I use mmol/L:

    oss << "SOLUTION 1\n";
    oss << "-units mmol/L " << "\n";
    oss << "-density 1 calc \n";
    oss << "-water 1 # kg " << "\n";

[dlparkhurst]

No, no, no. After the initial solution calculation, C is 12C in all cases; it is not the total carbon in solution; C does not include 13C or 14C. To get total carbon in solution, you would need to add TOT("C") + TOT("[13C]") + TOT("[14C]"). TOT("C") is the total 12C in solution. Sorry if it is confusing, but if [12C] were to be defined as an element, it would be (almost) exactly the same as the definitions for C. It is a convenience to use "C" to represent 12C in the isotopic calculations.

You should not subtract a contribution from 13C, HCO3 is by definition H 12C O3-.

R = {H[13]CO3-} / {HCO3-}

I'm not sure if that is your complaint, but if you give an example of a calculation that you think is incorrect, I will look at it.

[GeeqC]

I suppose this syntax is used for C as an argument, as in TOT("C")?

But if HCO3- stands under the SOLUTION_SPECIES keyword, does it still represent the whole species, right?
The same for HCO3- concentration on the module (c_species vector)?

[dlparkhurst]

HCO3- is the 12C version. H[13C]O3- is the 13C version. Total molality of bicarbonate is MOL("HCO3-") + MOL("H[13C]O3-").

TOT("C") is total 12C (after initial solution calculation).

After the initial solution calculation C represents only 12C. All 13C is represented as the "element" [13C].

In reaction calculations, whenever you see or use the element "C", think 12C.

[GeeqC]

If this is the case, this may have diverse consequences I cannot oversee right now.

For example, does this not cause problems with charge balance, if parts of the species are missing? My transport model makes sure the charge balance is maintained, but this may lead to a hypercorrection if disbalance is created in each time step.

Another issue might be the calculation of alkalinity? Would this not be affected by the missing contribution of 13C?

Even if the errors are only small, they could possibly amplify in a numerical solution.

[dlparkhurst]

[13C] is a separate element. It is just like any other element, say boron. It has alkalinity and charge associated with each of its species, in fact, just like C (which equals 12C). So, charge balance, alkalinity, mole balance, and other properties are preserved. There are no small, cumulative errors; mass balance is maintained for 12C and 13C, as well as charge balance.

The whole point to this approach to isotopes is that 13C can react slightly differently than 12C, which gives rise to fractionation. You must think of the two isotopes as two separate components, as much as Cl and Br are separate components. As a consequence, there is a complete set of species definitions for [13C] (just like any other element), and minerals become solid solutions; for example, calcite is a solid solution between CaCO3 and Ca[13C]O3. If you want surface complexation with carbon, you must define complexation for C (12C) and [13C] as two distinct reactions. Similarly, a gas phase would contain both CO2(g) and [13C]O2(g), and kinetic reactions would have to account for slightly different rates between 12C and 13C (which gives rise to fractionation); the rate of methanogenesis is different for 12C than 13C as evidenced by the fractionation factor for this process.

You must understand this approach to be doing the calculations you are attempting. If you want to use a more traditional fractionation factor approach, you must implement it yourself, or use another geochemical model.

[GeeqC]

Now I understand. H13CO3- and H12CO3- are treated like too separate species also on the module.

This means both isotopes have to be transported separately and both are included in charge balance. This is also how we calculated isotopes in previous studies. This makes things easier.

Then there is only one issue left: the I_R-value in the selected output does not match the delta calculated from isotope concentrations on the module. It would help to know how exactly the value in the selected output is calculated.

[dlparkhurst]

Please give an example. Be clear as to which values differ, but should be the same.

[dlparkhurst]

Perhaps this will help:

Code: [Select]

SOLUTION
pH 7 charge
Na 1
Ca 1
C 3
[13C] 10
END
MIX
1 1.0
SELECTED_OUTPUT
-reset false
-iso R(13C) R(13C)_HCO3-
USER_PUNCH
10 r_std = 0.0111802
110 r = CALC_VALUE("R(13C)_HCO3-")
120 r2 = [MOL("H[13C]O3-") + MOL("CaH[13C]O3+") + MOL("NaH[13C]O3")]/[MOL("HCO3-") + MOL("CaHCO3+") + MOL("NaHCO3")]
130 permil = [r/r_std - 1]*1000
150 permil2 = [r2/r_std - 1]*1000
160 print "R(13C)_HCO3-:    ", permil, permil2
210 r = CALC_VALUE("R(13C)_CO3-2")
220 r2 = [MOL("[13C]O3-2") + MOL("Ca[13C]O3") + MOL("Na[13C]O3-")]/[MOL("CO3-2") + MOL("CaCO3") + MOL("NaCO3-")]
230 permil = [r/r_std - 1]*1000
250 permil2 = [r2/r_std - 1]*1000
260 print "R(13C)_CO3-2:    ", permil, permil2
END



[GeeqC]

Sorry for being unclear. Here, I try to be more specific:

If setting in the selected output:

SELECTED_OUTPUT
...
-iso R(13C)_HCO3-

This yields a value I_R(13C)_HCO3- in the selected output.

I tried to add the user punch you sent me, but this did not change anything. The value I_R(13C)_HCO3- seems to be calculated by default if the -iso keyword is set.

The question is what is being calculated?
I tried to reproduce it by manually calculating (R/Rstd - 1)*1000 using the concentrations on the module (c_species), but the result doesn't match. I noticed that if I calculate the R including the ion pairs, the value gets closer.

In any case, it would help to know what is actually calculated in the selected output.

[dlparkhurst]

My previous post shows the calculation of R(13C)_HCO3- in a Ca, Mg, C system. What is not clear?

Yes, the HCO3- ion pairs are included. If you want something else, you can calculate it with CALCULATE _VALUES definitions.

[GeeqC]

I realize that the notation "sum_species("*H[13C]O3*","[13C]")" with asterisks in the database includes all the ion pairs. Now, the calculations match.