Modelling isotope fractionation during evaporation

[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: [Select]

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


[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: [Select]

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
-bad_step_max 500
END
INCREMENTAL_REACTIONS
USE solution 1
USE kinetics 1
USER_GRAPH 1
    -headings               day Solution_D Delta_kin_D Kin_D
    -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
    -headings               day Solution_18O Delta_kin_18O Kin_18O
    -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

[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

[dlparkhurst]

Line 20 in the two RATES definition are not the same. In Evap_HDO the ratio of D (2H) to 1H is calculated, while in Evap_H[18O], the ratio of 18O to 16O is calculated.

In line 30, GET(1) is the rate of H2O evaporation and alpha is the equilibrium fractionation for D to H, and 18O to O respectively. So, in Evap_HDO, the rate of loss of HDO relative to H2O is calculated by multiplying by 2 (there are two H in H2O that are lost, so multiplying by 2 arrives at the total atoms of H that are lost), R_D_soln is the current ratio of D to H in solution, and alpha accounts for the fractionation between solution and gas. For Evap_18O, the factors are analogous, although there is only one O in H2O (so no factor of 2).

Code: [Select]

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
Lines 30 calculates the ratio of the change in D to the change in H over the time step, which in turn is used to calculate the permil of this current increment in the kinetic reaction. Line 50 considers the integrated amount of the kinetic reaction, that is the total amount of D and H that have been evaporated since the beginning of the calculation. Lines 130 and 150 are analogous for 18O.

M and M0 are zero for convenience using the sign conventions for kinetic reactions. If the SAVEd amount is negative, then M increases. The product of the SAVEd amount and the coefficient in KINETICS (HDO 1) adds the reaction to solution if positive, and removes from solution if negative, so HDO is removed from solution. So, M starts from zero and increase; M is always the total amount of HDO that has evaporated up to the current point in the integration.

The evaporation rates, averaged over the last time step, are KIN_DELTA("Evap_HDO")/KIN_TIME and KIN_DELTA("Evap_[18O]")/KIN_DELTA.

Code: [Select]

USER_GRAPH 3
    -headings               permil Rate_HDO
    -axis_titles            "D solution, permil" "Evaporation HDO, mol/day" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
5 r_std = 155.76e-6               
20 GRAPH_X (CALC_VALUE("R(D)")/r_std - 1) * 1000
40 GRAPH_Y KIN_DELTA("Evap_HDO")/KIN_TIME * 86400
  -end


[evansmanu]

Dear David,

Thanks very much for your elaborate explanation. Now I understand the line of codes.

However, in line 30 under Evap_HDO, you wrote 'R_d_soln' I think this is a mistake because it should rather be 'r_d_soln' which is taken from line 20.

 

[dlparkhurst]

The Basic implemented in PHREEQC is case independent.