Craig-Gordon Isotope Fractionation During Evaporation - Delta Values Not Updati

[evansmanu]

Scientific Background and Objective

I am attempting to simulate Craig-Gordon evaporation fractionation of stable water isotopes (δ¹⁸O and δD) in a tropical river system at 28°C. The starting water represents an ocean vapour condensate (δ¹⁸O = -10.0‰, δD = -80.0‰ VSMOW) and the target is a measured river water composition (δ¹⁸O = -0.47‰, δD = +2.69‰ VSMOW). The goal is to run four relative humidity scenarios (h = 0.85, 0.70, 0.55, 0.40) and identify which scenario best reproduces the measured river composition, thereby constraining the evaporative conditions of the system.

The Problem
I have tried three different PHREEQC implementations using iso.dat and none of them correctly updates the delta values of the residual liquid. I am attaching the full input file. Here is a summary of each approach and what went wrong:

Approach 1 - KINETICS with CALC_VALUE
Used three parallel KINETICS reactions (Evap_H2O, Evap_HDO, Evap_H2[18O]) with rate equations that compute α_eff and remove isotopologues proportional to CALC_VALUE("R(D)") and CALC_VALUE("R(18O)"). Result: completely flat lines - CALC_VALUE("Delta_D") and CALC_VALUE("Delta_[18O]") did not change at all over 30 days regardless of humidity scenario.

Approach 2 - MIX with negative fractions

Pre-computed the Rayleigh trajectory analytically, created SOLUTION blocks for each vapour end-member, and used MIX with a negative fraction to remove water at each step, followed by SAVE solution + END. Result: also flat - a single point at 0.0‰ was plotted. CALC_VALUE("d18O_val") and CALC_VALUE("dD_val") returned 0 rather than the actual delta values after MIX, suggesting MIX does not propagate isotope tracers through mass balance.

Approach 3 - REACTION with ISOTOPE_RATIOS

Used REACTION to remove H2O and ISOTOPE_RATIOS to specify the delta value of the removed water in permil VSMOW (e.g. D = -157.60 H2O, [18O] = -21.25 H2O). Result: the following errors terminated the run:

ERROR: For ISOTOPE_RATIO D, did not find ISOTOPE definition for this isotope
ERROR: For ISOTOPE_RATIO D, did not find SOLUTION_MASTER_SPECIES for isotope
ERROR: For ISOTOPE_RATIOS D, did not find corresponding CALCULATE_VALUE definition
ERROR: For ISOTOPE_RATIO [18O], did not find ISOTOPE definition for this isotope

Questions

1. What is the correct PHREEQC/iso.dat method to implement a stepwise batch evaporation model where the vapour leaving the solution has a different isotopic composition from the residual liquid, and where CALC_VALUE("Delta_D") and CALC_VALUE("Delta_[18O]") correctly reflect the enriched residual liquid after each step?
2. Why does the KINETICS approach not update delta values? Is PHREEQC overwriting the isotope tracer concentrations from speciation after each kinetic step, making it impossible to drive fractionation this way?
3. Is ISOTOPE_RATIOS the appropriate keyword for this problem? If so, what is the correct syntax - does it accept permil delta values, absolute R values, or something else? The manual examples (e.g. Example 17) do not appear to demonstrate this use case directly.

Any guidance or a minimal working example would be greatly appreciated.

Code: [Select]

TITLE Craig-Gordon river water evaporation - 28C, four humidity scenarios
TITLE Starting water: d18O = -10 permil, dD = -80 permil (ocean vapour)
TITLE Target river water: d18O = -0.47 permil, dD = +2.69 permil
TITLE Scenarios: h = 0.85, 0.70, 0.55, 0.40
TITLE Reference: Horita and Wessolowski (1994), Merlivat and Jouzel (1979)

