REVIEW

Special Section: Tribute to Rien van Genuchten, Recipient of the 2023 Wolf Prize for Agriculture

Vadose Zone Journal

# Analysis of experimental and simulation data of evaporation-driven isotopic fractionation in unsaturated porous media

Jana Schneider<fnref n="1" />,*,

<figure><img src="image_1.png" /></figure>

Isabelle Braud<fnref n="3" />

Stefanie Kiemle<fnref n="2" />,*

<figure><img src="image_2.png" /></figure>

Rainer Helmig<fnref n="2" />

Katharina Heck<fnref n="2" />

<figure><img src="image_3.png" /></figure>

Jan Vanderborght<fnref n="1" />

<figure><img src="image_4.png" /></figure>

Youri Rothfuss<fnref n="1" />

<figure><img src="image_5.png" /></figure>

## Correspondence

Jan Vanderborght, Institute of Bio-, and Geosciences: Agrosphere (IBG-3), Forschungszentrum Jülich GmbH, Jülich 52425, Germany.

Email: j.vanderborght@fz-juelich.de

Assigned to Associate Editor Binayak Mohanty.

## Funding information

Deutsche Forschungsgemeinschaft, Grant/Award Numbers: 2531/14-1, 327154368

## Abstract

Stable water isotopologs can add valuable information to the understanding of evaporation processes. The identification of the evaporation front from isotopolog concentration depth profiles under very dry soil conditions is of particular interest. We compared two different models that describe isotopolog transport in a drying unsaturated porous medium: SiSPAT-Isotope and DuMuₓ. In DuMuₓ, the medium can dry out completely whereas in SiSPAT-Isotope, drying is limited to the residual water saturation. We evaluated the impact of residual water saturation on simulated isotopic concentration. For a low residual water saturation, both models simulated similar isotopolog concentrations. For high residual water saturation, SiSPAT-Isotope simulated considerably lower concentrations than DuMuₓ. This is attributed to the buffering of changes in isotopolog concentrations by the residual water in SiSPAT-Isotope and an additional enrichment due to evaporation of residual water in DuMuₓ. Additionally, we present a comparison between high-frequency experimental data and model simulations. We found that diffusive transport processes in the laminar boundary layer and in the dried-out surface soil layer need to be represented correctly to reproduce the observed downward movement of the evaporation front and the associated peak of isotopolog enrichment. Artificially increasing the boundary layer thickness to reproduce a decrease in evaporation rate leads to incorrect simulation of the location of the evaporation front and isotopolog concentration profile.

**Abbreviations:** air, dry air; ff, free flow (atmosphere); g, gaseous; H₂O, water; i, stable water isotopolog; l, liquid; pm, porous medium (soil); REV, representative elementary volume; s, solid; VG, van Genuchten; α, phase; κ, component.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

© 2024 The Author(s). *Vadose Zone Journal* published by Wiley Periodicals LLC on behalf of Soil Science Society of America.

<footnote n="1">¹Institute of Bio- and Geosciences: Agrosphere (IBG-3), Forschungszentrum Jülich GmbH, Jülich, Germany</footnote>

<footnote n="2">²Institute for Modelling Hydraulic and Environmental Systems, Department of Hydromechanics and Modelling of Hydrosystems, University of Stuttgart, Stuttgart, Germany</footnote>

<footnote n="3">³INRAE, Riverly, Villeurbanne, France</footnote>

<footnote>*Jana Schneider and Stefanie Kiemle contributed equally to this paper.</footnote>
<figure><img src="image_6.png" /></figure>

# 1 | INTRODUCTION

Evaporation is a key variable in the global water cycle— depending on the land surface type, soil evaporation makes up 10% to 87% of the total evapotranspiration (Rothfuss et al., 2021). Especially in dry regions, soil evaporation can account for up to 95% of the total evapotranspiration (Kool et al., 2014; Wilcox et al., 2003).

Soil evaporation is defined as the water loss from the soil surface that is not associated with carbon assimilation by the vegetation. Reducing the water loss caused by soil evaporation increases the amount of water that can be used for plant transpiration, which is linked with carbon uptake by vegetation. Consequently, this uptake leads to an increase in vegetation growth and crop growth, and finally to enhanced crop yields. Especially in dry regions and in terms of climate change, reducing soil evaporation is a key process to optimize water use management (Anderson et al., 2017; Boutraa et al., 2010; Ghumman et al., 2020; Kool et al., 2014; Milly et al., 2005; Oki & Kanae, 2006; Stoy et al., 2019).

Water-stable isotopologs ($^1H^2H^{16}O$ and $H_2^{18}O$) are powerful tracers for water flow in soil–plant–atmosphere systems (Dubbert & Werner, 2019; Braud et al., 2005; Rothfuss & Javaux, 2017; Rothfuss et al., 2021; Yakir & da SL Sternberg, 2000), as their transport and gas–liquid exchange behavior differs from ordinary water. This difference in transport behavior leads to a change in isotopolog concentrations in the pore water that is related to evaporation rates. Furthermore, they can be used to infer how incoming water with a different isotopic signature than the water present in the soil pores intermixes and displaces pore water (Braud et al., 2005; Herczeg & Leaney, 2011; Rothfuss et al., 2010, 2015). Consequently, the variation in the ratio of soil water isotopes ($\delta^2H$, $\delta^{18}O$) provides information about flow, mixing within the unsaturated zone, the evaporation front, and interactions between atmosphere and soil. Especially for evaporation, measurements of isotope concentration profiles can be used to assess evaporation rates and the depth of the evaporation front in the soil profile, as evaporation processes at the soil surface and within the soil lead to an enrichment of the heavier isotopologs in the liquid phase (Rothfuss et al., 2021).

Modeling the transport of isotopologs in the unsaturated zone can enhance the understanding of the influencing parameters of evaporation and can be used to assess water management strategies. There are several models available that can model isotopolog transport and fractionation for that purpose. State-of-the-art for water-stable isotopolog modeling is the Simple Soil–Plant–Atmosphere Transport model (SiSPAT–Isotope), which is a soil–vegetation–atmosphere transfer model including a module for isotopolog transport in the soil and within plants (Braud et al., 2005; Haverd & Cuntz, 2010). The liquid water flow and vapor transport model is based on the Philip and de Vries (1957) model. Another

## Core Ideas

* Isotopolog fractionation in the unsaturated zone was studied.

* A comparison between high-frequency experimental data of isotope enrichment and model simulations was carried out.

* Effects of drying below residual water content and of boundary conditions on isotope concentrations were analyzed.

* The evaporation front depth was estimated by using the depth of the concentration gradient maximum.

model based on the same physical processes is Soil–Litter–Iso (Haverd & Cuntz, 2010; Haverd et al., 2011), developed to be used on a regional scale based on Ross’ fast solution to the Richards equation (Ross, 2006). Additionally, the HYDRUS-1D model can simulate isotopolog transport while accounting for evaporation fractionation (Zhou et al., 2021). The fractionation process is based on the same equations as the SiSPAT–Isotope module. Recently, Kiemle et al. (2023) presented a model for isotopolog transport in a coupled soil–atmosphere model within the framework of the numerical simulator DuMu$^x$ (Koch et al., 2021).

As isotopolog fractionation mainly occurs at the evaporation front, that is, the interface between dry and wet soil, the model should be able to describe the drying of the unsaturated zone to dry soil conditions. Generally, when the soil dries out, capillary pressure rises. To describe that process on the continuum scale, commonly a relationship between capillary pressure and saturation is used (e.g., the van Genuchten relationship van Genuchten, 1980).

Under dry soil conditions, this relationship generally leads to large gradients and theoretically infinite capillary pressures, which makes the numerical modeling of these conditions particularly challenging. There are different possibilities to overcome these difficulties. One possibility is to regularize the capillary–pressure–saturation relationship after a certain threshold in water saturation is reached, which is the method that is used in DuMu$^x$. Another possibility, which is used in SiSPAT–Isotope, is to not allow complete drying of the soil but to keep a certain amount of water saturation immobile without further evaporation by setting a maximum capillary pressure that cannot be exceeded. If this maximum capillary pressure can be set sufficiently high so that the vapor pressure in equilibrium with this capillary pressure is lower than the minimum vapor pressure simulated or observed, then the fluxes from the porous medium to the soil surface can be matched to the fluxes in the free flow boundary layer, taking into account the flow and transport resistances in the porous medium and the boundary layer defined by their properties
and states. However, when this maximal capillary pressure is reached in the simulation and kept constant, the vapor pressures at the soil surface cannot further drop so that fluxes through the boundary layer and evaporation rates would stay constant whereas the fluxes from the drying porous medium toward the soil surface decrease further. To avoid this discrepancy in a coupled simulation, the boundary layer thickness should be increased. These differences in the handling of very dry soil conditions, not only change simulated evaporation rates but also have implications for isotopolog transport and fractionation. One goal of this paper is to compare the isotopolog transport described in SiSPAT-Isotope with the implementation in DuMu<sup>x</sup>, especially highlighting the differences that occur in dry soil with water saturation below residual water saturation. Further, we will present a comparison of the simulated isotope concentration profiles with the high-frequency experimental data from Rothfuss et al. (2015).

Since the experimental conditions in the laboratory experiments did not allow us to precisely define the surface boundary conditions, which are essential for modeling evaporation, we had to derive them from the measured evaporation rates. Therefore, the aim of our simulation study is not to simulate the evaporation (which we basically prescribe) but to simulate the evaporation-driven isotope fractionation, that is, the evolution in time and depth of the isotope concentrations in the soil profile.

A comparison of modeled isotopolog fractionation in the unsaturated zone with detailed experimental data is missing in the literature. Last, we will show that the evaporation front can be detected at the maximum gradient of isotopolog concentration.

## 2 | METHODS

In this section, we introduce the general mathematical concepts of DuMu<sup>x</sup> and SiSPAT-Isotope to describe mass and energy transport within a porous medium.

### 2.1 | Mathematical concept

Note, that we use the following abbreviations to describe our systems. We use the components κ: water (H₂O), dry air (air), and stable water isotopolog (i, for all isotopologs). The phase (α) is either referred to as liquid (l) or gaseous (g) or solid (s). The abbreviation (ff) is used for the free-flow domain (atmosphere domain) and the abbreviation (pm) is for the porous-medium domain (soil domain). All processes are described on the representative elementary volume (REV)-scale.

Note, a detailed description of the assumptions for the mass and energy balance concepts used in DuMu<sup>x</sup> and SiSPAT-isotope can be found in Vanderborght et al. (2017).

#### 2.1.1 | Mass and energy transport in the porous medium

**DuMu<sup>x</sup>—Mass balance**

To describe non-isothermal, multi-component, multi-phase flow within a porous medium, DuMu<sup>x</sup> solves a multi-phase Darcy's law in combination with a mass and energy balance. The component mass balance for all components κ can be described as follows:

$$ \sum_{\alpha \in \{l,g\}} \left( \phi \frac{\partial}{\partial t} (\rho_{\alpha} X_{\alpha}^{\kappa} S_{\alpha}) + \nabla \cdot (\mathbf{v}_{\alpha} \rho_{\alpha} X_{\alpha}^{\kappa}) - \nabla \cdot (D_{\alpha,pm}^{\kappa} \rho_{\alpha} \nabla X_{\alpha}^{\kappa}) \right) = q^{\kappa} \quad \forall \kappa \in \{\text{H}_2\text{O, air, } i\}. \quad (1) $$

