Unexplained isotope fractionation effects

[GeeqC]

While in the previous thread "Modifying database for isotope speciation" we could clarify the different outputs for the isotopologues, I have now started to evaluate the values. Fractionation between carbonate species (CO2, HCO3-, CO3-2) in the initial solution appeared unrealistically large (values ranging between -200 and +150 permille PDB).

As a test, I tried to set in the database all fractionation factors to: -ln_alpha1000   0.0
Activity parameters are only set for the species. For the isotopologues they should be implicitly defined via the fractionation factors.

Input is:

SOLUTION 1
    ...
    C  2.4
    [13C]  0

This should then result in zero fractionation, but the large fractionation effects still occur.
I do not understand what is happening.

[dlparkhurst]

Sorry, you'll have to provide a simple example--perhaps a NaHCO3 solution. Please use iso.dat and include any changes to the database in the input file.

[GeeqC]

Here is a basic example, which I ran in regular Phreeqc using iso.dat:

SOLUTION 1
    -temp     10
    -units    mmol/l
Alkalinity 2.3
C(4)  2.4
[13C]  0.0271  #should yield a d13C around 10 permille
Ca  10.3
Na  470
Cl  520

SELECTED_OUTPUT 1
    -totals       C(4)  [13C]
    -isotopes   R(13C)_CO2(aq) R(13C)_HCO3-

I also tested the input "-isotope  13C  10" but according to the user manual, this is only used for inverse modelling.

In the output, [13C] is zero, and the R's show some default value.

[dlparkhurst]

I'm still not sure what you are getting at. Here is a script that removes fractionation for 13C in a simple NaHCO3 solution.

You could also simply redefine a couple species as in the SOLUTION_SPECIES below. The NAMED_EXPRESSION is the better way to do it.

Code: [Select]

NAMED_EXPRESSIONS
Log_alpha_13C_CO2(aq)/CO2(g) # 1000ln(alpha(25C)) -0.84
Log_alpha_13C_CO2(g)/CO2(aq) # 1000ln(alpha(25C)) = 0.84
Log_alpha_13C_HCO3-/CO2(g)   # 1000ln(alpha(25C)) = 7.82
Log_alpha_13C_HCO3-/CO2(aq)  # 1000ln(alpha(25C)) = 8.7
Log_alpha_13C_CO3-2/CO2(g)   # 1000ln(alpha(25C))
Log_alpha_13C_CO3-2/CO2(aq)  # 1000ln(alpha(25C))

SOLUTION_SPECIES
#HCO3- + [13C]O2 = H[13C]O3- + CO2
#CO3-2 + [13C]O2 = [13C]O3-2 + CO2

SOLUTION
pH 7
Na 1 charge
C(4) 1
[13C](4) 10
END
RUN_CELLS
-cells 1
END

[GeeqC]

Attempting to provide a simple example, I have switched back to regular Phreeqc (not PhreeqcRM), trying to calculate d13C in a single solution.

As in your example, I have now defined:
C(4)  1
[13C](4)  10

At this point, I encounter several issues:

1.) The simulation does not return any output for the isotopes (using iso.dat with unmodified fractionation factors).

2.) Are you writing the NAMED_EXPRESSIONS directly into the input file? Do these expressions then overwrite the ones in the database?

3.) In your latest example, the solution species commented out are:
        #HCO3- + [13C]O2 = H[13C]O3- + CO2
        #CO3-2 + [13C]O2 = [13C]O3-2 + CO2
     But would this not mean that the isotopologues are not calculated?
     I would need the isotopologues calculated, just with zero fractionation.

[dlparkhurst]

1. The initial solution calculation does not include isotopologues; it simply calculates the amount of 13C in solution. You need to then use the solution in a reaction calculation to see the isotopologues, which is why I included the RUN_CELLS data block.

2. The NAMED_EXPRESSION data block in the input file overwrites any definitions from the database file.

3. Did you even run the file that I gave you? You will see in the RUN_CELLS calculation that there is no fractionation. By commenting out the SOLUTION_SPECIES, there is no change to the database for these species. However, the NAMED_EXPRESSIONS overwrite the database definitions, which provides zero fractionation.

If you were to use the SOLUTION_SPECIES definitions that are commented, and remove the NAMED_EXPRESSIONS, you would also see results with no fractionation. However, there would be some loose ends in the printout that still use the NAMED_EXPRESSIONS, so it is better to remove fractionation by redefining the NAMED_EXPRESSIONS.