# ================================================================
# EXPLANATION OF THE CORRECTED RATE FORMULATION
# ================================================================
# The previous version showed flat lines because removing HDO and
# H2[18O] in proportion to R(D) and R(18O) does NOT change the
# delta value -- it just scales with the bulk removal.
#
# The key insight: isotopic ENRICHMENT of the residual liquid only
# occurs because the vapour removes the HEAVY isotopologue at a
# rate FASTER than the bulk ratio by the factor alpha_eff.
#
# The Evap_H2O reaction already removes bulk water (and implicitly
# carries away isotopes at the natural ratio). The Evap_HDO and
# Evap_H2_18O reactions need to remove only the EXTRA amount
# caused by fractionation -- the DEFICIT above bulk proportional
# removal:
#
#   extra_HDO    = rate_H2O * 2 * R(D)   * (alpha_eff - 1) * TIME
#   extra_H2_18O = rate_H2O * R(18O) * (alpha_eff - 1) * TIME
#
# This extra removal enriches the residual liquid in heavy isotopes,
# producing the expected upward trend in delta values over time.
#
# alpha_eff deficit values at 28 deg C verified:
#   h=0.85: D deficit=0.07760,  18O deficit=0.01125
#   h=0.70: D deficit=0.07961,  18O deficit=0.01340
#   h=0.55: D deficit=0.08163,  18O deficit=0.01555
#   h=0.40: D deficit=0.08365,  18O deficit=0.01770
#
# Horita & Wessolowski (1994) give 1000*ln(alpha) -- DIVIDE BY 1000
# before EXP() to get alpha.
# TK is reserved in PHREEQC Basic -- use TEMP_K for your variable.
# ================================================================

RATES

# ----------------------------------------------------------------
# Evap_H2O: bulk water evaporation - reference rate
# PARM(1) = evaporation rate mol/s
# Removes bulk H2O and stores rate in slot 1 for isotope reactions
# ----------------------------------------------------------------
Evap_H2O
-start
10  rate  = PARM(1)
20  moles = rate * TIME
30  SAVE  -moles
40  PUT(rate, 1)
-end

# ----------------------------------------------------------------
# Evap_HDO: removes the EXTRA HDO above bulk proportional removal
# This drives isotopic enrichment of the residual liquid in D
# PARM(1) = relative humidity h
#
# extra_HDO = rate_H2O * 2 * R(D) * (alpha_eff - 1) * TIME
#
# Horita & Wessolowski (1994) for D:
#   1000*ln(alpha_eq) = 1158.8*(T^3/1e9) - 1620.1*(T^2/1e6)
#                       + 794.84*(T/1e3) - 161.04
#                       + 2.9992*(1e9/T^3)
# Merlivat & Jouzel (1979):
#   eps_k_D = 12.5*(1-h) permil
# ----------------------------------------------------------------
Evap_HDO
-start
10  h            = PARM(1)
20  TEMP_K       = TK

# Equilibrium fractionation - divide by 1000 before EXP
30  tln_aeq_D    = 1158.8*(TEMP_K^3/1e9) - 1620.1*(TEMP_K^2/1e6) + 794.84*(TEMP_K/1e3) - 161.04 + 2.9992*(1e9/TEMP_K^3)
40  alpha_eq     = EXP(tln_aeq_D / 1000)

# Kinetic fractionation
50  eps_k        = 12.5e-3 * (1 - h)
60  alpha_k      = 1 + eps_k

# Total effective fractionation
70  alpha_eff    = alpha_eq * alpha_k

# Current D/H ratio of solution
80  r_d_soln     = CALC_VALUE("R(D)")

# Remove only the EXTRA HDO above bulk proportional evaporation
# This is what drives enrichment of the residual liquid
90  rate_H2O     = GET(1)
100 moles        = rate_H2O * 2 * r_d_soln * (alpha_eff - 1) * TIME
110 SAVE -moles
-end

# ----------------------------------------------------------------
# Evap_H2_18O: removes the EXTRA H2[18O] above bulk removal
# This drives isotopic enrichment of the residual liquid in 18O
# PARM(1) = relative humidity h
#
# extra_H2_18O = rate_H2O * R(18O) * (alpha_eff - 1) * TIME
#
# Horita & Wessolowski (1994) for 18O:
#   1000*ln(alpha_eq) = -7.685 + 6.7123*(1e3/T)
#                       - 1.6664*(1e6/T^2)
#                       + 0.35041*(1e9/T^3)
# Merlivat & Jouzel (1979):
#   eps_k_18O = 14.2*(1-h) permil
# ----------------------------------------------------------------
Evap_H2_18O
-start
10  h            = PARM(1)
20  TEMP_K       = TK