Here, φ [-] is the porosity, ρ<sub>α</sub> (kg m<sup>-3</sup>) the mass density, X<sub>α</sub><sup>κ</sup> (kg kg<sup>-1</sup>) is the mass fraction, D<sub>α,pm</sub><sup>κ</sup> (m<sup>2</sup>s<sup>-1</sup>) is the effective diffusion coefficient, and q<sup>κ</sup> (kg m<sup>-3</sup> s<sup>-1</sup>) accounts for the source/sink term. As we are regarding all processes on the REV-scale, we introduce an effective parameter, namely the phase saturation S<sub>α</sub> (m<sup>3</sup>m<sup>-3</sup>). The sum of the saturation for all phases within a REV equals 1. The fluid-phase flux v<sub>α</sub> (m s<sup>-1</sup>) is described using the multi-phase Darcy's law:

$$ \mathbf{v}_{\alpha} = - \frac{k_{r,\alpha}}{\mu_{\alpha}} \mathbf{K} (\nabla p_{\alpha} - \rho_{\alpha} g), \quad (2) $$

with the relative permeability k<sub>r,α</sub> [-], the dynamic viscosity μ<sub>α</sub> [Pa s], the phase pressure p<sub>α</sub> [Pa], and the intrinsic permeability tensor K (m<sup>2</sup>). The gravity is described by the vector g (m s<sup>-2</sup>).

The effective diffusion coefficient is defined as

$$ D_{\alpha,pm}^{\kappa} := \phi S_{\alpha} D_{\alpha}^{\kappa,j} \tau. \quad (3) $$

Here, τ [-] is the tortuosity factor and D<sub>α</sub><sup>κ,j</sup> (m<sup>2</sup>s<sup>-1</sup>) is the binary diffusion coefficient (of components κ and j in phase α). Because of the low advection velocities and the high diffusion coefficient in the gaseous phase, dispersion is not considered in the model concepts.

**DuMu<sup>x</sup>—Energy balance**

For energy transport, we assume a local thermodynamic equilibrium between the phases. Thus, for all phases, the energy balance is written as follows:

$$ \sum_{\alpha \in \{l,g\}} \left( \phi \frac{\partial (\rho_{\alpha} S_{\alpha} u_{\alpha})}{\partial t} + \nabla \cdot (\rho_{\alpha} h_{\alpha} \mathbf{v}_{\alpha}) \right) + (1 - \phi) \frac{\partial (\rho_{s} c_{s} T)}{\partial t} - \nabla \cdot (\lambda_{pm} \nabla T) = 0. \quad (4) $$

The energy transport is described by the specific internal energy u<sub>α</sub> (J kg<sup>-1</sup>), the specific enthalpy h<sub>α</sub> (J kg<sup>-1</sup>), and the
temperature T [K]. Further, to account for the energy storage of the solid phase, we use the solid heat capacity $c_s$ (J kg$^{-1}$ K$^{-1}$) and the solid density $\rho_s$ (kg m$^3$).

The heat conductivity is defined by using an effective thermal conductivity $\lambda_{pm}$ (W m$^{-1}$ K$^{-1}$) to account for the changes in thermal conductivity based on the saturation of the porous medium (Somerton et al., 1974):

$$ \lambda_{pm} := \lambda_{\text{dry}} + \sqrt{(S_l)(\lambda_{\text{wet}} - \lambda_{\text{dry}})}, \quad (5) $$

with $\lambda_{wet} := \lambda_s^{(1-\phi)} \lambda_l^{\phi}$ and $\lambda_{dry} := \lambda_s^{(1-\phi)} \lambda_g^{\phi}$, where $\lambda_s$, $\lambda_g$, and $\lambda_l$ (W m$^{-1}$ K$^{-1}$) are the solid phase, gas phase, and liquid phase conductivities, respectively.

### SiSPAT-Isotope—Mass balance

The SiSPAT-Isotope model simulates the transport of heat, water, and its heavier stable isotopolog within the soil-plant-atmosphere continuum making a few simplifications compared to the full two-phase four-component approach used in DuMu$^x$. Coupled heat and water transfer equations within the soil, including the liquid and vapor transfer, are solved with the approach of Philip and de Vries (1957), modified by Milly (1982). SiSPAT-Isotope considers only vertical fluxes. The pressure in the gas phase is assumed to be constant and not influenced by flow in the gas or liquid phase so that the flow in the liquid phase is not coupled to flow in the gas phase. Considering only vertical fluxes, the gas velocity, $v_g$ (ms$^{-1}$) at a certain depth is calculated from the change in volumetric gas content per unit of time that is integrated over depth from that depth to the bottom of the simulation domain assuming a no-flow boundary condition for the gas flow at the bottom of the simulation domain. Further, the advection in the liquid phase is not influenced by the pressures of the gas phase, and gradients in the liquid phase density are assumed to be negligible. The approach does not further consider a coupling between the transport of the water components and the dry air component, which is based on the approximation that the dry air component in the liquid phase can be neglected $\sum_{K \neq a} X_l^K \approx 1$. Also, we assume that in Equation (6), all water isotopologs are combined as one component $X_\alpha^{all_{H_2O}}$. Assuming these simplifications for the multiphase flow equation, the following equation can be named as a 1-component, 1.5-phase equation, as water is regarded as the only component and (Vanderborght et al. 2017):

$$ \phi \frac{\partial \rho_g X_g^{all_{H_2O}} S_g}{\partial t} + \phi \frac{\partial \rho_l X_l^{all_{H_2O}} S_l}{\partial t} + \frac{\partial}{\partial z} (v_l \rho_l X_l^{all_{H_2O}}) + \frac{\partial}{\partial z} (v_g \rho_g X_g^{all_{H_2O}}) - \frac{\partial}{\partial z} \left[ \rho_g D_{g.pm}^{all_{H_2O}} \frac{\partial X_g^{all_{H_2O}}}{\partial z} \right] = 0. \quad (6) $$

### SiSPAT-Isotope—Energy balance

The basic heat equation used in SiSPAT-Isotope is taken from de Vries (1958):

$$ \frac{\partial}{\partial z} \left[ (D_{g.pm}^{H_2O} L_v) \frac{\partial \rho_g X_g^{H_2O}}{\partial z} + \lambda_{pm} \frac{\partial T}{\partial z} - v_g \rho_g c_g T - v_l \rho_l c_l T \right] - (\phi S_g c_g \rho_g + \phi S_l c_l \rho_l + (1-\phi) c_s \rho_s) \frac{\partial T}{\partial t} - L_v \frac{\partial \rho_g X_g^{H_2O}}{\partial t} = 0, \quad (7) $$

where $c_\alpha$ (J kg$^{-1}$ K$^{-1}$) is the heat capacity of phase $\alpha$ and $L_v$ (J kg$^{-1}$) is the latent heat of vaporization of liquid water. Compared to Equation (4), the heat used during gas phase expansion is neglected.

### 2.1.2 | Constitutive relations for modeling of drying within the porous medium domain

In the unsaturated zone at equilibrium conditions, a pressure difference—capillary pressure, $p_c$ [Pa]—between the fluid phase (wetting, $p_w$ [Pa]) and gas phase (non-wetting phase, $p_n$ [Pa]), is observed:

$$ p_c = p_n - p_w. \quad (8) $$

On the REV scale, capillary pressure depends on water saturation. A common method to describe the dependence of capillary pressure on saturation for two-phase systems is the van Genuchten model (van Genuchten, 1980). Here, the wetting phase saturation is expressed as a function of capillary pressure using empirically determined parameters:

$$ S_{we} := \left[ 1 + (\alpha_{vG} p_c)^{n_{vG}} \right]^{-m_{vG}}, \quad (9) $$

with $m_{vG} := 1 - \frac{1}{n_{vG}}$ and $S_{we} := \frac{S_w - S_{wr}}{1 - S_{wr} - S_{nr}}$. $\alpha_{vG}$, $m_{vG}$, and $n_{vG}$ are empirical parameters describing soil properties. Residual water saturation $S_{wr}$ and residual gas saturation $S_{nr}$ are empirical parameters expressing the saturation threshold where respectively the wetting or non-wetting phase becomes immobile. $S_{we}$ denotes the effective wetting saturation. In SiSPAT-Isotope, the van Genuchten model describes the volumetric water content $\theta := \phi S_w$, using the residual water content $\theta_r := \phi S_{wr}$ and the saturated water content $\theta_{sat} := \phi(1 - S_{nr})$.

The relative permeability $k_{r,\alpha}$ likewise can be modeled as a function of saturation (Mualem, 1976; van Genuchten, 1980) with
<figure><img src="image_7.png" /><figcaption>FIGURE 1 Capillary pressure–saturation relationship and its regularization (denoted with the red, dotted line) once residual water saturation is reached. The regularization in DuMu<sup>x</sup> is described by a linear relation between capillary pressure and water saturation between S<sub>l</sub> = 0 and S<sub>l</sub> = S<sub>wr</sub>.</figcaption></figure>

$$
\begin{align}
k_{r,w} &:= (S_{we})^{0.5} \left[ 1 - \left( 1 - S_{we}^{\frac{1}{m_{vG}}} \right)^{m_{vG}} \right]^2, \\
k_{r,n} &:= (1 - S_{we})^{0.5} \left[ 1 - S_{we}^{\frac{1}{m_{vG}}} \right]^{2m_{vG}}. \tag{10}
\end{align}
$$

As mentioned before, DuMu<sup>x</sup> and SiSPAT-Isotope handle very dry soil conditions differently. During drying, before the residual water saturation is reached, both models calculate the relationship between capillary pressure and saturation with the van Genuchten model. Assuming chemical equilibrium in the porous medium, the vapor concentration in the gas phase is related to the capillary pressure using Kelvin's equation p<sup>k</sup><sub>sat,Kelvin</sub> [Pa], and the saturated vapor pressure p<sup>k</sup><sub>sat</sub> (when p<sub>c</sub> = 0) is related to the temperature:

$$
P_{sat,Kelvin}^{K} := P_{sat}^{K} \exp\left(-\frac{p_{c}}{\rho_{l}RT}\right), \quad (11)
$$

where **R** (J mol⁻¹K⁻¹) describes the universal gas constant. If vapor transport is considered, liquid water can be lost via vapor transport and water saturation may decrease to values lower than the residual saturation. As a consequence, the capillary pressure–saturation relation can be extended between S<sub>l</sub> = S<sub>wr</sub> and S<sub>l</sub> = 0. Once the residual water saturation is reached, in DuMu<sup>x</sup>, the capillary pressure–saturation curve is regularized to avoid infinite capillary pressures. This allows the fact that the water saturation can fall below residual water saturation, as drying due to evaporation and transport of water vapor in the gas phase is still possible. Advective transport in the liquid phase is not possible anymore, due to the relative permeability becoming zero once residual water saturation is reached (see Figure 1). Using such regularization methods results in numerical stable solutions, which is discussed by Mosthaf et al. (2014).

In SiSPAT-Isotope, drying is limited by the residual water
saturation. An approach describing the physical behavior

