A Fully Coupled Numerical Solution of Water, Vapor, Heat, and Water Stable Isotope Transport in Soil

Abstract

Modeling water stable isotope transport in soil is crucial to sharpen our understanding of water cycles in terrestrial ecosystems. Although several models for soil water isotope transport have been developed, many rely on a semi‐coupled numerical approach, solving isotope transport only after obtaining solutions from water and heat transport equations. However, this approach may increase instability and errors of model. Here, we developed an algorithm that solves one‐dimensional water, heat, and isotope transport equations with a fully coupled method (MOIST). Our results showed that MOIST is more stable under various spatial and temporal discretization than semi‐coupled method and has good agreement with semi‐analytical solutions of isotope transport. We also validated MOIST with long‐term measurements from a lysimeter study under three scenarios with soil hydraulic parameters calibrated by HYDRUS‐1D in the first two scenarios and by MOIST in the last scenario. In scenario 1, MOIST showed an overall NSE, KGE, and MAE of simulated δ¹⁸O of 0.47, 0.58, and 0.92‰, respectively, compared to the 0.31, 0.60, and 1.00‰ from HYDRUS‐1D; In scenario 2, these indices of MOIST were 0.33, 0.52, and 1.04‰, respectively, compared to the 0.19, 0.58, and 1.15‰ from HYDRUS‐1D; In scenario 3, calibrated MOIST exhibited the highest NSE (0.48) and KGE (0.76), the smallest MAE (0.90) among all scenarios. These findings indicate MOIST has better performance in simulating water flow and isotope transport in simplified ecosystems than HYDRUS‐1D, suggesting the great potential of MOIST in furthering our understandings of ecohydrological processes in terrestrial ecosystems.

Introduction

Hydrogen and oxygen stable isotopes are powerful tools to reveal and document water cycles and ecohydrological processes in terrestrial ecosystems (Vereecken et al., 2016). Their applications encompass plant water sourcing (Brooks et al., 2010), evaporation estimation (Walker et al., 1988; Xiang et al., 2021), water age identification (Evaristo et al., 2019; Ferguson et al., 2023; Hoffmann et al., 2004; Jouzel et al., 2007; Wang et al., 2023; Zhang et al., 2019), and evapotranspiration partitioning from plots to continental scales (Evaristo et al., 2015; Vereecken et al., 2016). However, despite significant advances in experimental techniques and in-situ sampling of water isotopes in terrestrial ecosystems (Jasechko et al., 2013), the related modeling work is still limited (Zhou et al., 2021). This limits our ability to accurately simulate isotope transport in soil-vegetation-atmosphere systems and to sharpen our understanding of evaporation and transpiration on an ecosystem scale.

Several efforts have been undertaken to simulate isotope transport in soil across various scales. At the catchment scale, notable models include Ech2o-iso (Kuppel et al., 2018a), UMSWIM (Smith et al., 2016), JoFlow (Knighton et al., 2017), and StaRR (Ala-Aho et al., 2017) have successfully applied in addressing key ecohydrological issues such as plant water uptake (Li et al., 2023), water age tracking (Smith et al., 2016), and evapotranspiration separation (Bao et al., 2021; Dubbert et al., 2013). However, most of them employed a simple linear mixing theory to model isotopic signals in soil water. This simplification facilitates computations, making the models easily applicable at large scales, but may overlook the influence of physical mechanisms on isotope transport, such as phase changes and vapor diffusion.

Physical processes associated with isotope transport within soil receive greater attention at smaller scales, such as one-dimensional isotope movement. Based on Melayah et al. (1996), Braud et al. (2005) developed the “SiSPAT-Isotope” (Simple Soil-Plant-Atmosphere Transfer) model, which incorporated the resistance to isotope transport between the soil surface and atmosphere (Braud, et al., 2009a, 2009b). Subsequently, Haverd and Cuntz (2010) developed “Soil-Litter-Iso,” which is also a one-dimensional model for the transport of heat, water, and stable isotopes in soil containing surface litter and actively transpiring vegetative cover. This model extended the linearization method of Ross (2003) to vapor transport and emphasized the importance of isotope vapor transport, which is later supported by Sprenger et al. (2018) using a SWIS model (Mueller et al., 2014). Compared to SiSPAT-Isotope, Soil-Litter-Iso is more efficient for thicker soil layers and larger time steps because of the implementation of the special linearization scheme in the model framework. However, Soil-Litter-Iso does not consider liquid and vapor heat capacity variation and the changes in vapor volume during heat transport, which may substantially bias the estimation of heat flux transport within the soil (Satio et al., 2006), thereby resulting in biased isotope transport fluxes.

Another approach capitalizes on the capability of the well-known model, HYDRUS-1D, for simulating water flow and chemical transport. Stumpp et al. (2012) simulated oxygen-18 movement in soil by modifying a solute transport module of HYDRUS-1D. However, Stumpp et al. (2012) originally neglected fractionation and thus, their model was only applicable in situations where the slope of the isotopic evaporation line from soil water was close to the local meteoric water line (Stumpp et al., 2012). The latest version of HYDRUS-1D incorporated isotopic fractionation by integrating with the SiSPAT-Isotope model (Zhou et al., 2021).

A common problem of these models related to isotope simulation is that each of these modeling approaches considers water and isotope transport separately, where the isotope transport is solved after solving the water and heat transport at each time step (semi-coupled method). This is reasonable because isotope transport does not affect soil water and heat transport. However, numerical errors and instabilities from soil, water and heat transport equations can be transferred into isotope transport equations at each time step (step 2 and three in Figure 1). Thus, either a tiny spatial discretization at the soil surface (10⁻⁶ m, Braud et al., 2005) for extremely small mass balance errors (in the magnitude of 10⁻¹⁶, Zhou et al., 2021), or complicated discretization schemes (Haverd & Cuntz, 2010) were required to protect the isotope solutions from numerical oscillation. We hypothesize that steps 2 and 3 from the semi-coupled method can be incorporated into step 1 in a fully coupled method (Figure 1) and thus, the fully coupled method reduces the impact of mass balance errors from soil water and heat transport equations on the isotope transport equation (step 2, Figure 1) compared to the semi-coupled method. Consequently, the fully coupled method should be more accurate and less affected by numerical instability than semi-coupled methods.

[IMAGE OMITTED. SEE PDF]

Therefore, the objectives of this study are: (a) to develop a one-dimensional model which can solve fully coupled soil water, heat, vapor, and isotope transport simultaneously and (b) to validate MOIST through analytical simulations with specific boundary conditions and field measurements. MOIST was written in MATLAB programming language and it is expected to be more efficient for simulating isotope transport within soil under various spatial and temporal scales than existing models by solving fully coupled soil water, heat, vapor, and isotope transport equations.

Materials and Methods

Model Description

Soil Water and Heat Transport

One-dimensional water and heat transport within soil can be described by mass and energy conservation equations with a downward-positive coordinate system (Farthing & Ogden, 2017; Saito et al., 2006): 1 $$\frac{\mathit{\partial }\left(\theta +{\theta }_{v}\right)\,}{\mathit{\partial }t}=-\frac{\mathit{\partial }\left({q}_{l}+{q}_{v}\right)}{\mathit{\partial }z}-{S}_{p}$$ 2 $$\frac{\mathit{\partial }{C}_{soil}T}{\mathit{\partial }t}+{\lambda }_{E}\frac{\mathit{\partial }{\theta }_{v}}{\mathit{\partial }t}=\frac{\mathit{\partial }\left({K}_{H}\frac{\mathit{\partial }T}{\mathit{\partial }z}\right)}{\mathit{\partial }z}-{C}_{w}\frac{\mathit{\partial }\left({q}_{l}T\right)}{\mathit{\partial }z}-{C}_{vh}\frac{\mathit{\partial }\left(T{q}_{v}\right)}{\mathit{\partial }z}-{\rho \lambda }_{E}\frac{\mathit{\partial }{q}_{v}}{\mathit{\partial }z}-{C}_{w}{S}_{p}T$$ where θ and θ_v are the volumetric soil water (m³ m⁻³) and vapor content (m³ m⁻³); t is the time (s); q_l is the liquid flux (m s⁻¹) and q_v is the vapor flux (m s⁻¹); z is the spatial distance (m); S_p is the sink term (s⁻¹), which is zero when there is no sink during simulation (e.g., root water uptake); C_soil is the soil volumetric heat capacity (J m⁻³ K⁻¹); T is the soil temperature (K); λ_E is the latent heat of vapourization (J kg⁻¹); K_H is the soil thermal conductivity (W m⁻¹ K⁻¹); C_w and C_vh are the heat capacities of liquid water (J m⁻³ K⁻¹) and vapor (J m⁻³ K⁻¹); and ρ is the water density (kg m⁻³). Detailed explanations of these variables can be found in the Supporting Information S1.

The liquid flux, q_l, is calculated by Buckingham-Darcy's law (positive downwards) (Haverd & Cuntz, 2010): 3 $${q}_{l}=-K(h)\frac{dh}{dz}+K(h)$$ where K(h) is the unsaturated hydraulic conductivity (m s⁻¹) as a function of h. The vapor flux, q_v, is calculated by the Fick's law (Shokri et al., 2009): 4 $${q}_{v}=-{D}_{vs}\frac{\mathit{\partial }{c}_{v}}{\mathit{\partial }z}$$ where D_vs is the diffusivity of vapor within soil (m² s⁻¹). Note that $${c}_{v}$$ is the function of h and T. Thus, the thermal vapor conduction is implicitly considered. Detailed derivations can be found in the Supporting Information S1.

Note that Equation 1 is the mixed form of Richards' equation. However, the current study uses the head-based Richards' equation: 5 $$\left(1-{c}_{v,sat}{R}_{H}\right)\frac{\mathit{\partial }\theta }{\mathit{\partial }h}\frac{\mathit{\partial }h}{\mathit{\partial }t}+\left({\theta }_{sat}-\theta \right)\frac{\mathit{\partial }\left({c}_{v,sat}{R}_{H}\right)}{\mathit{\partial }t}=-\frac{\mathit{\partial }\left({q}_{l}+{q}_{v}\right)}{\mathit{\partial }z}-{S}_{p}$$ where $${c}_{v,sat}$$ is the saturated vapor concentration within soil air (m³_{liquid water} m⁻³_air); and R_H is the relative humidity in soil. Equations for c_v, D_vs, $${c}_{v,sat}$$, and R_H can be found in the Supporting Information S1.