# Equilibrium fractionation - divide by 1000 before EXP
30  tln_aeq_18O  = -7.685 + 6.7123*(1e3/TEMP_K) - 1.6664*(1e6/TEMP_K^2) + 0.35041*(1e9/TEMP_K^3)
40  alpha_eq     = EXP(tln_aeq_18O / 1000)

# Kinetic fractionation
50  eps_k        = 14.2e-3 * (1 - h)
60  alpha_k      = 1 + eps_k

# Total effective fractionation
70  alpha_eff    = alpha_eq * alpha_k

# Current 18O/16O ratio of solution
80  r_18O_soln   = CALC_VALUE("R(18O)")

# Remove only the EXTRA H2[18O] above bulk proportional evaporation
90  rate_H2O     = GET(1)
100 moles        = rate_H2O * r_18O_soln * (alpha_eff - 1) * TIME
110 SAVE -moles
-end

# ================================================================
# STARTING SOLUTION
# Ocean vapour condensate: d18O = -10 permil, dD = -80 permil
# temp = 28 deg C
# ================================================================

SOLUTION 1  Ocean vapour condensate - starting water
    temp      28
    pH        7.0   charge
    units     mg/L
    [18O]     -10.0
    D         -80.0
    O(0)      1     O2(g)  -0.7

SAVE solution 1
END

# ================================================================
# SCENARIO 1: h = 0.85  (wet season, high humidity)
# alpha_eff_D deficit   = 0.07760
# alpha_eff_18O deficit = 0.01125
# Slowest enrichment -- smallest kinetic fractionation
# ================================================================

TITLE Scenario 1 - h = 0.85 (wet season, slow evaporation)

USE solution 1

KINETICS 1
    Evap_H2O
        -formula   H2O      1
        -m         0
        -m0        0
        -parms     1e-5
        -tol       1e-10
    Evap_HDO
        -formula   HDO      1
        -m         0
        -m0        0
        -parms     0.85
        -tol       1e-10
    Evap_H2_18O
        -formula   H2[18O]  1
        -m         0
        -m0        0
        -parms     0.85
        -tol       1e-10
    -cvode
    -steps         2592000  in  150 steps
    -bad_step_max  1000

INCREMENTAL_REACTIONS true

SELECTED_OUTPUT 1
    -file          river_h085.csv
    -reset         false
    -time          true

USER_PUNCH 1
    -headings  Day  d18O_liquid  dD_liquid  d18O_vapour  dD_vapour  d_excess
    -start
    10  PUNCH TOTAL_TIME / 86400
    20  PUNCH CALC_VALUE("Delta_[18O]")
    30  PUNCH CALC_VALUE("Delta_D")
    # Instantaneous vapour delta = liquid delta minus equilibrium epsilon at 28 C
    # epsilon_18O_eq = 9.1 permil, epsilon_D_eq = 73.0 permil
    40  d18O_v = CALC_VALUE("Delta_[18O]") - 9.1
    50  dD_v   = CALC_VALUE("Delta_D") - 73.0
    60  PUNCH d18O_v
    70  PUNCH dD_v
    # d-excess = dD - 8 * d18O  (deviation below GMWL)
    80  PUNCH CALC_VALUE("Delta_D") - 8 * CALC_VALUE("Delta_[18O]")
    -end

USER_GRAPH 1
    -headings           Day  d18O_liq_h085  dD_liq_h085
    -axis_titles        "Days"  "Isotope delta (permil VSMOW)"  ""
    -initial_solutions  false
    -connect_simulations  true
    -plot_concentration_vs  x
    -start
    10  GRAPH_X  TOTAL_TIME / 86400
    20  GRAPH_Y  CALC_VALUE("Delta_[18O]")
    30  GRAPH_SY CALC_VALUE("Delta_D")
    -end

END

# ================================================================
# SCENARIO 2: h = 0.70  (moderate humidity)
# alpha_eff_D deficit   = 0.07961
# alpha_eff_18O deficit = 0.01340
# Slightly faster enrichment than h=0.85
# ================================================================

TITLE Scenario 2 - h = 0.70 (moderate humidity)

USE solution 1