of the capillary pressure–saturation relationship more accu-
rately was also implemented in SiSPAT-Isotope but was
not used in our simulations because of numerical insta-
bilities. In this approach, the capillary pressure–saturation
relation is described by a log-linear relation that was found
to describe empirical observations well (Schneider & Goss,
2012). Other approaches to describe the capillary pressure–
water saturation relation in dry soils can be found in, for
example, Campbell et al. (1992) and Fayer and Simmons
(1995).

This difference in handling of porous media with very
low water saturations not only leads to differences in the
description of evaporation but also has implications for the
fractionation of the isotopologs as demonstrated in the fol-
lowing after introducing the processes leading to fractionation
in general.

## 2.2 | Physical processes of evaporation-driven stable water isotopolog fractionation

2.2.1 | Stable water isotopolog fractionation in general—Equilibrium and kinetic fractionation

Natural stable water isotopologs move in the soil pore space both in the gaseous and liquid phases and exchange continuously between the two phases. Due to a difference in molecular weight of different isotopologs, the ratio of the thermodynamic equilibrium concentrations of an isotopolog in the two phases differs from the ratio of the water molecule concentrations, which leads to a higher abundance of heavier isotopologs in the liquid phase and lower abundance in the gas phase, called equilibrium fractionation. Equilibrium fractionation can be expressed using the fractionation factor α<sup>i</sup><sub>eq</sub> [-], which is equal to the ratio of the molecular isotopic ratios in the gas and liquid phases and related to the ratio of concentrations of the isotopologs in both phases as (Braud et al., 2005):
$$ \alpha_{eq}^{i} = \frac{x_{g}^{i}/x_{g}^{H_{2}O}}{x_{l}^{i}/x_{l}^{H_{2}O}} \approx \frac{\rho_{l}}{\rho_{g}} \frac{C_{g}^{i}}{C_{l}^{i}} := \exp\left(-\left(\frac{a}{T^{2}} + \frac{b}{T} + c\right)\right), \quad (12) $$

where $x_g^i$ are the mole fractions of $\kappa \in \{H_2O, i, air\}$ and $i \in \{^1H^2H^{16}O, ^1H_2^{18}O\}$ in the phase $\alpha \in \{l, g\}$, $\rho_l$ and $\rho_g^{H_2O}$ (kg m⁻³) are the volumetric mass densities of the liquid phase and of the water vapor in the gas phase, respectively, and $C_g^i$ and $C_l^i$ (kg m⁻³) are the concentrations of the isotopolog in the gas and liquid phases, respectively. The relation is approximate for the concentrations and is based on the assumption that the molar and mass fractions of the heavier isotopologs are negligible compared to the mass fraction of H₂O. The coefficients a, b, and c for the isotopologs $^1H^2H^{16}O$ and $^1H_2^{18}O$ are listed in Table A.1.

Due to the different transport behavior of isotopologs, that is, due to different binary molecular diffusion coefficients of isotopologs, vapor transport resulting from a gradient in vapor concentration and molecular diffusion will correspond with a greater flux of lighter than heavier isotopologs. The different binary molecular diffusion coefficients are listed in Table A.2.

This flux difference between lighter and heavier isotopes leads to the so-called kinetic fractionation. When vapor transport in the porous medium is described in the transport models, kinetic fractionation resulting from vapor fluxes in the porous medium emerges from the transport equations. The kinetic fractionation that results from the evaporation of water at the free-flow porous-medium interface is due to a slower molecular diffusive transfer of the heavier isotopologs through the boundary layer in the free flow (see Section 3 for more details).

### 2.2.2 | Simulation of isotopic fractionation in the drying porous medium

Isotopic fractionation in the porous medium is simulated by solving the component transport equations (or mass balance equations) in the different phases. In DuMu<sup>x</sup>, this corresponds to solving mass balances for four components (Equation 1): the components H₂O, air, $^1H_2^{18}O$, and $^1H^2H^{16}O$. In SiSPAT-Isotope, one mass balance equation for all water components in both liquid and gaseous phases is solved (Equation 6). Since H₂O is the main component, it is assumed that the total water balance can be described using transport coefficients (diffusion in gas and liquid phases and relations between capillary pressure and water vapor pressures) of H₂O and is not influenced by the transport of the heavier water-stable isotopologs water components.

The transport of heavier water-stable isotopologs is described in SiSPAT-Isotope using a mass conservation equation with isotopolog-specific transport coefficients (diffusion

coefficients and equilibrium fraction factors) (Braud et al., 2005; Bruckler et al., 1996):

$$ \frac{\partial}{\partial z} \left[ D_{l,g}^{i} \frac{\partial \rho_{l} X_{l}^{i}}{\partial z} - v_{l,g} \rho_{l} X_{l}^{i} \right] - \frac{\partial [\Theta^{i} \rho_{l} X_{l}^{i}]}{\partial t} = 0, \quad (13) $$

$$ \Theta^i = \phi[S_l + (1 - S_l)\beta^i], \quad (14) $$

$$ v_{l,g} = \left( v_l + \beta^i v_g - D_{g,pm}^i \frac{\partial \beta^i}{\partial z} \right), \quad (15) $$

$$ D_{l,g}^{i} = (D_{l,pm}^{i} + D_{g,pm}^{i} \beta^{i}), \quad (16) $$

where $\beta^i = \alpha_{eq}^i \frac{\rho_g^{H_2O}}{\rho_l}$, where $\alpha_{eq}^i$ is the equilibrium fractionation factor as defined in Equation (12) and $\Theta^i$ accounts for variations of $\beta^i$ in time. Equation (13) is the same as Equation (1).

As mentioned before, the difference in the handling of very dry porous media, when water saturation reaches residual water saturation, can have a strong influence on isotopic fractionation, as depicted in Figure 2.

In DuMu<sup>x</sup>, water saturation can drop below residual saturation, which means that more water can evaporate leading to more diffusive transport, kinetic fractionation, and more enrichment of the heavier isotopologs. On the other hand, when all liquid water is evaporated, the capacity term in the isotopolog transport equation (i.e., $\Theta^i$, Equation 14) becomes very small so that the isotopolog concentrations near the free-flow porous-medium interface can rapidly equilibrate with the lower isotopolog concentrations in the free flow.

## 3 | SIMULATION SCENARIO

This section describes the simulation scenario of isotopic fractionation during soil drying. The simulations were set up to reproduce a soil-drying experiment (Rothfuss et al., 2015) in which stable water isotopologs were monitored non-invasively and continuously at different depths in the soil column with a high frequency. The continuous and high-frequency isotopolog dataset allows us to test whether the simulation models can reproduce the observed dynamics of isotopic fractionation during soil drying. To our knowledge, our study is the first that compares numerical simulations of isotopic fractionation in drying soils against such high-frequency measurements.

In the following, we briefly present the experimental setup (Rothfuss et al., 2015) and how we set up the simulations accordingly (see Figure 3).
<figure><img src="image_8.png" /><figcaption>**FIGURE 2** Left: Fractionation processes in the porous medium. Right: Drying and fractionation of isotopologs in DuMuₓ and SiSPAT-Isotope under dry soil conditions.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**FIGURE 3** Schematic description of the experimental setup (left side) and its adaption for numerical simulations (right side). The three-dimensional experiment was reduced to a one-dimensional simulation setup for both simulators. The domain was discretized into 600 cells (0.001 m cell⁻¹) in DuMuₓ (finite volume method) and into 289 nodes in SiSPAT-Isotope (finite difference method). In SiSPAT-Isotope, the discretization grid was additionally refined toward the top and the bottom of the domain (min. 0.000115 m node⁻¹, max. 0.0201 m node⁻¹).</figcaption></figure>
<figure><img src="image_10.png" /></figure>

<table><caption>TABLE 1 Parameter for transport, hydraulic, and thermal properties.</caption><thead><tr><th>Parameter</th><th>Values</th></tr></thead><tbody><tr><td>Tortuosity, τ</td><td>0.67 [-]</td></tr><tr><td>Porosity, φ</td><td>0.44 [-]</td></tr><tr><td>Permeability, K</td><td>1.29 × 10<sup>-11</sup> [m<sup>2</sup>]</td></tr><tr><td>Residual water saturation, S<sub>wr</sub></td><td>0.1364 [-]</td></tr><tr><td>Residual gas saturation, S<sub>nr</sub></td><td>0.1911 [-]</td></tr><tr><td>van Genuchten parameter, n<sub>vG</sub></td><td>5.4 [-]</td></tr><tr><td>van Genuchten parameter, α<sub>vG</sub></td><td>0.00031 [Pa<sup>-1</sup>]</td></tr><tr><td>Solid density, ρ<sub>s</sub></td><td>2650 [kg m<sup>-3</sup>]</td></tr><tr><td>Solid thermal conductivity, λ<sub>s</sub></td><td>2.8 [W m<sup>-1</sup> K<sup>-1</sup>]</td></tr><tr><td>Solid heat capacity, c<sub>s</sub></td><td>507.53 [J kg<sup>-1</sup> K<sup>-1</sup>]</td></tr></tbody></table>

During the experiment, which was carried out in a laboratory, soil temperature, water content, and the concentrations of the isotopologs in the liquid phase (here: $H_2^{18}O$ and $^1H^2H^{16}O$) were measured at different depths (0.01, 0.03, 0.05, 0.07, 0.10, 0.20, 0.40, and 0.60 m) inside an initially water-saturated sand (FH31) column (inside diameter = 0.11 m, height = 0.6 m) over a period of 290 days. Soil temperature was measured in each depth with a type K thermocouple (Greisinger Electronic GmbH). For soil water content measurements, electric conductivity sensors (EC-5, Decagon Devices) were taken. Isotopic analysis of liquid water and water vapor was performed using a cavity ring-down spectrometer (L1102-i, Picarro, Inc.) (Rothfuss et al., 2015).

In addition to the soil measurements, the atmospheric conditions 2 m above the sand column (air temperature, relative humidity, and the concentration of the isotopologs in the gaseous phase) and the evaporation rate were monitored. Air relative humidity and temperature were monitored with a combined relative humidity and temperature sensor (RFT-2, UMS GmbH). For the analysis of isotopologs, laboratory air was sampled passively with a stainless steel tubing at 0.6 m above the soil surface (see Figure 3 “Atmospheric measurements”).

The simulation setup was designed to describe the flow and transport of the components $H_2O$, air, $H_2^{18}O$, and $^1H^2H^{16}O$ during evaporation from a column in one dimension (here: vertical direction, domain height = 0.6 m). The parameters of the transport, hydraulic, and thermal properties are listed in Table 1. The values for porosity, residual water saturation, and residual non-wetting saturation were derived from the experiment by Rothfuss et al. (2015). The other parameter values are taken from experiments conducted by Stingaciu et al. (2010).

The initial conditions were derived from the soil measurements. In Figure 3, the initial profiles for saturation, temperature, and isotopic concentration are displayed. While

temperature and isotopic concentration are initially given as block profiles, saturation is initially defined as a smooth profile. The gas pressure ($p_g = 10^5$ Pa) is kept initially constant over the column height. The boundary conditions are described by Neumann conditions. At the bottom, we apply a no-flow condition and at the top, we use solution-dependent flux conditions. The mass transfer across the soil-atmosphere interface is described by evaporation fluxes $E^\kappa$ (kg m$^{-2}$ s$^{-1}$):