[GeeqC]

Thank you for the explanations. The example now ran through successfully:

1.) The problem was indeed that RUN_CELLS was missing.

2.) Isotopologue concentrations and delta values are calculated correctly for zero fractionation.

3.) If "NAMED_EXPRESSIONS" are switched off, delta values are in the right range: d13C(CO2)= 1.25; d13C(HCO3-)=11.32, and d13C(CO3-2)=9.58

There was, however, a further complication: When defining Alkalinity = 2.3 and C(4) = 2.4, the iso.dat database defaults to pH7. But the pH should be calculated based on alkalinity and ionic composition. This is done correctly in the phreeqc.dat database.

[dlparkhurst]

>There was, however, a further complication: When defining Alkalinity = 2.3 and C(4) = 2.4, the iso.dat database defaults to pH7. But the pH should be calculated based on alkalinity and ionic composition. This is done correctly in the phreeqc.dat database.

Nah. Both phreeqc.dat and iso.dat give a pH of 7.669.

Of the quantities pH, C(4), and Alkalinity, you can only define two; the third will be a calculated value.

Code: [Select]

SOLUTION
Alkalinity 2.3
C(4) 2.4
END

[GeeqC]

Aha, this seems to be a problem of the iso.dat database in phreeqc-3.5.0-14000 for mac. If I use the iso.dat file from PhreeqcRM the pH is calculated correctly.

As the isotope calculation works in Phreeqc, I will now try to fix the problems in my cpp code. It seems the overall sequence is wrong:
- RunString
- Set initial conditions (ic1)
- Load initial conditions to the module: InitialPhreeqc2Module(ic1)
- Populate boundary condition vector: InitialPhreeqc2SpeciesConcentrations(bc_conc, bc1)
- RunCells

According to your previous comment, RunCells must come first. However, in order to run the cells, I probably would need already the initial conditions?

[dlparkhurst]

The "advect" examples in the directory Tests/ show the sequence of calls for using PhreeqcRM.

[GeeqC]

The problem with the example Advect_cpp.cpp is that it requires several other files to be included, and it contains a lot of stuff that I do not understand. Nevertheless, I was now able to run it in its original form.

Having managed this, I changed all references to the iso.dat database, also in the included files, except in Gas_cpp.cpp.

Then it was ready to run an example input file. I used the example advect.pqi and changed SOLUTION 0 as follows:

SOLUTION 0
        units  mmol/kgw
        temp  10
        Alkalinity  2.3
        C(4)  2.4
        [13C](4)  10
        Ca  10.3
        Mg  50
        Na  470
        Cl  520
        S(6)  28
END

RUN_CELLS
        cells  0
END

Isotopologues are calculated and appear in the right range. But the concentrations do not match the delta values given in the selected output. The concentrations of the totals in the selected output are truncated to the fourth digit, so that I could not calculate their isotope ratio.

[dlparkhurst]

Use SELECTED_OUTPUT; -high_precision true.

And I don't respond well to whining.

[GeeqC]

I added "-high_precision  true" to SELECTED_OUTPUT. But the values are still truncated:
           ...
           8 m_HCO3-(mol/kgw):     0.0017
           9 m_H[13C]O3-(mol/kgw):     0.0000
           ...

[dlparkhurst]

Here is a file that prints the ratios of 13C to 12C for CO2, HCO3- and CO3-2 to the selected output file with high precision. You can look at the CALCULATED_VALUES section of the iso.dat database file to see how the ratios are calculated.

I'm not sure what values you are looking at with lesser precision; I don't think they would be USER_PUNCH. In any of the Basic scripts (USER_PRINT or USER_PUNCH) you can control the format with STR_E$ and STR_F$ functions. I used STR_F$ to control the print in the USER_PRINT printout.

Code: [Select]

