The vadose zone (often called the unsaturated zone) is commonly defined as the geologic media spanning from the land surface to the groundwater table of the first unconfined aquifer (Stephens 2018). The vadose zone is known to play a critical role within the biosphere: (1) as a storage medium to supply water to the plants and atmosphere, and (2) as a controlling agent in the transmission of recharging water as well as contaminants from the land surface to groundwater (Nimmo 2005). Past interest in the vadose zone has largely been a result of public concern about the need to protect groundwater reserves for drinking and agricultural purposes. The main sources of groundwater contamination typically originate in the vadose zone as leaking underground storage tanks, municipal solids and hazardous waste landfills, waste management sites, unlined pits, ponds, and lagoons, household septic systems, pesticide application areas, or surface spills (LaGrega et al. 1994). More importantly, these contaminants are modified by interactions within the vadose zone. For example, Małecki and Matyjasik (2003) reported that the average total dissolved solids values for rainwater changed dramatically from 30.2 mg/L at the land surface to 318 mg/L in groundwater. They also noted a corresponding change in the water type from a SO4–Cl–Ca–NH4 to a HCO3–SO4–Ca–Mg type due to vadose zone hydrogeochemical interactions. This and other studies document the significant role of vadose zone in contaminant distribution and migration (e.g., Oostrom et al. 2016; Arora and Mohanty 2017; Wan et al. 2018). Research efforts have thus far focused on the vadose zone as a significant reservoir for the capture, storage, and release of contaminants.

Since vadose zone involves numerous coupled physical, geochemical, and microbial processes, the interplay of these processes is difficult to understand especially from a remediation perspective. Considerable attention has been given to characterize and quantify contaminant transport in subsurface heterogeneous media. For example, heterogeneity caused by macropores and fractures can lead to preferential flow and rapid infiltration, which, in turn, increases groundwater vulnerability to potential contamination. In contrast, clay lenses and layers can provide opportunities to sorb contaminants in place. More recently, a growing body of work has focused on the outsized impacts of nutrient cycling and higher reaction rates at ritical interfaces. Herein, critical interfaces are defined as the interacting boundaries between zones of distinct hydrological, biogeochemical and lithological properties (Li et al. 2017). Biogeochemical processes and biodiversity at these interfaces are often orders of magnitude higher than the surrounding matrix. For example, iron- and sulfate-reducing bacteria showed one to two orders of magnitude greater community numbers in a layered sand-over-loam column when compared to homogeneous sand and loam columns (Hansen et al. 2011). McGuire et al. (2005) studied the impact of a moderate-sized rainfall event on redox processes at a shallow, sandy aquifer contaminated with petroleum hydrocarbons and chlorinated solvents, and concluded that recharge effects on the progression of terminal electron accepting processes existed primarily at the interface between infiltrating water and the aquifer, and not at the average aquifer scale. Quantification of these critical interfaces and plume evolution are further complicated by transient flow and transport conditions in the vadose zone. Transient conditions, including seasonal, annual, and long-term variability in recharge as well as transient interactions between climatic, hydrologic and biogeochemical factors complicate the study of contaminant fate and transport in the vadose zone. For example, Arora et al. (2016) concluded that changes in concentrations of redox-sensitive chemicals in a shallow alluvial aquifer were directly related to rainfall and recharge events at bi-monthly and annual time scales. Han et al. (2001) determined that the wetting-drying moisture regime resulted in higher metal reactivity in arid regions, resulting in redistribution and fractionation of heavy metals such as Ni, Zn and Cu. Although these studies have clearly demonstrated that critical interfaces and transient conditions regulate biogeochemical reaction rates, it remains unclear how to characterize these interfaces and time periods, and quantitatively incorporate them in numerical models.

To capture complex vadose zone characteristics, significant advances have been made in subsurface characterization techniques which are non-invasive, have finer resolution and capture a large spatial extent (Wainwright et al. 2018). In particular, these include in situ water, solute and gas monitoring, geophysical and remotely sensed observations, as well as weather data, which all offer diverse sets of data. The data stream, however, comes with its own set of issues, such as sparse, heterogeneous and uncertain measurements, a variety of data sources, and differences in frequency of measurements. The major open questions related to these datasets include: (i) how to incorporate such diverse data types into reactive transport models; (ii) which methods to use to separate bulk properties from small scale heterogeneities and interfaces; and (ii) how to characterize error characteristics and incorporate this uncertainty into model predictions. The geoscience community is embracing the potential of additional data to improve model structure and therefore create opportunities for targeted sampling to better understand vadose zone dynamics (Mcmillan et al. 2012; Finsterle 2015).

To analyze complex environmental pollution problems, a quantitative study of water flow and contaminant transport in the vadose zone can be valuable and complementary to measurements of the heterogeneous and complex subsurface, as complicated by transient conditions. The past decade has seen a significant development of flow and reactive transport models as indispensable tools for studying vadose zone processes, including incorporating vadose zone complexities (e.g., sharp permeability interfaces) and transient conditions (e.g., climatic variability) to support decision making in such diverse areas as assessment of acid mine drainage, evaluation of the consequences of managed aquifer recharge, safety analysis of nuclear waste repositories, and analysis of global climate change. However, major gaps and challenges remain with regard to further conceptualization of vadose zone complexities, incorporation into and parameterization of models.

The focus of this chapter is to present a combination of modern approaches for vadose zone characterization and modeling. The key questions addressed here are how modeling activities can better serve quantification of vadose zone complexities, including heterogeneities and critical interfaces (e.g., Fig.1), while accounting for transient phenomena, and what areas and key challenges need to be addressed to improve the applicability and usefulness of current vadose zone models.

The chapter is organized as follows. The next section provides a review of field-based measurements to acquire diverse data streams relevant to the vadose zone. This is followed by a description of geochemical modeling of the vadose zone. At this point, the focus shifts to describing soil moisture studies from measurements to modeling. From here, the chapter proceeds to investigate some case studies which illuminate the power of geochemical modeling in deciphering vadose zone complexities and interactions. We then provide a brief discussion of current model limitations and future opportunities.

Vadose zone characterization and monitoring methods are commonly used for the development of a complete and accurate assessment of the inventory, distribution, and movement of contaminants in the unsaturated zone; development of improved predictive methods for liquid flow and contaminant transport; design of remediation systems (barrier systems, stabilization of buried wastes in situ, cover systems for waste isolation, in situ treatment barriers of dispersed contaminant plumes, bioreactive treatment methods of organic solvents in sediments and groundwater); design of chemical treatment technologies to destroy or immobilize highly concentrated contaminant sources (metals, radionuclides, explosive residues, and solvents) accumulated in the subsurface.

In the world of widely applied modern remote sensing and geophysical technologies for observations of vadose zone and groundwater processes, we need local scale measurements (Rubin et al. 1998; Faybishenko 2000), which is often called using a Latin word in situ—literally means on site or at the original place. Local hydro-geophysical methods are usually minimally invasive, conducted at small scales < 1 m, and are used to supplement and constrain larger scale geophysical and remote spatial-temporal measurements. The goal of this section is to provide a synoptic view and a basic guide for common in situ vadose zone and groundwater measurements, as well as large scale remotely sensed observations, which are essential for hydro-geophysical characterization of the near-surface hydrological and biogeochemical properties.


