Processes > Inverse modelling

Modelling isotope fractionation during evaporation

(1/2) > >>

evansmanu:
Dear all,

I encountered a new concept that I am unfamiliar with in phreeqc simulations. My idea is that I have surface water isotope ratio data on O-18 (-0.47) and deuterium (2.69), which I would want to numerically assess the influence of evaporation on the fractionation of these two stable isotopes. I plan to start with the O-18 isotope ratios of the ocean vapour (δ18O vapour = -10.)  and simulate the fractionation during evaporation at different rates. I looked at a similar process that considered fractionation during the cooling of water vapour but could not comprehend it.

I would appreciate it if anyone could share a script to perform what I have described.

Regards to all

dlparkhurst:
Let me say first that I do not know how you want to do the calculation.

Let's simplify and talk only of 18O. The PHREEQC approach treats 16O and 18O as separate elements, and uses absolute concentrations rather than ratios. There are species H2O(aq), H2[18O](aq), H2O(g), and H2[18O](g).  Further, there are Henry's law equations K(16O) = H2O/H2O(g) and  K(18O) = H2[18O]/H2[18O](g). The two Ks are similar, but not the same, and the difference explains equilibrium fractionation between gas and solution. Fractionation during the evaporation of water involves the loss of H2O(aq) and H2[18O](aq) from solution. The relative rates of the loss of these two molecules determine the fractionation process.

It is possible to calculate equilibrium between gas and liquid. The following example calculates what the solution 18O permil would be in equilbrium with a gas phase of -10 permil. First it is calculated in USER_PRINT using fractionation factors from iso.dat. The water composition is -0.72 permil. A calculation is made with 1 L of water at -0.72 permil and a small empty gas phase, which, when equilibrated, produces the expected -10 permil gas phase. This example shows how to manipulate isotope ratios, permil, and absolute concentrations of isotopes.

You could assume that the water vapor that evaporates is the equilibrium water vapor (-10 permil for a solution that is -0.72 permil). You could then calculate the evolution of a given volume of solution as the water evaporated. Or, you could specify differing rates (RATES and KINETICS) of loss of H2O and H2[18O] and calculate how the liquid and gas evolve. You will have to decide how you want to proceed.

--- Code: ---USER_PRINT
10 alpha = 10^LK_NAMED("Log_alpha_18O_H2O(g)/H2O(l)")
20 REM calculate 18O ratio for -10 permil
30 r_std = 2005.2e-6
40 REM (r_smp/r_std - 1)*1000 = permil
50 permil = -10
60 r_gas = (permil/1000+1)*r_std
70 r_soln = r_gas / alpha
80 PRINT "                R        Permil"
90 soln_permil = (r_soln/r_std - 1)*1000
100 gas_permil = (r_gas/r_std - 1)*1000
110 PRINT  "18O(aq): ", r_soln, soln_permil
120 PRINT "18O(g):  ", r_gas, gas_permil
SOLUTION 1
[18O] -0.72
END
PRINT
-user_print false
END
USE solution 1
GAS_PHASE 1
-fixed_volume
H2O(g) 0
H2[18O](g) 0
-vol 1e-3
END

--- End code ---

dlparkhurst:
Here is a KINETICS calculation that approximates equilibrium evaporation from a pure water solution.  1 mol of H2O(g) is removed per day (arbitrary) from the solution along with HDO and H2[18O] as prescribed by equilibrium between gas phase and liquid phase. Essentially, water continuously evaporates into a vacuum. The incremental addition to the gas phase (delta) and the average isotopic composition of all gas removed are plotted for D and 18O. Ultimately, the average gas phase (kinetic reactant amounts) isotopic composition must equal the initial liquid isotopic composition as the system approaches pure gas.

--- Code: ---RATES
Evap_H2O
-start
10 rate = 1 / 86400
20 moles = rate * TIME
30 SAVE -moles
40 PUT(rate, 1)
-end
Evap_HDO
-start
10 alpha = 10^LK_NAMED("Log_alpha_D_H2O(g)/H2O(l)")
20 r_d_soln = CALC_VALUE("R(D)")
30 moles = GET(1)*2*R_d_soln*alpha*TIME
40 SAVE -moles
-end
Evap_H2[18O]
-start
10 alpha = 10^LK_NAMED("Log_alpha_18O_H2O(g)/H2O(l)")
20 r_18O_soln = CALC_VALUE("R(18O)")
30 moles = GET(1)*r_18O_soln*alpha*TIME
40 SAVE -moles
-end
SOLUTION 1
pH 7 charge
[18O] -0.47
D     2.69
O(0)  1 O2(g) -0.7 10
END
KINETICS 1
-cvode
Evap_H2O
-formula  H2O  1
-m        0
-m0       0
-tol      1e-08
Evap_HDO
-formula  HDO  1
-m        0
-m0       0
-tol      1e-08
Evap_H2[18O]
-formula  H2[18O]  1
-m        0
-m0       0
-tol      1e-08
-steps       54*86400 #  54 moles of H2O removed
END
INCREMENTAL_REACTIONS
USE solution 1
USE kinetics 1
USER_GRAPH 1
-axis_titles            "Days" "D, permil" ""
-initial_solutions      false
-connect_simulations    true
-plot_concentration_vs  x
-start
5 r_std = 155.76e-6
10 GRAPH_X TOTAL_TIME / 86400
20 GRAPH_Y (CALC_VALUE("R(D)")/r_std - 1) * 1000
30 r = KIN_DELTA("Evap_HDO")/(2*KIN_DELTA("Evap_H2O"))
40 GRAPH_Y (r/r_std - 1) * 1000
50 r = (KIN("Evap_HDO"))/(KIN("Evap_H2O")*2)
60 GRAPH_Y (r/r_std - 1) * 1000
-end
USER_GRAPH 2
-axis_titles            "Days" "18O permil" ""
-initial_solutions      false
-connect_simulations    true
-plot_concentration_vs  x
-start
10 GRAPH_X TOTAL_TIME / 86400
100 r_std = 2005.2e-6
110 GRAPH_SY (CALC_VALUE("R(18O)")/r_std - 1) * 1000
130 r = KIN_DELTA("Evap_H2[18O]")/KIN_DELTA("Evap_H2O")
140 GRAPH_SY (r/r_std - 1) * 1000
150 r = (KIN("Evap_H2[18O]"))/(KIN("Evap_H2O"))
160 GRAPH_SY (r/R_std - 1) * 1000
-end
END
--- End code ---

evansmanu:
Dear David,

Thanks very much for this code; that is exactly what I was looking for. I am going to implement the code right away.

Regards

evansmanu:
Dear David,
I have come again with a few questions about your code. I have been struggling to understand some of the lines. Could you please elaborate on the below lines on your code?
1.   Numbers 20 and 30 under Evap_HDO are also the same for Evap_H2[18O]
2.   Lines 30, 50, 130 and 150 under the two USER_GRAPHS
Also, why are the -m and mo all assigned 0

Is it possible to plot the corresponding evaporation rates on the y or x-axis? for instance, the evaporation rates against the ratios (O18, D).

Regards