Both the head-based and the mixed forms of Richards' equation are commonly used for modeling water flow in porous media. The current study utilizes the head-based form of Richards' equation because the head is continuous across the soil interface under typical natural conditions, especially for the layered soils (Zha et al., 2019). However, $$\frac{\mathit{\partial }\theta }{\mathit{\partial }h}\frac{\mathit{\partial }h}{\mathit{\partial }t}$$ may not be numerically equal to $$\frac{\mathit{\partial }\theta }{\mathit{\partial }t}$$ (Celia et al., 1990; Clark et al., 2021) because $$\frac{\mathit{\partial }\theta }{\mathit{\partial }h}$$ could introduce mass balance errors during linearization. Nevertheless, the mass balance of head-based Richards' equation can be significantly improved by a second-order approximation to the time derivative (Celia et al., 1990) and effectively controlled by adaptive time-stepping schemes (Ireson et al., 2023). MATLAB provides well-tested ordinary differential equations solvers, which implement adaptive time stepping schemes. As such, these solvers effectively meet the criteria.

Soil Water Isotope Transport

The abundance of oxygen and hydrogen stable isotopes is conventionally expressed as δ values, in units of per mil (‰). However, for convenience, the abundances in MOIST were presented as concentrations (kg m⁻³), $${c}_{i}$$, where the relationship between isotopic ratio ($${R}_{i}$$) and concentration ($${c}_{i}$$) can be expressed as (Braud et al., 2005; Melayah et al., 1996): 6 $${c}_{i}=\frac{{M}_{i}}{{M}_{w}}{R}_{i}\rho $$ where M_i and R_i are the molar mass (kg mol⁻¹) and the isotopic ratio of isotopic species i, respectively, for a given isotope species; R_i can be calculated by: 7 $${R}_{i}=\delta \frac{{R}_{ref}}{1000}+{R}_{ref}$$ where $${R}_{ref}$$ is the reference value for the isotopic ratio, which is 155.76 × 10⁻⁶ and 2,005.2 × 10⁻⁶ for HDO and H₂¹⁸O, respectively (Gonfiantini, 1978).

The isotope mass conservation equation for both liquid and vapor phases is: 8 $$\frac{\mathit{\partial }\left({c}_{i}^{l}\theta +{c}_{i}^{v}\left({\theta }_{sat}-\theta \right)\right)}{\mathit{\partial }t}=-\frac{\mathit{\partial }{q}_{i}}{\mathit{\partial }z}-{c}_{i}^{l}{S}_{p}$$ where $${c}_{i}^{l}$$ and $${c}_{i}^{v}$$ are the concentrations of isotope species i in liquid (kg m⁻³) and vapor phases (kg m⁻³). Assuming an instantaneous equilibrium between liquid and vapor phases (Vanderborght et al., 2017), a relationship between liquid and vapor isotopic concentration can be expressed as (Braud et al., 2005; Haverd & Cuntz, 2010): 9 $${c}_{i}^{v}={c}_{v}\,{\alpha }_{i}^{\ast }{c}_{i}^{l}$$ where $${\alpha }_{i}^{\ast }$$ is the equilibrium fractionation coefficient of isotope species i, which can be written as (Braud et al., 2005): 10 $${\alpha }_{i}^{\ast }={e}^{\left(\mathit{-}\tfrac{{a}^{\ast }}{{(T+273.15)}^{2}}-\tfrac{{b}^{\ast }}{T+273.15}\mathit{-}{c}^{\ast }\right)}$$ where $${a}^{\ast }$$, $${b}^{\ast }$$, $${c}^{\ast }$$ have values of 24,844, −76.248, and 0.052612 for ²H, and 1,137, −0.4156, and −0.0020667 for ¹⁸O.

The total isotopic flux, q_i (positive downwards), consists of liquid isotopic flux, $${q}_{i}^{l}$$, and vapor isotopic flux, $${q}_{i}^{v}$$: 11 $${q}_{i}^{l}={c}_{i}^{l}{q}_{l}-{D}_{i,s}^{l}\frac{\mathit{\partial }{c}_{i}^{l}}{\mathit{\partial }z}$$ 12 $${q}_{i}^{v}=-{D}_{i,s}^{v}\frac{\mathit{\partial }{c}_{i}^{v}}{\mathit{\partial }z}$$ where $${D}_{i,s}^{l}$$ and $${D}_{i,s}^{v}$$ are the liquid (m² s⁻¹) and vapor (m² s⁻¹) diffusivity in soil for isotope species i, respectively. Detailed equations of $${D}_{i,s}^{l}$$ and $${D}_{i,s}^{v}$$ can be found in the Supporting Information S1.

Combing Equations 4 and 9 (Equation S3–10 in Supporting Information S1), and (Equation S3–11 in Supporting Information S1) with Equation 12 leads to (Haverd & Cuntz, 2010): 13 $${q}_{i}^{v}=-{\left({\alpha }_{diff}\right)}^{nD}\hspace*{.5em}{\alpha }_{i}^{\ast }{c}_{i}^{l}{q}_{v}-{D}_{i,s}^{v}{c}_{v}\left({\alpha }_{i}^{\ast }\frac{\mathit{\partial }{c}_{i}^{l}}{\mathit{\partial }z}+{c}_{i}^{l}\frac{\mathit{\partial }{\alpha }_{i}^{\ast }}{\mathit{\partial }z}\right)$$

Equation 13 illustrates that the isotope transport in vapor phase consists of convection and diffusion terms, where the latter can be expressed as a concentration gradient in the liquid phase (first term within the brackets), and an equilibrium fractionation coefficient gradient (second term within the brackets). Detailed derivations can be found in Supporting Information S1.

Root Water Uptake

Root water uptake is modeled as a sink term in Equations 1 and 8 (Li et al., 2001): 14 $${S}_{p}=\frac{{\Omega }{F}_{i}{P}_{t}}{{\Delta }z}$$ where $${\Omega }$$ is the efficiency coefficient (between 0 and 1), which can be obtained from a prescribed stress function; Δz is the thickness of the considered layer (m) and P_t is the potential transpiration rate (m s⁻¹), which can be obtained from the Penman-Monteith equation (Allen et al., 1999) and the leaf area index. The current version of MOIST does not account for isotope fractionation during root water uptake. Root water uptake in MOIST is described by the Feddes' model (Feddes et al., 1978). Essentially, potential evapotranspiration (pET) is calculated using the Penman-Monteith equation (Allen et al., 1999). Following this, pET partitioned into potential evaporation and potential transpiration based on leaf area index (Šimůnek et al., 2013). The potential transpiration is then distributed across individual soil layers by multiplying root proportion (F_i) of each layer. Ultimately, the actual root water uptake flux from each layer is determined by multiplying the allocated potential root water uptake flux by the stress factor (Equations 14 and 15, Figure 2).

[IMAGE OMITTED. SEE PDF]