Conceptualization, or in other words—development of a conceptual model, is an approach to select why, what, and how to design and conduct measurements. Conceptualization of vadose zone systems is needed for the integrated qualitative and quantitative characterization of unsaturated flow and transport processes affected by natural behavior and man-induced changes. Development of appropriate conceptual models of water flow and chemical transport in the vadose zone is critical for developing adequate predictive modeling methods and designing cost-effective remediation techniques. Conceptual models of complex and heterogeneous vadose zone must take into account the processes of preferential and fast water seepage and contaminant transport toward the underlying aquifer. Such processes are enhanced under episodic natural precipitation, snowmelt, and extreme chemistry of waste leaks from tanks, cribs, and other surface sources. However, until recently, the effects of episodic infiltration and preferential flow on a field scale have not been taken into account when predicting flow and transport and developing remediation procedures. The pronounced temporal and spatial structure of water seepage and contaminant transport, which is difficult to detect, poses unique and difficult problems for characterization, monitoring, modeling, engineering of containment, and remediation of contaminants. Lack of understanding in this area has led to severely erroneous predictions of contaminant transport and incorrect remediation actions.

Because hydrological and biogeochemical cycles differ under conditions of the natural environment, for example, at the watershed scale, irrigation of agricultural lands, at contaminated sites, and urban areas, different types of conceptual models need to be developed, and different types of local scale measurements are needed for different environments.

A typical feature of the near-surface hydrological processes in different environments is preferential flow phenomenon (Allaire et al. 2009; Robinson et al. 2008). Preferential flow, often called fingering phenomena, is caused by the near-surface heterogeneity, cracks, and local surface depressions. Water moves faster through the preferential flow zones. Preferential flow is one of the main reasons of the spatial heterogeneity of the moisture content in soil. Due to the fingering effect, one would reasonably expect significant spatial and temporal variations of the soil moisture content and water potential in the near-surface zone. Note that more significant changes occur down to about 6–8 m depths.

To understand the vadose zone processes, one needs to measure both the soil moisture content and water pressure. The relationship between the moisture content and water pressure is fundamentally important to understanding soil water behavior. The graph of this relationship is commonly called the water retention curve.

In situ measurements in the vadose zone

Collecting soil samples.

Soil investigations involve collecting soil samples—either disturbed or undisturbed soil samples, called monoliths, using special rings. Sensors are installed in boreholes—near-surface drilling can be done using hand held augers, usually to 2–3 m depths. Deeper drilling is carried out using mechanical augers, which are often integrated with electrical and gamma logging tools. The main requirement is to minimally disturb the soil during drilling of boreholes and sample collection.

Soil moisture measurements.

Soil moisture is widely recognized as a key parameter in the mass and energy balance between the land surface and the atmosphere. The importance of soil moisture can be evident from the fact that researchers around the world are collecting soil moisture data, and nearly four decades of soil moisture data are now available from the International soil moisture network (Dorigo et al. 2011). Here, we describe only in situ measurements of soil moisture. A detailed description of remotely sensed observations to obtain soil moisture content and incorporate it into modeling is described in later sections.

Methods of measurements and collection of water samples are different depending on the moisture content. For example, one can easily collect water samples in the presence of so called gravitational water, but need to apply vacuum to withdraw water when the moisture content is below the residual water content. Deep lysimeters are used for collecting water samples below the depth of ~ 5 m.

The standard reference method for determining the soil moisture content is to oven dry soil samples at 105 oC. For organic soils and gypsiferous soils, the temperature should not exceed 70 oC (because of the problem with organic matter volatilization). Disturbed soil samples are commonly used to determine the gravimetric moisture content of soil. The moisture content is often expressed in volumetric units—gravimetric moisture content times the soil dry density will give you a volumetric moisture content.

Neutron moisture meter.

A neutron moisture meter is based on the principle of utilizing neutron scattering around the borehole. Due to its hydrogen content, water is an effective neutron moderator, slowing high-energy neutrons. The technique is non-destructive, and is sensitive to moisture in the bulk of the target material.

Time Domain Reflectometry (TDR).

TDR is a relatively new method for measurement of soil water content, first reported by Topp et al. (1980). It is based on the standard electromagnetic (EM) method for determining the moisture content. An electromagnetic pulse is generated by the TDR cable tester, and the propagation velocity the electromagnetic wave is measured along a transmission line (wave-guide) embedded in the soil. The support volume of TDR measurements depends on the probe design (Ferré et al. 1998; Skierucha et al. 2012). The main advantages of the TDR method over other methods for repetitive soil water content measurement (e.g., neutron probe) are: accuracy to within 1 or 2% of volumetric water content; calibration requirements are minimal (in many cases soil-specific calibration is not needed); averts radiation hazards associated with neutron probe or gamma-attenuation techniques; and capable of providing continuous soil water measurements through automation and multiplexing.

Capacitance/frequency domain technology.

ECH2O EC-5 is a simple sensor with excellent accuracy based on measuring the dielectric constant of soil based on the capacitance/frequency domain technology. Advantages are: minimum salinity and textural effects, factory calibration, small area of influence, and high resolution.

Figure 2 shows a schematic of the operational ranges of field and laboratory methods used in monitoring the matric suction in the vadose zone for soil water physical processes (Gee and Ward 1999).

Water pressure measurements.

A tensiometer is used to directly measure the soil water pressure. The tensiometer consists of a porous tip/cup, a water-filled tube and a vacuum gauge. The porous cup of the tensiometer is buried in the soil. As water is pulled out of the tensiometer, the vacuum inside the tube increases. Soil suction creates the vacuum inside the tensiometer, which is measured by the gauge. When the vacuum created in the tensiometer corresponds to the soil suction, the water outflow from the tensiometer is ceased. The soil water suction is also called the soil water pressure. In conjunction with a water retention curve, tensiometers can be used to determine the soil moisture content. Tensiometers are used in irrigation scheduling to help farmers and other irrigation managers to determine when to water.

There are different types of porous cups and different designs. Theoretically, tensiometers with porous cups can be used in the range of water pressure from 0 (that is the reference value indicating the atmospheric pressure) to −1 bar, but practically to only −0.7–0.8 bars. There were recently developed the probes to measure the water pressure to −5 bars.

Water sampling.

Water sampling is based on the idea of collecting soil water into the cup by applying vacuum inside the porous cup to exceed the capillary pressure holding water in soil. For shallow depths to ~ 5–7 m, you can bring the collected water to the surface using vacuum. For deeper depths, we need to first apply a vacuum to draw water inside the lysimeter and to push it to the surface by applying positive air/gas pressure.

Soil gas measurements.