$$E^\kappa = \rho_g D_g^{\kappa, \text{air}} \frac{[X_g^\kappa]^{pm} - [X_g^\kappa]^{ff}}{\delta_z} \quad \kappa \in \{\text{H}_2\text{O}, i\}, \quad (17)$$

and the energy transfer across the soil-atmosphere interface is described by a heat flux $H_g$ (J m$^{-2}$ s$^{-1}$):

$$H_g = \lambda_g \frac{[T]^pm - [T]^{ff}}{\delta_z}. \quad (18)$$

To account for the correct impact of the atmosphere on the isotopologs in the soil profile, the atmospheric conditions (here: $[X_g^\kappa]^{ff}$ and $[T]^{ff}$) are prescribed through the experimental data. The boundary layer thickness $\delta_z$ [m] is kept constant in DuMu$^x$ ($\delta_z = 0.005$ m) and varies with time in SiSPAT-Isotope. In both cases, the measured cumulative evaporation rate fits this value. The derivation of the boundary layer thickness in both simulators is explained in Section A.2 in the Appendix. At the interface, we assume that only molecular diffusion occurs. Therefore, the kinetic fractionation is considered through the isotope-specific binary diffusion coefficient, as presented in Barnes and Allison (1983). The fractionation processes within the porous medium are defined through the component description of the isotopologs. The binary diffusion coefficients and the vapor pressure of the isotopologs (the change in vapor pressure can be also expressed by the equilibrium fractionation factor, which is explained in Section 2.2) have been adjusted using parametrizations from the literature. A detailed overview of these parameters can be found in Section A.1 in the Appendix.

Please note, for the following analysis, the concentration of the isotopologs is evaluated using the commonly used $\delta$-notation that relates the ratio of isotopologs to ordinary water $R_\alpha^i$ to a standard value $R_{V-SMOW}$: $\delta_\alpha^i := \frac{R_\alpha^i - R_{V-SMOW}}{R_{V-SMOW}} \times 1000[\%]$ with $R_\alpha^i = \frac{N_\alpha^i}{N_{H_2O}^i}$ (here: $N^\kappa$ [-] describes the particle number) and $R_{V-SMOW}$ the standard mean ocean water (Gonfiantini (1978)) (here: $R_{V-SMOW}^{1H^2HO} = 0.00015576$ and $R_{V-SMOW}^{H_2^{18}O} = 0.00200520)$. Concerning the $\delta$-notation, the superscript $i$ describes only the heavier atom of the isotopolog instead of the entire molecule (here: $H_2^{18}O$ is expressed by $\delta_\alpha^{18O}$ and $^1H^2H^{16}O$ by $\delta_\alpha^{2H}$).
# 4 RESULTS

## 4.1 Study of the dry zone

Despite the different treatment of the dynamics of the boundary layer thickness and related to it, the aerodynamic resistance in the boundary condition, both models simulated similar cumulative evaporation curves (see Figure A.1b). The impact of the definition of the boundary layer thickness gets more significant regarding the non-cumulative evaporation rates (see Figure A.1a). The increase in the aerodynamic resistance (see Figure A.2) directly after the start of the numerical simulation results in a rapid decrease in the simulated evaporation in SiSPAT-Isotope. When using a constant boundary layer thickness and thus a constant aerodynamic resistance, as it is used in DuMu<sup>x</sup>, a nearly constant evaporation rate is modeled in the first 25 days of the simulation, which corresponds to the so-called phase-I evaporation or energy-limited evaporation. Later, due to transport limitations in soil, the evaporation rate simulated by DuMu<sup>x</sup> drops as well and transits into the so-called stage-II evaporation or diffusion-limited evaporation. In SiSPAT-Isotope, this transport limitation is a result of the increasing aerodynamic resistance. Besides the different treatments of the aerodynamic resistance in the boundary conditions, the models use different approaches dealing with the drying of the soil (see Figure 4).

The water saturation modeled by SiSPAT-Isotope reduces over time until it reaches a constant value that is equal to the residual water saturation. DuMu<sup>x</sup> simulates a complete drying of the top-soil layers. Once the residual water saturation is reached, the saturation degree suddenly drops to zero saturation, for example, at -0.01 m depth after 40 days and at -0.05 m depth after 220 days. In DuMu<sup>x</sup>, we assume that: (i) the hydraulic conductivity can drop to zero, which leads to a stoppage of the liquid flow, and (ii) that the water retention curve is linearly regularized in the regime below the residual water saturation. In the range where the saturation and capillary pressure are linearized, the capillary pressure is still quite low, which leads to a vapor concentration that is close to the saturated vapor concentration. Thus, the remaining residual liquid water can be removed via vapor diffusion. The regularization of the saturation-capillary pressure curve and assuming that no fluid flow could take place, that is, film flow was neglected, led to a rapid decrease in liquid saturation after reaching the residual saturation. Comparing the numerical results with the experimental data, the measured water saturation is higher than the simulated one. Only at -0.01 and -0.05 m depths, the water saturation gradually decreases toward the residual water saturation after 75 and 230 days of the experiment. The interpretation of the sensor measurements should be considered with caution. Since the simulations reproduce the cumulative evaporation well and the simulations close the soil water balance, the measurements seem to overestimate the

<figure><img src="image_11.png" /><figcaption>**FIGURE 4** Water saturation of $S_i$ over time at different soil column depths. The simulated values (DuMu<sup>x</sup> [blue] and SiSPAT-Isotope [orange]) are plotted against experimental observations (gray).</figcaption></figure>

water saturation. Figure 5 shows the simulated and measured ¹H²H¹⁶O concentrations at different depths over time.

In Section A.3 in the Appendix, a similar plot for $H_2^{18}O$ is shown. We focus the discussion on the ¹H²H¹⁶O plots since the same conclusions can be drawn from the $H_2^{18}O$ plots. In the drying soil, an enrichment of the isotopologs toward the soil surface was observed. After reaching a maximum, the isotopolog concentrations decreased slowly, which reflects the intrusion of isotopic-depleted atmospheric air into the soil's top layers. This behavior is modeled by both, SiSPAT-Isotope and DuMu<sup>x</sup>.

Focusing on the concentrations at -0.01 m, DuMu<sup>x</sup> simulates a first local maximum directly after the start of the simulation, which was not observed in the experimental data and not modeled by SiSPAT-Isotope. This local maximum occurred during stage-I evaporation and was also found in simulations described in Kiemle et al. (2023). We discuss the underlying mechanisms that caused this first peak in a sensitivity analysis later. For further discussions, we will refer
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>FIGURE 5 Isotopolog concentrations of ¹H²H¹⁶O in pore water, expressed as δ values δᵢ²H, over time at different soil column depths. The simulated values (DuMuₓ [blue] and SiSPAT-Isotope [orange]) are plotted against experimental observations (dark gray). In depths near the surface, the fluctuations of the atmospheric isotopolog concentrations, δg²H, are shown (note the different y-axis for the atmospheric concentrations).</figcaption></figure>

to the first maximum as the stage-I peak and define the second peak as the stage-II peak.

To demonstrate the impact of the isotopolog concentrations in the atmosphere and their fluctuations on the concentrations in the soil, the measured concentrations in the atmosphere are included in the plots representing the near soil surface. High-frequency fluctuations in the atmospheric isotopolog concentrations were observed. These fluctuations were transmitted to the concentrations in the soil and the coherence between the measured concentration fluctuations in the soil and the atmosphere was strong after 120 days for both isotopologs (see also Figure A.4). DuMuₓ reproduced these measured fluctuations in the soil after 120 days remarkably well, but they were smoothed in the SiSPAT-Isotope simulations. Since the top-soil layer in the SiSPAT-isotope simulations kept a residual water content, the buffer capacity

of this layer to buffer isotopolog concentration fluctuations in the atmosphere and the gas phase at the soil surface was larger than the buffer capacity of the completely dry surface layer in the DuMuₓ simulations. The comparison study shows that both models are able to reproduce the general isotopic behavior from the experiments. However, the models use different assumptions with respect to water transport in the dry top-soil layer and the definition of the boundary layer at the porous medium-free flow interface. The increase of the boundary layer thickness to describe the reduction of the evaporation rate leads to a numerically stable solution. Nevertheless, its use should be avoided when simulating isotopolog enrichment since it leads to an overestimation of the isotopolog enrichment near the soil surface and does not reproduce isotopolog concentration profiles correctly.
<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>FIGURE 6 Impact of residual water saturation $S_{wr}$ on the isotopolog concentrations of (a) $H_2^{18}O$ and (b) $^1H^2H^{16}O$, (c) the water saturation, and (d) the evaporation rate due to the different regularization methods of the drying on soils near the surface (–0.01 m depth) in SiSPAT-Isotope (orange) and DuMu$^x$ (blue). The green lines indicate the time of occurrence of the different isotopic peaks (Case: $S_{wr} = 0$): stage-1 peak (straight; day 16), lowest isotopolog concentration (dashed; day 111 [SiSPAT-Isotope]; day 116 [DuMu$^x$]), and stage-II peak (dash-dotted; day 131 [SiSPAT-Isotope]; day 137 [DuMu$^x$]). Note, that this plot contains values only after 35 h of simulation time and no initial values. This illustration was chosen for better readability due to rapid and strong changes compared to the initial values.</figcaption></figure>

Under field conditions, the evaporative demand is much higher than in laboratory settings leading to comparable fast dry soil conditions after a wetting event. In combination, as laser-based measurement techniques are becoming state-of-the-art in isotopolog monitoring, we anticipate that describing the processes near the surface will be relevant to interpret the observation data. Our simulation showed a sensitivity of the isotopic concentration once the soil was dried below residual saturation, which can be seen in the enrichment behavior and in response to the atmospheric fluctuations. To further study and highlight the impact of this difference in the handling of the dry zone, we next compare and analyze the sensitivity of the isotopolog concentrations to the residual water saturation $S_{wr}$. In Figure 6, a comparison of DuMu$^x$ and SiSPAT-Isotope can be seen for varying values of $S_{wr}$. To highlight the impact of that single parameter, the previous simulation setup was simplified by using constant boundary conditions. The values used for the constant setup are listed in Table A.3. With the change to stable atmospheric boundary conditions, the boundary layer thickness in SiSPAT-Isotope was kept constant. Thus, SiSPAT-Isotope and DuMu$^x$ are using the same boundary conditions. To prevent numerical

instabilities caused by the constant boundary layer thickness, the van Genuchten parameter $n_{VG}$ was changed to $n_{VG} = 1.6$.

With this setup, both models simulate stage-I evaporation before transferring to stage-II (Figure 6d), which further results in the formation of a stage-I peak in both models (Figure 6a,b). As the evaporation rate is transferring into stage-II evaporation, both models simulate a minimum in the isotopolog concentrations followed by an enrichment of the isotopologs until the residual water saturation (SiSPAT-Isotope) or zero saturation (DuMu$^x$) is reached (Figure 6c). While DuMu$^x$ allows further drying until reaching a saturation equal to zero, SiSPAT-Isotope limits its minimum saturation equal to the residual water saturation value resulting in a constant wet phase in soil. The choice of the residual water saturation and the choice of the method for describing the dry zone have an effect on the evaporation rate and on the isotopolog concentration.