KINETICS 2
    Evap_H2O
        -formula   H2O      1
        -m         0
        -m0        0
        -parms     1e-5
        -tol       1e-10
    Evap_HDO
        -formula   HDO      1
        -m         0
        -m0        0
        -parms     0.70
        -tol       1e-10
    Evap_H2_18O
        -formula   H2[18O]  1
        -m         0
        -m0        0
        -parms     0.70
        -tol       1e-10
    -cvode
    -steps         2592000  in  150 steps
    -bad_step_max  1000

INCREMENTAL_REACTIONS true

SELECTED_OUTPUT 2
    -file          river_h070.csv
    -reset         false
    -time          true

USER_PUNCH 2
    -headings  Day  d18O_liquid  dD_liquid  d18O_vapour  dD_vapour  d_excess
    -start
    10  PUNCH TOTAL_TIME / 86400
    20  PUNCH CALC_VALUE("Delta_[18O]")
    30  PUNCH CALC_VALUE("Delta_D")
    40  d18O_v = CALC_VALUE("Delta_[18O]") - 9.1
    50  dD_v   = CALC_VALUE("Delta_D") - 73.0
    60  PUNCH d18O_v
    70  PUNCH dD_v
    80  PUNCH CALC_VALUE("Delta_D") - 8 * CALC_VALUE("Delta_[18O]")
    -end

USER_GRAPH 2
    -headings           Day  d18O_liq_h070  dD_liq_h070
    -axis_titles        "Days"  "Isotope delta (permil VSMOW)"  ""
    -initial_solutions  false
    -connect_simulations  true
    -plot_concentration_vs  x
    -start
    10  GRAPH_X  TOTAL_TIME / 86400
    20  GRAPH_Y  CALC_VALUE("Delta_[18O]")
    30  GRAPH_SY CALC_VALUE("Delta_D")
    -end

END

# ================================================================
# SCENARIO 3: h = 0.55  (dry season)
# alpha_eff_D deficit   = 0.08163
# alpha_eff_18O deficit = 0.01555
# Stronger enrichment -- larger kinetic fractionation
# ================================================================

TITLE Scenario 3 - h = 0.55 (dry season, faster evaporation)

USE solution 1

KINETICS 3
    Evap_H2O
        -formula   H2O      1
        -m         0
        -m0        0
        -parms     1e-5
        -tol       1e-10
    Evap_HDO
        -formula   HDO      1
        -m         0
        -m0        0
        -parms     0.55
        -tol       1e-10
    Evap_H2_18O
        -formula   H2[18O]  1
        -m         0
        -m0        0
        -parms     0.55
        -tol       1e-10
    -cvode
    -steps         2592000  in  150 steps
    -bad_step_max  1000

INCREMENTAL_REACTIONS true

SELECTED_OUTPUT 3
    -file          river_h055.csv
    -reset         false
    -time          true

USER_PUNCH 3
    -headings  Day  d18O_liquid  dD_liquid  d18O_vapour  dD_vapour  d_excess
    -start
    10  PUNCH TOTAL_TIME / 86400
    20  PUNCH CALC_VALUE("Delta_[18O]")
    30  PUNCH CALC_VALUE("Delta_D")
    40  d18O_v = CALC_VALUE("Delta_[18O]") - 9.1
    50  dD_v   = CALC_VALUE("Delta_D") - 73.0
    60  PUNCH d18O_v
    70  PUNCH dD_v
    80  PUNCH CALC_VALUE("Delta_D") - 8 * CALC_VALUE("Delta_[18O]")
    -end

USER_GRAPH 3
    -headings           Day  d18O_liq_h055  dD_liq_h055
    -axis_titles        "Days"  "Isotope delta (permil VSMOW)"  ""
    -initial_solutions  false
    -connect_simulations  true
    -plot_concentration_vs  x
    -start
    10  GRAPH_X  TOTAL_TIME / 86400
    20  GRAPH_Y  CALC_VALUE("Delta_[18O]")
    30  GRAPH_SY CALC_VALUE("Delta_D")
    -end

END

# ================================================================
# SCENARIO 4: h = 0.40  (very dry, maximum kinetic fractionation)
# alpha_eff_D deficit   = 0.08365
# alpha_eff_18O deficit = 0.01770
# Fastest enrichment -- strongest kinetic fractionation
# ================================================================

