Isotopes in reaction-transport

[GeeqC]

I am now trying to implement stable isotopes in PhreeqcRM. Therefore, equilibrium fractionation factors are available from the newiso.dat database. However, the isotope concentrations (e.g. [13C]) do not appear in the c_species variable, which would be essential to transfer them to the transport module.

Isotope concentrations can be transferred via c_species if they are defined as decoupled elements. But this gets complicated if trying to equilibrate isotopes among the diverse aqueous species. I wonder if there might be a better way to transfer isotope concentrations to the transport module?

[dlparkhurst]

I assume by the use of c_species you are using GetSpeciesConcentrations? If so, the total [C13] will not be transferred in the returned array, rather concentration of every aqueous [C13] species will be included in the list--[13C]O2, H[13C]O3-, CaH[C13]O3+, and possibly many, many others.

If you want to retrieve the total [C13] concentration in your specified units (say, mol/L), then you must use GetConcentrations, or retrieve the value through BMI GetValue or selected output.

If you are calculating individual ion diffusion, then, in fact, you will need to define a diffusion coefficient for every [C13] aqueous species, use GetSpeciesConcentrations, and calculate the diffusion of every aqueous species--a formidable task. If you are modeling advective-dispersive transport, then you should use the total element concentrations returned  by GetConcentrations; you will have many, many fewer transport calculations.

[GeeqC]

Yes, I am indeed using GetSpeciesConcentrations to get the full list of species. But none of the isotopologues shows up.
I define the totals in the input:

oss << "SOLUTION 1\n";
...
oss << "C " << 2.4 << "\n";
oss << "[13C] " << c13dic << "\n";

Printing the input shows that the values are correctly parsed to PhreeqcRM.

A problem could be that I copied the isotope section from the newiso.dat to the phreeqc.dat database, even if I thoroughly checked that all relevant information is copied. I also tried switching the master species from [13C]O2 to [13C]O3-2 to be compatible with the phreeqc.dat database, but it didn't help. The isotopes do not show up in c_species.

[dlparkhurst]

No, you got me. That is a bug. Isotopes initially are handled differently than other elements because they are input in different units. It is necessary to do a calculation to convert the initial units (often permil) to moles. The way the lists of elements and species were created did not take that calculation into account, and the lists did not include the isotopic species.

Here is a code change for PhreeqcRM.cpp about line 2930. I have included a few lines for context, but the marked line is the one that is needed.

Code: [Select]

for (i = 0; i < components.size(); i++)
{
if (components[i] == "H") continue;
if (components[i] == "O") continue;
if (components[i] == "H2O") continue;
if (components[i] == "Charge") continue;
in << components[i] << " 1e-6\n";
}
// Added line for isotopes
in << "END; MIX; " << next << " 1.0; END\n";   /////////////////
int status = phast_iphreeqc_worker->RunString(in.str().c_str());
if (status != 0)
{
this->ErrorMessage(phast_iphreeqc_worker->GetErrorString());
throw PhreeqcRMStop();
}

I will make the correction in the next release.

[GeeqC]

Thanks very much for the update, but there still seems to be a problem.

I have added the indicated line in PhreeqcRM.cpp, then recompiled the program. The date in the finder confirms that the .dylib library was updated. That the compilation was successful is confirmed as my code runs through as before. But the isotopologues still do not appear in c_species.

Could it be that I did something wrong with the input? Or the way I copied the isotope section into phreeqc.dat?

[dlparkhurst]

Can't tell you the problem. Could be you or me.

Do you use VisualStudio? Can you run the project TestRM?

I can tell you how to modify the Species test case (C++, Fortran, or C) to include isotopes. I'm not sure where you got newiso.dat; the database iso.dat is distributed with the codes, and I would use it for testing.

[GeeqC]

The Species_cpp.cpp seems to be an all-inclusive test file. Please understand that it is difficult for me to understand this architecture which seems to refer to many functions in other files.

More helpful would be some kind of minimal code, only testing a simple input solution and returning the c_species vector. If this successfully returns the isotope values I can take it from there to search for the error.

P.S. To answer your questions:
- I use Xcode on mac for compilation.
- Yes, I also found the iso.dat database. I only copied the data relevant for carbon isotopes into the modified phreeqc.dat database with redox decoupled gases. It is possible that some error occurred in this process.

[dlparkhurst]

I think you are on your own. The test cases are a known quantity that we could use to debug, but I don't use Macs and I could not guide you through the couple of changes necessary to test isotopes. It would require changing the database to iso.dat, and add isotopes to advect.pqi.

I'm assuming you have some PhreeqcRM code that functions. All you need to test isotopes is an input file that includes a SOLUTION with C and [13C], and I would use the iso.dat as the database to avoid any mistakes you may have made in copying.  You need to link with the library compiled with the change I posted. You must use RunFile (or RunString) to define the building blocks for your system, then FindComponents and then GetSpeciesNames. The return from GetSpeciesNames should have isotopic aqueous species of [13C].

[GeeqC]

You were right. I missed the whole speciation section of the iso.dat database. Now it works, in principle, i.e. c_species returns the isotope values.

Comparing phreeqc.dat and iso.dat I noticed some fundamental differences: the phreeqc.dat contains much more data, most importantly the enthalpies and molar volumes. Also, the master species are sometimes defined differently (e.g. H+ vs. H3O+). Now I have a formidable task to bring the information together and, in addition, decouple the redox states.

[dlparkhurst]

I strongly urge you to try the thermodynamic approach to kinetics first.

You apparently have the non-isotopic decoupled calculations. Are the calculations adding CH2O kinetically that different than those calculations? At least work through the non-isotopic calculations before you launch into a major effort with isotopes that may not be worth it.

[GeeqC]

I can see that this raises several follow-up questions.

The main problem of this thread, that isotope values did not appear in the c_species vector, has been resolved. For clarity, it may be better discuss any further question in a different thread.