Stage-I evaporation rate and the stage-I peak are barely affected by a higher residual water saturation (examined residual saturation range: 0–0.1364). With ongoing drying and the transition into stage-II, the impact of the residual water saturation value becomes relevant. The residual water
saturation determines the availability of mobile water so that a higher residual saturation value corresponds with an earlier end of stage-I evaporation and with an earlier minimum in the isotopolog concentration. The stage I peak and subsequent decline of isotopolog concentrations are a consequence of the very high evaporation rates at the beginning of the simulation that decline rapidly afterward until they reach a steady stage-I evaporation. This decline in evaporation rate is related to the equilibration of the soil-surface temperatures. Initially, the soil column and soil-surface temperatures are equal to the free-flow temperature. The consumption of heat for evaporation at the soil surface leads initially to a decrease in surface temperature, surface soil vapor concentration, and evaporation rate. The lower soil-surface temperatures generate a sensible heat flux from the free flow to the porous medium-free flow interface. Constant surface temperature and evaporation rates are reached when the sensible heat flux compensates for the latent heat flux. The initially high but rapidly decreasing evaporation rates correspond with an initially strong source of heavier isotopologs due to high kinetic fractionation. The decreasing source in combination with drying leads to an increasing air phase diffusion coefficient and more back diffusion of the isotopologs into the soil and to the advection of more depleted air into the column generating a decrease after a first strong rapid increase of the isotopolog concentrations near the soil surface during stage-I evaporation. At the onset of stage-II evaporation, the liquid flow with more isotopically depleted water toward the evaporation surface decreases as water becomes immobile whereas the evaporation and the kinetic fractionation source continue. This leads to an increase in isotopolog concentrations until a stage-II peak is formed. For SiSPAT-Isotope, this corresponds with reaching the residual water saturation level, and in DuMu<sup>x</sup>, with a complete dried-out soil (S<sub>w</sub> = 0). However, there is a marked difference in the dependence of the peak value on S<sub>wr</sub> between the two models. In SiSPAT-Isotope simulations, the peak values decrease more with increasing S<sub>wr</sub> than in DuMu<sup>x</sup> simulations. There are two reasons that can explain this difference. When S<sub>wr</sub> is larger than zero, there is always liquid water present in SiSPAT-Isotope. The liquid phase functions as a buffer to changes in isotope concentrations that are generated by kinetic fractionation. The larger this buffer, the smaller the enrichment will be. Second, in DuMu<sup>x</sup>, the residual water can evaporate. This is an additional source of evaporation at the observation depth, which leads to a stronger isotope enrichment. The slightly higher peaks that are simulated by DuMu<sup>x</sup> for a lower S<sub>wr</sub> are due to the longer stage-I evaporation and the stronger drying with more enrichment. After the stage-II peak is reached at the observation depth, there is no evaporation happening at this depth anymore, the evaporation front moved deeper in the soil column, and isotopolog concentrations decrease due to intrusion of isotope-depleted atmospheric water vapor.

## 4.2 | Study of evaporation front

The localization of the evaporation front in drying soils is one aspect of analyzing stable water isotopologs (Walker et al. (1988)). A small difference in the chemical potential between the isotopologs leads to an enrichment of the heavier isotopologs at the evaporation front. This can be seen in the isotopolog concentration profiles in Figure 7. In Figure 7, the isotopolog concentration profiles that are simulated by DuMu<sup>x</sup> and SiSPAT-Isotope during a drying period are plotted together with the experimental data points. The different stages of the drying period are represented through the liquid phase saturation. In Figure A.5, a detailed description of the isotopolog concentration profiles simulated by DuMu<sup>x</sup> and how they are affected by the different states of the evaporation rate is given.

The results presented in Figure A.5 show three different stages of the drying soil: At day 1, the soil is still partly saturated in all soil layers. The atmospheric water demand can be fulfilled from the upper soil layers. Thus, the evaporation front and the maximum isotopolog concentration are located at the top of the domain. Both numerical models show this behavior.

Day 140 represents the middle of the drying scenario. The soil has dried out further, and the evaporation rate has switched to a stage-II evaporation behavior. In the experimental data and the numerical results derived by DuMu<sup>x</sup>, the evaporation front indicated through a maximum in the isotopolog concentrations is well-defined. In SiSPAT-Isotope, the enrichment in soils is modeled, but not the intrusion of isotope-depleted air in the dried-out zone. Here, the top boundary condition is defined by increasing the aerodynamic resistance. The higher atmospheric resistance ensures that the gas vapor concentration stays near the saturated concentration at the surface. The evaporation is simulated to occur mainly at the soil surface resulting in a maximal isotopolog concentration at the top. Switching to a constant aerodynamic resistance in the top boundary conditions, evaporation is simulated at the moving evaporation front, and SiSPAT-Isotope simulates a similar characteristic isotope profile as simulated by DuMu<sup>x</sup> (see Figure A.3).

Day 289 states the end of the simulations and experiments. The findings from this state are similar to the ones on Day 140 (Figure A.5). However, the offset between the measurements and simulations is larger. The offset can have multiple reasons, for example, long run-time of experiments, correct description of drying soil.

Our results show the importance of defining the boundary conditions to simulate the correct isotopic behavior. If a decreasing evaporation rate is simulated by an increasing boundary layer thickness, the model can simulate the isotopic response in each layer, as demonstrated in Figure 5, but it fails to model the characteristic isotopic profile as depicted in Figure 7.
<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /><figcaption>**FIGURE 7** Comparison of simulated isotopolog concentration profiles (left: $\delta_i^{2H}$, middle: $\delta_i^{18O}$, right: $S_i$) against experimental measurements after (a) 1 day, (b) 140 days, and (c) 289 days of drying using DuMu$^x$ (blue) and SiSPAT-Isotope (orange).</figcaption></figure>

From an experimental perspective, a common measure to evaluate the evaporation front in soils is to equate the evaporation front with the height of the maximum isotopolog concentration (Rothfuss & Javaux, 2017). However, Rothfuss et al. (2015) stated that this measure can be used only at the isotopic steady state. For the transient state, they suggested equating the height of the maximum gradient of the isotopolog concentration with the evaporation front. We used the highly resolved numerical models to analyze the suggestion made by Rothfuss et al. (2015). The high resolution of the numerical simulations allows us to analyze both methods on a detailed scale, but also on a coarser scale within an order of the sampling points of the experiments. Note, as we used the experimental setup on which this statement has been formulated, the numerical analysis is only performed with DuMu$^x$.

The procedure of both methods is quite similar. For a given time step, the sampling point where the maximal enrichment is observed, is used to determine the location of the evaporation front in the common approach. Sampling points were at: -0.01, -0.03, -0.05, -0.07, -0.10, -0.20, -0.40, and -0.60 m. In the gradient approach, the average gradient of isotopolog concentrations between each pair of neighboring observation points is calculated. The depth of the evaporation front is assigned to the layer above the interval where the maximal gradient is observed: -0.01 to -0.03 m: 0 m, -0.03 to -0.05 m: -0.02 m, -0.05 to -0.07 m: -0.04 m, -0.07 to -0.1 m: -0.06 m. A more detailed explanation can be found in Rothfuss et al. (2015).

The comparison of both methods using the numerical simulation results using DuMu$^x$ is displayed in Figure 8. To check on the performance of both methods, we compare the solutions against a refined solution where the resolution matches the resolution of the discretization grid. The analysis of the evaporation front with the fine resolution was carried out for both methods, the common and the alternative. If using a fine resolution, both methods show identical results. As the methods show identical results for the fine resolution, we compared the solutions of the coarse scale to the one of the fine scale. In 53.89% of all data points, the alternative approach is closer to the refined solution than the common approach. The error analysis can be found in Kiemle et al. (2024). Based on these results, we cannot make a suggestion whose approach is more accurate, but we recommend further investigations on this topic to analyze under which conditions which approach is more suitable.

5 | CONCLUSION

In this study, we present a comparison study of two numerical simulators, SiSPAT-Isotope and DuMu$^x$. The description of water vapor movement and energy transport differs slightly
<figure><img src="image_19.png" /><figcaption>FIGURE 8 Evaporation front propagation in soil simulated by DuMu<sup>x</sup> for (a) <sup>1</sup>H<sup>2</sup>H<sup>16</sup>O and (b) H<sub>2</sub><sup>18</sup>O. Two approximation methods for the localization of the evaporation front, (i) alternative approach (blue line) and (ii) common approach (light-blue line), are compared against a refined solution (dark-blue line).</figcaption></figure>

between the models. Especially under very dry soil conditions, for example, occurring at the evaporation front, there are notable differences: While drying and before the soil water saturation reaches residual saturation, both models describe the relationship between capillary pressure and water saturation with the van Genuchten relationship. Once residual water saturation is reached, both models differ in their approach: In DuMu<sup>x</sup>, the capillary pressure–saturation curve is regularized to avoid infinite capillary pressures, while in SiSPAT–Isotope, drying below residual water saturation is not possible.

Another difference between the two models is the treatment of the upper boundary conditions. In DuMu<sup>x</sup>, a constant boundary layer thickness with a constant resistance was considered, whereas in SiSPAT-Isotope, the resistance in the boundary layer was increased to simulate a reduction in evaporation rate. We analyzed how these different methods affect isotopolog fractionation and compared the simulations with high-frequency experimental data.

The high-frequency variations of isotopolog concentrations that were measured in the dried-out soil layer in response to variations in the free air isotopolog concentrations were reproduced very accurately by DuMu<sup>x</sup>, which simulated a complete drying out of the surface layer. The reduction in water content below residual saturation was simulated to occur relatively fast and corresponded with an increase in drying rate that was observed at a certain depth. This is due to the rapid decrease in liquid water content with increasing capillary pressure whereas the vapor concentration remains close to saturation over the range of capillary pressures where the capillary pressure–saturation relation was regularized. Such high-frequency variations in isotopolog concentrations were not simulated by SiSPAT-Isotope in which the soil did not dry out below the residual water content. Neither was an increase

in drying rate at a certain depth simulated but the drying rate decreased over time.

We analyzed the influence of the value of the residual water saturation on the peak enrichment and showed that for lower residual water saturations the simulation results of both simulators match more closely. For higher residual water saturations, the simulated peak concentrations by SiSPAT-Isotope were considerably smaller than the peaks simulated by DuMu<sup>x</sup>. This could be attributed to the presence of liquid water that functions as a buffer to changes in the isotopic concentrations, and to evaporation of residual water that is simulated in DuMu<sup>x</sup>, which leads to a higher enrichment of the isotopologs. However, the comparison with experimental data did not show that a larger simulated enrichment improved the prediction of the measurements considerably. A rapid and complete drying out of the soil surface layer was not observed by the soil moisture sensors in the experiment. However, an increase in drying rate was observed at –0.05 m at the time when the simulations predicted rapid drying. Such an increase in drying rate was also measured at –0.01 m depth but later and much smaller than in the DuMu<sup>x</sup> simulation. To further analyze the impact of the water content (or its absence) in the dried-out soil layer on isotopolog concentrations and their dynamics, alternative descriptions of the capillary pressure–water saturation relation in dry soil and the implementation of relationships that have been developed for this range (e.g., Campbell et al., 1992; Fayer & Simmons, 1995; Schneider & Goss, 2012 should be considered).