TITLE Scenario 4 - h = 0.40 (very dry, maximum kinetic fractionation)

USE solution 1

KINETICS 4
    Evap_H2O
        -formula   H2O      1
        -m         0
        -m0        0
        -parms     1e-5
        -tol       1e-10
    Evap_HDO
        -formula   HDO      1
        -m         0
        -m0        0
        -parms     0.40
        -tol       1e-10
    Evap_H2_18O
        -formula   H2[18O]  1
        -m         0
        -m0        0
        -parms     0.40
        -tol       1e-10
    -cvode
    -steps         2592000  in  150 steps
    -bad_step_max  1000

INCREMENTAL_REACTIONS true

SELECTED_OUTPUT 4
    -file          river_h040.csv
    -reset         false
    -time          true

USER_PUNCH 4
    -headings  Day  d18O_liquid  dD_liquid  d18O_vapour  dD_vapour  d_excess
    -start
    10  PUNCH TOTAL_TIME / 86400
    20  PUNCH CALC_VALUE("Delta_[18O]")
    30  PUNCH CALC_VALUE("Delta_D")
    40  d18O_v = CALC_VALUE("Delta_[18O]") - 9.1
    50  dD_v   = CALC_VALUE("Delta_D") - 73.0
    60  PUNCH d18O_v
    70  PUNCH dD_v
    80  PUNCH CALC_VALUE("Delta_D") - 8 * CALC_VALUE("Delta_[18O]")
    -end

USER_GRAPH 4
    -headings           Day  d18O_liq_h040  dD_liq_h040
    -axis_titles        "Days"  "Isotope delta (permil VSMOW)"  ""
    -initial_solutions  false
    -connect_simulations  true
    -plot_concentration_vs  x
    -start
    10  GRAPH_X  TOTAL_TIME / 86400
    20  GRAPH_Y  CALC_VALUE("Delta_[18O]")
    30  GRAPH_SY CALC_VALUE("Delta_D")
    -end

END


[dlparkhurst]

Here is an example that removes H2O, HDO, and D2O kinetically. The rate equations use a rate dependent on 1/sqrt(mass) for the water isotopomers and the number of moles of each species in the solution. The absolute rates and the time of the simulation are arbitrary to show the evolution.

The key is that H2O evaporates faster than HDO so that the ratio of HDO to H2O increases--that is, the permil value of D gets  higher (heavier.) To get to 2 permil with these rates requires removal of about 96 percent of the water.

Now for relative humidity, the rate is probably dependent on the difference between the calculated partial pressures (fugacities) of H2O(g), HDO(g), and D2O(g) for the solution and the partial pressures in the gas phase. But you can consider how you want to define the kinetic reactions.

Code: [Select]

RATES
H2O_evap
-start
10 k = 1/sqrt(GFW("H2O"))
20 rate = k * MOL("H2O")*TOT("Water")
30 moles = -rate*TIME
40 SAVE moles
-end
HDO_evap
-start
10 k = 1/sqrt(GFW("HDO"))
20 rate = k * MOL("HDO")*TOT("Water")
30 moles = -rate*TIME
40 SAVE moles
-end
D2O_evap
-start
10 k = 1/sqrt(GFW("D2O"))
20 rate = k * MOL("D2O")*TOT("Water")
30 moles = -rate*TIME
40 SAVE moles
-end
END
SOLUTION 1  Ocean vapour condensate - starting water
    temp      28
    pH        7.0   charge
    units     mg/L
    D         -80.0
END
KINETICS
H2O_evap
-formula H2O 1
      -m 0
HDO_evap
-formula HDO 1
      -m 0
D2O_evap
-formula D2O 1
      -m 0
-step 25 in 25
-cvode
-step_divide 10
END
USER_GRAPH 1
    -headings               fraction D_permil H2O HDO D2O
    -axis_titles            "Fraction of H2O removed" "D, permil" "Log moles"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X 1 - TOT("water")
20 GRAPH_Y ISO("D")
30 GRAPH_SY log10[MOL("H2O")*TOT("water")], log10[MOL("HDO")*TOT("water")], log10[MOL("D2O")*TOT("water")]
  -end
    -active                 true
END
INCREMENTAL_REACTIONS
USE solution 1
USE kinetics 1
END