Unsaturated zone travel time, Central Valley, California

An important aspect of diffuse (nonpoint) pollution transport in the subsurface is the movement of contaminants through the unsaturated zone.

In general, the simulation of flow and contaminant transport in the unsaturated zone is a highly nonlinear and complex problem due to the often highly heterogeneous structure of the unsaturated zone. At the groundwater basin scale a highly detailed, fully three-dimensional simulation of unsaturated zone flow using tools like HYDRUS is computationally not feasible.  Therefore, simpler approaches are typically sought. Among the most efficient physical estimates of unsaturated zone flow is the vertical piston flow approach. It assumes that water moves vertically through the unsaturated zone at an annualized average constant flow rate, R.  The travel time t that water takes to travel from the land surface to the water table can then be calculated as

t = θ Dgw / R

where θ[-] is the mobile water content, Dgw[L] is the depth to groundwater and R is the average long-term vertical flow rate [L/T] in the unsaturated zone.

Unsaturated Zone Flow Rate

Under the piston flow assumption, water moves vertically through the unsaturated (vadose) zone in a uniform, plug-like manner. Within this framework, the average annual vertical flow rate is assumed to be equal to the groundwater recharge rate. This provides a practical and physically consistent way to estimate travel times from the land surface to the water table.

Groundwater recharge—along with its spatial distribution across the Central Valley and its temporal variability—has been extensively evaluated using large-scale regional groundwater models, including C2VSim or CVHM. These models integrate hydrologic processes such as precipitation, irrigation return flow, and soil properties to produce spatially distributed estimates of deep percolation and recharge.

In CV-NPSAT framework, two independent recharge estimates are used to represent uncertainty and variability in long-term conditions:

• CVHM-based estimate:
  The average deep percolation rates simulated by CVHM over the period 2009–2019.
• C2VSim-based estimate:
  The average deep percolation rates simulated by C2VSim over the period 2000–2015.

By incorporating both datasets, CV-NPSAT captures a range of plausible recharge conditions, allowing for more robust evaluation of unsaturated zone travel times and their impact on contaminant transport.

Depth to groundwater

Accurate estimation of depth to groundwater is a critical component in modeling subsurface transport processes. In the Central Valley, this information is available from two primary sources, each with distinct advantages and limitations.

• Measured Data: Depth to groundwater is directly measured at a set of wells distributed across the Central Valley. These measurements provide high accuracy at specific locations, reflecting true local conditions. However, the spatial distribution of wells is irregular, and the data are not easily interpolated to produce continuous regional maps.
• Simulated Data: Regional groundwater models such as C2VSim and CVHM provide spatially continuous estimates of water table depth over large areas and long time periods, capturing key processes such as pumping, surface water interactions, and geologic variability. However, these simulated results are subject to uncertainty and may not fully match observed conditions.

Here we employ a hybrid approach using both data sources to take advantage of the precision of measured data on one hand, and the continuity of simulated data on the other hand.  The approach preserves the undulating nature of the water table surface, estimated from the simulation model, as much as possible, while forcing the water table surface to coincide with the measured depth to the water table at the exact locations of the measured wells.

Available data

Measured groundwater depth data were obtained from the Department of Water Resources, while simulated water table depths are available from the C2VSim and CVHM models for the Central Valley. For consistency, the simulated depth estimates use the same time periods as those adopted for groundwater recharge.
we define as

• Measured data at well locations: Xmeas, Ymeas, Dmeas (x,y coordinates of the measured locations and measured water table depth)
• Simulated data at the finite difference cell or finite element of the simulation grid: Xsim, Ysim, Dsim (x,y coordinates of the simulated water table depth locations, usually one to several points per square mile throughout the Central Valley)

In the following we denote as F(x, y, X, Y, D) ageneral interpolation operator (e.g., kriging) that uses the data at X,Y locations to interpolate D at the x,y locations.

Methodology

The steps to condition the simulated data with the measured data are as follows:

1. Simulated Values at Well Locations (Dsim_meas): Extract the simulated depth to the water table by interpolating simulated water table depth to each location at which measurement data are available (well locations):
   Dsim_meas (Xmeas, Ymeas) = F(Xmeas, Ymeas, Xsim, Ysim, Dsim)
   Given the high resolution of the simulation grid, interpolation error at these points is assumed negligible.
2. Interpolation of Simulated Well Values (Dsim_interp): Using only the simulated data, Dsim_meas, at the measured well locations, we perform an interpolation of those simulated well measurements, which yields a continuous, interpolated map of an estimated depth to the  (simulated) water table across the entire Central Valley. This interpolated map of the simulated water table depth is different from the simulated map of water table depth, except at the simulated well measurement locations, which are also the locations of real well measurements.  The interpolated water table depth across the finite difference cells/finite elements of the simulation grid, based on simulated well measurements is
   Dsim_interp(Xsim, Ysim) = F(Xsim, Ysim, Xmeas, Ymeas, Dsim_meas)
3. Estimation of Interpolation Error (Derror): We now consider the difference between the original simulated map of water table depth, Dsim, and the estimated map of simulated water table depth based on simulated well measurement data, Dsim_interp.  This difference represents the exact error introduced by the interpolation method (e.g., kriging) relative to the physically simulated water table depth at locations away from the well locations with (simulated) measurements. The error at any one location will be a function of the density and spatial distribution of the measured data. The estimation error is: 
   Derror(Xsim, Ysim) = Dsim - Dsim_interp
   Note that Derror = 0 at the measured well locations and, hence, simulated well measurement locations.
4. Interpolation of Measured Data (Dmeas_interp): We now repeat the exact same interpolation of step 2, but using the actual measured data instead to obtain an estimated depth of water table at the simulation grid cell locations:
   Dmeas_interp(Xsim, Ysim) = F(Xsim, Ysim, Xmeas, Ymeas, Dmeas)
5. Final Conditioned Groundwater Depth (Dgw): Finally, we obtain a conditionally interpolated depth to groundwater, Dgw, by adding the estimation error computed in step 3 to the interpolated depth to groundwater obtained from measured data:
   Dgw(Xsim, Ysim) =  Dmeas_interp + Derror
   This step preserves, as closely as possible, the undulations in the water table captured by the simulation, here reflected in the magnitude of Derror, while preserving the measured depth to groundwater at the well locations. Derror = 0 at the locations of these wells due to having selected those locations when we obtained Dsim_meas.

CV-UNSAT Earth Engine Application

To support interactive exploration of the unsaturated flow model used in CV-NPSAT online tool, we provide an Earth Engine application that allows users to visually compare, side by side, the key datasets used in estimating unsaturated zone travel times.