Based on the comparison between the experimental data and the simulation results, we can see that DuMu<sup>x</sup> can simulate the vertical isotopolog profiles better than SiSPAT-Isotope. SiSPAT-Isotope did not simulate a downward movement of the isotopolog concentration peak since the reduction in evaporation during soil drying was modeled by
increasing the thickness of the laminar boundary layer. This
kept the evaporation front near the soil surface. To simulate
a downward movement of the evaporation front and of the
isotopolog concentration peak, the reduction in evaporation
should be simulated due to an increased resistance to water
transfer through a drying surface layer.

SiSPAT-Isotope can simulate such effects as well, but the
numerical solution may become unstable, especially in soils
with a very narrow pore size distribution (large *n* value). This
problem especially occurs in uniform sandy materials that
are often used in lab experiments. For field applications, this
might be less of a problem since field soils typically have
much more heterogeneous pore size distributions. Addition-
ally, we showed that either using the maximum isotopolog
concentration or using the maximum isotopic gradient can
determine the location of the evaporation front in soils.

The comparison shows both, the complexity of describing
dry soil conditions in numerical models and the sensitivity
of isotopolog concentrations and fractionation to transport
processes in the dry soil layer and the free flow bound-
ary layer. This sensitivity implies that high-frequency in situ
measurements of isotopolog concentrations can help to gain
more insight into different transport processes in the free
flow boundary layer and upper soil layer controlling the
evaporation process but with different impacts on isotopolog
fractionation. The dependency of isotopolog fractionation
on diffusive transport, that is, kinetic fractionation, could
therefore be used to evaluate the role of airflow fluctuations
(turbulent pumping) and diffusive transport in the dry soil
surface layer and the free flow boundary layer and how they
are influenced by the pore network structure of the porous
medium and the roughness of the porous medium-free flow
interface. Advective processes are expected to play a more
important role under field conditions than under the labora-
tory conditions of the experiment that were analyzed in this
study. Both numerical models should be extended to include
the impact of airflow fluctuations in the dry surface soil layer
and across the porous medium-free flow interface to describe
these processes and to infer their impact on evaporation from
the isotopic behavior at the soil-atmosphere interface. In a
further step, this can be expanded to include a comparison
with data from larger scales and under outdoor conditions, for
example, lysimeters.

ACKNOWLEDGMENTS

This work was funded by the Deutsche Forschungsgemein-
schaft (DFG, German Research Foundation)- Project Number
Satin 2531/14-1. We thank the Deutsche Forschungsgemein-
schaft (DFG, German Research Foundation) for supporting
this work by funding SFB 1313, Project Number 327154368.

CONFLICT OF INTEREST STATEMENT

The authors declare no conflicts of interest.

AUTHOR CONTRIBUTIONS

**Jana Schneider:** Conceptualization; data curation; formal analysis; investigation; methodology; resources; software; validation; visualization; writing—original draft; writing—review and editing. **Stefanie Kiemle:** Conceptualization; formal analysis; investigation; methodology; software; validation; visualization; writing—original draft; writing—review and editing. **Katharina Heck:** Formal analysis; investigation; methodology; software; validation; visualization; writing—original draft; writing—review and editing. **Youri Rothfuss:** Data curation; investigation; methodology; resources; supervision; writing—review and editing. **Isabelle Braud:** Data curation; investigation; methodology; resources; software; supervision; validation; visualization; writing—review and editing. **Rainer Helmig:** Investigation; methodology; software; supervision; validation; writing—review and editing. **Jan Vanderborght:** Conceptualization; investigation; methodology; resources; supervision; validation; writing—original draft; writing—review and editing.

DATA AVAILABILITY STATEMENT