Soil gas is collected by applying vacuum through simple tubes installed in a borehole at different depths. Gas samples are collected into an impermeable container. Soil gas sampling is used to detect the source and locations of contaminated soils or groundwater. A soil gas probe can be connected to a preliminary vacuumed canister called Summa canisters (

Evaporation and evapotranspiration.

Evaporation is a vapor flux from the open water surface. Evaporometers for direct measurements of evaporation are installed at many meteorological stations around the world. Evapotranspiration is a term to define a cumulative vapor flux from both the surface and transpiration by plants. The most advanced apparatus to measure evapotranspiration is a weighing lysimeter ( This is a tension controlled soil column with Tensiometers and soil moisture sensors to create the water and temperature conditions that are similar to those in the field. The lysimeter is installed on the scale to monitor the weight to calculate the total moisture content in the soil. For example, in Germany, there is a national lysimeter network for observing soil evapotranspiration called Tereno (Pütz et al. 2018). The observatory sites are located in different climatic conditions, characterized by different precipitation, elevation, and temperature (

Eddy Covariance Stations.

Measurements of gas exchange between terrain and the atmosphere are conducted at the Eddy Covariance Stations, which are part of the AmeriFlux monitoring network. There are 395 registered AmeriFlux sites. Measurements include monitoring of CO2, vapor, temperature, wind, humidity, and precipitation. The main idea is measurements of high frequency wind and concentrations, which are then used to calculate fluxes. The detailed description of the Eddy Covariance stations can be found at or

Surface runoff and soil loss due to erosion.

One of the most important but difficult types of measurements is the evaluation of the surface runoff and associated soil loss due to erosion, for example, due to flooding. Measurements can be conducted using specially designed sites with different slope gradients and lengths, where water and sediments are collected at the exit from site.

Stable isotopes.

One of the most effective methods to assess evapotranspiration is the use of measurements of stable isotopes—deuterium and oxygen-18 in soil water (Böhnke et al. 2002). The idea is that the stable hydrogen and oxygen isotopic compositions of fresh waters are highly correlated according to the ‘global meteoric water line.’ At some locations there is a small departure from this line, and the departure is even more increased due to evapotranspiration. So the departure from the local meteoritic water line can be used to assess evapotranspiration.

In situ measurements in the saturated zone

Groundwater level measurements.

Groundwater level measurements are conducted in boreholes, which are usually equipped with the metal or plastic casing. Water enters the borehole through the screen being installed near the bottom of the well. Water level measurements are often accompanied by the barometric pressure measurements to make corrections on the effect of barometric pressure. The water level can be measured with a tape or using a water pressure transducer connected to the data acquisition system. A multi-parameter monitoring system, consisting of several sensors, can also be installed in a borehole to make measurements of different parameters.

Water sampling from boreholes.

The chemical composition of a groundwater sample collected from a monitoring borehole is dependent on two factors: the length of the screened interval in the borehole, and the variability in permeability of the adjacent strata. During water sampling, the cone of depression of the water table is developed around the borehole. The radius of influence of pumping is largest at the depth of high permeability layer. In this sense collected samples can contain mixture of concentrations from different formations. Inflows from more permeable strata or fissures will dominate and bias the sample in the borehole and produce a flow-weighted average sample (FWA). In this case, an unrepresentative water sample is collected. Mixing mechanisms within boreholes can be very complex and research demonstrates that boreholes with longer screened intervals lead to greater bias, and therefore uncertainty, in the sampling result. Inflows from more permeable strata or fissures will dominate and bias the sample in the borehole and produce a flow-weighted average sample (FWA). Mixing mechanisms within boreholes can be very complex and research demonstrates that boreholes with longer screened intervals lead to greater bias, and therefore uncertainty, in the sampling result.

Remote sensing and geophysical observations

A suite of properties are expected to influence the vadose zone flow and biogeochemistry. Examples include the density, type, and distribution of plants and their dynamics; geomorphology and microtopography; the texture, density of near-surface soils; pore water chemistry and connectivity. Many of these properties can be measured through in situ sampling, but with the caveat that it is difficult to distribute these properties over space. Recently, there have been a variety of new technologies developed to capture the heterogeneity of the surface and subsurface properties—relevant to the vadose zone processes—by taking advantage of the state-of-the-art spatially integrating techniques—such as fiber optics sensing, geophysics, and unmanned aerial vehicle (UAVs).

Geophysical methods.

Geophysical methods—including electrical resistivity, seismic, and radar—have been increasingly used to characterize the vadose zone in a noninvasive manner (e.g., Rubin and Hubbard 2005; Binley et al. 2015). They can be used to image soil moisture and biogeochemical properties (e.g., Dafflon et al. 2011, 2013; Johnson et al. 2012, 2010; Wainwright et al. 2014). Autonomous electrical resistivity and phase tomography (ERT) monitoring, in particular, has the potential to achieve rapid and automated detection and identification of water infiltration within the vadose zone (Tran et al. 2016). In addition, recent developments in hydrogeophysical inversion methods allow for near-real-time characterization of vadose zone infiltration (e.g., Johnson et al. 2015). The inversion method also enable the estimation of soil hydraulic properties (Tran et al. 2016), which are critical for model predictions. These geophysical approaches can bridge the gap in sparse wellbore locations, by providing high-resolution and spatially and temporally extensive information in a minimally invasive manner (e.g., Johnson et al. 2015).

Remote sensing observations.

Multiple platforms exist to remotely sense soil moisture and other landscape features for characterizing vadose zone processes. These include satellite, airplane-based and drone-based approaches. Among these, the UAV technologies are the latest growing area in the applications for environmental monitoring. The high-resolution surface elevation data from Light Detection and Ranging (LiDAR) or Photo Detection and Ranging (PhoDAR), for example, provides the surface elevation in a centimeter resolution which enables to map surface runoff patterns and infiltration. A recent study combined the UAV-based surface elevation mapping with isotopic analysis, and was successful in identifying surface flow accumulation and fast flow paths into the vadose zone (Christensen et al. 2018). In addition, thermal cameras are increasingly being used to map the locations of groundwater seep zones and river-groundwater interactions. Moreover, coupled spectral and structural information from UAV images can map the heterogeneity of plant species, photosynthetic activity as well as evapotranspiration, which are key components of vadose zone dynamics.

Fiber optics.

Fiber optics technology uses optical fibers for measuring environmentally relevant parameters including temperature (Suarez et al. 2011), strain, soil moisture, acoustic waves (Bao and Chen 2012; Cox et al. 2012), and some chemical signatures (e.g., Potyrailo and Hieftje 1998). Each pulse of light samples the state of the fiber at all locations, yielding property measurements along the entire length at a fine lateral resolution (~ 0.25–1 m). In contrast to autonomous point sensors, the fiber optics technology can acquire real-time high-frequency datasets combining large extent and fine space/time sampling.

Field vadose zone characterization and monitoring methods are commonly used to constrain the results of remote geophysical and air-borne measurements, and to provide scientific and technical background to predict contaminant migration. By incorporating a combination of geophysical cross-borehole tomography and local-type measurements, numerical models can be developed to simulate the movement of organic and radioactive contaminants and their exchange between the zones of preferential flow and the soil or rock matrix. For example, the zonation approach (described below) serves as an important framework to include vadose zone complexities into models.

Over the last 20 years, there has been a growing reliance on reactive transport modeling (RTM) to address some of the most compelling issues facing our planet: nuclear waste management, contaminant remediation, and pollution prevention. While these issues are motivating the development of new and improved capabilities for subsurface environmental modeling using RTM, there remain longstanding challenges in characterizing the natural variability of hydrological, biological, and geochemical properties in subsurface environments and limited success in transferring models between sites and across scales. To achieve these ambitious objectives for subsurface reactive transport simulation, the subsurface science and engineering community is being driven to provide accurate assessments of engineering performance and risk for important issues with far-reaching consequences. In this regard, we present an overview of vadose zone complexities and how to incorporate these into models.


To develop predictive capabilities of the vadose zone processes, a conceptual framework is necessary. Conceptual models typically describe the main physical processes taking place in the vadose zone, and how these processes interact or dominate the system. In general, the conceptual physical framework includes various processes and features such as a hydrologeologic setting describing the hydro-stratigraphic and structural details and relevant hydrological and biogeochemical processes. Relevant hydrogeological processes include surface infiltration rates (e.g., precipitation and evapotranspiration rates), flow rates (lateral and vertical) and details of the geologic formation. At the same time, relevant biogeochemical details include the presence of major reactive minerals, sediment characterization, pore water chemistry.

Because the presence of fast and slow flow paths significantly impact hydrological and biogeochemical processes in the vadose zone, it is important to capture the spatial and temporal signatures of vadose zone dynamics in the conceptual and consequently numerical frameworks. In particular, incorporation of preferential flow paths, clay lenses, and hydrostratigraphic layers with different properties need to be adequately represented in the numerical model. Several numerical solvers (e.g., TOUGHReact, PFLOTRAN, CrunchTope) exist that can be used to describe the hydrologic and biogeochemical implications of these vadose zone complexities. In the following section, we will briefly describe how to translate a conceptual framework into numerical models using appropriate parametrization depending on the level of complexity.

Description of the hydrological system

Variably saturated flow can be described using a set of equations for multiphase flow, explicitly considering liquid and air in the vadose zone (Steefel et al. 2015). Alternatively, the Richards equation can also be used for modeling of variably saturated flow in the vadose zone, which assumes the liquid phase is active while the gas phase is passive (Richards 1931). The Richards equation can be formulated either on the basis of aqueous phase saturation or fluid pressure. The Richards equation has been widely used in the literature for representing vadose zone processes (e.g., Yabusaki et al. 2017; Dwivedi et al. 2018a,b). However, a wide variety of problems such as nuclear waste disposal, acid rock drainage, and geological CO2 sequestration may require multi-phase formulation. In the multi-phase formulation, the mass conservation equation requires the solution of each component present in different phases.

Simulations of hydrological and biogeochemical processes in heterogeneous soil and fractured-porous media require an application of the multi-continuum class models. In general, the multi-continuum approach includes different interacting regions with different hydrogeologic properties. For example, the dual-continuum approach includes two interacting regions—soil matrix and fracture domains; soil matrix is characterized with the less permeable intra-aggregate pore region, while the fracture domain is associated with the inter-aggregate pore region. Some examples of multi-continuum class of models include (1) the equivalent continuum model (ECM), which can incorporate the effects of fractures (Faybishenko et al. 2000); (2) dual permeability model (DPM) and dual, multi-porosity model, and multiple interacting continua approach (MINC)—representing fractures as regions of high permeability-low porosity in heterogeneous porous media (Barenblatt et al. 1960; Warren and Root 1963; Pruess 1985); and (3) discrete fracture and matrix model, which can represent a single fracture or an infinite number of equally spaced fractures in the porous media (Snow 1965). More details on different multi-continuum approaches including exchange functions and governing equations can be found in the literature (Faybishenko et al. 2000; Šimůnek et al. 2003;Arora et al. 2011).

One of the main parameters included in numerical models is the unsaturated hydraulic conductivity that is a nonlinear function of capillary pressure and liquid saturation; permeability values typically decrease as the water content and/or liquid pressure decrease (Mualem 1976; van Genuchten 1980). Brooks-Corey and van Genuchten formulations have been extensively used in the literature to relate the relationship between the capillary pressure and aqueous phase saturation. Additionally, the Mualem and Burdine formulations are used for computing the relative permeability as a function of the effective aqueous saturation (Brooks and Corey 1964; Mualem 1976; van Genuchten 1980).

Description of the geochemical system

The key geochemical processes relevant in the vadose zone are aqueous chemistry, aqueous speciation, mineral precipitation and dissolution, sorption, ion exchange and redox reactions. To evaluate biogeochemical implications in the vadose zone, solute transport and reactions are coupled in the continuity equation (Dwivedi et al. 2016; Steefel et al. 2005). These reactions represent the transformation of reactant to product species. Reactions can either be classified as homogeneous or heterogeneous depending upon whether they occur within the same phase (e.g., aqueous) or across different phases (e.g., mineral reactions). Nevertheless, these reactions involve multiple geochemical species including major cations, anions, and trace elements. Correspondingly, a set of primary species is defined that determine the geochemical system. At the same time, a set of secondary species is defined on the basis of different combinations of primary species. The reactive transport codes can sweep through the database to pick secondary species, or secondary species can be explicitly defined in the input files as well. To compute the concentration of the secondary species, a geochemical equilibrium between aqueous species is assumed.

Biogeochemical reactions can either be kinetically controlled or equilibrium based. In general, mineral precipitation and dissolution, homogeneous reactions (aqueous phase), as well as microbially mediated reactions are kinetically controlled. On the contrary, sorption or surface complexation are equilibrium-based reactions. A kinetic treatment is more general; however, the type of treatment (kinetic vs. equilibrium) is contingent upon the relative time scales of the processes involved (Steefel and Lasaga 1994). For example, if the rate of mineral dissolution is significantly faster than the transport rate of a reactant, then it is reasonable to assume equilibrium-based mineral dissolution. An in-depth review on kinetically controlled and equilibrium based reactions is available elsewhere (Steefel et al. 2005; Steefel and Maher 2009; Li et al. 2010).

Model domain and discretization

Temporal discretization.

Most of the modern reactive transport models have adaptive time steps. However, the maximum time step can be chosen for simulations based on the Courant–Friedrichs–Lewy (CFL) condition. The CFL condition ensures that the distance between mesh elements exceeds the distance that any particle travels during the time step (Sonnenthal et al. 2014).

Spatial discretization.

At the watershed-, basin- or catchment-scales, the spatial discretization of the grid can be as fine as the digital elevation model (DEM). DEM can be derived from a high-resolution LiDAR dataset. DEM data are typically high-resolution (sub-meter) and can describe microtopographic features as well. However, it may be numerically intractable to run such simulations at these large scales. In that case, variable resolution meshes are preferred. In addition, mesh can be structured or unstructured. There are various tools that can create these meshes such as LaGrit, CUBIT, and Meshmaker (Hammond et al. 2014; Sonnenthal et al. 2014; Coon et al. 2016).

Material properties.

As suggested above, hydrostratigraphic details are essential for appropriately representing the conceptual framework in numerical models. For example, details such as the thickness of the alluvial deposits and fluvial sediments including debris-flow can be spatially variable, and the soil properties can also vary in texture from sandy, silty to gravelly. These variations have a significant impact on flow and transport properties. Hence, different hydrostratigraphic layers need to be incorporated in the mesh by assigning different permeability values. The average hydrological properties (e.g., van Genuchten parameters) including permeability values can be estimated using pedotransfer functions from textural data of sediments from the site. Pedotransfer functions are described in more detail in the section.

An adequate representation of heterogeneity in model domain is also important to capture spatial and temporal signatures of vadose zone processes. In this regard, geophysical data and zonation approach are useful to map hydrogeological and geochemical properties such as permeability, percent of fines, and mineral ratio in the subsurface environment. In addition, reactive minerals can be characterized using XRD or digestion methods on soil samples. These properties can be assigned on the mesh by identifying appropriate regions. These regions include a set of control volumes having different material properties including soil types (e.g., van Genuchten parameters), permeability values, and minerals.

Finally, it is important to realize that the geochemical transformation of compounds in particular mineral reactions can alter the porosity and permeability of the porous media. Small changes in porosity can lead to substantial changes in permeability, and subsequently, these alteration in material properties are likely to modify the unsaturated flow properties of the porous media (Vaughan 1987). Most reactive transport codes can simulate these changes in material properties by computing volume changes of the matrix and fractures (Steefel et al. 2015; Xie et al. 2015).

Initial and boundary conditions

To evaluate vadose zone dynamics, reactive transport models require appropriate boundary conditions incorporating spatio-temporal variability of hydrologic and biogeochemical signatures. Hydrologic boundary conditions can be simulated using Dirichlet (or first-type) boundary condition, also known as hydrostatic or seepage face boundary conditions as shown in Figure 3.

For example, a river can be simulated using a dynamic seepage face boundary condition by applying observed transient river-stage measurements. A seepage face boundary is defined as an interface at atmospheric pressure (Patm), where the water exits the porous medium (Bear 1975). On the contrary, the hydrostatic boundary condition is applied at the interface of the water table in the vadose zone, where negative pressure exists along the face of the model domain (Hammond et al. 2014; Sonnenthal et al. 2014).

It is essential to understand that redox or geochemical conditions change in response to perturbations such as precipitation events. Usually, recharge varies in both in space and time due to the topographic effects that cause substantial variability. Spatio-temporal recharge is typically applied as Neumann (or second-type) boundary condition (Hammond et al. 2014).

Pore water chemistry also changes with time; therefore, it is important to consider variations in the pore water composition to define geochemical boundary conditions in the model. It is possible to apply time-varying geochemical boundary conditions by assigning different pore water chemistry as a function of time in modern reactive transport codes. Moreover, it is important that models are appropriately initialized and run with proper meteorological, hydrological, and geochemical inputs. Typically, models should be run long enough so that model results are independent of initial conditions.

As described above, soil moisture is a key variable controlling water, energy and nutrient fluxes from the land surface to the atmosphere. Although small in proportion (~0.15% of global freshwater), it is still an influential store of water in the hydrologic cycle (Dingman 1994; Western et al. 2002). Several studies have documented the role of soil moisture in modulating climate interactions, rainfall-runoff processes and plant growth. In the context of biogeochemical transport through the vadose zone, soil moisture is often the main variable controlling contaminant flow and transport. However, at the scale of the global hydrologic cycle, the processes governing water flow and solute transport through the vadose zone are the least understood (Wu and Li 2006). This is because a number of complexities such as, heterogeneity in soil texture, structure and composition, as well as transient environmental conditions can change soil moisture status in the vadose zone (Jana et al. 2007; Jana and Mohanty 2012a; Vereecken et al. 2014). Given the importance of the effect of soil moisture on the hydrologic cycle and related processes such as contaminant transport, agricultural runoff, crop production and flood control, it is obvious that concerted efforts must be made to measure, understand, and model soil moisture dynamics. Several studies have demonstrated that a better depiction of the spatio-temporal heterogeneity of soil moisture results in more accurate estimates from hydrologic and related process models (Koster et al. 2004; McCabe et al. 2005; van den Hurk et al. 1997). Here, we describe issues and complexities related to measurements and modeling of soil moisture within the vadose zone.

Soil moisture measurements

Soil moisture measurements, made in situ or by remote sensing, typically involve some trade-off between the scale triplet—support, spacing, and extent (Blöschl and Sivapalan 1995). Herein, support is the area (or time) over which the measurement is valid; spacing is the interval at which measurements are made; and extent is the overall coverage of the domain where measurements are made (Skøien and Blöschl 2006). In situ measurements offer support dimensions in the order of centimeters. However, it becomes impractical to have spacing of measurements in the same order, and to have such measurements over large extents. Intermediate support resolution (of the order of a few hundred meters to a kilometer) can be obtained using sensors mounted on airplanes or UAVs, but the extent of their coverage is typically limited to a few tens of kilometers. Satellite-based remote sensors, on the other hand, provide data with abundant (mostly global) coverage, but their support dimensions are of the order of hundreds and thousands of meters for soil moisture. Such poor resolution results in smoothing out much of the soil moisture heterogeneity in the footprint, which results in the omission of critical information from a remediation perspective. Recently, certain satellite-based sensors such as NASA's Soil Moisture Active Passive (SMAP) instrument were launched with the aim of providing soil moisture information at finer resolution (~3 km) (Entekhabi et al. 2010). However, these active (radar) sensors are no longer active due to technical issues. Other possible sources of relatively fine resolution remotely sensed soil moisture products can be obtained from Sentinel-1 and forthcoming missions such as NISAR and Tandem-L. These new efforts are described in more detail elsewhere (Mohanty et al. 2017).

In addition to spatial scales, temporal resolutions of the data also vary based on measurement platforms. The quickest return time of satellite remote sensors is of the order of 1-2 days. This means that information regarding intra-day variations in soil moisture are lost and only discrete snapshots of the moisture state are available. As a result, the immediate effects of certain low-intensity, short-duration precipitation events could go unobserved. In comparison, in situ sensors can provide almost real-time data as a continuous stream. Based on the type of application for soil moisture, different frequency of time series data may be preferred.

Time-stable locations.

Large spatial extent, coarse resolution soil moisture measurements from satellite-based sensors are often calibrated and validated using in situ measurements. However, due to the mismatch in support scales and the inherent complexity of the vadose zone (such as, variations in topography, vegetation, soil physical properties and meteorology), a comprehensive validation becomes difficult to achieve (Cosh et al. 2004). Further, due to practical considerations, the number of in situ measurements is often limited in terms of both spatial as well as temporal spread. Such a situation necessitates the identification of locations within a remote-sensing footprint that can generally describe the mean soil moisture behavior within the footprint that is invariant over a sufficient period of time (Martínez-Fernández and Ceballos 2005). First proposed by Vachaud et al. (1985), such locations are referred to as time stable or temporally stable locations. Identification of such locations can help greatly reduce the number of locations to be sampled, and also enable permanent installation of in situ sensors for long-term studies relating measurements at the two disparate scales of measurement. It must be noted, however, that information regarding the soil moisture variability across the watershed is lost when measuring at only time stable locations. In order to overcome this drawback, newer non-linear model reduction techniques such as mode decomposition and discrete empirical interpolation method have been applied (Chaturantabut and Sorensen 2009; Hu and Si 2013). In particular, these newer approaches extend the concept of time stable locations to arrive at locations that provide not only the average soil moisture values for the area of interest, but also those that can help recover the dynamics across all locations in the domain. Using these techniques, it is possible to determine the least number of good candidate observation locations that would be required for a given accuracy of prediction across the footprint.

Soil moisture modeling

A major issue with quantifying soil moisture in the vadose zone is that the scale of measurement of soil moisture does not match that of the process/model (Jana and Mohanty 2012b). A further shortcoming of remotely sensed soil moisture data is that they are representative of only the near-surface layers (0-5cm). Due to the nature of the microwave remote sensing, deeper penetration of signals below this depth is not feasible. While the remote sensing footprint can provide excellent coverage (extent) and possibly intermediate level resolution in some cases (support and spacing), it is limited to only the skin (surface) and no vertical profile information is available. However, agricultural or contaminant transport studies may require and unfortunately rely on deeper soil moisture profiles. Consequently, extrapolation of remotely sensed soil moisture by modeling or scaling is necessary to obtain deeper layer moisture profiles. While we describe scaling in more detail below, assimilation of vadose zone model outputs with remotely sensed data has been performed in several studies. For example, González-Zamora et al. (2016) used a soil water index with a single parameter that was characteristic of the water travel time. This index served as a proxy for root zone soil moisture. Similarly, Pollacco and Mohanty (2012) used a joint parameter inversion scheme using remotely sensed soil moisture and evapotranspiration data to obtain root zone soil moisture profiles.

Scaling of soil moisture

The wide range of support and extent of available soil moisture data makes it necessary to bridge the gap between scales at which soil moisture data are available and the scale at which process representation is required. Two widely used approaches to bridge this disconnect between measurement and application scales are pedo-transfer functions (PTFs) and scaling. PTFs involve obtaining the required soil hydraulic parameters from available or easily measurable soil properties such as the texture and structure (Pachepsky and Rawls 2010). Scaling involves conversion of available measured data at one resolution to another resolution by aggregation (upscaling) or disaggregation (downscaling).

Pedo-transfer functions.

The pedo-transfer function (PTF) approach is used to estimate the soil water retention relationship based on several types of basic soil properties, such as particle size distribution, i.e., sand, silt, and clay content, as well as organic carbon) (Tietje and Tapkenhinrichs 1993). The point regression PTF approach is used to predict the soil water content values at specific soil water potentials using empirically derived regression equations (e.g., Rawls et al. 1982; Ahuja et al. 1985). Function parameter PTFs such as those developed by Vereecken et al. (1989) predict the parameters for the entire water retention and/or hydraulic conductivity functions. Both approaches have been widely used for various soil databases. However, it has been demonstrated that point regression PTFs generally perform better (Tomasella et al. 2003). As the relationship between basic soil properties and soil water retention parameters is rather complicated, the variability in the retention parameters is often driven by different subsets of soil physical properties at distinct ranges of soil water pressure (Jana et al. 2007). Some attempts have been made in the past to apply PTFs across spatial resolutions in order to obtain scaled soil hydraulic parameters (e.g., Jana et al. 2007, 2008, 2012).

Scaling algorithms.

The essential aim of scaling is to describe vadose zone complexity and heterogeneous features, and more importantly, the integrated effect of this complexity on flow or transport processes (Western et al. 2002). Scaling algorithms as applied in the vadose zone literature can be broadly classified into two categories. The first, and most common, approach is to derive parameters of the constituent flow equations (e.g., the Richards' equation) that effectively produce the same fluxes at the scale of interest as the scale of their measurement. This approach is, logically, termed as parametric scaling. The other approach, termed process-based scaling, attempts to modify the forms of the equations describing the soil moisture dynamics based on key governing variables at each scale of interest. The reasoning behind this approach is that different heterogeneous features and vadose zone complexities will become relevant at different scales, so the crucial question is to determine the key governing variables (Vereecken et al. 2016). Because there is significant knowledge gap about the nature of the relationships between these governing variables and the soil moisture dynamics, process scaling is less favored in practice.

Upscaling is the process of supplanting a heterogeneous domain of soil properties with an effective, homogeneous one. Most upscaling efforts for soil hydraulic parameters are based on the assumption that the movement of water occurs only in the vertical plane, and that no lateral movements take place (Khaleel et al. 2002; Vereecken et al. 2007). However, as would be expected, this assumption is valid only at relatively finer scales of around a few hundred meters (Vereecken et al. 2014). At larger footprints, it is no longer safe to assume that complex soil features such as topography, land cover, and meteorology do not exert significant influence on soil moisture dynamics. Hence, it becomes necessary to incorporate these complex features into the scaling methodology. One of the early steps in this direction was taken by Jana and Mohanty (2012a) who incorporated the influence of topography in an upscaling algorithm. A power averaging approach was used where the upscaled parameter, P*, was defined as




and the scale parameter, η, is given by


In the above system of equations, Sup(pi; pj) is the support from parameter pi to parameter pj, i and j are indices for locations where the parameters are measured, x, y, and z are Cartesian coordinates of the measurement location, and S is the scale to which the parameters are being aggregated. As Eq. (4) suggests, the scale parameter is a product of the normalized difference in elevation between the two locations i and j, and the linear distance between measurement values, normalized by the scale dimension. When upscaling soil hydraulic parameters to resolutions of a kilometer or beyond, topography-based scaling outperformed other algorithms which ignored the effect of such complex features (Jana and Mohanty 2012a).

Downscaling is the process of disaggregating large homogeneous footprints of soil moisture into a heterogeneous field of finer resolution values. Several downscaling approaches have been applied in the past including data fusion (e.g., Das et al. 2011; Montzka et al. 2016), statistical modeling (e.g., Ines et al. 2013; Verhoest et al. 2015) and data assimilation (e.g., Sahoo et al. 2013). A comprehensive review of a number of spatial downscaling algorithms is provided by Peng et al. (2017).


A recent survey of the hydrological sciences community showed the need for better integration of field work and modeling across a range of spatial and temporal scales in order to advance understanding in this area (Blume et al. 2016). The study pointed out that one of the greatest challenges to this advancement is the lack of improved hydrological theories. In research related to soil moisture, it has been argued that concepts such as the Richards' equation, valid at the continuum scale, are routinely applied across much larger extents wherein assumptions made at the continuum scale are no longer valid (e.g., Blöschl and Sivapalan 1995). More than two decades after this observation, little has changed.

We suggest that future efforts be focused on developing comprehensive scaling algorithms that are capable of refactoring the role of vadose zone complexities, especially under varied saturation conditions and transient phenomenon (Fig. 4). Further, a holistic approach to characterize and quantify soil moisture, beyond efforts from the hydrological community, will be useful. For example, statistical data mining approaches can prescribe important approaches like identification of time stable locations.

To gain predictive understanding of heterogeneous and temporally variable vadose zone processes, it is critical to integrate multi-type multi-scale data encompassing both surface and the subsurface. In particular, spatially extensive characterization technologies do not often directly measure the properties of interest; such as permeability, or they have a large uncertainty (such as petrophysical relationships for geophysics). It is cost prohibitive to measure all the properties across many locations. Similarly, temporal vadose zone patterns can be difficult to interpret because their variability can depend on several factors, including temperature, precipitation, vegetation and stream discharge. To tackle this challenge, statistical approaches combined with geophysical and geochemical data have been developed recently to capture spatial and temporal variability in field scale environments (e.g., Arora et al. 2013, 2019; Sassen et al. 2012; Wainwright et al. 2014). This section describes these zonation-based and time variable approaches and its application at two field studies—the Rifle floodplain and the Savannah River site. A brief description of these sites is also included in the following section.

Identifying critical interfaces through zonation

The zonation approach is based on the suppositions that: (1) many subsurface and surface hydrogeomorphlogical, vegetation, and physiochemical properties are likely to be correlated with each other; (2) some properties may exert a dominant control on the vadose zone processes (such as clay content, soil moisture, microtopography); and (3) some spatially extensive measurements (such as geophysical and remote sensing) may provide reasonable proxies for controlling variables, thus offering an approach to characterize key properties over large regions and in a minimally invasive manner. This is the extension of lithofacies and hydrofacies concepts (e.g., Fogg et al. 1998; Klingbei et al. 1999; Weissmann et al. 2002; Heinz et al. 2003; Yabusaki et al. 2011), in which geological units represents a set of common hydrological or geochemical properties. The zonation approach made these concepts more general and quantitative such that the zones or units are discovered and defined by a suite of spatially extensive datasets such as geophysics and remote sensing as well as point measurements. The zonation approach works by using a statistical data-mining approach to integrate spatially extensive datasets for identifying zones that have distinct characteristics relative to neighboring zones. For example, Sassen et al. (2012) and Wainwright et al. (2014) used geophysical data and zonation approach to map hydrogeological and geochemical properties (i.e., permeability, percent of fines, mineral ratio) in the subsurface at a uranium contaminated site. Sassen et al. (2012) showed that the correlation between hydraulic and geochemical properties would need to be taken into account to describe the flow processes. Hubbard et al. (2013) and Wainwright et al. (2015) applied the zonation-based approach to characterize the spatial organization of the Arctic tundra ecosystem and permafrost environment. They showed that the zonation approach could capture the heterogeneity of near-surface soil characteristics and greenhouse gas fluxes, and that it enables mapping of soil properties and fluxes across the site based on sparse point-scale measurements. In this manner, the zonation approach provides a tractable and effective way to capture the spatial heterogeneity of co-varied properties and employs varied datasets across scales. We further describe the zonation approach through two case studies, wherein the approach was successful in describing and quantifying contaminant fate and transport.

Identifying critical time periods through wavelet and entropy approaches

It is important to address time invariance in vadose zone dynamics because biogeochemical interactions are inherently non-linear and complex, so much so that changes to measured water chemistry parameters (such as pH, SO4) and gas fluxes (e.g., CO2) can indicate the influence of multiple processes simultaneously. This is problematic when developing conceptual models and identifying governing controls on vadose zone processes that need to be incorporated into numerical models. Wavelet and entropy-based approaches offer the opportunity to systematically interrogate complex, multivariate datasets to characterize natural variability and identify the governing processes that exert control over vadose zone patterns at different time scales. For example, Arora et al. (2013) used wavelets to associate variations in geochemical concentrations at a landfill site to water table and precipitation dynamics, which showed dominance at 8 month scales. They further documented that temporal variability in sulfate concentrations was significantly impacted by FeS cycling and uptake by vegetation, which would need to be incorporated into models to accurately represent redox dynamics at the site. Similarly, Arora et al. (2019) applied the entropy-based approach to characterize temporal variability of the Arctic tundra ecosystem and permafrost environment. They showed that the entropy approach could capture the variable nature of relationships between gas fluxes and near surface soil characteristics (e.g., soil moisture, temperature), and thereby document the impact of changing environmental conditions (e.g., warming) on these patterns in a future climate. In this manner, wavelet and entropy-based techniques provide a simplistic yet tractable approach to identify temporal patterns and capture the influence of co-variability of hydrological and biogeochemical processes, as well as climatic controls on vadose zone dynamics and conditions.

Case studies

We describe the zonation approach through two case studies, wherein spatially heterogeneous clusters were successful in describing and quantifying contaminant fate and transport. We further describe a wavelet-entropy approach, whose characterization of temporal patterns was successfully embedded into a reactive transport model of the Rifle site.

Rifle case study.

The Rifle site is located in the town of Rifle, CO, adjacent to the Colorado River. The Department of Energy (DOE) Rifle site is a former uranium mill tailings site, where soil and groundwater have been contaminated with uranium above the background level. Extensive characterization of the site has occurred over the past ten years to identify and quantify biogeochemical cycling in the Rifle floodplain. In particular, studies have found regions near riverbanks—known as naturally reduced zones (NRZs)—that have relatively high organic carbon, reducing conditions and elevated uranium concentration (Campbell et al. 2012). Historically, NRZs have been considered important as diffusion-limited interfaces rich in uranium and pyrite. A zonation-based approach was used for the probabilistic mapping of NRZs within a three-dimensional subsurface domain using a geophysical method of induced polarization imaging, since it identifies the electrical chargeability of subsurface materials and hence is sensitive to certain minerals (i.e., reduced sulfides) associated with NRZs (Wainwright et al. 2016).

Furthermore, a novel wavelet-entropy approach was applied to identify temporal variability in biogeochemical processes (Arora et al. 2016a). In particular, this approach works by interrogating complex, multivariate geochemical datasets to identify transient phenomenon within the floodplain environment. Results indicated that seasonal and annual hydrologic perturbations associated with rainfall and snowmelt events drove biogeochemical cycling at the site.

A 2-D reactive transport model was developed to investigate the impact of these NRZs and transient phenomenon—i.e., water table variations and temperature gradients—on subsurface carbon fluxes in the Rifle floodplain (Arora et al. 2016b). Results from model simulations ignoring NRZs underpredicted atmospheric CO2 fluxes by almost 230% compared to simulations with NRZ relevant microbial reaction pathways. Results further indicated that subsurface carbon exports from the Rifle site were underestimated by almost 170% (to 3.3 g m−2 d−1) when transient conditions were ignored in the simulations. This modeling study concluded that spatially (e.g., naturally reduced zones) and temporally discrete conditions (e.g., temperature fluctuations) can significantly contribute to carbon cycling in floodplains and need to be appropriately represented in model simulations. Other modeling efforts at the site have therefore explicitly used 3D map of NRZs to quantify oxygen and nitrogen fluxes (Fig. 5) (e.g., Yabusaki et al. 2017; Dwivedi et al. 2018a).

Savannah case study.

The DOE's Savannah River Site F-Area includes three unlined seepage basins, where low-level radioactive waste solutions were disposed (Bea et al. 2013; Wainwright et al. 2014). The site located near Aiken, South Carolina is a former nuclear weapons production site. Groundwater at the site has been contaminated with uranium and other radionuclides associated with plutonium separation for over 30 years. The site has been extensively characterized in terms of subsurface structure (Sassen et al. 2012; Wainwright et al. 2014) and uranium geochemistry (Dong et al. 2012). These studies found that the hydrogeological and geochemical parameters (such as permeability, percent of fines, mineral ratio) are tightly correlated to each other, and distinct geological units—reactive facies—exist, each of which has a unique distribution of co-varied properties. Thus, reactive facies based zonation approach has been found to provide a relevant description of physiochemical controls at plume relevant scales. In particular, a seismic technique was used to identify the spatial distribution of reactive facies along the uranium plume center line (Wainwright et al. 2014).

A modeling study found that incorporating reactive facies based properties and parameters improved the accuracy of model predictions, which are important to understand long-term uranium mobility and plume evolution at the site (Sassen et al. 2012). Moreover, a reactive transport model including the saturated and unsaturated (vadose) zones, U(VI) and H+ adsorption (surface complexation) onto sediments, dissolution and precipitation of Al and Fe minerals, and key hydrodynamic processes was developed for the site (Bea et al. 2013). Herein, uncertainty quantification analyses were conducted to assess the sensitivity of model results to various input parameters affecting the migration of the acidic-U(VI) plume at the F-Area. This study concluded that model (and parameter) sensitivity evolves in space and time, and recommended that uncertainty-based analyses be conducted to assess the temporal efficiency of remediation strategies in contaminated sites. A recent study further compared two surface complexation models in describing U(VI) plume evolution at two different time periods of the site operation—the basin infiltration and post-closure times (Arora et al. 2018). The study concluded that there was uncertainty in describing the pre-closure behavior of U(VI) from the F-Area because of the different model approaches. This suggests that a careful evaluation and buildup of the reactive transport models is necessary for use in evaluating remediation strategies, nuclear cleanup, and other opportunities.

Advancement in computational methods and measurement techniques have expanded the role of RTMs to address a wide range of spatial and temporal scales—from pore to catchment scale and across 100 to 103 years (Vereecken et al. 2016; Li et al. 2017). With this advancement, however, there is a growing need to evaluate the robustness and performance of RTMs and describe model limitations. In this regard, we present the hindrances to model application, formal treatment of uncertainty, the role of benchmarking activities, and future opportunities therein.

Model validation

Geochemical model validation is important if model results are being used to provide accurate assessments of engineering performance, risk analysis, and support policy formulations, which have far-reaching consequences. However, in contrast to hydrologic models, the extent to which geochemical models can be validated has been questioned (Tsang 1987). In this regard, model validation has been primarily related to careful evaluation of the conceptual framework and support for better model calibration (e.g., through the incorporation of spatial and temporal heterogeneities) (Arora et al. 2018).

Although geochemical models have been successful in describing processes within the heterogeneous vadose zone, as shown through case studies above, several aspects need to be determined before developing such a model. In particular, a careful evaluation of the input data for the species of interest, determination of the accuracy required of the modeling results, and whether the conditions of the system being modeled match the range of conditions for which the model is valid is needed (Herbert 1996). For example, Um et al. (2008) concluded that a nonelectrostatic surface complexation model developed using variable pH measurements on synthetic ferrihydrite resulted in inaccurate predictions of U(VI) mobility in subsurface environments. This and other studies demonstrate the importance of obtaining experimental data over a sufficiently wide range of chemical conditions when developing a surface complexation model (Fox et al. 2006, 2012; Hyun et al. 2009). Similarly, a critique of the thermodynamic database would be necessary to evaluate to what extent input parameters affect model predictions.

Inverse estimation.

A reactive transport model generally consists of a number of flow, transport and kinetic parameters, whose values are not necessarily known a priori. Many of these parameters are determined through experimental or laboratory studies, and yet others need to be estimated through inverse modeling. Further, the choice of data types for the inverse analyses has a significant impact on the resulting model predictions and model performance. For example, Arora et al. (2012) concluded that lumped macropore-matrix observations of water content, cumulative outflow, and effluent concentration increased errors in inversely estimating model parameters representing preferential flow and bromide transport at the column scale in comparison to using separate measurements of macorpore and soil matrix.

Uncertainty quantification

Uncertainty in vadose zone models can arise from a number of different sources including inappropriate conceptual framework and assumptions, uncertainty in the input and thermodynamic parameters, and the inherent stochastic nature of the subsurface system. Conceptual model uncertainties for vadose zone models have been quantified using Bayesian model averaging (e.g., Wöhling et al. 2008; Gupta et al. 2012). This Bayesian approach evaluates conceptual uncertainty by using multimodel ensemble simulations.

Uncertainties in input model parameters can arise from measurement errors and inverse estimation. Inverse parameter estimation is challenging because of nonidentifiability of the solution set, ill-posedness of the inverse problem, and non-unique parameter set (Carrera and Neuman 1986). Moreover, heterogeneous vadose zone models that incorporate chemical and mineral heterogeneities or preferential flow typically suffer from increased interdependence of model parameters, which increases the risk of reaching local minima in the parameter set (e.g., Arora et al. 2012). Thus, uncertainty estimation seems important from a geochemical perspective. Bayesian methods present a powerful solution to the challenges of model nonlinearity and identifying the maximum likelihood of parameters.

Development of benchmarks

Despite the availability of a large number of models available to simulate vadose zone flow and reactive transport, substantial differences exist among them—some that are related to the conceptual framework, process couplings, and yet others that ascribe numerical differences. Considering these differences, benchmarking activities are critical to building confidence in models and providing measures of model performance. Several benchmarking studies have been conducted in the past that compared thermodynamic database and speciation (e.g., Nordstrom and Archer 2003), and others that compared conceptual and numerical capabilities of models (e.g., Mayer et al. 2015). These are described in more detail elsewhere (Arora et al. 2015; Steefel et al. 2015; Dwivedi et al. 2016).

Future opportunities

There is a significant need for the further development of monitoring and predictive tools for allocation and management of water resources, decision making for environmental protection, and increasing cost saving for remediation efforts. New tools and approaches are needed that are open source and easy to use, with capabilities ranging from inverse parameter estimation and uncertainty quantification, to launching and monitoring simulations. Herein, we describe some opportunities for future growth.

Data-worth analysis.

Addressing remediation issues in a data-driven, cost-effective manner is of increasing significance, especially because of their economic and environmental implications. It is therefore essential to determine the influential vadose zone properties with sufficient accuracy so that relevant predictions can be made with confidence. Data-worth analysis is one such metric that measures the extent to which the collected data can reduce the uncertainty of target predictions that are critical for decision-making. The data-worth analysis works by ranking the contribution that each (potential or existing) data point makes to the solution of an inverse problem (used for model calibration and parameter estimation) and a subsequent predictive simulation (Finsterle 2015).


One of the biggest challenges for vadose zone modeling is to apply laboratory experimental results to natural field systems. Geoscientists have often delved into the rate conundrum particularly dealing with significant discrepancies between rate constants derived from experimental studies to application at the field scale (e.g., Maher et al. 2006). While several scaling algorithms have focused on soil moisture estimation in the vadose zone, relatively fewer approaches have considered scaling that have biogeochemical implications.

Comprehensive vs parsimonious approach.

Some authors have argued that rather than developing a geochemical model with a large number of parameters that are poorly defined at a larger extent, a parsimonious approach may work better (e.g., Basu et al. 2010). The latter approach can still describe contaminant distribution as a function of landscape characteristics while being more numerically efficient than a reactive transport model. However, the former approach has several advantages in teasing out individual components and processes especially to describe the stochastic and non-linear nature of the vadose zone (e.g., Fatichi et al. 2016). In particular, reactive transport models may be better poised to forecast the spatiotemporal evolution of plume or non-linear feedbacks with climate than empirically based models.

A recent review paper (Vereecken et al. 2016) suggests that a major challenge in developing a numerical model for the vadose zone is the coherent integration of all the information on (i) the multiscale architecture (including the heterogeneous vadose zone architecture); (ii) the process formulation for the chosen range of scales; (iii) the system's coupling to the environment, which is typically represented as an external forcing but should also include the feedbacks to the atmosphere and/or groundwater; and (iv) the available data, which often need to be transferred into the chosen range of scales. Thus, future opportunities in the development of vadose zone models should focus on holistic approaches that honor the new data types and streams, but also include all aspects of the critical zone—from the bedrock to the top of the canopy.

Vadose zone hydrologic and biogeochemical processes play a significant role in the capture, storage and distribution of contaminants between the land surface and groundwater. One major issue facing geoscientists in dealing with investigations of the unsaturated zone flow and transport processes is the evaluation of heterogeneity of subsurface media. This chapter presents a summary of approaches for monitoring and modeling of vadose zone dynamics in the presence of heterogeneities and complex features, as well as incorporating transient conditions. Modeling results can then be used to provide early warning of soil and groundwater contamination before problems arise, provide scientific and regulatory credibility to environmental management decision-making process to enhance protection of human health and the environment.

We recommend that future studies target the use of RTMs to identify and quantify critical interfaces that control large-scale biogeochemical reaction rates and ecosystem functioning. Improvements also need to be made in devising scaling approaches to reduce the disconnect between measured data and the scale at which processes occur.

This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Biological and Environmental Research (as part of the Watershed Function Scientific Focus Area), and Office of Science, Office of Advanced Scientific Computing (as part of the project “Deduce: Distributed Dynamic Data Analytics Infrastructure for Collaborative Environments”) under Contract No. DE-AC02-05CH11231.

Open access: Article available to all readers online.