SOLUTION
pH 7
Na 1 charge
C(4) 1
[13C](4) 10
END
SELECTED_OUTPUT
-reset false
-high_precision true
-isotopes R(13C)_CO2(aq) R(13C)_HCO3- R(13C)_CO3-2
USER_PUNCH
-headings UP_R_13CO2 UP_R_H13CO3- UP_R_H13CO3-2
10 PUNCH CALC_VALUE("R(13C)_CO2(aq)")
20 PUNCH  CALC_VALUE("R(13C)_HCO3-")
30 PUNCH CALC_VALUE("R(13C)_CO3-2")
USER_PRINT
10 PRINT "R(CO2):   ", STR_F$(CALC_VALUE("R(13C)_CO2(aq)"), 10, 8)
20 PRINT "R(HCO3-): ", STR_F$(CALC_VALUE("R(13C)_HCO3-"), 10, 8)
30 PRINT "R(CO3-2): ", STR_F$(CALC_VALUE("R(13C)_CO3-2") , 10, 8)
RUN_CELLS
-cell 1
END

[GeeqC]

Thank you, this was insightful. As a result, I could track down several issues:

1.) Running your input file:
Sorry, I have to come back to the architecture of Advect_cpp.cpp: In order to run your file, I had to modify and embed it in advect.pqi. I am not sure I was supposed to do this, but it worked.

2.) Selected output:
The selected output -isotopes R(13C)_CO2(aq) ... always worked fine.
The problem was with the total concentrations, where the values are truncated after the fourth digit. You can see this issue if running -totals in SELECTED_OUTPUT in the default advect.pqi. As a result, also the UP's were not calculated correctly.

3.) The isotopologue concentrations
I am looking at the output in the console: Under "Species:" the lefthand column should represent the concentrations. In principle, the isotopologue concentrations yielded the correct delta values. However, for the seawater solution (given above) the delta values deviate between different ion pairs and free species. In the iso.dat database no fractionation between ion pairs is indicated (except for CaHCO3+ a gamma is defined).

[dlparkhurst]

Sorry, I am slow to get into your issues, but you have basically worn me out. I have no enthusiasm for your problems.

I was misleading in that PhreeqcRM does not automatically write any selected output files and -high_precision is thus not applicable. PhreeqcRM has methods that give you the capability to retrieve the selected output values, possibly for multiple SELECTED_OUTPUT definitions (that have unique user numbers). There are a whole series of methods to select which selected output definition is to be retrieved, and the method GetSelectedOutput that retrieves all the values. The method GetSelectedOutputHeading can be used to get the headings of the selected output file. The retrieved values are full double precision, so there is no issue with precision, but it is up to you what you want to do with the values.    The following allows retrieving the selected output values at the last time step.      

bool print_selected_output_on = (steps == nsteps - 1) ? true : false;

You can use USER_PRINT if you write out typical PHREEQC output-file data. The advection example turns on the output file for the last time step.

bool print_chemistry_on = (steps == nsteps - 1) ? true : false;

You can use USER_PRINT to do the isotope ratio calculations and it will appear in the output file. You can use MOL, TOT, CALC_VALUE, and other Basic functions in USER_PRINT to convince yourself that the calculations are correct. I am willing to help you with USER_PRINT, but not SELECTED_OUTPUT. Of course, it is much easier to use PHREEQC, rather than PhreeqcRM, to investigate the basic isotope calculations.

 

[GeeqC]

I am really sorry that I have worn you out.
It may not be necessary to further discuss the issues with SELECTED_OUTPUT. Anyway, this seems to be a problem specific for the example run Advect_cpp.cpp. Because in my own code, the selected output appears in full "double" precision.

To avoid all the output issues, we can look directly at the values from the species vector: c[j * nxyz + i]
The concentrations show fractionation between ion pairs of the same species. Why is this? No fractionation is defined for the ion pairs in the database, so that the fractionation should be zero.

If we could clarify this for the example run, I can take it from there and track the issue to my own code. I think solving these issues may be of interest for a broader community, in order to apply isotope fractionation in PhreeqcRM. Indeed, PhreeqcRM would be a powerful tool to calculate isotope reaction and transport!

[dlparkhurst]

OK, you say there is unexpected fractionation. You should be able to create a PHREEQC input file that demonstrates it. In fact, it seems to me that any SOLUTION definition with the same set of elements should show similar fractionation. So, send me a PHREEQC input file with SOLUTION and RUN_CELLS that shows fractionation and tell me explicitly which ratios of species show this fractionation. Use iso.dat for the database.

[GeeqC]

It seems I could locate the source of the problem:
The whole time I was looking at the output on the console. But there the values are formatted and therefore rounded. The inconsistencies in the isotope ratios calculated from isotopologue concentrations are due to rounding.

In my own simulation I export all values of the c_species and so vectors to a .csv file. There appears the full precision.