Concerning DuMu<sup>x</sup>, all code relevant to obtaining the numer-
ical examples is implemented in DuMu<sup>x</sup> (Koch et al., 2021
and can be found under Gitlab (https://git.iws.uni-stuttgart.de/
dumux-pub/kiemle2023a). Additionally, the code is available
in the data repository of the University of Stuttgart (DaRUS)
(Schneider et al., 2020) and can be accessed via the following
source: Kiemle and Heck (2023a).

Concerning SiSPAT-Isotope, the relevant code to repro-
duce the numerical examples is also available in the
data repository of the University of Stuttgart (DaRUS)
and can be accessed via the following source: Schneider
(2023).

The numerical data sets generated by DuMu<sup>x</sup> and SiSPAT-
Isotope and the experimental data from Rothfuss et al. (2015);
plus, the scripts to reproduce the displayed figures are also in
the data repository of the University of Stuttgart (DaRUS) and
can be accessed via the following source: Kiemle et al. (2024).

ORCID

Jana Schneider https://orcid.org/0000-0002-3644-5146

Jan Vanderborght https://orcid.org/0000-0001-7381-3211

REFERENCES

Anderson, R. G., Zhang, X., & Skaggs, T. H. (2017). Measurement and partitioning of evapotranspiration for application to vadose zone studies. *Vadose Zone Journal*, *16*(13), 1–9. https://doi.org/10.2136/ vzj2017.08.0155

Barnes, C., & Allison, G. (1983). The distribution of deuterium and 18O in dry soils: 1. Theory. *Journal of Hydrology*, *60*(1-4), 141–156.
https://doi.org/10.1016/0022-1694(83)90078-1

Boutraa, T., Akhkha, A., Al-Shoaibi, A., & Alhejeli, A. (2010). Effect of water stress on growth and water use efficiency (WUE) of some wheat cultivars (*Triticum durum*) grown in Saudi Arabia. *Journal*
<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

of Taibah University for Science, 3, 39–48. https://doi.org/10.1016/
S1658-3655(12)60019-3

Braud, I., Bariac, T., Biron, P., & Vauclin, M. (2009). Isotopic com-
position of bare soil evaporated water vapor. Part II: Modeling of
RUBIC IV experimental results. *Journal of Hydrology*, **369**(1-2),
17-29. https://doi.org/10.1016/j.jhydrol.2009.01.038

Braud, I., Bariac, T., Gaudet, J. P., & Vauclin, M. (2005). SiSPAT-Isotope, a coupled heat, water and stable isotope (HDO and H218O) transport model for bare soil. Part I. Model description and first verifications. *Journal of Hydrology*, **309**(1-4), 277–300. https://doi.org/10.1016/j.jhydrol.2004.12.013

Bruckler, L., Bariac, T., & Melayah, A. (1996). Modeling the transport of water stable isotopes in unsaturated soils under natural conditions. 1. Theory. *Water Resources Research*, 32(7), 2047–2054. https://doi.org/10.1029/96WR00674

Campbell, G., & Shiozawa, S. (1992). Prediction of hydraulic proper-
ties of soils using particle-size distribution and bulk density data. In
*Proceedings of International Workshop on Indirect Methods for Esti-
mating the Hydraulic Properties of Unsaturated Soils* (pp. 317–328).
University of California.

de Vries, D. (1958). Simultaneous transfer of heat and moisture in porous media. *Eos, Transactions American Geophysical Union*, *39*(5), 909–916. https://doi.org/10.1029/TR039i005p00909

Dubbert, M., & Werner, C. (2019). Water fluxes mediated by vegeta-
tion: Emerging isotopic insights at the soil and atmosphere interfaces.
*New Phytologist*, **221**(4), 1754–1763. https://doi.org/10.1111/nph.
15547

Fayer, M. J., & Simmons, C. S. (1995). Modified soil water retention functions for all matric suctions. *Water Resources Research*, *31*(5), 1233–1238.

Ghumman, A. R., Ghazaw, Y. M., Alodah, A., Rauf, A., Shafiquzzaman,
M., & Haider, H. (2020). Identification of parameters of evaporation
equations using an optimization technique based on pan evaporation.
*Water*, **12**(1), 228. https://doi.org/10.3390/w12010228

Gonfiantini, R. (1978). Standards for stable isotope measurements in
natural compounds. *Nature*, **271**(5645), 534–536. https://doi.org/10.
1038/271534a0

Haverd, V., & Cuntz, M. (2010). Soil–Litter–Iso: A one-dimensional model for coupled transport of heat, water and stable isotopes in soil with a litter layer and root extraction. *Journal of Hydrology*, **388**(3-4), 438–455. https://doi.org/10.1016/j.jhydrol.2010.05.029

Haverd, V., Cuntz, M., Griffith, D., Keitel, C., Tadros, C., & Twining, J. (2011). Measured deuterium in water vapour concentration does not improve the constraint on the partitioning of evapotranspiration in a tall forest canopy, as estimated using a soil vegetation atmosphere transfer model. *Agricultural and Forest Meteorology*, **151**(6), 645–654. https://doi.org/10.1016/j.agrformet.2011.02.005

Herczeg, A. L., & Leaney, F. (2011). Environmental tracers in arid-zone
hydrology. *Hydrogeology Journal*, **19**(1), 17–29. https://doi.org/10.
1007/s10040-010-0652-7

Hernández-López, M. F., Braud, I., Gironás, J., Suárez, F., & Muñoz, J. F. (2016). Modelling evaporation processes in soils from the Huasco salt flat basin, Chile. *Hydrological Processes*, **30**(25), 4704–4719. https://doi.org/10.1002/hyp.10987

Kiemle, S., & Heck, K. (2023). Dumux code for modelling stable water isotopologue transport within soils using fractionation parameterizations [Data set]. DaRUS. https://doi.org/10.18419/darus-3330

Kiemle, S., Schneider, J., & Heck, K. (2024). Replication data for
analyzing stable water isotopologue transport within soils using frac-

tionation parameterizations [Data set]. DaRUS. https://doi.org/10.
18419/darus-3572

Kiemle, S., Heck, K., Coltman, E., & Helmig, R. (2023). Stable water isotopologue fractionation during soil-water evaporation: Analysis using a coupled soil-atmosphere model. *Water Resources Research*, 59(2), e2022WR032385. https://doi.org/10.1029/2022WR032385

Koch, T., Gläser, D., Weishaupt, K., Ackermann, S., Beck, M., Becker,
B., Burbulla, S., Class, H., Coltman, E., Emmert, S., Fetzer,
T., Grüninger, C., Heck, K., Hommel, J., Kurz, T., Lipp, M.,
Mohammadi, F., Scherrer, S., Schneider, M., ... Flemisch, B. (2021).
DuMux 3—An open-source simulator for solving flow and transport
problems in porous media with a focus on model coupling. *Comput-
ers & Mathematics with Applications*, 81, 423–443. https://doi.org/
10.1016/j.camwa.2020.02.012

Kool, D., Ben-Gal, A., Agam, N., Šimůnek, J., Heitman, J., Sauer, T., & Lazarovitch, N. (2014). Spatial and diurnal below canopy evaporation in a desert vineyard: Measurements and modeling. *Water Resources Research*, 50(8), 7035–7049. https://doi.org/10.1002/2014WR015409

Mathieu, R., & Bariac, T. (1996). An isotopic study (2H and 18O) of water movements in clayey soils under a semiarid climate. *Water Resources Research*, 32(4), 779–789. https://doi.org/10.1029/96WR00074

Merlivat, L. (1978). The dependence of bulk evaporation coefficients on air-water interfacial conditions as determined by the isotopic method. *Journal of Geophysical Research: Oceans*, 83(C6), 2977–2980. https://doi.org/10.1029/JC083iC06p02977

Milly, P. C., Dunne, K. A., & Vecchia, A. V. (2005). Global pattern of trends in streamflow and water availability in a changing climate. *Nature*, **438**(7066), 347–350. https://doi.org/10.1038/nature04312

Milly, P. C. D. (1982). Moisture and heat transport in hysteretic, inhomogeneous porous media: A matric head-based formulation and a numerical model. *Water Resources Research*, 18(3), 489–498. https://doi.org/10.1029/WR018i003p00489

Mosthaf, K., Helmig, R., & Or, D. (2014). Modeling and analysis of evaporation processes from porous media on the REV scale. *Water Resources Research*, 50(2), 1059–1079.

Mualem, Y. (1976). A new model for predicting the hydraulic conduc-
tivity of unsaturated porous media. *Water Resources Research*, **12**(3),
513–522. https://doi.org/10.1029/WR012i003p00513

Oki, T., & Kanae, S. (2006). Global hydrological cycles and world water resources. *Science*, **313**(5790), 1068–1072. https://doi.org/10.1126/science.1128845

Philip, J., & de Vries, D. (1957). Moisture movement in porous materials under temperature gradients. *Eos, Transactions American Geophysical Union*, **38**(2), 222–232. https://doi.org/10.1029/TR038i002p00222

Ross, P. J. (2006). *Fast solution of richards' equation for flexible soil hydraulic property descriptions* (Land and Water Technical Report, vol. 39(06)) CSIRO. https://doi.org/10.4225/08/5859741868a90

Rothfuss, Y., Biron, P., Braud, I., Canale, L., Durand, J.-L., Gaudet, J.-P., Richard, P., Vauclin, M., & Bariac, T. (2010). Partitioning evap-
otranspiration fluxes into soil evaporation and plant transpiration
using water stable isotopes under controlled conditions. *Hydrological
processes*, **24**(22), 3177–3194. https://doi.org/10.1002/hyp.7743

Rothfuss, Y., & Javaux, M. (2017). Reviews and syntheses: Isotopic approaches to quantify root water uptake: A review and comparison of methods. *Biogeosciences*, **14**(8), 2199–2224. https://doi.org/10.5194/
bg-14-2199-2017
Rothfuss, Y., Merz, S., Vanderborght, J., Hermes, N., Weuthen, A., Pohlmeier, A., Vereecken, H., & Brüggemann, N. (2015). Long-term and high-frequency non-destructive monitoring of water stable isotope profiles in an evaporating soil column. *Hydrology and Earth System Sciences*, 19(10), 4067–4080. https://doi.org/10.5194/hess-19-4067-2015

Rothfuss, Y., Quade, M., Brüggemann, N., Graf, A., Vereecken, H., & Dubbert, M. (2021). Reviews and syntheses: Gaining insights into evapotranspiration partitioning with novel isotopic monitoring methods. *Biogeosciences Discussions*, 18(2), 3701–3732. https://doi.org/10.5194/bg-18-3701-2021

Schneider, J. (2023). SiSPAT-isotope code for modelling stable water isotopologue transport within soils. https://doi.org/10.18419/darus-3595

Schneider, M., Flemisch, B., Frey, S., Hermann, S., Iglezakis, D., Ruf, M., Schembera, B., Seeland, A., & Steeb, H. (2020). Daten-management im SFB 1313. *Bausteine Forschungsdatenmanagement*, 2020(1), 28–38. https://doi.org/10.17192/bfdm.2020.1.8085

Schneider, M., & Goss, K.-U. (2012). Prediction of the water sorption isotherm in air dry soils. *Geoderma*, 170, 64–69. https://doi.org/10.1016/j.geoderma.2011.10.008

Somerton, W., Keese, J., & Chu, S. (1974). Thermal behavior of uncon- solidated oil sands. *Society of Petroleum Engineers Journal*, 14(05), 513–521. https://doi.org/10.2118/4506-PA

Stingaciu, L., Weihermüller, L., Haber-Pohlmeier, S., Stapf, S., Vereecken, H., & Pohlmeier, A. (2010). Determination of pore size distribution and hydraulic properties using nuclear magnetic res- onance relaxometry: A comparative study of laboratory methods. *Water Resources Research*, 46(11), W11510. https://doi.org/10.1029/2009WR008686

Stoy, P. C., El-Madany, T. S., Fisher, J. B., Gentine, P., Gerken, T., Good, S. P., Klosterhalfen, A., Liu, S., Miralles, D. G., Perez-Priego, O., & Rigden, A. J. (2019). Reviews and syntheses: Turning the challenges of partitioning ecosystem evaporation and transpiration into opportu- nities. *Biogeosciences*, 16(19), 3747–3775. https://doi.org/10.5194/bg-16-3747-2019

van Genuchten, M. T. (1980). A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. *Soil Science Society of America Journal*, 44(5), 892–898. https://doi.org/10.2136/sssaj1980.03615995004400050002x

Vanderborght, J., Fetzer, T., Mosthaf, K., Smits, K. M., & Helmig, R. (2017). Heat and water transport in soils and across the soil- atmosphere interface: 1. theory and different model concepts. *Water Resources Research*, 53(2), 1057–1079. https://doi.org/10.1002/2016WR019982

Walker, G., Hughes, M., Allison, G., & Barnes, C. (1988). The move- ment of isotopes of water during evaporation from a bare soil surface. *Journal of Hydrology*, 97(3), 181–197. https://doi.org/10.1016/0022-1694(88)90114-X

Wilcox, B. P., Breshears, D. D., Seyfried, M. S., & Trimble, S. (2003). Water balance on rangelands. In *Encyclopedia of water science* (pp. 791–794). CRC Press. https://doi.org/10.1081/E-EWS120010097

Yakir, D., & da SL Sternberg, L. (2000). The use of stable isotopes to study ecosystem gas exchange. *Oecologia*, 123(3), 297–311. https://doi.org/10.1007/s004420051016

Zhou, T., Šimůnek, J., & Braud, I. (2021). Adapting HYDRUS-1D to simulate the transport of soil water isotopes with evaporation fraction- ation. *Environmental Modelling & Software*, 143, 105–118. https://doi.org/10.1016/j.envsoft.2021.105118

**How to cite this article:** Schneider, J., Kiemle, S., Heck, K., Rothfuss, Y., Braud, I., Helmig, R., & Vanderborght, J. (2024). Analysis of experimental and simulation data of evaporation-driven isotopic fractionation in unsaturated porous media. *Vadose Zone Journal*, e20363. https://doi.org/10.1002/vzj2.20363

APPENDIX

A.1 | Component-specific factors

As mentioned before, the component description and, thereof, the fractionation of the isotopologs are dependent on parametrizations. In Table A.1, the values for the equilib- rium fractionation factor are listed. The kinetic fractionation is driven by diffusion in the liquid and the gas phase. The respec- tive diffusion coefficients and the parametrization values are listed in Table A.2.

A.2 | Derivation of boundary layer thickness for the upper boundary flux

The boundary layer thickness is an unknown degree of freedom in the upper boundary conditions. As the total evapo- ration rate is measured in the experiments, the boundary layer thickness can be derived from the measured total evaporation rate and used to define the component evaporation fluxes. The

<table><caption>TABLE A.1 Parametrization values for the equilibrium fractionation factor for isotopologs ¹H²H¹⁶O and H₂¹⁸O.</caption><thead><tr><th>Coefficient</th><th>a</th><th>b</th><th>c</th></tr></thead><tbody><tr><td><sup>1</sup>H<sup>2</sup>H<sup>16</sup>O</td><td>24,844</td><td>-76.248</td><td>0.052612</td></tr><tr><td>H<sub>2</sub><sup>18</sup>O</td><td>1137</td><td>-0.4156</td><td>-0.00206677</td></tr></tbody></table>
<figure><img src="image_22.png" /><figcaption>FIGURE A.1 Resulting total evaporation fluxes for DuMu<sup>x</sup> (blue) and SiSPAT-Isotope (orange) in comparison to experimental data (gray).</figcaption></figure>

<table><caption>TABLE A.2 Binary diffusion coefficients for <sup>1</sup>H<sup>2</sup>H<sup>16</sup>O, H<sub>2</sub><sup>18</sup>O and H<sub>2</sub><sup>16</sup>O.</caption><thead><tr><th>Definition</th><th>Parameterization value</th></tr></thead><tbody><tr><td>Liquid diffusion coefficients</td><td></td></tr><tr><td>Self-diffusion coefficient of H<sub>2</sub>O in the liquid phase</td><td></td></tr><tr><td>D<sub>l</sub><sup>H<sub>2</sub>O,self</sup> = 10<sup>-9</sup>exp(-<sup>535400</sup>/<sub>T<sup>2</sup></sub> + <sup>1393.3</sup>/<sub>T</sub> + 2.1876)</td><td></td></tr><tr><td>Binary diffusion coefficient of H<sub>2</sub>O - i in the liquid phase</td><td>(Mathieu & Bariac, 1996)</td></tr><tr><td>D<sub>i</sub><sup>H<sub>2</sub>O,i</sup> = a<sup>i</sup>D<sub>l</sub><sup>H<sub>2</sub>O,self</sup><br>i ∈ {<sup>1</sup>H<sup>2</sup>H<sup>16</sup>O,<sup>1</sup>H<sub>2</sub><sup>18</sup>O}</td><td>a<sup>H<sub>2</sub><sup>18</sup>O</sup> = 0.9669<br>a<sup>1H<sup>2</sup>H<sup>16</sup>O</sup> = 0.9833</td></tr><tr><td>Gas diffusion coefficients</td><td></td></tr><tr><td>Binary diffusion coefficient of H<sub>2</sub>O - air in the gas phase</td><td></td></tr><tr><td>D<sub>g</sub><sup>air,H<sub>2</sub>O</sup> = 2.17 · 10<sup>-7</sup> (<sup>T</sup>/<sub>273.15</sub>)<sup>1.88</sup></td><td></td></tr><tr><td>Binary diffusion coefficient of i - air in the gas phase</td><td>(Merlivat, 1978)</td></tr><tr><td>D<sub>g</sub><sup>air,i</sup> = <sup>1</sup>/<sub>b<sub>i</sub></sub>D<sub>g</sub><sup>air,H<sub>2</sub>O</sup>,<br>i ∈ {<sup>1</sup>H<sup>2</sup>H<sup>16</sup>O,<sup>1</sup>H<sub>2</sub><sup>18</sup>O}</td><td>b<sup>1H<sup>2</sup>H<sup>16</sup>O</sup> = 1.0251<br>b<sup>H<sub>2</sub><sup>18</sup>O</sup> = 1.0285</td></tr></tbody></table>

derivation of the boundary layer thickness was different for DuMu<sup>x</sup> and SiSPAT-Isotope. The calculated total evaporation rate with the derived boundary layer thickness is compared with the measured evaporation rate in Figure A.1.

In DuMu<sup>x</sup>, a constant boundary layer thickness is set for the entire simulation period. The boundary layer thickness is derived through an interactive process until the deviation of the calculated evaporation rate to the measured evaporation rate is minor.

<figure><img src="image_23.png" /><figcaption>FIGURE A.2 Derivation of the aerodynamic resistance in SiSPAT-Isotope over time. The aerodynamic resistance is either imposed by using the measured evaporation flux (gray lines) or further modified by step-wise linearizing the computed aerodynamic resistance (orange line).</figcaption></figure>

In SiSPAT-Isotope, it is assumed that the boundary layer thickness changed over time. The boundary layer thickness is derived by introducing the aerodynamic resistance, $r_a$ [s m<sup>-1</sup>], which is defined as $r_a := \frac{\delta_z}{D_g^{H_2O,air}}$. The aerodynamic resistance is then derived in two steps. In the first step, the daily measured evaporation flux is imposed to SiSPAT-Isotope, so the aerodynamic resistance is inversely calculated (Figure A.2). In the second step, the aerodynamic resistance over time is described using a step-wise
<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /><figcaption>FIGURE A.3 Effect of different boundary layer conditions in SiSPAT-Isotope on: (a) relative humidity at soil surface level, (b) the $δ^2H$, and (c) $δ^{18}O$ profiles on day 140. Constant (straight line) and fitted (dashed line) boundary layer thickness have been chosen as variants for the boundary layer conditions.</figcaption></figure>

linear function. A detailed description of the optimization of the aerodynamic resistance can be found in Braud et al. (2009) and Hernández-López et al. (2016). Figure A.3 shows how different descriptions of the boundary layer thickness (here: constant and fitted) affect the relative humidity and the isotopolog concentrations in SiSPAT-Isotope (see also Section 4.1).

### A.3 | Supplementary material about the description of the isotopolog concentrations profiles

Figure A.4 supplements the findings in Figure 5. Figure A.5 describes in detail the impact of the phase status of the evaporation rate on the isotopolog concentration profile in soil. Here, we highlighted the enrichment of isotopologs during stage-I evaporation and the transition into stage-II.

### A.4 | Input and boundary condition values for setup change used for analyzing the dry zone

In Table A.3, we summarize the conditions used for analyzing the variations of the residual water saturation $S_{wr}$. Concerning the transport, hydraulic, and thermal property parameters, the values remained as listed in Table 1 beside the van Genuchten parameter $n_{VG}$.

<figure><img src="image_26.png" /><figcaption>FIGURE A.4 Isotopolog concentrations $δ^{18}O$ over time at different depths. The simulated values (DuMu$^x$ [blue] and SiSPAT-Isotope [orange]) are plotted against experimental observations (gray).</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>FIGURE A.5 Simulated isotopolog concentrations (b) $\delta_l^{2H}$ and (c) $\delta_l^{18O}$ by DuMu$^x$ at different phases of (a) the evaporation rate.</figcaption></figure>

<table><caption>TABLE A.3 Initial and boundary conditions for analyzing the impact of the residual water saturation on the dry zone.</caption><thead><tr><th>Parameter</th><th>Values</th></tr></thead><tbody><tr><td colspan="2"><strong>Soil conditions (initial)</strong></td></tr><tr><td>Water saturation, S<sub>l</sub></td><td>0.74 [m<sup>3</sup>m<sup>-3</sup>]</td></tr><tr><td>Temperature, T</td><td>289 [K]</td></tr><tr><td>Isotopolog concentration, δ<sub>l</sub><sup>2H</sup></td><td>-52 [%c]</td></tr><tr><td>Isotopolog concentration, δ<sub>l</sub><sup>18O</sup></td><td>-8 [%c]</td></tr><tr><td>Gas pressure, p<sub>g</sub></td><td>10<sup>5</sup> [Pa]</td></tr><tr><td colspan="2"><strong>Atmospheric conditions (initial and boundary)</strong></td></tr><tr><td>Relative humidity, S<sub>l</sub></td><td>0.4 [-]</td></tr><tr><td>Temperature, T</td><td>289 [K]</td></tr><tr><td>Isotopolog concentration, δ<sub>g</sub><sup>2H</sup></td><td>-100 [%c]</td></tr><tr><td>Isotopolog concentration, δ<sub>g</sub><sup>18O</sup></td><td>-15 [%c]</td></tr><tr><td>Gas pressure, p<sub>g</sub></td><td>10<sup>5</sup> [Pa]</td></tr><tr><td>Boundary layer thickness, δ<sub>z</sub></td><td>0.005 [m]</td></tr><tr><td colspan="2"><strong>Hydraulic parameters</strong></td></tr><tr><td>van Genuchten parameter, n<sub>VG</sub></td><td>1.6 [-]</td></tr></tbody></table>
## A.5 | Nomenclature

### A.5.1 | Latin symbols

<table><thead><tr><th>Symbol</th><th>Unit</th><th>Definition</th></tr></thead><tbody><tr><td>c<sub>α</sub></td><td>J kg<sup>-1</sup> K<sup>-1</sup></td><td>Heat capacity</td></tr><tr><td>C</td><td>kgm<sup>-3</sup></td><td>Concentration</td></tr><tr><td>D<sup>K</sup><sub>a.pm</sub></td><td>m<sup>2</sup>s<sup>-1</sup></td><td>Effective diffusion coefficient</td></tr><tr><td>D<sup>K,j</sup><sub>α</sub></td><td>m<sup>2</sup>s<sup>-1</sup></td><td>Binary diffusion coefficient</td></tr><tr><td>E<sup>κ</sup></td><td>kgm<sup>2</sup>s<sup>-1</sup> or mmh<sup>-1</sup></td><td>Evaporation flux</td></tr><tr><td>g</td><td>ms<sup>-2</sup></td><td>Gravity vector</td></tr><tr><td>h<sub>α</sub></td><td>Jkg<sup>-1</sup></td><td>Specific enthalpy</td></tr><tr><td>H<sub>g</sub></td><td>Jm<sup>-2</sup>s<sup>-1</sup></td><td>Heat flux</td></tr><tr><td>K</td><td>m<sup>2</sup></td><td>Intrinsic permeability tensor</td></tr><tr><td>k<sub>rα</sub></td><td>–</td><td>Relative permeability</td></tr><tr><td>L<sub>v</sub></td><td>Jkg<sup>-1</sup></td><td>Latent heat vaporization</td></tr><tr><td>m</td><td>–</td><td>Empirical parameter</td></tr><tr><td>n</td><td>–</td><td>Empirical parameter</td></tr><tr><td>N<sup>κ</sup></td><td>–</td><td>Particle number</td></tr><tr><td>p<sub>α</sub></td><td>Pa</td><td>Phase pressure</td></tr><tr><td>p<sub>c</sub></td><td>Pa</td><td>Capillary pressure</td></tr><tr><td>p<sub>sat,Kelvin</sub></td><td>Pa</td><td>Kelvin pressure</td></tr><tr><td>q<sup>K</sup></td><td>ms<sup>-1</sup></td><td>Source/sink term</td></tr><tr><td>R</td><td>–</td><td>Isotopic ratio</td></tr><tr><td>R</td><td>Jmol<sup>-1</sup>K<sup>-1</sup></td><td>Universal gas content</td></tr><tr><td>r<sub>α</sub></td><td>sm<sup>-1</sup></td><td>Aerodynamic resistance factor</td></tr><tr><td>S<sub>α</sub></td><td>m<sup>3</sup>m<sup>-3</sup></td><td>Phase saturation</td></tr><tr><td>S<sub>we</sub></td><td>m<sup>3</sup>m<sup>-3</sup></td><td>Effective wetting saturation</td></tr><tr><td>S<sub>nr</sub></td><td>m<sup>3</sup>m<sup>-3</sup></td><td>Residual gas saturation</td></tr><tr><td>S<sub>wr</sub></td><td>m<sup>3</sup>m<sup>-3</sup></td><td>Residual water saturation</td></tr><tr><td>T</td><td>K</td><td>Temperature</td></tr><tr><td>t</td><td>s</td><td>Time</td></tr><tr><td>u<sub>α</sub></td><td>Jkg<sup>-1</sup></td><td>Specific internal energy</td></tr><tr><td>vG</td><td>–</td><td>Van Genuchten</td></tr><tr><td>vGM</td><td>–</td><td>Van Genuchten-Mualem</td></tr><tr><td>v<sub>α</sub></td><td>ms<sup>-1</sup></td><td>Fluid-phase velocity</td></tr><tr><td>v<sub>g</sub></td><td>ms<sup>-1</sup></td><td>Gas flux</td></tr><tr><td>v<sub>l</sub></td><td>ms<sup>-1</sup></td><td>Liquid water flux</td></tr><tr><td>X</td><td>kgkg<sup>-1</sup></td><td>Mass fraction</td></tr><tr><td>x</td><td>molmol<sup>-1</sup></td><td>Mole fraction</td></tr><tr><td>z</td><td>m</td><td>Vertical coordinate</td></tr></tbody></table>

### A.5.2 | Greek symbols

<table><thead><tr><th>Symbol</th><th>Unit</th><th>Definition</th></tr></thead><tbody><tr><td>α<sub>vG</sub></td><td>m<sup>-1</sup></td><td>Empirical parameter</td></tr><tr><td>α<sup>i</sup><sub>eq</sub></td><td>–</td><td>Equilibrium fractionation factor</td></tr><tr><td>β<sup>i</sup></td><td>–</td><td>Factor relating liquid and Vapor isotope concentration</td></tr><tr><td>δ<sub>l</sub></td><td>per mil</td><td>Isotope concentration (delta notation)</td></tr><tr><td>δ<sub>z</sub></td><td>m</td><td>Boundary layer thickness</td></tr><tr><td>λ<sub>pm</sub></td><td>Wm<sup>-1</sup>K<sup>-1</sup></td><td>Effective thermal conductivity</td></tr><tr><td>μ<sub>α</sub></td><td>Pas<sup>-1</sup></td><td>Dynamic viscosity</td></tr><tr><td>φ</td><td>–</td><td>Porosity</td></tr><tr><td>ρ<sub>α</sub></td><td>kgm<sup>-3</sup></td><td>Mass density</td></tr><tr><td>θ</td><td>m<sup>3</sup>m<sup>-3</sup></td><td>Volumetric water content</td></tr><tr><td>θ<sub>sat</sub></td><td>m<sup>3</sup>m<sup>-3</sup></td><td>Saturated volumetric water content</td></tr><tr><td>θ<sub>r</sub></td><td>m<sup>3</sup>m<sup>-3</sup></td><td>Residual volumetric water content</td></tr><tr><td>τ</td><td>–</td><td>Tortuosity factor</td></tr></tbody></table>

### A.5.3 | Indices and exponents

<table><thead><tr><th>Air</th><th>Dry air</th></tr></thead><tbody><tr><td>α</td><td>Phase</td></tr><tr><td>ff</td><td>Free-flow</td></tr><tr><td>g</td><td>Gas phase</td></tr><tr><td>H<sub>2</sub>O</td><td>Water</td></tr><tr><td>i</td><td>Isotopolog</td></tr><tr><td>κ</td><td>Component</td></tr><tr><td>l</td><td>Liquid phase</td></tr><tr><td>nw</td><td>Non wetting phase</td></tr><tr><td>pm</td><td>Porous-medium domain</td></tr><tr><td>s</td><td>Solid</td></tr><tr><td>sat</td><td>Saturated</td></tr><tr><td>w</td><td>Wetting phase</td></tr></tbody></table>