As demonstrated by Figure 2, $${\Omega }$$ is a function of pressure head, h, and can be modeled with four critical pressure heads (h₁, h₂, h₃, h₄) that may vary among different plant species: 15 $${\Omega }=\left\{\begin{array}{@{}c@{}}0,h\,< \,{h}_{4}\,\text{or}\,h\, > \,{h}_{1}\\ \frac{h-{h}_{4}}{{h}_{3}-{h}_{4}},{h}_{4}\mathit{\le }h\mathit{\le }{h}_{3}\\ \frac{h-{h}_{1}}{{h}_{2}-{h}_{1}},{h}_{2}\mathit{\le }h\mathit{\le }{h}_{1}\\ 1,{h}_{3}\mathit{\le }h\mathit{\le }{h}_{2}\end{array}\right.$$

F_i is the fraction of root length in the depth interval of z to z + ∆z, defined as: 16 $${F}_{i}=\frac{\int \nolimits_{z}^{z+\mathit{{\increment}}z}R(x)dx}{\int \nolimits_{0}^{Zr}R(x)dx}$$ where z is the depth (m) of the upper boundary of layer i; z + $$\mathit{{\increment}}$$ z is the depth of the lower boundary (m); Z_r is the maximum depth of the soil column (m); and R(x) is a predefined root length density distribution function which can be obtained by fitting to fine root length density measurements. For convenience, here we define R(z) as (Dwyer et al., 1988): 17 $$R(z)=R{D}_{0}\left[1-\frac{1}{1+{e}^{-bz}}+\frac{{e}^{-bz}}{2}\right]$$ where $$R{D}_{0}$$ is root length density at z = 0 and it is canceled out when calculating F_i; b is an empirical root distribution parameter (m⁻¹), and the equation was given by Li et al. (2001).

Note that R(z) can be static or treated as dynamic with $${Z}_{r}$$ being described as the classical Verhulst-Pearl logistic growth function (Kot, 2001): 18 $${Z}_{r}(t)={L}_{m}\frac{{L}_{0}}{{L}_{0}+\left({L}_{m}-{L}_{0}\right){e}^{-rt}}$$ where L₀ is the initial value of the rooting depth at the beginning of the growing season (m); L_m is the maximum rooting depth during the growing season (m); r is the growth rate (day⁻¹) (Šimůnek et al., 2013).

Boundary Conditions

Boundary Conditions for Soil Water Flow

According to the Soil-Litter-Iso model (Haverd & Cuntz, 2010), the coupled energy and water mass conservation equations were solved at the soil-air interface: 19 $$\frac{1}{{r}_{bw}}\left[{c}_{vs}\left({T}_{s},{R}_{Hs}\right)-{c}_{va}\right]=-{D}_{v}\frac{\mathit{\partial }{{c}_{v,sat}R}_{H}}{\mathit{\partial }z}-K\frac{\mathit{\partial }h}{\mathit{\partial }z}+K$$ 20 $${R}_{net}=\frac{{C}_{p}}{{r}_{bh}}\left({T}_{s}-{T}_{a}\right)+\frac{\rho {\lambda }_{E}}{{r}_{bw}}\left[{c}_{vs}\left({T}_{s},{R}_{Hs}\right)-{c}_{va}\right]-{K}_{H}\frac{\mathit{\partial }T}{\mathit{\partial }z}$$ where r_bw and r_bh are the resistance to vapor (s m⁻¹) and heat (s m⁻¹) transfer at the soil-air interface, respectively; $${c}_{vs}$$ and $${c}_{va}$$ are the vapor concentrations (m³_water m⁻³_air) at the soil surface and in the atmosphere, respectively; R_net is the net radiation (W m⁻² K⁻¹); C_p is the volumetric heat capacity of air at constant pressure (J m⁻³ K⁻¹); and T_a is the air temperature (K).

Equations 19 and 20 were used only in the top half of the soil surface layer for obtaining the unknown soil surface temperature (T_s) and surface relative humidity (R_Hs) by using a Jacobian iteration method at the beginning of each time step. T_s and R_Hs are required to calculate surface evaporation flux, surface liquid and vapor fluxes, surface sensible heat flux, surface latent heat flux and heat flux into the soil (Haverd & Cuntz, 2010).

The lower boundary condition for soil water flow can be set to free drainage (zero gradient of soil water pressure head at the lower boundary), seepage surface (constant soil water pressure head at the lower boundary), or zero water flux. The lower boundary condition for heat transport can be set to zero temperature gradient.

Boundary Conditions for Isotope Transport

Isotopic evaporation fluxes through the soil surface to the atmosphere were calculated using the Craig-Gordon (1965) model: 21 $${E}_{i}=\frac{{\alpha }_{k}}{{r}_{bw}}\left[{c}_{vs}{c}_{i,s}^{l}{\alpha }_{i}^{\ast }\left({T}_{s}\right)-{c}_{i}^{va}\right]$$ where $${\alpha }_{k}$$ is the kinetic fractionation coefficient; $${c}_{i}^{va}$$ is the isotopic concentration of the atmospheric water vapor (kg m⁻³); $${c}_{i,s}^{l}$$ is the isotopic concentration at the soil surface (kg m⁻³); and $${E}_{i}$$ is the surface soil isotopic flux (m s⁻¹) and can be specified as: 22 $${E}_{i}=-{D}_{i}^{v}\frac{\mathit{\partial }\left({c}_{v}{c}_{i}^{l}{\alpha }_{i}^{\ast }T\right)}{\mathit{\partial }z}-{D}_{i}^{l}\frac{\mathit{\partial }{c}_{i}^{l}}{\mathit{\partial }z}+{q}_{ls}{c}_{i,s}^{l}$$ where $${q}_{ls}$$ is the liquid flux at the soil surface (m s⁻¹).

The kinetic fractionation coefficient, $${\alpha }_{k}$$, can be written as (Haverd & Cuntz, 2010; Mathieu & Bariac, 1996): 23 $${\alpha }_{k}={\left({\alpha }_{diff}\right)}^{nk}$$ with 24 $$nk=\frac{\left({\theta }_{s}-{\theta }_{res}\right){n}_{a}+\left({\theta }_{sat}-{\theta }_{s}\right){n}_{s}}{{\theta }_{sat}-{\theta }_{res}}$$ where n_a and n_s are coefficients with values of 0.5 and 1.0, respectively, and θ_s is the soil surface water content (m³ m⁻³). There are several equations available to calculate $${\alpha }_{k}$$. In the current study, a sensitivity analysis was conducted on these equations (Supporting Information S1), which revealed that MOIST was not sensitive to available formulations. Thus, we chose Equation 21 for its relatively simple formulation, and it considers the influence of surface soil water content on isotope fractionation.

Combining Equations 21 and 22 leads to the final expression for surface soil isotopic concentration, $${c}_{i,s}^{l}$$ (kg m⁻³): 25 $$\frac{{\alpha }_{k}}{{r}_{bw}}\left[{c}_{vs}{c}_{i,s}^{l}{\alpha }_{k}\left({T}_{s}\right)-{c}_{i}^{va}\right]=-{D}_{i}^{v}\frac{\mathit{\partial }\left({c}_{v}{c}_{i}^{l}{\alpha }_{i}^{\ast }T\right)}{\mathit{\partial }z}-{D}_{i}^{l}\frac{\mathit{\partial }{c}_{i}^{l}}{\mathit{\partial }z}+{q}_{ls}{c}_{i,s}^{l}$$ Using this relationship, $${c}_{i,s}^{l}$$ can be readily solved as it is the only unknown quantity.

It is important to note that energy, water, and isotope mass conservation equations (Equations 19, 20 and 25) at the air-soil interface are identical to that of Haverd and Cuntz (2010), except that the influence of litter layer is not considered. This means that the fluxes from the center of the top-soil layer to the soil-air interface are balanced by the fluxes from the soil-air interface to the atmosphere (Equations 19, 20 and 25). The lower boundary conditions for isotope transport are determined by soil water flux or can be customized. Normally, the lower boundary condition for isotope transport is zero gradient or zero flux.

Numerical Implementations

MOIST adopts a cell-centered spatial discretization scheme along with the finite volume method (FVM). While the finite element method (FEM) is another commonly used approach for solving partial differential equations, it is particularly advantageous for complex problems such as two- and three-dimensional water transport in intricate geometries with complex boundary conditions. However, this study focuses on one-dimensional soil water movement, where the simplicity of the scenario reduces performance differences between FVM and FEM, making the choice of method less impactful on the performance of the fully coupled scheme. Furthermore, compared to FEM, FVM is inherently suited to conservation law problems, as it conserves fluxes both locally and globally, and it also offers simpler mathematical implementation. Therefore, FVM is used in MOIST.

To solve Equations 1, 2 and 8 simultaneously, we expand them into: 26 $$\left[\left(1-{c}_{v,sat}{R}_{H}\right)\frac{\mathit{\partial }\theta }{\mathit{\partial }h}+\left({\theta }_{sat}-\theta \right){c}_{v,sat}\frac{\mathit{\partial }{R}_{H}}{\mathit{\partial }h}\right]\frac{\mathit{\partial }h}{\mathit{\partial }t}+\left[\left({\theta }_{sat}-\theta \right){c}_{v,sat}\frac{\mathit{\partial }{R}_{H}}{\mathit{\partial }T}+\left({\theta }_{sat}-\theta \right){R}_{H}\frac{\mathit{\partial }{c}_{v,sat}}{\mathit{\partial }T}\right]\frac{\mathit{\partial }T}{\mathit{\partial }t}=-\frac{\mathit{\partial }\left({q}_{l}+{q}_{v}\right)}{\mathit{\partial }z}-{S}_{p}$$ 27 $$\left[\rho {\lambda }_{E}{c}_{v,sat}\left(\left({\theta }_{sat}-\theta \right)\frac{\mathit{\partial }{R}_{H}}{\mathit{\partial }h}-{R}_{H}\frac{\mathit{\partial }\theta }{\mathit{\partial }h}\right)\right]\frac{\mathit{\partial }h}{\mathit{\partial }t}+\left[\rho {\lambda }_{E}\left({\theta }_{sat}-\theta \right)\left({c}_{v,sat}\frac{\mathit{\partial }{R}_{H}}{\mathit{\partial }T}+{R}_{H}\frac{\mathit{\partial }{c}_{v,sat}}{\mathit{\partial }T}\right)+{C}_{soil}\right]\frac{\mathit{\partial }T}{\mathit{\partial }t}=-\frac{\mathit{\partial }{q}_{T}}{\mathit{\partial }z}$$ 28 $$\left[\left({c}_{i}^{l}-{c}_{i}^{l}{\alpha }_{i}^{\ast }{c}_{v,sat}{R}_{H}\right)\frac{\mathit{\partial }\theta }{\mathit{\partial }h}+{c}_{i}^{l}\left({\theta }_{sat}-\theta \right){c}_{v,sat}{\alpha }_{i}^{\ast }\frac{\mathit{\partial }{R}_{H}}{\mathit{\partial }h}\right]\frac{\mathit{\partial }h}{\mathit{\partial }t}+{c}_{i}^{l}\left({\theta }_{sat}-\theta \right)\left({c}_{v,sat}{R}_{H}\frac{\mathit{\partial }{\alpha }_{i}^{\ast }}{\mathit{\partial }T}+{\alpha }_{i}^{\ast }{R}_{H}\frac{\mathit{\partial }{c}_{v,sat}}{\mathit{\partial }T}+{c}_{v,sat}{\alpha }_{i}^{\ast }\frac{\mathit{\partial }{R}_{H}}{\mathit{\partial }T}\right)\frac{\mathit{\partial }T}{\mathit{\partial }t}+\left(\theta +{\alpha }_{i}^{\ast }{c}_{v,sat}{R}_{H}{\theta }_{sat}-\theta {\alpha }_{i}^{\ast }{c}_{v,sat}{R}_{H}\right)\frac{\mathit{\partial }{c}_{i}^{l}}{\mathit{\partial }t}=-\frac{\mathit{\partial }{q}_{i}}{\mathit{\partial }z}-{c}_{i}^{l}{S}_{p}$$ $$\frac{\mathit{\partial }h}{\mathit{\partial }t}$$, $$\frac{\mathit{\partial }T}{\mathit{\partial }t}$$, and $$\frac{\mathit{\partial }{c}_{i}^{l}}{\mathit{\partial }t}$$ were obtained by solving Equations 26–28 for h, T, and $${c}_{i}^{l}$$ at each time step simultaneously. Since c_v,sat is a function of T, $${R}_{H}$$ is a function of h and T, and $${\alpha }_{i}^{\ast }$$ is a function of T (Equation 10), the analytical expressions of $$\frac{\mathit{\partial }h}{\mathit{\partial }t}$$, $$\frac{\mathit{\partial }T}{\mathit{\partial }t}$$, and $$\frac{\mathit{\partial }{c}_{i}^{l}}{\mathit{\partial }t}$$ can be written as: 29 $$\frac{\mathit{\partial }h}{\mathit{\partial }t}=\frac{1}{A}\left[-\frac{\mathit{\partial }{q}_{T}}{\mathit{\partial }z}-B\frac{\left(-\frac{A}{C}\left(\frac{\mathit{\partial }q}{\mathit{\partial }z}+{S}_{p}\right)+\frac{\mathit{\partial }{q}_{T}}{\mathit{\partial }z}\right)}{\frac{D}{C}A-B}\right]$$ 30 $$\frac{\mathit{\partial }T}{\mathit{\partial }t}=\frac{-\frac{A}{C}\left(\frac{\mathit{\partial }q}{\mathit{\partial }z}+{S}_{p}\right)+\frac{\mathit{\partial }{q}_{T}}{\mathit{\partial }z}}{\frac{D}{C}A-B}$$ 31 $$\frac{\mathit{\partial }{c}_{i}^{l}}{\mathit{\partial }t}=\frac{F}{E}$$ with coefficients A to F: 32 $$\mathrm{A}=\rho {\lambda }_{E}\left[\left({\theta }_{sat}-\theta \right)\,{c}_{v,sat}\frac{\mathit{\partial }{h}_{r}}{\mathit{\partial }h}-{c}_{v,sat}{R}_{H}\frac{\mathit{\partial }\theta }{\mathit{\partial }h}\right]$$ 33 $$B={\,C}_{soil}+\rho {\lambda }_{E}\left[\left({\theta }_{sat}-\theta \right){c}_{v,sat}\frac{\mathit{\partial }{h}_{r}}{\mathit{\partial }T}+\left({\theta }_{sat}-\theta \right){R}_{H}\frac{\mathit{\partial }{c}_{v,sat}}{\mathit{\partial }T}\right]$$ 34 $$\mathrm{C}=\left(1-{c}_{v,sat}{R}_{H}\right)\frac{\mathit{\partial }\theta }{\mathit{\partial }h}+\left({\theta }_{sat}-\theta \right){c}_{v,sat}\frac{\mathit{\partial }{R}_{H}}{\mathit{\partial }h}$$ 35 $$\mathrm{D}=\left({\theta }_{sat}-\theta \right){c}_{v,sat}\frac{\mathit{\partial }{R}_{H}}{\mathit{\partial }T}+\left({\theta }_{sat}-\theta \right)\,{R}_{H}\frac{\mathit{\partial }{c}_{v,sat}}{\mathit{\partial }T}$$ 36 $$\mathrm{E}=\theta +{\alpha }_{i}^{\ast }{c}_{v,sat}{R}_{H}\left({\theta }_{sat}-\theta \right)$$ 37 $$\mathrm{F}=-\frac{\mathit{\partial }{q}_{i}}{\mathit{\partial }z}-{c}_{i}^{l}{S}_{p}-{c}_{i}^{l}\frac{\mathit{\partial }\theta }{\mathit{\partial }h}\frac{\mathit{\partial }h}{\mathit{\partial }t}-{c}_{i}^{l}{\theta }_{sat}\left[{c}_{v,sat}{R}_{H}\frac{\mathit{\partial }{\alpha }_{i}^{\ast }}{\mathit{\partial }T}\frac{\mathit{\partial }T}{\mathit{\partial }t}+{\alpha }_{i}^{\ast }{R}_{H}\frac{\mathit{\partial }{c}_{v,sat}}{\mathit{\partial }T}\frac{\mathit{\partial }T}{\mathit{\partial }t}+{\alpha }_{i}^{\ast }{c}_{v,sat}\left(\frac{\mathit{\partial }{R}_{H}}{\mathit{\partial }h}\frac{\mathit{\partial }h}{\mathit{\partial }t}+\frac{\mathit{\partial }{R}_{H}}{\mathit{\partial }T}\frac{\mathit{\partial }T}{\mathit{\partial }t}\right)\right]+{c}_{i}^{l}\left[{c}_{v,sat}{R}_{H}\theta \frac{\mathit{\partial }{\alpha }_{i}^{\ast }}{\mathit{\partial }T}\frac{\mathit{\partial }T}{\mathit{\partial }t}+{\alpha }_{i}^{\ast }{R}_{H}\theta \,\frac{\mathit{\partial }{c}_{v,sat}}{\mathit{\partial }T}\frac{\mathit{\partial }T}{\mathit{\partial }t}+{\alpha }_{i}^{\ast }{c}_{v,sat}\theta \left(\frac{\mathit{\partial }{R}_{H}}{\mathit{\partial }h}\frac{\mathit{\partial }h}{\mathit{\partial }t}+\frac{\mathit{\partial }{R}_{H}}{\mathit{\partial }T}\frac{\mathit{\partial }T}{\mathit{\partial }t}\right)+{\alpha }_{i}^{\ast }{c}_{v,sat}{R}_{H}\frac{\mathit{\partial }\theta }{\mathit{\partial }h}\frac{\mathit{\partial }h}{\mathit{\partial }t}\right]$$

Equations 26–28 are expressed as a system of coupled ordinary differential equations of Equations 29–37. This system of equations ($$\frac{\mathit{\partial }h}{\mathit{\partial }t}$$, $$\frac{\mathit{\partial }T}{\mathit{\partial }t}$$, and $$\frac{\mathit{\partial }{c}_{i}^{l}}{\mathit{\partial }t}$$) are solved at the center of each layer (cell-centered discretization) by MATLAB solver ode113 (or ode23tb). Furthermore, the semi-coupled scheme (Figure 1) is also incorporated in MOIST, where the derivative vector of the soil water and heat transport equations (Equations 29 and 30) are constructed by Equations 32–35 and passed to the solver first. Then, Equation 31 is rewritten as: 38 $$\frac{\mathit{\partial }{C}_{i}^{l}}{\mathit{\partial }t}=\frac{-\frac{\mathit{\partial }{q}_{i}}{\mathit{\partial }z}-H}{I}$$ with 39 $$H={\,c}_{i}^{l}\left[\frac{\mathit{\partial }\theta }{\mathit{\partial }t}-{c}_{v,sat}{R}_{H}\,{\alpha }_{i}^{\ast }\frac{\mathit{\partial }\theta }{\mathit{\partial }t}+\frac{\mathit{\partial }\left({c}_{v,sat}{R}_{H}\,{\alpha }_{i}^{\ast }\,{\theta }_{sat}\right)}{\mathit{\partial }t}-\theta \frac{\mathit{\partial }\left({c}_{v,sat}{R}_{H}\,{\alpha }_{i}^{\ast }\,\right)}{\mathit{\partial }t}\right]$$ 40 $$I=\theta +{c}_{v,sat}{R}_{H}{\alpha }_{i}^{\ast }\left({\theta }_{sat}-\theta \right)$$

The solutions from the coupled soil water and heat transport equation (h and T) at new time point are used to update parameter values ($${c}_{v,sat}$$, $${R}_{H}$$, $${\alpha }_{i}^{\ast }$$, $$\theta $$) in Equations 39 and 40. Subsequently, the derivative vector of isotope transport can be constructed and solved (Equations 38–40).

MATLAB Solver ode113

The ode113 solver uses an adaptive, variable-order, variable-step-size (VOVS) method (Shampine, 2002). This is implemented with a variable order Adams-Bashforth-Moulton (ABM) method, which is a combination of an explicit type of the Adams-Bashforth (AB) and an implicit type of Adams-Moulton (AM) methods. Specifically, the AB method is used to obtain the solution at the new time step by taking multiple previous time steps into account, while the AM method is used to make corrections.

The ode113 selects automatically between the first and twelfth order approximation during the computation based on the estimation errors. This is helpful for minimizing the estimated errors and for achieving high efficiency in time. Moreover, the time step size is adjusted according to the estimation error. Therefore, ode113 can handle a wide range of problems with high accuracy and efficiency. The ode113 is well suited to transport in relatively uniform media but is susceptible to numerical oscillation when hydraulic conductivities between layers differ greatly.

MATLAB Solver ode23tb

Ode23tb is a solver specifically designed for solving ordinary differential equations with highly oscillatory solutions, such as those arising from heterogeneity in hydraulic conductivities between soil layers. The solver combines a trapezoidal rule (sometimes referred as the second order AM method) with a second order backward differentiation formula (BDF), hence the “tb” suffix. As an implicit solver, ode23tb is generally more computationally expensive than other solvers that adopt explicit numerical schemes. However, because it adopts the trapezoidal BDF method, it is more efficient than other types of implicit methods, such as the fully implicit Euler method or the backward Euler method. Therefore, ode23tb is better suited than ode113 when the soil physical properties differed greatly between layers. Like ode113, ode23tb can adjust the step size automatically based on the oscillatory behavior of the solution. Thus, ode23tb is an efficient and accurate solver for stiff systems, making it less likely to have numerical instability.

Modeling Efficiency

To evaluate the model performance quantitatively, the Nash-Sutcliffe efficiency (NSE) (Nash & Sutcliffe, 1970) is used: 41 $$NSE=1-\frac{\sum\limits _{t=1}^{{T}_{ol}\,}{\left[{M}_{0}(t)-{M}_{m}\right]}^{2}}{\sum\limits _{t=1}^{{T}_{ol}}{\left[{M}_{0}(t)-\overline{{M}_{0}}\right]}^{2}\,}$$ where M₀ and M_m are observations (measurements) and modeling values (simulations) respectively. $$\overline{{M}_{0}}$$ is the average of the observations over time, t is the temporal point, and T_ol is the total number of temporal steps. NSE ranges between negative infinity to 1. NSE of one is indicative of excellent performance of the model in predicting the temporal variations of variables, while an NSE of 0 suggests the model can only reflect average values. A negative NSE implies poor performance of the model in regenerating the temporal variations of variables.

Moreover, Kling-Gupta efficiency (KGE, Gupta et al., 2009) has become popular in recent years to evaluate the model performance (Knoben et al., 2019; Zhou et al., 2022). KGE considers correlation, bias, and variability between simulations and measurements, providing a more comprehensive assessment of the model performance and can be written as: 42 $$KGE=1-\sqrt{{(\text{cor}-1)}^{2}+{\left(\frac{{\sigma }_{m}}{{\sigma }_{o}}-1\right)}^{2}+{\left(\frac{{\mu }_{m}}{{\mu }_{o}}-1\right)}^{2}}$$ where cor is the linear correlation between observations and simulations, $${\sigma }_{m}$$ and $${\sigma }_{o}$$ are standard deviations of simulations and observations, respectively; $${\mu }_{m}$$ and $${\mu }_{o}$$ are mean of simulations and observations, respectively. Like NSE, KGE = 1 indicates perfect agreement between simulations and observations. A model had a negative KGE is generally considered as “not satisfactory” (Andersson et al., 2017).

Sometimes, NSE and KGE may disagree with each other during model comparisons (Knoben et al., 2019). A third index, the absolute error (MAE), is used to evaluate the model performance. Here, we decided to use MAE rather than the root mean square error (RMSE) because the residuals between simulations and measurements are non-normally distributed (Chai & Draxler, 2014; Sprenger et al., 2018). MAE is given as: 43 $$MAE=\frac{1}{{T}_{ol}\,}\sum\limits _{t=1}^{{T}_{ol}}\,\vert {M}_{0}(t)-{M}_{m}(t)\vert $$

Model Validation

To validate MOIST, we conducted sensitivity analysis on various temporal and spatial discretization, through theoretical, semi-analytical, and long-term field tests. However, we mainly reported the findings of sensitivity under various temporal and spatial discretization, semi-analytical tests under unsaturated and non-isothermal conditions, and a long-term field test. This is because the results from other semi-coupled models are readily available for comparison. Detailed information about the theoretical tests, and the semi-analytical test under saturated isothermal conditions, can be found in Supporting Information S1.

Comparison Between the Fully Coupled and the Semi-Coupled Method Under Various Temporal and Spatial Discretization

The feature of MOIST is its use of the fully coupled method to simultaneously solve all governing equations. To demonstrate the difference between fully coupled and semi-coupled methods, both versions of MOIST employing fully and semi-coupled methods were utilized to calculate the δ²H transient distribution under unsaturated non-isothermal conditions with various temporal and spatial discretization. We expect that this comparison can effectively reflect the differences between the semi-coupled and fully coupled methods.

The current test assumes a 1 m soil column containing Yolo light clay. The relationships between soil water content, pressure head, and unsaturated hydraulic conductivity for the Yolo light clay were described by the Brooks-Corey (Brooks & Corey, 1964) model: 44 $$S=\frac{\theta -{\theta }_{res}}{{\theta }_{sat}\mathit{-}{\theta }_{res}}=\left\{\begin{array}{@{}c@{}}{\left(\frac{h}{{h}_{e}}\right)}^{-\lambda },h\mathit{\le }he\\ 1\quad h\mathit{\ge }he\end{array}\right.$$ 45 $${S}^{\eta }=\left\{\begin{array}{@{}c@{}}\frac{K}{{K}_{sat}},h\mathit{\le }he\\ 1\ h\mathit{\ge }he\end{array}\right.$$ where S is the effective saturation; h_e is the air-entry value (m); and λ and η are the shape coefficients, where η = 2/λ + 3 (Table 1).

Table 1 Hydraulic Properties of the Simulated Soil Used in Theoretical Tests

he (m)	λ	Ksat (m s−1)	θsat (m3 m−3)	θres (m3 m−3)
−0.193	0.22	1.23 × 10−7	0.35	0.01

The initial conditions assume the soil water content as 70% of its saturated value and the net radiation as a constant value of 200 W m⁻² at the soil surface. The air temperature and relative humidity during the simulation period remained at 30°C and 0.2, respectively. Water in the soil column can escape only through evaporation from the top of the column. The upper boundary conditions for soil water and isotope transport were calculated by Equations 19, 20 and 25. For the lower boundary, the rate of water supply from the bottom of the profile is equivalent to the evaporation rate at each time step. The simulation length was 700 hr.

Three spatial steps, Δz = 0.01 m, Δz = 0.02 m, and Δz = 0.005 m were used. The fixed initial temporal step was 25 s. By contrast, for the various temporal steps, three initial steps, Δt = 100 s, Δt = 50 s, and Δt = 25 s, were used and the spatial step was fixed at 0.01 m. Note that if a large temporal step is used, the result may be affected by the time-adaption, which could reduce the difference between semi-coupled and fully coupled methods. Thus, we used small temporal steps to approximately transform the time-adaptive solver into a fixed time step solver. This is because the solver can obtain a satisfactory solution within the provided initial time step without the need for time adjustments when boundary conditions are stable.

Semi-Analytical Test Under Unsaturated Non-Isothermal Conditions

Barnes and Allison (1984) developed a semi-analytical solution to predict δ²H and δ¹⁸O profiles under unsaturated, non-isothermal conditions as: 46 $$\frac{d{\delta }_{i}}{dz}+\frac{{\delta }_{i}-{\delta }_{i,\mathit{sup}}}{{z}_{l}+{R}_{H}{z}_{v}}=\frac{{R}_{H}{z}_{v}\left({\alpha }_{k}-{\alpha }_{i}^{\ast }\right)}{{z}_{l}+{R}_{H}{z}_{v}}\frac{d\left[\mathit{ln}\left(\rho {R}_{H}{c}_{v,sat}\left({\alpha }_{k}-{\alpha }_{i}^{\ast }\right)\right)\right]}{dz}$$ 47 $${z}_{l}=\frac{{D}_{i,s}^{l}}{{q}_{evap}}$$ 48 $${z}_{v}=\frac{{D}_{i,s}^{v}{c}_{v,sat}}{{q}_{evap}}$$ where $${z}_{l}$$ and $${z}_{v}$$ are the liquid and vapor characteristic lengths (m), respectively. Note that all modeling conditions were the same as that described in Section 2.2.1, except the 1 m soil column was divided uniformly at a step of 0.01 m and the total simulation length was 250 days to ensure that steady state is achieved.

Long-Term Simulations at HBLFA Raumberg-Gumpenstein, Austria

Site Description

Precipitation and seepage water from lysimeters were collected from May 2002 to February 2007 by Stumpp et al. (2012) at the HBLFA Raumberg-Gumpenstein, Austria. During the experiment, the air temperature exhibited a sinuous variation, ranging between −15°C and 27°C (Figure 3a), with a mean of 8.2°C. Similarly, the atmospheric relative humidity showed seasonal fluctuations and varied between 0.30 and 0.99 (Figure 3b), with a mean of 0.89. In addition, most precipitation events occurred during the summer (Figure 3c), with a daily mean rainfall of 2.8 mm day⁻¹.

[IMAGE OMITTED. SEE PDF]

Five lysimeters were used to investigate the influence of land cover and fertilization on soil water and solute transport by Stumpp et al. (2012). For simplicity, the current study used only lysimeter-3 for the comparison of numerical simulations. Lysimeter-3 had a surface area of 1 m², depth of 1.5 m, and was filled with three soil horizons (0–0.25 m, 0.25–1.0 m, 1.0–1.5 m) of undisturbed Dystric Cambisol (Stumpp et al., 2012) (Figure 4).

[IMAGE OMITTED. SEE PDF]

Each year the lysimeter was planted with winter rye, which had a maximum rooting depth of 1 m. Weekly precipitation and drainage water samples from the bottom of Lysimeter-3 were collected between May 2002 and February 2007. Isotopic compositions were analyzed by using dual-inlet mass spectrometry. Further information about the site and experimental procedures can be obtained from Stumpp et al. (2012).

Model Setup

The MOIST simulations adopted the solute transport and soil hydraulic parameter values provided by Stumpp et al. (2012). The soil hydraulic properties were described by the van Genuchten (1980) model. Because the measured saturated hydraulic conductivities varied greatly within different soil horizons (Stumpp et al., 2012), water flux at the interface of different soil layers may vary drastically. Therefore, to minimize the oscillation of numerical solutions, the ode23tb solver, which is designed for stiff problems, was used in this simulation.

The initial soil water content and δ¹⁸O profiles were provided by Stumpp et al. (2012). The upper boundaries of soil water and heat transport were calculated by Equations 14 and 16, while the upper boundary of isotope transport was calculated by Equation 25. The lower boundary condition for soil water flow was defined as seepage surface. Zero temperature and zero isotopic concentration gradients were set as the lower boundary conditions of heat and isotope transport, respectively.

The root distribution of winter rye at the site varied during the growing season and has been described in HYDRUS-1D (Stumpp et al., 2012), along with the water stress function described by the Feddes model (Feddes et al., 1978). Following this approach, the water stress function parameters h₁, h₂, h₃, and h₄ were set as 0, 0.01, 5, and 160 m, respectively (Stumpp et al., 2012).

Environmental parameters such as air temperature, air relative humidity, and precipitation were assumed to be constant within each hour. Regarding the isotopic compositions of atmospheric vapor, $${\delta }_{a}$$, it is expressed as following equation when rainfall occurs (Skrzypek et al., 2015): 49 $${\delta }_{a}=\frac{{\delta }_{rain\,}-\left({\alpha }^{+}+1\right)1000}{{\alpha }^{+}}$$ where $${\alpha }^{+}$$ is the equilibrium fractionation factor.

Obviously, Equation 49 cannot be used when there is a rain-free period. However, the field measurements of $${\delta }_{a}$$ (oxygen-18) in Austria (where the lysimeter study conducted) ranged between −27‰ and −13‰ annually (Kurita et al., 2012), with an average of −20‰ (Galewsky et al., 2016). Therefore, we used −20‰ as the value of $${\delta }_{a}$$ for the rain-free period. Values −27‰ and −13‰ were used to explore the influence of $${\delta }_{a}$$ values on the simulated outflow isotopic signals.

Soil hydraulic parameters were calibrated by both Stumpp et al. (2012) and Zhou et al. (2022) with this data set (Table 2). We obtained a new set of hydraulic parameters by calibrating MOIST using the same data set. Consequently, simulations were conducted in three scenarios: Scenario 1, calibrated parameters from Stumpp et al. (2012) were utilized in MOIST; Scenario 2, calibrated parameters from Zhou et al. (2022) were employed in MOIST; Scenario 3, parameters calibrated by MOIST were used. In the first two scenarios, we aimed to highlight the advantages of the coupling scheme because the parameters are independent of MOIST, allowing for a comparison that reflects model differences. In Scenario 3, we aimed to demonstrate that model calibration can further enhance the performance of MOIST.

Table 2 Calibrated Soil Hydraulic Parameters by Stumpp et al. (2012) and Zhou et al. (2022)

	Layers (m)	α (m−1)	n	θsat (m3 m−3)	θres (m3 m−3)	Ksat (m s−1)
Stumpp et al. (2012)	0–0.30	2.3	1.14	0.30	0	1.27 × 10−5
0.31–0.90	7.6	1.07	0.32	0	6.94 × 10−4
0.91–1.50	1.6	1.07	0.32	0	1.27 × 10−5
Zhou et al. (2022)	0–0.30	2.0	1.15	0.30	0	2.55 × 10−5
0.31–0.90	30.0	1.11	0.41	0	3.32 × 10−5
0.91–1.50	8.2	1.10	0.30	0	2.55 × 10−5

The soil column was divided into 150 layers with a spatial step of 0.01 m in the model simulation. The simulation length was 1,736 days, and the initial temporal step was 86,400 s. The minimum time step was 1 × 10⁻⁸ s and the maximum time step was 86,400 s. The time step can be automatically adjusted within the range between minimum and maximum. Lastly, the lysimeter simulation (including Sections 2.2.1 and 2.2.2) was completed by MATLAB R2022a on a personal computer with an Intel processor (i7-10700k).

Model Calibration

Three soil horizons were defined in this lysimeter study (Stumpp et al., 2012). Limited by the computation time of MOIST, we did not calibrate all 13 parameters (4 parameters of van Genuchten model of each horizon plus the longitude dispersivity). Instead, we calibrated the saturated hydraulic conductivity of each horizon because K_sat are generally prone to large uncertainties (Zhang & Schaap, 2019) and highly sensitive to the choice of models (Table 2). Moreover, the longitude dispersivity was also calibrated because it directly affects the isotope transport (Zhou et al., 2022). One dispersivity for three horizons can be estimated because only the isotopic composition of the drainage water was measured in the Stumpp et al. (2012) data set. The remaining soil hydraulic parameters, such as α, n, θ_sat, and θ_res, were identical to that from Stumpp et al. (2012).

The calibrated saturated hydraulic conductivities from Stumpp et al. (2012) served as base values, and each was assigned a coefficient ranging between 0.1 and 10 for scaling purposes. These coefficients were uniformly sampled through a Monte Carlo simulation. Then, the Nash-Sutcliffe Efficiency (NSE) of simulated outflow isotope signals, along with the disparity between modeled and measured total outflow amount, was calculated. In cases where MOIST failed to produce results under specific coefficients, a negative NSE was assigned. The set of coefficients yielding the highest NSE and the smallest difference between modeled and measured total outflow was selected. The calibration was conducted on Cedar, a heterogeneous high-performance cluster (HPC) within the Digital Research Alliance of Canada. Note that the validation is simplified, but it provided valuable insights into the potential impact of calibration on the performance of MOIST when applied to real-world situations.

Results

MOIST was validated against semi-coupled solutions under various temporal and spatial discretization, theoretical solutions, semi-analytical solutions under saturated isothermal conditions, semi-analytical solutions under unsaturated non-isothermal conditions, and long-term field measurements. Here we focused primarily on the validation of MOIST against semi-coupled solutions under various temporal and spatial discretization, semi-analytical tests under unsaturated non-isothermal conditions and the long-term field measurements. Validations against theoretical tests and semi-analytical tests under saturated isothermal conditions can be found in the Supporting Information S1.

Comparison Between Fully Coupled and Semi-Coupled Methods Under Various Spatial and Temporal Discretization

The feasibility and stability of fully coupled and semi-coupled methods were assessed. The initial δ²H profile was uniformly distributed with a value of 0‰. Given that water (δ²H = 0‰) was continuously supplied from the bottom of the profile, and water escaped the column through evaporation only at the soil surface, δ²H of soil water in the column at any point in space and time is unlikely to be ≤0‰, especially at its lower boundary. However, considering numerical errors, a negative value of −0.3‰ was set as the threshold for the largest acceptable deviation from 0‰ (δ²H of the bottom-supplied water) because the largest error of simulated δ²H from MOIST in the semi-analytical test is 0.31‰ (Table 4). Additionally, a 0.1‰ threshold was used as the maximum acceptable absolute deviation between the fully coupled and semi-coupled methods, recognizing that the two methods cannot produce identical numerical solutions.

Various Temporal Discretization

Under different temporal discretization, the fully coupled method consistently displayed δ²H values greater than −0.3‰ (see Figures 5a–5c). By contrast, the semi-coupled method exhibited unreasonable negative δ²H values (depicted by dark blue areas in Figures 5d–5f) across all temporal step lengths. This discrepancy can be attributed to errors originating from soil water and heat equations being incorporated into the isotope transport equation. Consequently, noticeable disparities in the estimated transient δ²H between the two methods emerged (highlighted by pink areas in Figures 5g–5i), and these differences became more pronounced with increasing temporal steps (Figures 5g–5i). This highlights the robustness of the fully coupled method, which is less sensitive to temporal discretization and capable of accommodating larger time steps compared to the semi-coupled method.

[IMAGE OMITTED. SEE PDF]

Solver CPU times of both fully coupled and semi-coupled methods were presented in Figure 6. For the fully coupled method, solver CPU times were 14.4 s, 7.2 s, and 3.5 s under temporal discretization of 100 s, 50 s, and 25 s, respectively (Figure 6). By contrast, the corresponding times for the semi-coupled method were 29.2 s, 15.2 s, and 7.7 s, respectively (Figure 6). Notably, the fully coupled method reduced solver CPU time by approximately 50% compared to the semi-coupled method. This gained efficiency can be attributed to the fact that all equations were solved simultaneously in the fully coupled method, requiring the solver to be called only once within a single iteration. On the other hand, the semi-coupled method involved solving water and heat transport equations before addressing the isotope equation, leading to the solver being called twice within one iteration (Figure 6). Additionally, the solver CPU time approximately decreased linearly with the initial step increased (Figure 6), which supports that employing a small initial time step can effectively alleviate the need for a time-adaptive solver.

[IMAGE OMITTED. SEE PDF]

Various Spatial Discretization

The comparison between fully coupled and semi-coupled method was also conducted under various spatial discretization (Figures 7a–7f). The semi-coupled method exhibited similar δ²H transit profiles as the fully coupled method. The difference between these two methods was diminished as spatial step decreased (Figures 7g–7i). However, the unreasonable negative δ²H values can be observed again at the bottom from semi-coupled method (Figures 7d–7f), as appeared under various temporal discretization (Figure 5). This illustrated that the semi-coupled method is more sensitive to spatial discretization than the fully coupled method.

[IMAGE OMITTED. SEE PDF]

Again, the solver CPU times for both fully coupled and semi-coupled methods under various spatial discretization were presented in Figure 8. Specifically, with spatial discretization of Δz = 0.005 m, 0.01 m, and 0.05 m, the solver CPU times for the fully coupled method were 21.0 s, 14.4 s, and 14.5 s, respectively (Figure 8). By contrast, the corresponding times for the semi-coupled method were 29.9 s, 29.2 s, and 28.9 s under the same conditions. The consistently lower CPU times from the fully coupled method than the semi-coupled method confirms the enhanced efficiency of the fully coupled approach.

[IMAGE OMITTED. SEE PDF]

Semi-Analytical Tests Under Unsaturated Non-Isothermal Conditions

Under unsaturated conditions, a drying layer appeared at the soil surface (Figure 9a) and water flow in this region was dominated by upward vapor diffusion (Figure 9b). The total water flux within the column was uniform with depth (Figure 9b) because a steady state flow was achieved. δ²H and δ¹⁸O increased sharply with depth until reaching their peak values at approximately 0.02 m below the soil surface, which is the location of the bottom of the drying layer (Figures 9c and 9d). This differs from the saturated conditions (Supporting Information S1), where the maximum δ²H and δ¹⁸O appeared at the soil surface. This is because in the unsaturated system the soil water content at the air-soil interface was close to residual soil water content: when a drying layer forms, the air invades into the drying layer, resulting in a downward shift of the isotope peak values as compared to the saturated system.

[IMAGE OMITTED. SEE PDF]

Visually, the analytical solution and the estimates from MOIST agreed very well across the entire soil profile (Figures 9c and 9d). If numeric models are biased, they tend to be biased at locations where sharp changes of isotope occur. Therefore, we compared the simulated peak concentrations obtained from the MOIST and the analytical solutions. Quantitatively, the approaches agreed well at the maximum δ²H and δ¹⁸O (Table 3). They also agreed in terms of the slopes of HDO/H₂¹⁸O in the liquid and vapor dominated regions (Table 3). This demonstrates that MOIST is robust with respect to the accuracy of calculating isotope transport within soil.

Table 3 Comparison of Maximum δ²H, δ¹⁸O, and the Slope of ²H/¹⁸O in Liquid Dominated Regions Estimated From MOIST and the Analytical Solution Under Unsaturated and Non-Isothermal Conditions

	MOIST	Semi-analytical solution
Maximum δ2H (‰)	40.34	40.65
Maximum δ18O (‰)	21.15	22.13
Slope of HDO/H218O in vapor dominated region	3.55	3.31
Slope of HDO/H218O in liquid dominated region	1.91	1.90

Zhou et al. (2021) presented the simulated maximum δ²H and δ¹⁸O from a revised HYDRUS-1D model (rHZ for short) under unsaturated non-isothermal conditions with fine (Δz = 1 × 10⁻⁶ m), medium (Δz = 0.005 m), and coarse (Δz = 0.01 m) spatial discretization. However, their spatial intervals were not uniformly distributed, except in the coarse spatial discretization scenario. Because MOIST currently does not support adaptive spatial steps (which will be addressed in the future), its results were compared with those of rHZ under the coarse scenario.

Table 4 showed that the largest differences in δ²H and δ¹⁸O between MOIST and the analytical solution were 0.31 and 0.98‰, respectively. By contrast, these values were 34.68 and 24.88‰ from rHZ. Compared to rHZ, MOIST significantly reduced the numerical error of the maximum isotopic compositions between numerical and analytical solutions (Table 4), illustrating that MOIST offered better accuracy under the same spatial discretization scheme as compared to rHZ. Many factors could lead to the superior performance of MOIST. The most striking difference is that MOIST is based on the fully coupled method while rHZ is based on the semi-coupled method. In the semi-coupled method, mass balance errors from soil water and heat transport can be inflated by coarse spatial resolutions and transferred into isotope transport at each time step. This is reflected in the fact that fine spatial discretization is needed for rHZ to pass this analytical test. By contrast, there was little error transfer from the soil water and heat equations to the isotope equation in MOIST, because all the equations were integrated by the fully coupled method. These findings illustrated that the fully coupled method has the potential to improve numerical accuracy and stability under large spatial steps.

Table 4 Difference of Maximum Isotopic Composition Between Numerical and Analytical Solutions With a Layer Thickness of 0.01 m

	δ2H (‰)	δ18O (‰)
Revised HYDRUS-1D (Zhou et al., 2021)	34.68	24.88
MOIST	0.31	0.98

Validation by the Long-Term Experiment at HBLFA Raumber-Gumpenstein

In the long-term validation, MOIST modeled the isotopic signals of drainage water under three scenarios: Soil hydraulic parameters calibrated by Stumpp et al. (2012) (Scenario 1, Figure 10a), Zhou et al. (2022) (Scenario 2, Figure 10b), and MOIST (Scenario 3, Figure 10c).

[IMAGE OMITTED. SEE PDF]

In scenario 1, the MOIST model, as well as the HYDRUS-1D revised by Stumpp et al. (2012) (rHS for short), reproduced temporal variations of measured δ¹⁸O in drainage water (Figure 10). However, the NSE, KGE, and MAE of rHS were 0.35, 0.60, and 1.00‰, respectively. By contrast, these indices from MOIST were 0.47, 0.58, and 0.92‰, respectively, suggesting that MOIST slightly outperformed rHS when using parameters calibrated by rHS. Interestingly, although fractionation would typically result in more enriched isotopic compositions, this was not observed between the curves generated by rHS (without fractionation, represented by the blue solid line in Figure 10) and MOIST (with fractionation, represented by the red solid line in Figure 10). This discrepancy may be attributed to the possibility that the calibrated parameters do not precisely match MOIST. However, when the fractionation process was removed from MOIST, as expected, the temporal isotopic compositions became more depleted (red dash line in Figure 10). This underscores the significance of incorporating fractionation when modeling isotope transport within soil. In scenario 2, parameters calibrated by another version of revised HYDRUS-1D (Zhou et al., 2022, rHZ for short) were employed in MOIST for comparison. It can be observed that rHZ and MOIST generated comparable temporal patterns of drainage water isotopic compositions, especially after 1,500 days (Figure 10b). Quantitively, NSE, KGE, and MAE of rHZ were 0.19, 0.58, and 1.15‰, respectively. By contrast, these indices from MOIST were 0.33, 0.52, and 1.04‰, respectively. The smaller errors (MAE) from MOIST supported the fully coupled method outperformed the semi-coupled method. Moreover, simulated isotopic signals of outflow from MOIST under various $${\delta }_{a}$$ values illustrated a marginal influence of $${\delta }_{a}$$ on the time series of outflow isotopic signals (Figure 10b).

Results from scenario 3 demonstrated that accuracy of predictions from MOIST benefited significantly from calibration (Figure 10c). The temporal trend of δ¹⁸O of drainage water fitted well to the measurements, especially from the first 600 days. Moreover, the highest NSE (0.48), KGE (0.76), and the smallest MAE (0.90) among the three scenarios not only confirmed the outperformance of MOIST, but also underscored the importance of calibration when applying a hydrological model into real conditions.

Calibrated Soil Hydraulic Conductivities and the Longitude Dispersivity

Only saturated hydraulic conductivity (K_sat) and longitudinal dispersivity were calibrated in MOIST (Table 5) and the remaining parameters were maintained as calibrated by Stumpp et al. (2012). This is because MOIST exhibited better accuracy with Stumpp et al. (2012) calibrated parameters compared to those calibrated by Zhou et al. (2022) (Figure 10).

Table 5 Initial Longitude Dispersivity and Saturated Hydraulic Conductivities From Stumpp et al., (2012) and Calibrated by MOIST

	Layers (m)	Longitude dispersivity (m)	Ksat (m s−1)	MAE of K (m s−1)	Overall MAE of K (m s−1)	Objective functions
Stumpp et al. (2012)	0–0.30	4.70 × 10−2	1.27 × 10−5	1.83 × 10−6	1.80 × 10−5	SW; δBF; K; RC; BF-t
0.31–0.90	6.94 × 10−4	4.70 × 10−5
0.91–1.50	1.27 × 10−5	3.47 × 10−6
MOIST	0–0.30	6.08 × 10−2	4.20 × 10−6	5.37 × 10−7	1.48 × 10−5	δBF; BF-T
0.31–0.90	5.80 × 10−4	3.92 × 10−5
0.91–1.50	5.72 × 10−6	3.89 × 10−7

MOIST yielded 6.08 × 10⁻² m for longitudinal dispersivity, which is close to initial values (4.70 × 10⁻² m, Table 5). For saturated hydraulic conductivities (K_sat), initial values (from rHS) were 1.27 × 10⁻⁵, 6.94 × 10⁻⁴, and 1.27 × 10⁻⁵ m s⁻¹ for three horizons, with corresponding MAE values of estimated K were 1.83 × 10⁻⁶, 4.7 × 10⁻⁵, and 3.47 × 10⁻⁶ m s⁻¹, and the overall MAE is 1.8 × 10⁻⁵ m s⁻¹. By contrast, MOIST calibrated K_sat values for three horizons to 4.20 × 10⁻⁶, 5.80 × 10⁻⁴, and 5.72 × 10⁻⁶ m s⁻¹, with respective MAE values of estimated K were 5.37 × 10⁻⁷, 3.92 × 10⁻⁵, and 3.89 × 10⁻⁷ m s⁻¹, and the overall MAE is 1.48 × 10⁻⁵ m s⁻¹. In comparison to initial values from rHS, MOIST exhibited smaller errors of K that calculated from calibrated K_sat. This suggested that the results from MOIST were acceptable, although the calibration of MOIST was simplified (fewer objective functions than rHS).

Discussion

Difference Between the Fully Coupled Method and the Semi-Coupled Method

The semi-coupled variation of MOIST exhibited unfeasible δ²H values across various spatial and temporal discretization, a problem not encountered by the fully coupled version of MOIST (Figures 5 and 7). This discrepancy can be attributed to the fact that mass balance errors arising from solving soil water and heat transport equations are introduced into the isotope transport equation within the semi-coupled method. By contrast, the fully coupled method avoids this issue by analytically integrating the soil water and heat transport equations into the isotope transport equation (Equations 36 and 37). This aligns with the findings of Pimenta and Alves (2019), who also observed that the fully coupled method excels in capturing coupled physical processes, mitigating errors, and delivering more precise outcomes when compared to the semi-coupled approach.

Moreover, it is worth noting that the disparities between the fully coupled and semi-coupled methods diminished as finer temporal and spatial steps were employed (Figure 5). This phenomenon can be attributed to the fact that the finer spatial and temporal discretization leads to reduced mass balance errors in solving the soil water and heat transport equations (Braud et al., 2005; Huang et al., 2018; Zhou et al., 2021). However, a surprising trend can be observed that as the layer thickness decreased in the semi-coupled simulations, a longer dark blue band appeared at the bottom of the column (Figures 7d–7f). This phenomenon can be attributed to numerical instability in the differential equations. To ensure numerical stability, adherence to CFL (Courant-Friedrichs-Lewy) condition is essential (Courant et al., 1967): 50 $$\vert \frac{\varsigma {\Delta }t}{{\Delta }z}\vert \,< \ {\varsigma }_{\mathit{max}}$$ where $${\varsigma }_{\mathit{max}}$$ is typically set to 1. This equation stipulates that the spatial step (Δz) should not be excessively small compared to the temporal step (Δt), as excessively small Δz values can lead to numerical oscillations. In the case of the semi-coupled method, where Δt was 25 s, decreasing Δz resulted in a larger $$\vert \frac{\varsigma {\Delta }t}{{\Delta }z}\vert $$ value, increasing the likelihood of instabilities. However, in the fully coupled method, as shown in Figures 5 and 7, the relationship between instability and spatial step size was less evident. This could be due to the scaling down of $$\varsigma $$ by coefficients from the soil water and temperature equations (Equations 26 and 27), which may reduce the value of $$\vert \frac{\varsigma {\Delta }t}{{\Delta }z}\vert $$. Consequently, compared to the semi-coupled method, the fully coupled method demonstrated greater stability across a range of spatial step sizes, highlighting its superior reliability in predicting temporal variations in the isotopic compositions of soil water.

While our tests were conducted under ideal boundary conditions, the results aligned with the findings presented in Panday and Huyakorn (2004). They further underscore the fact that the fully coupled method effectively minimizes errors and displays reduced sensitivity to variations in spatial and temporal step sizes. This attribute enhances the accuracy of one-dimensional isotope transport simulations. However, the semi-coupled method has advantages in steady-state systems. By assuming linear changes in the water and heat solution within the time step, the isotope equation in the semi-coupled method can use a smaller time step, improving the stability of the isotope solution, although at the cost of increased computation time. In our simulations (Sections 2.2.1), to ensure consistency, the time step for the isotope equation in the semi-coupled method was kept the same as that for the water and heat equations, but the time consumption of semi-coupled method was doubled compared to the fully coupled method, as shown in Figure 6. This contrasts with the theoretical expectation that the semi-coupled method should generally be faster due to its lower memory and computational requirements per iteration (Gee & Gracie, 2021; Pimenta & Alves, 2019). Several factors can contribute to this discrepancy.

First, the choice between using the semi-coupled method and the fully coupled method depends on the degree of coupling exhibited by the variables involved (Pimenta & Alves, 2019). In our study, all variables are tightly coupled due to intricate interactions between water and heat transport, which concurrently influence the movement of isotopes. Solving a coupled system in a segregated manner may lead to increased computational time, as it may require more iterations for convergence to occur (Ammara & Masson, 2004). Notably, the semi-coupled method consumed more time because it necessitates two solver calls within a single time step, whereas the fully coupled method requires only one call. This observation is supported by the solver CPU time (Figures 6 and 8).

Second, it is essential to consider the scale of the problem. In our tests, we examined a relatively small-scale problem, with the minimum spatial discretization (Δz) is 0.005 m. This means that the semi-coupled method processed two matrices with sizes of 400 × 400 and 200 × 200, respectively (Figure 11a). By contrast, the fully coupled method handled a single matrix with a size of 600 × 600 (Figure 11a). Both matrix size and quantity are critical factors influencing solver CPU times. With a smaller matrix size (resulted from a larger Δz), the fully coupled method can be faster than the semi-coupled method because both methods require similar CPU time per solver call, but the semi-coupled method calls the solver twice. However, for a larger matrix size (resulted from a smaller Δz), the efficiency of the fully coupled method may decrease, as processing a single large matrix can significantly increase RAM usage and processing time compared to handling two sub-matrices in the semi-coupled method (Figure 11a). For example, the solver CPU times of the semi-coupled method were 2 times larger than that of the fully coupled method when Δz = 0.02 and 0.01 m, but 1.4 times larger when Δz = 0.005 m (Figure 8). This suggested that the efficiency of fully coupled method depends on the scale of the problem (e.g., the number of soil layers) (Khorrami & Banihashemi, 2019).

[IMAGE OMITTED. SEE PDF]

To better illustrate the influence of problem scale (Δz) on the solver CPU time, we attempted finer simulations, such as Δz = 0.001 m, but non-convergence of the upper boundary solutions occurred in MOIST. To resolve this, Δt must be reduced but reducing the initial Δt alters the simulation conditions, making the solver CPU time non-comparable. Although the current version of MOIST cannot explicitly explore larger simulation problems without modifying initial Δt, we developed a simplified example (codes are available from the Supporting Information S1) with six different Δz levels (0.0005 m, 0.001 m, 0.005 m, 0.01 m, 0.05 m, and 0.1 m) to compare time consumption of the fully coupled method and semi-coupled method across different problem scales (Δz). It can be found that the solver CPU times for the fully coupled and semi-coupled methods followed the same trend as seen in Figure 8 when Δz > 0.005 m (Figure 11b). For smaller Δz values, the solver CPU time for both methods increased sharply due to larger matrix sizes. When Δz was 0.0005 m, the time consumption of both methods was similar (Figure 11b). Consequently, it is expected that for even finer discretization, such as Δz = 0.00005 m, the fully coupled method may spend more time than the semi-coupled method (Gee & Gracie, 2021; Pimenta & Alves, 2019). However, such a fine spatial discretization is uncommon in practice, and the resulting matrix size is extremely large (60,000 × 60,000 for the fully coupled method, 40,000 × 40,000 and 20,000 × 20,000 for the semi-coupled method), potentially causing an “out of memory” error. Presently, research on isotope transport in soil and plant root water uptake typically focuses on depths around 2 m (Li et al., 2020; Nehemy et al., 2021; Ogle et al., 2004; Rothfuss & Javaux, 2017) with a spatial interval of Δz = 0.01 m (e.g., the lysimeter simulation in this study). Under these conditions, the matrix to be solved typically has a size of 600 × 600. Therefore, the fully coupled method offers advantages in terms of computation time when compared to the semi-coupled method.

Outperformance of MOIST Under the Semi-Analytical Test and the Long-Term Lysimeter Simulation

The fully coupled approach empowers MOIST to exhibit superior accuracy in theoretical tests when compared to the revised HYDRUS-1D (rHZ, Zhou et al., 2021). There are several distinctions between MOIST and rHZ, including model structure and algorithms. However, theoretical testing scenarios with fixed boundary conditions and parameters, as well as known true values, provide a controlled environment that isolates model differences and underscores the impact of the fully coupled and semi-coupled methodologies on the results.

It is worth noting that numerical errors originating from the soil water flow equation can be propagated to the isotope transport equations, potentially leading to the observed numerical oscillations (Figures 5 and 7). This issue has prompted previous models like SiSPAT-Isotope (Braud et al., 2005) and the revised HYDRUS-1D (Zhou et al., 2021) to necessitate an extremely small thickness for the first soil layer (e.g., 1 × 10⁻⁶ m) to pass theoretical tests (Supporting Information S1).

However, in semi-analytical tests where spatial discretization is 0.01 m, error control of rHZ may not be as effective as it is under a finer spatial discretization of 1 × 10⁻⁶ m. This discrepancy resulted in a significant difference between rHZ and the semi-analytical solution. Therefore, it is crucial for a semi-coupled method to rigorously manage mass balance errors (Zha et al., 2019). By contrast, when equations are solved using the fully coupled method, fewer mass balance errors are transferred between equations (Ammara & Masson, 2004; Panday & Huyakorn, 2004; Pascau et al., 1996), a fact that has been verified through temporal and spatial discretization tests (Figures 5 and 7). Consequently, MOIST outperformed rHZ under a coarse spatial discretization in the semi-analytical test (Table 3).

Under the long-term lysimeter simulation, the fully coupled method outperforms the semi-coupled method in simulating isotope transport. Both MOIST and rHZ consider fractionation. Notably, MOIST utilizes parameters that have been calibrated by rHZ, and it exhibits a slight but discernible advantage in performance over rHZ (Figure 10). Moreover, the calibrated parameters from Stumpp et al. (2012) were applied to rHZ, the resulting Nash-Sutcliffe Efficiency (NSE) was 0.19 (Zhou et al., 2021), while MOIST achieved a notably higher NSE of 0.47. This discrepancy, with rHZ yielding a smaller NSE compared to MOIST, supports the idea that the semi-coupled method may introduce more mass balance errors into the process of isotope transport when compared with the fully coupled method. Additionally, as demonstrated by the theoretical test, under the spatial step of 0.01 m (used in the lysimeter simulation), the maximum errors between the semi-analytical solution and the semi-coupled method from rHZ exceeded those from MOIST (Table 4). This evidence, along with the above-mentioned discrepancy underscores the robustness and effectiveness of the fully coupled approach as adopted in MOIST.

Applications for MOIST

As a one-dimensional isotope transport simulator, MOIST contributes to enhancing the precision of calibrated model parameters, such as saturated hydraulic conductivities and isotopic longitudinal dispersivity, as demonstrated in this study (Table 5). This improvement is attributed to the fact that the fully coupled method employed in MOIST reduces errors of simulated isotopic signals. Previous studies confirmed that incorporating isotopic and soil water information can enhance the accuracy of model calibration (Stadnyk & Holmes, 2020; Zhou et al., 2022). However, model accuracy is a fundamental requirement. When the model itself generates numerical errors, these errors are compensated by calibrated parameters, potentially leading to biased parameter values. For instance, Zhou et al. (2022) obtained the calibrated soil saturated hydraulic conductivity was two orders of magnitude smaller than that reported by Stumpp et al. (2012). The greater K_sat from Stumpp et al. (2012) is a result of ignoring fractionation. As is known, fractionation reduces the evaporative isotopic flux. To achieve a satisfactory fit between observations and predictions of outflow isotopic signals, neglecting fractionation forces longitudinal dispersivity to increased and downward convection to decreased (Zhou et al., 2022). However, it is worth noting that the spatial step used in Zhou et al. (2022) was 0.01 m, which could produce significant numerical errors (Table 3). Such errors can lead to overly enriched isotopic signals, thereby amplifying the influence of fractionation (Smith et al., 2016). Consequently, it remains unclear whether the impact of fractionation on model calibration results from fractionation itself or from compensation for numerical errors. By contrast, the fully coupling scheme of MOIST leads to more accurate isotope simulations and consequently has a smaller impact on parameter calibration compared to the semi-coupled method (Table 5). This not only enables MOIST to derive more accurate model parameters but also potentially offers a better understanding of the impact of evaporative fractionation on model calibration.

Moreover, MOIST can be used to test assumptions and validate methods in water transit time estimations (Benettin et al., 2022) and evaporation estimation studies (Sprenger et al., 2017; Wang et al., 2021; Xiang et al., 2021). This is because MOIST can construct high-precision and high-resolution spatial-temporal distributions of isotopic species, which is crucial for studying water transit time (Benettin et al., 2022). Understanding water transit time, also known as water age distribution, is essential for comprehending the connectivity between sources of transpiration, streamflow, groundwater recharge (Nehemy et al., 2022), and the water cycles in the critical zone (Smith et al., 2022). However, challenges persist in obtaining high-precision soil and plant water isotope data (Smith et al., 2022), despite related research efforts (Rode et al., 2016; Von Freyberg et al., 2017). Additionally, various methods for water age estimation based on discrete isotopic measurements often exhibit discrepancies and large uncertainties (Benettin et al., 2022; Borriero et al., 2023). For instance, in the Bruntland Burn catchment, Kuppel et al. (2018a) simulated a watershed water age of approximately 2 years (using the instantons age mixing method), contrasting with 1.2 years (using flux tracing) reported by Soulsby et al. (2015) and approximately 1 year (using StorAge Selection) reported by Benettin et al. (2017). Due to the absence of true values, the accuracy and effectiveness of different methods remain unknown. Addressing this issue, MOIST can generate high-resolution spatial-temporal isotopic profiles, enabling the derivation of transit time from these profiles to serve as ground truth, especially incorporating with recent developed continuous plant water sourcing models, such as CrisPy (Fu et al., 2024) and PRIME (Neil et al., 2024). While other semi-coupled method-based numerical models can perform similar tasks, the accuracy of isotopic profiles may be affected by numerical oscillations (Figures 5 and 7), leading to increased uncertainties. Although MOIST is a one-dimensional model and water transit time studies are primarily conducted at the watershed scale (Benettin et al., 2022), smaller-scale experiments coupled with accurate modeling serve as the foundation for upscaling (Benettin et al., 2021; Douinot et al., 2019). Therefore, MOIST holds promise for transit time studies and cross-validating various water transit time estimation methods.

Regarding the evaporation estimation, MOIST can validate methods for calculating short-term and long-term evaporation based on measured isotopic compositions of soil water. For example, Sprenger et al. (2017) determined evaporation rates from measured soil water isotopic compositions, while Xiang et al. (2021) computed long-term evaporation rates from deep soil isotopic compositions. Both approaches rely on the assumption that sampled soil isotopic compositions at specific depths (shallow or deep soil) retain the evaporation signal. Understanding how much evaporation signal is preserved by the isotopic composition at a given depth is crucial for assessing the accuracy of calculated evaporation rates from methods such as those of Springer et al. and Xiang et al. This assessment necessitates high-resolution (both temporal and spatial) data on soil evaporation rates and soil water isotopic compositions, which can be challenging to obtain. MOIST can bridge this gap through simulations.

Limitations and Future Work

First, in the long-term lysimeter study simulation (Figure 10), a notable observation emerged regarding the simulated δ¹⁸O values. Regardless of whether the model is calibrated or not, or whether it considers fractionation, there was a consistent underestimation of δ¹⁸O between the 1,200th and 1,500th days. This discrepancy can likely be attributed to the phenomenon of preferential flow. During this period, there was a significant influx of precipitation that was relatively enriched in δ¹⁸O, occurring between the 1,050th and 1,300th days. This enriched rainfall contributed to the δ¹⁸O enrichment of water reaching the bottom of the lysimeter. Consequently, a portion of this enriched water was transported to the bottom of the soil column via preferential flow pathways, resulting in an overall increase in δ¹⁸O of drainage water. It is worth noting that in the MOIST, rHZ, and rHS models, preferential flow is not considered for isotope transportation. This omission leads to the systematic underestimation of δ¹⁸O.

Previous research has demonstrated that different modes of soil water mobility, such as preferential and piston flow, as well as mobile and immobile soil water, can significantly influence the isotopic composition of bulk soil water (Finkenbiner et al., 2022; Sprenger et al., 2018). This effect becomes particularly pronounced when substantial rainfall infiltration occurs, potentially giving rise to significant preferential flow events that closely mirror the isotopic signals of the input water. By contrast, the current version of MOIST does not account for preferential flow, resulting in incomplete simulations. To rectify this limitation and better understand hydraulic processes, it is imperative to incorporate preferential flow and its impact on isotopic variations in bulk soil water. Neglecting these factors in isotope simulations can lead to biased interpretations of runoff generation and evaporation (Good et al., 2015). Therefore, future iterations of MOIST should consider implementing a two-pore domain model to provide a more comprehensive simulation of isotope transport processes within the soil. Such enhancements are likely to result in improved simulation accuracy (Sprenger et al., 2018).

Second, the current version of MOIST relies on field measurements from Kurita et al. (2012) to represent $${\delta }_{a}$$ during rain-free period when simulating the long term lysimeter study. By contrast, the approach presented by Zhou et al. (2021) for estimating $${\delta }_{a}$$ is locally suitable. However, this local method was unable to be applied in MOIST because in Zhou et al. (2021), the isotope value of topmost soil layer can be computed from the solution of each time step and used in this local method directly. By contrast, MOIST operates on a cell-centered scheme, resulting in the surface soil water isotopic compositions being unknown at the commencement of each time step. Deriving the isotope concentrations of surface soil water in MOIST necessitates computation based on isotope mass balance between the top-half layer and the isotope value in the atmosphere (refer to Equation 25). Consequently, the formula proposed by Zhou et al. (2021) is inapplicable in MOIST as one equation cannot resolve two unknowns simultaneously. In the future iterations of MOIST, we aim to incorporate the equations for local $${\delta }_{a}$$ to better describe isotope transport processes.

Lastly, the primary objective of this study is to assess the performance of MOIST under various conditions and demonstrate the performance of the fully coupled method. Consequently, during the model comparison, we refrained from conducting an exhaustive calibration for MOIST. Notably, with the inclusion of calibrated saturated hydraulic conductivities, our model exhibited smaller errors compared to other models. Furthermore, the calibrated ranges of saturated hydraulic conductivities were smaller than the initial values (as shown in Table 5), aligning with the findings of Zhou et al. (2022) under the same data set. This alignment underscores the acceptability of the calibrated saturated hydraulic conductivities and longitudinal dispersivity. It instills confidence that a more comprehensive calibration process could further enhance the capacity of MOIST for accurate simulations.

We concur with the notion that for practical model applications, comprehensive calibration using independent data is indispensable (Kuppel et al., 2018b; Tsai et al., 2021). This necessity arises due to the inherent uncertainties associated with many parameters. Calibration serves the purpose of refining the model by imposing constraints, resulting in a more accurate representation of potential hydrological processes. Therefore, in forthcoming iterations, we intend to implement comprehensive calibration as an integral step toward practical applications.

Conclusions

We developed MOIST, a novel soil water and isotope transport model using MATLAB programming language. MOIST is unique in that it solves water, vapor, heat, and isotope transport simultaneously with a fully coupled method and with a cell-centered numerical approximation to the derivatives. MOIST successfully passed well-known theoretical tests, and semi-analytical solutions. Due to its adoption of fully coupled solution method, MOIST can predict isotope profiles under various temporal and spatial discretization with greater accuracy and stability than the existing semi-coupled models. We also tested MOIST against well-controlled long-term lysimeter studies at HBLFA Raumberg-Gumpenstein in Austria. Results showed good agreement between measured and predicted values and outperformed HYDRUS-1D. Given the evidence accrued in this study, it can be concluded that MOIST can be a robust tool for simulating one-dimensional isotope transport within soil. Moreover, the model is open-source and thus can be customized according to different requirements. As such, the adopted equations and parameters can be easily updated as our knowledge about isotope fractionation and transport continues to expand, making MOIST a suitable tool for current and future exploration of isotope transport in the soil-vegetation-atmosphere continuum.

Acknowledgments

This research was supported by the China Scholarship Council (Grant number 201906300103), and the Shandong (China)-Israel Cooperation Program (2023KJHZ007). We thank Dr. Andrew Ireson for the numerical simulation lectures. We also acknowledge editors and reviewers for improving this manuscript.

Data Availability Statement

Source codes of MOIST (Fu & Si, 2023) used for semi-analytical tests, theoretical tests, calibration, and the long term lysimeter study tests were developed by MATLAB R2022a and available via Creative Commons Attribution 4.0 International license through . The long term lysimeter measurements are available in Stumpp et al. (2012).