Processes > Inverse modelling

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: ---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

-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

--- 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

Navigation

[0] Message Index

[#] Next page

Go to full version