Multiple questions regarding data entry for a surface complexation problem

Hi there,

I am working on a coal combustion product (CCP) leachate project. The goal is to understand the natural attenuation and reversibility for some remediation drivers at CCP sites, including but not limited to As, Sb, Co, Mo, Cr, and Se.


I have two sets of chemistry data, one for the background water and one for the leachate. I wonder if anyone can help with the following questions (I am not a geochemist, so my apologies if they sound basic):


1-      I need to have background water as “initial water” and the leachate as “boundary water” for my simulation. In the background water, the HCO3- concentration is known, while in the leachate, CO3-- is known. The Toughreact (TR) thermodynamic database only has HCO3- as a primary specie (I know I can switch CO3-- with HCO3- in the database using KSWITCH). Hence, I can’t include the CO3-- concentration for the leachate. How can I include the alkalinity for the leachate when I don’t know the HCO3- concentration? And, what if I have both carbonate and bicarbonate concentrations?

2-      The redox value for both background water and leachate are known. We also have the O2 concentration for the leachate, but the O2 concentration for the background water is unknown. In TR, it’s recommended to use O2 concentration for redox. When I convert redox to O2 concentration for the leachate, the resulting value is not close to the leachate known O2 concentration. Hence, I am not sure if using the same methodology is valid for the background water. What is the best approach to include the redox as O2 concentration in TR?

3-      The equivalent for some chemical constituents is not available in TR thermodynamic database. For example, I have the Boron concentration, but in the TR database, we only have b(oh)3(aq) as a primary specie. Should I use PHREEQC to find the fractionation of each chemical constituents and then only use the concentration for the ones that are available in the TR thermodynamic database?

4-      For some species, TR thermodynamic database has more than one oxidation or reduction primary species. For example, for Selenium, the database has HSe-, SeO3--, and SeO4-- as primary species. From the leachate chemistry, I know the Selenium concentration, but I don’t have the concentration of relevant Se constituents in the TR thermodynamic database. Like #3, should I use PHREEQC first, and then use the resulting concentrations for the primary species of interest in TR initial and boundary waters? I think the same condition applies to Fe vs Fe++ and Fe+++.

5-      I need to include the surface complexes in my simulation too. Like sample problem 10, I need to include hfo_soh and hfo_woh as primary species. I usually run a 1-cell chemical initialization for a significant time (e.g., 1 million year) to reach pseudo steady state, and then use the resulting water chemistry and mineral fractions in the big mesh simulation. So, should I use the recommended “1.0E-07, 1.0E+00” for cguess and ctot, respectively, in my chemical initialization, and then use the resulting values for hfo_soh and hfo_woh for the big mesh simulation?

6-      Following #5, what’s the best cguess and ctot values for hfo_soh and hfo_woh for the boundary water? From a trial and error, TR runs faster when I use the resulting values of hfo_soh and hfo_woh from chemical initialization of the initial water.


Thanks for your help and time!

Babak (aka Bobby)


Background (initial) water

    temp      22

    pH        6.8317

    redox     124.3

    units     mg/l

    Al        0.0001

    HCO3-     54.15

    B         0.0328

    Ba        0.0163

    Ca        16.8357

    Cl        12.9073

    Cr        0.00137

    Cu        0.001182

    Fe        0.1784

    K         0.7212

    Mg        4.7639

    Na        8.6965

    SO4       8.1399

    Se        0.00229

    Si        19


Leachate (boundary) water

    temp      30.6

    pH        11.65

    redox     124.48

    units     mg/l

    Ag        0.0001

    Al        44

    CO3--     26.71

    As        0.06864

    B         109

    Ba        0.14

    Be        0.0005

    Ca        15.2853

    Cd        0.06469

    Cl        54.744

    Co        0.00152

    Cr        0.00376

    Cu        0.0206

    Fe        1.53

    K         228.702

    Li        0.0599

    Mg        19.742

    Mn        0.008

    Mo        39.6

    Na        731.445

    Ni        0.1283

    O2        2.9%

    Pb        0.00458

    SO4       910.438

    Sb        0.002354

    Se        0.193

    Si        19

    Sr        1.2

    Tl        0.000287

    U         0.000677

    Zn        0.1305

