Processes > Inverse modelling

Modelling isotope fractionation during evaporation

<< < (2/2)

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

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

--- End code ---

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.

Navigation

[0] Message Index

[*] Previous page