Equilibrium fractionation occurs in closed systems. Equilibrium fractionation of water during phase changes is mainly controlled by temperature. The initial isotopic ratio

$$R_A = 2005.2 \cdot 10^{-6}$$

of seawater is changed by a fractionation fator $\alpha$. This fractionation factor can be calculated using the empirical relationship

$$ 1000 \: ln(\alpha) = D \cdot \frac{10^6}{T^2} + E \cdot \frac{10^3}{T} + F$$

where $D$, $E$ and $F$ are constants that are specific for each phase change, e.g. for water to vapour. The direction of the phase change is significant. The fractionation factor is defined as

$$\alpha = \frac{R_A}{R_B}$$

In this case as

$$\alpha_{evaporation} = \frac{R_{water}}{R_{vapour}}$$ The term

$$1000 \cdot ln(\alpha)$$

is chosen, to produce values that can be read and interpreted more easily and as $ln(\alpha)$ is quite close (within the error margin of the measurement) to $(alpha-1)\cdot(1000)$ for small deviations of the fractionation factor $\alpha$ from 1.0 (or for small fractionations – which is very often the case in hydrology). Once the term $1000 \cdot ln(\alpha)$ has been calculated for a given phase change, it can be converted to $\alpha$. Finally, the isotope ratio of the resulting phase B can be calculated using the relationship $\alpha=R_A/R_B$

In hydrology, isotope ratios are given in permil as

$$\delta { }^{18}O = \frac{R_{Standard}-R_{Probe}}{R_{Standard}} \cdot 1000 $$

Hence, once the isotope ratio $R_{Probe}$ has been calculated, it can be converted to a $\delta$ value using that relationship. Many $\alpha$ values can be found at alphadelta

Fractionation.R

library(ggplot2)
library(latex2exp)
 
# Rs Seawater is the Standard for stable Isotope analysis of water
# The ratios are given with their errors
R18O = 2005.2E-6
R18Oe= 0.45E-6
R2H  = 155.76E-6
R2He= 0.05E-6
d18Oe = 0.1
d2He  = 1.0
 
# Majoube 18O/16O, equilibrium coefficients for water - vapour system
D18O = 1.137
E18O = -0.42
F18O = -2.07
 
# 2H/H
D2H =  24.844
E2H = -76.25
F2H =  52.61
 
# Temperature C & K - define a temperature range
TC  = seq(0,35,5)
 
# Calculate equilibrium fractionation
equilibrium <- function(TC){
  # convert TC to TK
  TK <- TC*0
  TK  = TC+273.15
  lnta18O = D18O*10^6/TK^2+E18O*10^3/TK+F18O 
  lnta2H  = D2H *10^6/TK^2+ E2H*10^3/TK+F2H 
  lntas   = data.frame(TC,TK,lnta18O,lnta2H)
 
  # ln(alpha)*1000 -> alpha
  alpha18O    = exp(lnta18O/1000)
  alpha2H     = exp(lnta2H/1000)
  lntas$alpha18O <- alpha18O   
  lntas$alpha2H  <- alpha2H
 
  # R(Start) -> R(Sample)
  S18O = R18O * alpha18O
  S2H  = R2H  * alpha2H
  lntas["S18O"] <- S18O                 # Add new column to data
  lntas["S2H"]  <- S2H                  # Add new column to data
 
  # Rs -> delta
  d18O = (R18O- S18O) /R18O*1000
  d2H  = (R2H - S2H)  / R2H*1000
 
  lnta  <- lntas                        # Replicate example data
  lntas <- cbind(lnta, d18O = d18O)     # Add new column to data
  lnta  <- lntas                        # Replicate example data
  lntas <- cbind(lnta, d2H  = d2H)     # Add new column to data
  return(lntas)
}
 
lntas <- equilibrium(TC)
attach(lntas)
 
ggplot(data=lntas, mapping=aes(x=TC, y=lnta18O)) +
  # ggtitle("Equilibrium fractionation of water") +
  geom_line(linetype="dashed",  size = 0.5, color="blue") +
  geom_point(shape=21, fill="white", color="blue", size=3) +
  xlab("Temperature (Kelvin)") +
  ylab(TeX("$1000 \\cdot ln(\\alpha)$")) 
 
ggplot(data=lntas, mapping=aes(x=d18O, y=d2H)) +
  # ggtitle("Equilibrium fractionation of water") +
  geom_line(linetype="dashed",  size = 0.5, color="blue") +
  geom_point(shape=21, fill="white", color="blue", size=3) +
  geom_errorbar(aes(ymin=(d2H-d2He), ymax=(d2H+d2He)), color="grey") +
  geom_errorbarh(aes(xmin = d18O-d18Oe, xmax = d18O+d18Oe), color="grey") +
  xlab(TeX("\\delta 18O $")) +
  ylab(TeX("\\delta 2H  $")) +
  geom_text(label=TC, nudge_x=0.15, nudge_y=0.1, check_overlap=T)

If you do not have R or your computer or do not want to install it now, you can use the window below to run the code online. Just copy the R code in the (grey) block and insert it into the window. The example code can be deleted.

Diffuse fractionation occurs, if a gas is moving through a vaccum or through another gas. A simpler case is given, if the gase travels through another inert gas. A more complex fractionation occurs, if the gas moves through another gase that contains the moving gas and that gets enriched with the moving gas during the process. This is the case, for example, when evaporation occurs over a water surface and changes the vapour content of the boundary layer during the process as a result of evaporation.