5 replies

    • Mikey_Hannon
    • 1 yr ago
    • Reported - view

    Eric Sonnenthal, any chance we can get your help with this?

    • TOUGHREACT Developer
    • Eric_Sonnenthal
    • 1 yr ago
    • Reported - view

    Hi Babak,


    I'll try to answer some of these now, and get back to you on the rest.

    1. When setting up the chemistry, the concentration you input is the total concentration, rather than the species concentration, so it includes all carbonate species represented as either HCO3- or CO3-2 (for example). So you should make sure you are injecting the correct total concentration. To make sure you have the correct CO3-2 concentration when using the total as HCO3-2, just run a 1-grid block speciation calculation and look at the actual concentrations in chdump.out. If you have the correct pH, temperature, etc. then all the individual species concentrations should be what you expect. Even if CO3-2 is dominant, the total is still given as HCO3-2, but in chdump.out you will see CO3-2 with a higher molality than HCO3-. You can also write out the individual species molalities in addition to the total concentrations by adding those species in the solute.inp.

    2. For an unknown O2(aq), it is usually best to fix it to something you know, like the atmosphere if the leachate was in equilibrium with that, or a mineral (e.g., hematite) if you expect it to be in  equilibrium. You may have to make some repeated estimates for the guess, since it could vary by 20 orders of magnitude or more and be difficult to converge initially. If the system is highly reducing, it may be better to use HS- and SO4-2 as the redox couple. Sometimes, though, parts of the system are in equilibrium with the atmosphere and parts highly reducing, so in that case I usually use O2(aq).

    more later,

    cheers, Eric

      • Mikey_Hannon
      • 1 yr ago
      • Reported - view

      Thanks Eric Sonnenthal , this was very helpful! We look forward to hearing the rest of your responses.

    • TOUGHREACT Developer
    • Eric_Sonnenthal
    • 1 yr ago
    • Reported - view

    Hi Bobby,

    Regarding questions 3 onwards:

    3. As for carbonate species, B(OH)3 is entered as a total concentration, so you just need to convert B concentration to B(OH3) concentration. In moles there is one mole of B in one mole of B(OH)3, so the molality is the same. Note that the actual concentration of the species B(OH)3 will be calculated during the speciation along with any other B-species, and may be much lower than the total concentration expressed in terms of B(OH)3.

    4. You don't need to know the individual concentrations of the species. TOUGHREACT will calculate that during the speciation of the initial water, similar to PHREEQC, and it will be written out in chump.out, as well as in the the other output files if you specify those species in solute.inp. You just need to enter the total concentration of Se in the form of the primary species the database is expressed in.

    5 & 6. I'm not exactly sure what you are asking here, but yes, running the chemical initialization (speciation) first and using that as input is the right approach and will help in the convergence. If you run a pseudo-steady-state with minerals, then that water will be closer to equilibrium with that mineral assemblage, and you can use the output in chump.out (in the form to copy into chemical.inp) as a new water, unless you use the inchem which will then ignore chemical.inp initial water compositions (but not injection waters). I haven't used surface complexes in my own simulations for many years (other than the test problems) since they are not so important in many geothermal problems, so the test problem set-up is probably a good way to start. You may want to test some problems others have run in the literature either with TOUGHREACT or PHREEQC just to make sure the results you get on a single block make sense. You can also look at some of the papers by Sengor, Spycher, and Zheng who have done a lot of previous work using surface complexation models in environmental and nuclear waste applications.

    • Babak_shabani
    • 1 yr ago
    • Reported - view

    Hi Eric Sonnenthal ,

    Thanks a lot! This was incredibly helpful.

    I have one more question and appreciate it if you clear it for me.

    I have a 2D mesh with 102 x 10 (X, Z) cells. I run a leachate injection simulation for 40 years, and then run 100 years post-injection without injection in restart mode. There is a natural flow from upgradient (left boundary) to downgradient (right boundary). I want to do a sensitivity analysis for post-injection on some parameters such as pH and Eh. The idea is to change these parameters for upgradient (left boundary) and see their effect on sorption. Because the post-injection needs to be run in restart mode, I need to change the relevant elements parameters in "inchem" manually, which could be time consuming and a source of error. I was wondering if there is an easier way to change primary values of some elements in a restart run?

    Thanks again for your time and support.



Content aside

  • 1 yr agoLast active
  • 5Replies
  • 91Views
  • 3 Following