Tuning diffusion coefficients for water isotopes
Hi Eric Sonnenthal
I have a few questions about modeling diffusion for water isotopologues in an evaporating soil. Basically, the setup (following Singleton et al., 2004) has a static block at the bottom to act as the water table (persistent water source), a set of soil blocks above it, and another static block at the top to act as the atmosphere (with relative humidity set by air mass fraction in single-phase conditions). After spinning the model up under constant rainfall, I turn off rainfall and just let the system evolve under pure evaporation conditions.
Gas diffusion is an important part of the transport here. The TOUGHREACT formulation for calculating gas species diffusion coefficients is a great approximation, but not precise enough for isotope work and I’m wondering how to appropriately tune the model to best represent reality.
1. Is tuning diffusion coefficients via the molecular diameter a valid approach?
For example, my first approach was to tune the diffusion coefficients (D) of the water vapor isotopologues by varying the molecular diameter. I forced H2-16O to match the diffusion coefficient of water (D_water = D_H2-16O). Then, for other isotopologues, like H2-18O, I used an experimentally-derived diffusion fractionation factor (alpha-diffusion = D_H2-18O/ D_H2-16O) to calculate the diffusion coefficient. Is this a valid approach?
2. How accurate should I expect diffusion transport to be? Is it a valid approach to tune the diffusion coefficient to make model output match observations?
The output I’m after from this model is the distribution of d18O-soil water with depth (see figure). From theory and observations, I expect that the maximum value d18O-soil water produced in the soil profile should be functionally independent of all model conditions (e.g., porosity, permeability, tortuosity, temperature) except atmospheric relative humidity. My simulation correctly produces the qualitative shape of the d18O-soil water profile and correctly produces a static maximum d18O-soil water value regardless of soil conditions.
However, the quantitative d18O-soil water values are incorrect - the maximum values of d18O-soil water are too small. My model setup matches Singleton et al. (2004), so I’m confident the inputs are correct. I’m wondering if this might instead point to an issue with how diffusive transport is occurring for isotopologues with similar diffusion coefficients.
I can force the model to match theoretical expectations by tuning the diffusion coefficients, but this leads to diffusion coefficients that are “significantly” different from reality. I mean “significantly” in the context of isotope space as changes to the diffusion coefficient values are occurring in the first decimal place – in other words, the diffusion coefficients are still realistic for water vapor in air but are completely unphysical for isotope chemistry (i.e., well outside experimental/theoretical estimates). Assuming all else in the model setup is correctly, is it valid to take this approach?
Attached models
I’ve attached two setups that show the difference between the Question1 and Question2 approaches. They each have a ‘flow only spin up’ and ‘soil evap simulation’ folder for the model spin up and evaporation-only run, respectively. The only difference between them is how water isotopologue coefficients are tuned.
Question1 approach – Diffusion coefficients are tuned so that the relative difference between isotopologues follows experimental observations (Merlivat, 1978).
Question2 approach – Diffusion coefficients are tuned so that the d18O-soil water output matches theoretical expectations. All coefficients are within reason for water vapor, but the relative difference between isotopologue coefficients is allowed to vary freely (and unphysically) to be whatever is needed to make d18O-soil water correct.
11 replies
-
Hi Tyler,
Very interesting problem! Thanks for sending the input files. I will take a look at the cases. I think the first thing is that the molecular diffusion coefficients for the isotopologues in the gas phase are correct. I do not think the diffusion coefficients should be tuned. One should always start with the correct molecular diffusivities, and then work out the why there are different water values. Then there is the equilibration between water and gas phase for the species (correct log K values). I'm not sure why you say that the d18O in the soil water should be independent of porosity, tortuosity, and liquid saturation? Diffusion in partially saturated media and the equilibration is dependent strongly on the tortuosity (e.g., Millington-Quirk, linear, etc.) and liquid saturation relationship. Maybe I didn't understand your point here. In any case I will try running the cases. Recently, I've done several cases of matching gas diffusion (O2 and CO2) and reactive-transport for acid-rock-drainage benchmark problems (Mayer et al., 2015) and they were the same as other reactive transport codes (PFlotran, MIN3P, CRUNCH), so I am confident gas-phase diffusion and water-gas equilibration is mathematically correct in TOUGHREACT. However, isotopes sometimes require much more stringent convergence criteria, timestepping, etc., and small errors can build up, especially once ratios are calculated, so I'm always skeptical of my runs until I've done a lot of testing. I'll test this on our latest version 4, since the gas species transport and equilibration is more accurate than in 3.0, however I don't think it is different enough to account for the discrepancies.
cheers,
Eric
-
Hi Tyler,
I forgot to answer your question about the approach to match diffusion coefficients by modifying the molecular diameter. This seems perfectly valid. Fyi, I ran V4 (also, testing it with smaller time steps and a smaller tolch) and it gives the same results you obtained with V3. I also ran it using Millington-Quirk and the max del18O is a little higher and deeper, but not as high as you had in Q2 (plot attached). Regarding the expected del18O value in water -- what are you expecting to get? I think I found the issue (which I think I figured out and then forgot many years ago!). The del18O values for evaporation are not correct using this approach without making some modifications. If you look at the top soil grid block, the liquid saturation decreases by a factor of 0.25245. However, the sum of the water isotopes decreased by only a factor of 0.44738, because they are treated as trace species and were concentrated in the solution, where water cannot become more concentrated. So this method works as long as there is no evaporation. It should be possible by modifying the log K's, but a change in the transport equation for water isotopes might be another approach. I'll think more about what is possible without a code modification.
Eric
-
I think the first thing is that the molecular diffusion coefficients for the isotopologues in the gas phase are correct. I do not think the diffusion coefficients should be tuned. One should always start with the correct molecular diffusivities, and then work out the why there are different water values. Then there is the equilibration between water and gas phase for the species (correct log K values).
I think I am with you here but didn’t make my point clearly. I agree that the molecular diffusivity of water vapor (H2-16O) in air as calculated within TOUGHREACT gives a good value – and the value for H2-18O is also, to first-order, correct. What I’m wondering about is the relative values of the isotopologues. Experimental work suggests that the ratio of diffusion coefficients for H2-16O and H2-18O should be around 0.9723 (as H2-18O/H2-16O; Merlivat, 1978) and the TOUGHREACT formulation gives 0.9486. This is a small absolute difference, but very meaningful on the per mil level of natural variability in isotope ratios.
So if we want TOUGHREACT to provide accurate isotopologue diffusion, shouldn’t the relative difference between the diffusion coefficients be required to match observations? I’m thinking here in analogy to how we can use observed the saturation vapor pressure of H2-16O and the equilibrium fractionation factor for 18O/16O to calculate log K values that have the correct relative differences (which will be maintained even if the absolute log K values may be slightly incorrect).
I'm not sure why you say that the d18O in the soil water should be independent of porosity, tortuosity, and liquid saturation? Diffusion in partially saturated media and the equilibration is dependent strongly on the tortuosity (e.g., Millington-Quirk, linear, etc.) and liquid saturation relationship.
What I’m referring to is the maximum value of d18O that develops in the soil, nominally at the evaporation front. The soil conditions will greatly affect the shape of the profile, like in Fig. 3 of Singleton et al. (2004) (attached). However, regardless of the overall d18O vs. depth profile, the maximum value of d18O developed in any soil (under quasi-steady state) will be controlled primarily by atmospheric relative humidity – all three soils from sandy-to-clayey produce the same maximum value of d18O-soil water. This is why I was thinking of using the “maximum d18O-soil water” as a metric for tuning the relative values of the diffusion coefficients. It’s a reasonably well described
Recently, I've done several cases of matching gas diffusion (O2 and CO2) and reactive-transport for acid-rock-drainage benchmark problems (Mayer et al., 2015) and they were the same as other reactive transport codes (PFlotran, MIN3P, CRUNCH), so I am confident gas-phase diffusion and water-gas equilibration is mathematically correct in TOUGHREACT. However, isotopes sometimes require much more stringent convergence criteria, timestepping, etc., and small errors can build up, especially once ratios are calculated, so I'm always skeptical of my runs until I've done a lot of testing. I'll test this on our latest version 4, since the gas species transport and equilibration is more accurate than in 3.0, however I don't think it is different enough to account for the discrepancies.
Thanks for the reference – it’s great to hear that the codes all converge on similar results. My background is just experimental+observational work, so I’m obviously relying on those for benchmarks to make sure I have the simulation set up appropriately. I haven’t been able to produce noticeably different results by changing the convergence criteria, but I will try reducing the timestepping (I’ve been using 100-1000s seconds as the maximum).
-
Hi Tyler,
Regarding the factor, I was comparing the saturation in the top soil block at time=0 to time=1 year, and then the same for the sum of the isotopes. Thinking about it a little more, the total concentration of water isotopes (H216O+H217O+H218O) should not change during evaporation, unless there is a very strong increase in salinity which is not a factor in this problem. Basically, the concentration of water in water needs to stay the same. The approach of treating them as trace species and equilibrating the gas and liquid phases is fine as long as there isn't evaporation. I haven't had a chance to go through the equations to see how this can be done. This might only be possible by modifying the transport equation which is where the effect of saturation changes owing to evaporation results in a concentration increase. I'll try just hardwiring the code for those isotopes and see what happens.
Yes, I see your point about the diffusion coefficients calculated directly using the simple model. If the molecular diameters need to be modified to get the correct fractionation factor, that is the way to go. The model in TOUGHREACT (from Lasagna) is really just ok for gas species having significantly different diameters.
I think your timestepping is fine -- I only saw very slight differences at a smaller time step.
More later...
cheers,
Eric
-
Hi Tyler,
I set up a modified 2-block problem that only involves evaporation of the top soil block and diffusion into the atmosphere fixed at 0.4 RH. The grid and hydrologic properties are a little different from yours, but the process is the same. What you see here is that the top block dries to the residual liquid saturation (Slr=0.01) after 0.9 years, and once down to the Slr the RH drops to 0.4. It is dropping slightly before then which is driving evaporation. Summing up the isotopes, the sum drops more smoothly but also to 0.4. So water exists (because of vapor pressure lowering) even at an RH of less than 1.0. Basically the sum of the isotopes is reflecting their activity in equilibrium with a lower vapor pressure (RH= 0.4). As far as the water activity that makes sense, however, the total concentration of water isotopes should still be 1 mol/kgH2O. So I think the mass balance for the isotopes is likely off when there is strong evaporation at an RH less than 1.0, even though the equilibrium partitioning is correct at first. Anyways, I think I'll need to work on this some more. Probably the answer may be just using the TOUGH-determined RH as a fugacity coefficient, but this will take a little more working through the equations and the code. Some of these effects are all lumped into kinetic fractionation factors, etc., however because there is transport of water and gas as well as changing liquid saturations, the mass balance is an issue. I also add a plot showing the isotopic ratios over time for this case where the liquid isotopic transport is skipped and another case where it isn't (this one). For this problem the only difference is that the change in liquid saturation during dryout doesn't increase the isotopic concentration. I'm not yet sure how to deal with that issue, but I think maybe it is ok, because the equilibrium constant should account for it. I'll know once I modify the water isotope fugacities.
Eric
-
Hi Tyler,
I haven't had a chance to get back to this issue, but it is really important to figure out. I started putting in some code specific for water isotopes, since those are not like other trace isotopes, but had to get back to some other things. I think this may have to be done through the coupling parameters like we do for water or CO2 in minerals, and the mass balance in the multiphase flow. I will start looking at it soon, since it would be a very useful capabilities for solving evapotranspiration problems.
regards,
Eric