If the system is open, as during evaporation from a pond or during rainfall from a cloud, non-equilibrium fractionation occurs. There are different types of non-equilbrium fractionation. The first type is called diffuce fractionation. The basis for the calculation of diffuse fractionation effects is the difference in velocities of moving molecules. The water molecules with different isotopes having different masses, their velocities differ according to :

$$v_i = \sqrt{\frac{k \cdot T}{2 \cdot pi \cdot m}}$$

with $v$ the velocity of the molecule(s) in $m/s$, $k$ the Boltzmann Konstant $1.380649 \cdot 10^{-23} \: \frac{J}{K}$, $T$ the temperature in Kelvin and $m$ the mass of the molecule in kg. The fractionation can be derived from the ratio of velocities of different isotopes and $\alpha$ is in proportion to

$$ \alpha_{vacuum} = \frac{v_L}{v_H} = \frac{\sqrt{\frac{k \cdot T}{2 \cdot pi \cdot m_L}}}{ \sqrt{\frac{k \cdot T}{2 \cdot pi \cdot m_H}}}=\sqrt{\frac{m_H}{m_L}}=\sqrt{\frac{20}{18}}=1.054093$$ for vapour with a heavy isotope of oxygen and respectively for a heavy isotope of hydrogen:

$$ \alpha_{vacuum} = \frac{v_L}{v_H} = \frac{\sqrt{\frac{k \cdot T}{2 \cdot pi \cdot m_L}}}{ \sqrt{\frac{k \cdot T}{2 \cdot pi \cdot m_H}}}=\sqrt{\frac{m_H}{m_L}}=\sqrt{\frac{19}{18}}=1.027402$$

It can be shown that for vapour moving through an inert gas such as e.g. dry air with a mas of 28.8 g per mol, the diffuse fractionation is in proportion to

$$ \alpha = \sqrt{\frac{m_H \cdot (m_L+28.8)}{m_L \cdot (m_H+28.8)}}= \sqrt{\frac{20 \cdot (18+28.8)}{18 \cdot (20+28.8)}}=1.032266$$

in the case of water vapour containing a heavy isotope of oxygen and for vapour containing a heavy isotope of hydrogen instead:

$$ \alpha = \sqrt{\frac{m_H \cdot (m_L+28.8)}{m_L \cdot (m_H+28.8)}}= \sqrt{\frac{19 \cdot (18+28.8)}{18 \cdot (19+28.8)}}=1.016599$$

This calculation shows, that oxygen isotopes fractionate more than hydrogen isotopes, a phenomenon commonly observed in evaporation of water from free surfaces.

However, as air usually contains vapour that has an isotopic composition as well and as the vapour in the air through which vapour diffuses interacts and exchanges with the vapour as well in terms of their isotopic composition, the fractionation of moist air is lower. Gonfiantini (1978) has suggested the empirical relationship:

$$\epsilon = (1-h) \cdot ( \frac{\rho_i}{\rho}-1)$$

Here, $\epsilon$ is the isotopic change resulting from the fractionation process with a fractionation factor of $\alpha$, hence

$$\epsilon_{A \rightarrow B} = \left[\frac{R_A}{R_B}-1\right] \cdot 1000 = \left[ \alpha_{A \rightarrow B}-1\right] \cdot 1000 $$ Inserting the densities for the different isotopic species of vapour with a heavy oxygen isotope and with a heavy hydrogen isotope, the empirical formulas for the isotopic change in permil read

$$\epsilon_{{}^{18} O}=14.2 \cdot (1-h) $$

and

$$\epsilon_{{}^{2} H}=12.5 \cdot (1-h) $$

In the case of evaporation from an open water surface, both effects combine and add up. The total fractionation at a temperature $T$ and at a humidity $h$ can be calculated as:

$$\delta_{final}=\delta_{initial}+\epsilon_{equilibrium,\, T}+\epsilon_{diffuse,\, h}$$ An application of this stepwise calculation indicates that the evaporation slope approaches the theoretical value for equilibration with increasing humidity. At the lowest possible humidity, the slope of the evaporation line is between 4.4 at lower temperatures and 3.4 at high temperatures in dry ambient air.

Diffuse.R

  hum   = 0.8
  diffusefrac <- function(humidity = 0.8){
    df18O = 14.2
    df2H  = 12.5
    eps18 = df18O*(1-humidity)
    eps2H = df2H *(1-humidity)
    epsdf = data.frame(humidity, eps18, eps2H)
    return(epsdf)
  }
  epsdf <- diffusefrac(hum)

Both processes can be combined to obtain the isotopic fractionation during evaporation. At a given temperature a phase equilibrium between water and vapour will develop according to the equilibrium fractionation process. In addition, diffuse fractionation will conctribute additional fractionation by diffusion. The result of both processes results in isotope ratios of heavy and light oxygen and hydrogen isotopes in water that define the slope of the evaporation curve. For a given temperature and humidity, the combined effect of equilibrium and diffuse fractionation can now be calculated (still provided that sufficient water is available to prevent any fractionation effects due to enrichment in the reservoir of evaporating water).

Rayleigh fractionation occurs, if a reservoir is depleted significantly. The initial istopic ratio $R_0$ changes towards the actual ratio $R$ at a certain degree of depletion $f$ if the fractionation factor for the depleting process is $\alpha$:

$$R = R_0 \cdot f^{\left(\alpha-1\right)}$$

Code for Rayleigh fractionation. Please not that it requires the function equilibrium that is provided above. To run this code, you can add it to the first block of code at the bottom and execute the entire block - the function is then available.

Rayleigh.R

alpha18O <- equilibrium(10)$alpha18O
fR      <- seq(0.9,1,0.001)
R0      <- 2005.2E-6
R       <-  R0*fR^(alpha18O-1)
plot(fR,R, type="l")