Question on modeling water isotopes – after Singleton et al. (2004)
I have been playing around with your water isotopes model of Singleton et al. (2004) and have a question about how fractionation factors are implemented in TOUGHREACT. I’m able to reproduce your work, so I’m confident the model setup is correct. I see how equilibrium fractionation factors are introduced via the thermodynamic database but do not follow how kinetic fractionation via diffusion is produced.
In short, my questions are:
Why do fractionation factors calculated from diffusion coefficients (which come from equation A.1 in the manual) not match the apparent values output by TOUGHREACT? Starting from equation A.1 and using H2-16O and H2-18O as the isotopologues of interest, I get 18alpha-kinetic is 1.0541, but the apparent value from the Singleton output is 1.0173.
Why do the apparent TOUGHREACT kinetic fractionation factors differ from experimental ones? The Singleton apparent 18alpha-kinetic is 1.0173 while, for example, Barkan and Luz (2007) find a value of 1.0283 in line with other studies.
My reasoning and an excel sheet with calculations are also attached:
As I understand, TOUGHREACT uses equation A.1 in the manual to calculate tracer diffusion coefficients after Lasaga (1998). Following the equations in the Singleton paper, which follows Criss (1999), the kinetic (functionally diffusion in the Singleton model) fractionation factor is calculated as:
alpha-kinetic = (D/D’)^n
Where D and D’ are the diffusion coefficients for the common and rare isotopologues (from equation A.1) and n is an exponent equal to 1 in the case of pure diffusion, as in the Singleton setup. I’ll use H2-16O and H2-18O as an example, which gives:
18alpha-kinetic = 1.0541
In turn, the equilibrium and diffusion fractionation factors can be used to calculate the overall fractionation associated with evaporation for a given relative humidity (RH). The Singleton discussion notes that RH=0 is a special case that can then be used to calculate theoretical properties of evaporating soils at non-zero relative humidity:
alpha-evaporation(RH=0) = alpha-equilibrium * alpha-kinetic
18alpha-evaporation(RH=0) = 1.0649
However, these calculated values are significantly higher than the apparent values output by Toughreact (e.g., your Fig. 4 would suggest: 18alpha-kinetic ≈ 1.0173 and 18alpha-evaporation(RH=0) ≈ 1.0278) and as found in experimental work (e.g., Barkan and Luz 2005, 2007 find: 18alpha-kinetic ≈ 1.0283 and 18alpha-evaporation(RH=0) = 1.0388).
Do you have any insight on why these discrepancies occur? I have not found any issues with tuning the model setup – in other words changing grid size, the number of digits output, or the convergence tolerances does not change the model behavior.
I'll have to try and set this up a again, since 2003 is more than a few years ago! There are a few factors that could be causing the discrepancy. Timestepping -- Make sure to run it with a very small maximum time step. Tortuosity, porosity, and gas saturation -- the final term in the transport equation is D*porosity*tortuosity*S_gas. Diffusion in the liquid phase. Equilibration during diffusion.
I'll let you know what I find.
Thanks for your reply. I just halved the max timestep (Singleton was 1000 s) and it did not change the result. I will check out the other terms as well.
Michael gave me most of his input files to modify. The attached ones are for the conditions in Figure 3 (medium sand, atmospheric RH = 40 %, gener rate of water = 0 - evap only scenario, water reservoir d18O = -16.5 per mil, atmospheric vapor d18O = -21 per mil). You'll see I've also added in H2-17O - this project is geared towards understanding triple oxygen isotopes in soil profiles, which is why I want to make sure I understand how diffusion acts.
Thanks for sending the files. It is much easier to figure this out if I can run it. I did run it with our latest unreleased version. Can you tell me which code version you used? We have released so many versions over the years, I want to see if there are any differences in the results. Also, if you can just send me the time.out file, then I can compare your results to those I am getting.
Here is the time.out file as well
flow.out gives my version as: TOUGHREACT V3.0-OMP (May 2014)
I also ran simulations with variable tortuosity, porosity, and gas saturation (i.e., sand vs. clay), which do not affect the "d18O-max" value (the maximum value of d18O that develops in an evaporation-only soil) and do not change the apparent fractionation factors. I think this makes sense as these parameters should theoretically only depend on the atmospheric relative humidity, d18O-atmosphere, and d18O-water reservoir - as you noted in the 2004 paper.
I wanted to check in and see if you had a chance to re-investigate the evaporating soil profile. Do you think I'm visualizing how TOUGHREACT implements the kinetic, or diffusion, coefficient correctly? Since the shape of the profiles are similar to what we would expect theoretically, I'm wondering if my hand calculations are using an incorrect (i.e., different from what TOUGHREACT actually sees) kinetic fractionation factor and that would explain my confusion.