3.2.1 Geophysical data acquisition

To characterize subsurface heterogeneity, the study site was surveyed with a towed transient electromagnetic (tTEM) system in fall 2017 (Behroozmand et al., 2019). The tTEM system can rapidly map meter-scale variations in subsurface resistivity over large areas. Given the link between resistivity and sediment texture (Knight and Endres, 2005), tTEM is a natural choice for recharge site assessment. In total, 46 line kilometers of tTEM data were acquired in the almond orchard and an adjacent field (Fig. 2). Within the orchard, tTEM data were collected in rows ∼6.8 m apart and, after data processing, resulted in a sounding every ∼10 m along each row. In the adjacent field, tTEM data were collected in lines 25–40 m apart with a sounding every 5 m along each line. These data were inverted in Aarhus Workbench (Aarhus GeoSoftware) using a smooth, 1D spatially constrained approach (Auken et al., 2015). For each sounding, a global sensitivity threshold (Christiansen and Auken, 2012) was used to determine the depth of investigation (30–80 m), below which reliable estimates of resistivity cannot be recovered. The resulting resistivity model (Fig. 3a) contains 24 layers. Layer thicknesses increase from 1 m at the surface to 9.7 m at the maximum depth of investigation (81 m). The resistivity model was used to estimate the depth to the water table at the site (45 m; Goebel and Knight, 2021), which is in agreement with nearby water table measurements (44 m).

Figure 2An aerial image of the study site and a map of California showing the location of the site within the Central Valley (light gray area). Individual tTEM soundings are shown as blue circles, while the location of each cone penetration test (CPT) is represented as a yellow square. The extent of the hydrologic model domain (dashed yellow line) is 200 m larger in every direction than the extent of the almond orchard (dashed orange line). Satellite imagery © Google Earth 2023.

3.2.2 Interpolation and extrapolation of geophysical data

Following geophysical inversion, the resistivity model consists of irregularly spaced point data. To transform these data onto a rectilinear grid for use in a hydrologic model, we use sequential Gaussian simulation (SGSIM). SGSIM is a common technique for populating a grid with a Gaussian random variable without the smoothing effect of deterministic methods (e.g., ordinary kriging). Because the hydrologic model domain extends beyond the bounds of the resistivity model, this procedure includes both interpolation (between any soundings > 10 m apart) and extrapolation (vertically from the depth of investigation to 90 m and laterally up to 200 m away from the orchard boundaries). SGSIM accounts for the uncertainty in this process by generating multiple equally probable realizations of the variable of interest. The full workflow is shown in the Supplement (Fig. S1) and described briefly here.

First, resistivity values were transformed to a Gaussian distribution using a normal score transformation. An implicit assumption of SGSIM is stationarity, meaning that statistical properties of the random variable (e.g., the mean and variance) do not vary spatially throughout the domain. However, experimental variograms of the normal score values displayed a trend with depth. To account for this, we decomposed each normal score value (z(u), where u is a location in the domain) into trend (m(u)) and residual (r(u)) components (Wackernagel, 2003):

$$\begin{array}{}\text{(1)}& z\left(\mathit{u}\right)=m\left(\mathit{u}\right)+r\left(\mathit{u}\right).\end{array}$$

We modeled the trend using a radial basis function with a Gaussian kernel parameterized to capture the smooth, large-scale changes in normal score values. To extrapolate the trend laterally, for each geostatistical realization, 15 random 1D sounding profiles were copied from the tTEM survey to random locations outside the survey footprint. These 15 profiles are enough to ensure the extrapolated m(u) values exhibit a similar trend beyond the orchard as within it. This is a logical assumption given that the characteristic trend within the study area is a decrease in resistivity with depth, which is consistent with regional surveys of resistivity (Knight et al., 2018) and sediment type (Mid-Kaweah GSA, 2019). Next, we transformed the residuals to a Gaussian distribution and generated 600 realizations of z(u) on a $\mathrm{1200}\times \mathrm{800}\times \mathrm{90}$ m (width × length × depth) uniform grid using SGSIM. This number of realizations was chosen so as to satisfy the approximate number of simulations required for global sensitivity analyses, described in Sect. 3.4.2. For horizontal cell sizes, we used the approximate spacing between tTEM soundings (10 m). For the depth direction, we used the finest vertical cell size in the tTEM-derived resistivity model (1 m). Following SGSIM, each realization of z(u) was added to a unique realization of the extrapolated trend, m(u). Finally, the combined trend and residual values were back-transformed to resistivity units and downsampled to the vertical resolution of the tTEM-derived resistivity model. An example realization is displayed in Fig. 3b.

Figure 3One example realization at each step of the workflow for parameterizing a hydrologic model with tTEM data. Inversion of the tTEM soundings produces a resistivity model of irregularly spaced points (a). The geostatistical workflow (Sect. 3.2.2) then interpolates and extrapolates the resistivity model onto a rectilinear grid (b). A rock physics transform is then applied to create realizations of coarse fraction (c). Finally, after clustering each coarse-fraction realization into distinct facies, each facies is assigned specific parameter values, such as horizontal saturated hydraulic conductivity (d). All plots are shown with 10× vertical exaggeration.

Download

3.2.3 Transforming resistivity to coarse fraction

To determine sediment properties from resistivity data, 1D tTEM resistivity profiles were compared with sediment types determined from collocated cone penetration test (CPT) logs. This procedure is described in detail in Goebel and Knight (2021) and briefly summarized below. At this site, CPT logs were acquired at five locations extending from the surface to depths of 19–32 m, with a sampling interval of 10 cm (Fig. 2). CPT logs were then used to determine the soil behavior type and horizontal hydraulic conductivity of each sampling interval (Robertson et al., 1992; Lunne et al., 1997; Robertson, 2009, 2016). Profiles of soil behavior type at these five locations were then compared to profiles of electrical resistivity recorded during CPT logging (Fig. S2a), and k-medoids clustering was used to classify soil behavior types with similar resistivity values (Goebel and Knight, 2021). This classification scheme resulted in a coarse-dominated cluster corresponding to high resistivity values (“clean sand to silty sand” and “gravelly sand to sand”) and a fine-dominated cluster corresponding to low resistivity values (all other soil behavior types). Clustering results were then used to convert soil behavior type logs to binary coarse–fine sediment logs (Fig. S2b).

The CPT data were then aggregated to the same resolution as the tTEM-derived resistivity model by calculating the fraction of coarse-classified sediment layers over a given depth interval, a statistic referred to as coarse fraction. Aggregating the CPT data resulted in 62 depth intervals with paired coarse-fraction and resistivity information. Though the logs do not extend to the full depth of investigation of the tTEM survey, the range of resistivity values in the overlapping profiles covered 97 % of all resistivity values in the tTEM resistivity model (Fig. 3a). The eddy currents induced in the subsurface during a TEM survey are primarily oriented horizontally. If we assume that sediment layers within a single tTEM cell are also oriented horizontally, we can model the electrical current flow through each cell as current flow through a set of resistors in parallel. Algebraically, the resistivity of a tTEM cell (ρ_tTEM) is described by the harmonic mean of the resistivity of the fine- and coarse-grained layers:

$$\begin{array}{}\text{(2)}& {\mathit{\rho }}_{\mathrm{tTEM}}={\left(C_{\mathrm{F}}{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{coarse}}}}+\left(\mathrm{1}-C_{\mathrm{F}}\right){\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{fine}}}}\right)}^{-\mathrm{1}}\end{array}$$

where C_F is the fraction of coarse layers within a cell, ρ_coarse is the resistivity of the coarse-grained layers, and ρ_fine is the resistivity of the fine-grained layers. Using the procedure of Knight et al. (2018), we then generate 1000 bootstrapped samples of 62 paired ρ_tTEM and C_F values and use non-negative least squares to solve for ρ_coarse and ρ_fine for each sample. This technique produces a distribution of the expected resistivity values for coarse-dominated and fine-dominated end-members (Fig. S3). The range of resistivity values for a given end-member represents both uncertainty in the field data and the variability in resistivity that exists for a given sediment type due to other factors (water content, salinity, etc.).

Typically, sediments exhibit different resistivity–lithology relationships above and below the water table (Knight et al., 2018; Knight and Endres, 2005). The procedure above was used to build resistivity–sediment-type distributions for the vadose zone but could not be used for the saturated zone as none of the CPT logs extend below the water table. Instead, we examined comparable resistivity–sediment-type distributions from airborne electromagnetic data over the site and nearby well logs (Knight et al., 2018). We then calculated the offset, in Ω m, between the vadose-zone and saturated-zone distributions for a given sediment type. We calculated separate scaling factors for coarse-dominated (10^−0.118 Ω m) and fine-dominated (10^−0.110 Ω m) sediments; the scaling factor for coarse sediments is expected to be of larger magnitude because the sediments typically exhibit a larger difference in saturation above and below the water table than do fine-grained sediments. We then applied this scaling factor to the vadose-zone CPT-derived resistivity–sediment-type distribution to generate equivalent distributions for the saturated zone (Fig. S3a).

Next, we used the vadose-zone and saturated-zone resistivity–sediment-type distributions to transform resistivity to coarse-fraction values between 0 and 1. Using increments of 0.001, for each coarse-fraction value in this range we calculated ρ_tTEM 10 000 times using Eq. (2), with ρ_coarse and ρ_fine randomly sampled from the distributions in Fig. S3a. Paired coarse-fraction and ρ_tTEM values (Fig. S3b) were then aggregated into bins by resistivity and used to calculate the distribution of coarse-fraction values within each bin. Finally, we used the binned distributions to transform each geostatistical realization of resistivity (Sect. 3.2.2) to coarse fraction. For each geostatistical realization, we iterated through all resistivity bins and sampled from each bin's coarse-fraction distribution. To vary the transform between each realization, we chose random quantiles at which to evaluate the coarse-fraction distributions for the lowest and highest resistivity values observed in that realization (lowres_q and highres_q, respectively). We then used linear interpolation between lowres_q and highres_q to determine quantiles at which to evaluate coarse-fraction distributions for all resistivity bins between the two end points. This random sampling is advantageous over other methods (e.g., choosing the maximum likelihood value from each distribution) because it accounts for the uncertainty in the distribution of coarse-fraction values corresponding to a single value of resistivity. The end result is a suite of 600 realizations of coarse fraction, several of which are presented in Figs. 3c and S4.

3.3.1 Boundary conditions and discretization

The hydrologic modeling domain measures 1200 m × 800 m × 150 m in width, length and depth, respectively, with uniform, 10 m horizontal cell sizes. In the upper 90 m, cell thicknesses in the depth direction match the resolution of the tTEM cells (1–9.7 m). Below that depth, we add a single layer of 60 m thick cells to account for the estimated thickness of the unconfined aquifer below the tTEM depth of investigation (Mid-Kaweah GSA, 2019). In total, the domain contains 240 000 grid cells ($\mathrm{120}\times \mathrm{80}\times \mathrm{25}$ cells). The impact of this discretization scheme on model output is discussed further in the Supplement (Fig. S7 and Table S1).

We impose a no-flow boundary on the bottom of the domain, which approximates the depth to the Corcoran Clay, an aquitard. On the sides, we apply a constant-head boundary equal to the estimated depth of the water table (45 m). During model development, we adjusted the lateral boundaries as far as 1000 m from the edge of the almond orchard and found that constant-head boundaries 200 m away (Fig. 2) are sufficiently far to have minimal impact on model results (Table S2). The boundary conditions used at the top of the model vary depending on the model stage and are described below.

3.3.2 Model parameterization

The geophysical and geostatistical workflow above (Sect. 3.2.1–3.2.3) results in 600 unique realizations of one sediment property, coarse fraction, below the site. To estimate hydrologic properties of the subsurface (hydraulic conductivity, porosity, etc.) from this metric, we assume that there are only two end-member sediment types at the site – coarse and fine – and that each cell is a mixture of the two. We then define parameters for each end-member and an averaging procedure to calculate properties of cells that are a mix of the two end-members. This approach has been applied to Central Valley groundwater flow models in the past (e.g., Phillips and Belitz, 1991; Belitz and Phillips, 1995; Faunt, 2009) and is supported by the observation that sediment texture and hydraulic properties are spatially correlated in geologic media (Russo and Bouton, 1992). Nonetheless, this assumption is a simplification of the real groundwater system; its implications are discussed more fully in Sect. 6.

Due to a paucity of site-specific data, we compile a range of potential parameter values and averaging schemes using literature data from the Central Valley (rows 1–12 in Table 1), with a preference towards published research in the vicinity of the site. For most parameters, we approximate prior uncertainty using a uniform or log-uniform distribution, with lower and upper bounds as the minimum and maximum of the range of values reported in the literature. The one exception is horizontal hydraulic conductivity, described below. Once each parameter's prior uncertainty has been defined, we use Latin hypercube sampling to sample from each distribution, jointly varying all parameters between every simulation. Latin hypercube sampling has some statistical advantages over random sampling (Saltelli, 2008) and ensures that all regions of the parameter space are sampled.

Table 1Parameter ranges sampled during Monte Carlo simulations.

^a Principal component scores from principal component analysis of all coarse-fraction realizations are also included in global sensitivity analyses but are not shown in the table. ^b Lunne et al. (1997) and Robertson et al. (1992). ^c Only used for averaging vertical conductivity. p=0 indicates geometric mean and $p=-\mathrm{1}$ indicates harmonic mean. ^d Freeze and Cherry (1979), Fogg et al. (2000), and Faunt (2009). ^e Botros et al. (2009). ^f Zhu and Mohanty (2002). ^g These are the quantiles used to evaluate the rock physics transform (RPT) for the lowest and highest bins of resistivity values. See Sect. 3.2.3. ^h Bachand et al. (2019), Ganot and Dahlke (2021), Murphy et al. (2021), and Ma et al. (2022).

Download Print Version | Download XLSX

For hydraulic conductivity, we approximate prior uncertainty as the distribution of conductivity values calculated from cone penetration test logs, with separate distributions for coarse- and fine-dominated sediments. Previous work has shown that cone penetration tests can effectively estimate hydraulic conductivity in alluvial sediments (Gribb et al., 1998; Voyiadjis and Song, 2003; Shen et al., 2015) and that the range of values in each distribution agrees with previous conductivity estimates within the southern Central Valley (Botros et al., 2009; Faunt, 2009; Belitz and Phillips, 1995; Bertoldi et al., 1991; Johnson et al., 1968; Zamora, 2008). To calculate hydraulic conductivity in cells that are a mix of fine and coarse sediments, we assume a horizontally layered system and use the arithmetic mean for horizontal hydraulic conductivity (flow parallel to layers). A harmonic mean is typically used for flow perpendicular to layers (Freeze and Cherry, 1979), though the geometric mean may be more applicable in systems with interconnected coarse facies, like the Central Valley (Fogg et al., 2000; Faunt, 2009). Given this uncertainty, we use a power mean function to calculate vertical hydraulic conductivity and include the power mean exponent as an additional parameter during Monte Carlo simulations (row 3 in Table 1), with values between −1 (harmonic mean) and 0 (geometric mean). Because we use different horizontal and vertical hydraulic conductivity averaging schemes, conductivity is anisotropic in cells with coarse fraction between 0 and 1.

Parameter ranges for porosity, residual saturation, van Genuchten α and van Genuchten n (van Genuchten, 1980) are derived from a database of 97 sediment cores collected on the east side of the southern Central Valley (Botros et al., 2009). We use the arithmetic mean when averaging porosity and residual saturation across sediment types. To average van Genuchten parameters, we use the arithmetic mean for n and a power mean averaging scheme for α (Zhu and Mohanty, 2002), with exponent values ranging from 0 (geometric mean) to 1 (arithmetic mean).

Before applying each parameter set to a realization of coarse fraction, we classify each domain into discrete hydrogeologic units. First, we bin coarse-fraction values in each realization into five groups using k-means clustering. Cluster centers vary between each realization given that each realization contains a slightly different distribution of coarse-fraction values. We then iterate through each cluster and replace all coarse-fraction values within that bin with the bin mean. We define a sixth hydrofacies for the bottom layer of model cells (90–150 m depth) and assume that coarse fraction in this facies is equal to the mean coarse fraction in the rest of the domain. Previous work modeling Central Valley recharge has relied on a similar facies designation for deep portions of an aquifer with limited geologic information (Maples et al., 2019).

3.3.3 Model stages

Recharge modeling at the site consists of three separate stages:

1. spin up to steady state using constant forcing;

2. spin up for 1 year with meteorological forcing using CLM;

3. followed by

   • a.

     for flood-MAR simulations, inundation, then a relaxation period with meteorological forcing or

   • b.

     for reference simulations, 2 years of meteorological forcing without any inundation.

To initialize pressure head throughout the model domain, each simulation is spun up to steady state using a specified flux of 0.15 m yr⁻¹ applied at the surface. This value represents the net recharge rate (precipitation + irrigation − evapotranspiration) of a hypothetical almond orchard in the southern San Joaquin Valley, calculated using precipitation and reference evapotranspiration data from the California Irrigation Management Information System (CIMIS; http://www.cimis.water.ca.gov/, last access: 4 August 2022) and localized irrigation data from the California Department of Water Resources (CA DWR, 2018a). Reference evapotranspiration was converted to crop evapotranspiration using monthly crop coefficients for almonds within the Central Valley (Doll and Shackel, 2015). Given the spatial variability of evapotranspiration and irrigation, we calculated separate recharge rates for 15 CIMIS stations between 1998 and 2015 (Fig. S5); the final rate used for model spin-up (0.15 m yr⁻¹) is a weighted average of these 15 individual rates, with stations weighted according to their proximity to the study site. This value constitutes the best estimate of modern recharge rates, though it is of the same order of magnitude as the long-term-average deep-percolation rate across the entirety of the Central Valley (0.27 m yr⁻¹; Faunt, 2009). We apply this forcing to each simulation until the annual change in subsurface storage is <0.01 % (Ajami et al., 2014), which takes between 8 and 131 years depending on the simulation.

Next, each model is run for 1 year using meteorological forcing (Xia et al., 2012; Rodell et al., 2004) covering 1 February 2010 to 31 January 2011. We chose this time period so that inundation occurs in spring 2011, when there was surplus surface water available in the southern Central Valley for flood-MAR (Bachand et al., 2014). The timing of inundation is important, given that antecedent soil moisture may have an impact on recharge rates. Though realistic, springtime inundation may introduce increased crop risk, as recharge operations may overlap with the growing season. During this step, we also simulate crop transpiration and irrigation using CLM. Irrigation is applied between 07:00 and 17:00 LT from March through October, with daily irrigation rates set to 105 % of daily crop evapotranspiration (Doll and Shackel, 2015).

Table 2Inundation scenarios used for flood-MAR.

Download Print Version | Download XLSX

The pressure distribution at the end of the CLM spin-up period is then used as the initial condition for two separate simulations, run concurrently: a flood-MAR simulation and a reference simulation without MAR. Among Central Valley almond orchards, on-farm recharge varies both in the amount of water applied during a season (0.065–0.76 m) and in the frequency of recharge events, ranging from weekly to once per season (Ganot and Dahlke, 2021; Bachand et al., 2014; Murphy et al., 2021). To investigate the effects of flood frequency on each outcome of interest, we vary the frequency of inundation as an additional parameter (Table 1). In all cases, a total of 256 000 m³ (0.8 m³ m⁻²) of water is applied to the orchard between February and May 2011, though we vary the number and frequency of events across four scenarios (Table 2).

To implement each recharge event in ParFlow, we apply a specified head at the upper boundary to cells within the almond orchard and a no-flow boundary condition to the cells outside it. At the start of an inundation step, the specified head over the orchard is equal to the total amount of water applied for that event (e.g., 0.20 m for monthly inundation). For each subsequent time step, the specified head is decreased by the amount of water that infiltrated below ground in the previous time step. This loop is then repeated until no water is left on the surface. This scheme is a computationally efficient approximation of a falling-head infiltration that allows for heterogeneous infiltration with ponded water that flows laterally towards preferential flow paths. Water applied to the orchard cannot move laterally over the surface outside the bounds of the orchard, which mimics an on-farm recharge operation that uses check dams to isolate the recharge area.

Between subsequent recharge events and following the final recharge event, we implement a relaxation phase in which meteorological and irrigation forcing is applied at the surface using CLM. Each flood-MAR simulation consists of a single season of recharge in spring 2011 followed by a relaxation period until 1 February 2013, 2 years after initial inundation. To investigate changes in recharge efficiency beyond 2 years, 10 randomly selected simulations were run for an additional 8 years (ending on 1 February 2021) using meteorological and irrigation forcing without any additional recharge events.

Reference simulations cover the same time period as flood-MAR simulations but do not include any inundation. Meteorological and irrigation forcing is applied for the entire 2-year period using CLM. Comparing reference simulations to simulations with flood-MAR permits differentiating the effects of recharge from other model processes that occur naturally over the same time period, such as an increase in vadose-zone storage due to above-average precipitation.

3.4.1 Evaluating recharge outcomes

As discussed in Sect. 2.1, for each simulation, we evaluate three recharge outcomes of interest.

1. Infiltration rate at the surface. To create maps of infiltration rate across the almond orchard, for each simulation we compute the time-varying volumetric flux of water across the surface into the domain. This is a spatially resolved flux, with variable infiltration rates from one surface cell to the next. We then tabulate the cumulative volume of water that flowed into each individual surface cell across all recharge events. To more easily compare infiltration maps with ponded depth, we divide cumulative volume by surface area and represent the infiltrated volume as a length.

2. Root zone residence time. To evaluate the duration of root zone saturation, we analyze the time series of water content within the root zone for each simulation. Although maximum rooting depth of almond trees can extend to 2 m, effective rooting depth typically ranges from 0.3–1.0 m (Andreu et al., 1997; Koumanov et al., 2006), so we define the root zone as the upper 1 m of the subsurface. Depending on the variety of rootstock, almonds can tolerate 2–14 d of anoxia before bud break (December–early February) and 48 h after bud break (late February–November) without decreased plant growth or crop yield (O'Geen et al., 2015). In our simulations, most recharge events occur after bud break, so we assume that trees can only withstand 48 h of anoxia. During MAR, soil anoxia can develop at water saturation values less than 100 % (Ganot and Dahlke, 2021). Thus, for each root zone cell, we calculate the percentage of simulations in which the cell exceeded 48 h of >90 % saturation over the 2-year inundation–relaxation cycle. Note that this 48 h threshold does not directly account for the risk of tree mortality during high-wind events – which can occur even during brief periods of saturation – though simulations that exceed the threshold have higher risk than those that do not.

3. Saturated-zone recharge efficiency. We define the final outcome of interest, recharge efficiency, as the increase in saturated-zone storage below a site plus the increase in lateral saturated-zone flow away from the site, which increases storage in the broader aquifer system. The two components of recharge are conceptualized in Fig. 1 and contribute to the recharge efficiency, R_eff, at time t as

   $$\begin{array}{}\text{(3)}& R_{\mathrm{eff}}\left(t\right)={\displaystyle \frac{\sum _{i=\mathrm{1}}^t\left(Q_i^{\mathrm{MAR}}-Q_i^{\mathrm{Ref}}\right)\mathrm{\Delta }t+S_t^{\mathrm{MAR}}-S_t^{\mathrm{Ref}}}{V_{\mathrm{app}}}}\times \mathrm{100},\end{array}$$

   

   where $Q_i^{\mathrm{MAR}}$ is the saturated-zone flux out of the domain at time i for the simulation with recharge, $Q_i^{\mathrm{Ref}}$ is the lateral saturated-zone flux at time i for the simulation without MAR, Δt is the time step, $S_t^{\mathrm{MAR}}$ is total saturated-zone storage at time t for the MAR simulation and $S_t^{\mathrm{Ref}}$ is the saturated-zone storage at time t for the reference simulation. Thus, flux recharge ($Q_i^{\mathrm{MAR}}-Q_i^{\mathrm{Ref}})$ is cumulative over time, while storage recharge ($S_t^{\mathrm{MAR}}-S_t^{\mathrm{Ref}}$) is transient. Note that we divide the volumetric recharge at time t by the total volume of water applied (V_app) and multiply it by 100 to report it as a percentage.

To facilitate comparison with previous work, we also quantify deep percolation induced by MAR. Deep percolation is calculated using a modified version of Eq. (3) in which we replace saturated-zone storage with storage below the root zone (below 2 m) and saturated-zone flux with lateral flux out of the domain below the root zone. As with saturated-zone recharge efficiency, we then normalize this volume by the volume of water applied at the surface and express it as a percentage.

3.4.2 Global sensitivity analyses

To identify key parameters and processes affecting recharge, we perform a global sensitivity analysis for each outcome of interest. We utilize distance-based generalized sensitivity analysis (DGSA; Fenwick et al., 2014; Park et al., 2016), a computationally efficient technique that has been successfully applied to groundwater flow and transport models in the past (Hermans et al., 2018; Hoffmann et al., 2019; Perzan et al., 2021). We chose this form of sensitivity analysis because DGSA can take discrete parameters as input (e.g., inundation frequency) and can analyze high-dimensional model responses (e.g., an infiltration map), but it only requires a limited number of model simulations compared to other techniques, like variance-based methods.

In DGSA, simulations are clustered into groups with similar model responses (e.g., similar recharge efficiency) using the pairwise distance between responses. The distribution of an input parameter (e.g., porosity of the fine-grained end-member) is then compared between each cluster; if the prior distribution of a parameter is significantly different from the parameter's distribution within a cluster, then the model is said to be sensitive to that parameter. The area between empirical cumulative distribution functions (CDFs), referred to as the CDF distance, is used to quantify the degree of similarity between parameter distributions. Large CDF distances could result either from model sensitivity or from random model variability due to the limited sample size, so a bootstrapping procedure is used to test the statistical significance of observed CDF distances. If the observed distance is larger than the bootstrapped distance at a chosen significance level, then that parameter is said to be influential (sensitivity metric > 1). If the observed distance is smaller than the bootstrapped distance, the parameter is insignificant (sensitivity metric < 1).

In the present study, we cluster simulations on the pairwise Euclidean distance between model responses with k-medoids clustering, using a different distance for each outcome of interest. For saturated-zone recharge efficiency, we compute distances directly on each simulation's R_eff value at 2 years. For infiltration and root zone residence time, we first reduce the dimension of each model response (80 × 40 pixel images) using a form of principal component analysis known as eigenimage decomposition. We then compute the distance between each simulation's eigenimage scores, using the minimum number of scores required to explain 98 % of variance (typically, 40–50 components). Finally, we cluster simulations using these distances and perform the bootstrapping procedure described above to calculate sensitivity metrics.

In each sensitivity analysis, we include all model input parameters listed in Table 1. For a given simulation, a change in model response can result either from explicit perturbation of these input parameters or from the particular spatial configuration of hydrofacies for that simulation. For example, one stochastic realization of the model domain may have a higher percentage of fine-grained facies than another, limiting recharge. To account for this spatial variability in sensitivity analyses, we perform principal component analysis across all realizations of coarse fraction (Fig. 3c) and include each simulation's principal component scores as additional input parameters to DGSA.

Figure 4Maps of the first two outcomes of interest, aggregated across all 284 simulations. All simulations received 0.8 m of inundation, though the median volume of water that infiltrated into each surface cell (a) ranges from 0.10 to 1.40 m³ m⁻². The percent of simulations in which each surface cell exceeded the 48 h root zone saturation threshold is shown in panel (b), with the median coarse fraction of each surface cell shown in panel (c). Note that each map only shows the portion of the model domain immediately below the almond orchard.

Download

4.2.1 Surface infiltration rate

Across all simulations, the average infiltration rate during flood-MAR is 0.15±0.29 m h⁻¹ (mean ± standard deviation). Infiltration rates are typically higher for simulations with smaller, more frequent inundation events: a mean of 0.06 m h⁻¹ for simulations with a single 80 cm recharge event as compared to a mean of 0.24 m h⁻¹ for simulations with weekly 5 cm recharge events. Within each recharge event, infiltration rates slow down over time; among simulations with a single 80 cm recharge event, on average, the infiltration rate during the last time step is 17 % of the initial rate.

Infiltration rates also vary with sediment type; on average, cells with coarse fraction > 0.6 account for 91 % of the volume of water that infiltrates into the subsurface even though they only make up 67 % of the orchard surface area. Maps of infiltration below the orchard (Figs. 4a and S8) show that water preferentially infiltrates through coarse-grained cells, with 1.8–458 (median 10.7) times more water infiltrating into 100 % coarse cells than into 100 % fine-grained cells. Eigenimage decomposition and k-medoids clustering of all infiltration maps reveals two distinct groups of simulations. Simulations in both clusters exhibit preferential infiltration through coarse facies, though the extent of preferential infiltration is greater in one cluster than the other (Fig. S8).

For a small subset of simulations (1.8 %), the initial infiltration rate through fine-grained facies is faster than through coarse-grained facies. However, this relationship quickly switches once sediments become saturated (<1 h), implying that coarse facies briefly acted as capillary barriers. Within each simulation, the greatest infiltration rates typically occur within coarse-grained cells on the edge of the orchard. These cells have higher infiltration capacity because – once below the surface – water can spread laterally into adjacent, dry cells beyond the orchard boundary.

4.2.2 Root zone saturation

In 65 % of simulations, at least one root zone model cell remains saturated longer than 48 h, which is the estimated time almond trees can endure soil anoxia without diminished growth or yield. Of the cells that exceed this threshold, 81 % are predominantly fine-grained, with coarse fraction < 0.6 (Fig. S9). Conversely, in 8 % of simulations ($\mathrm{23}/\mathrm{284}$), at least one 100 % coarse-grained cell exceeds the saturation threshold. Aggregating results across simulations reveals a similar trend; contiguous blocks of fine-grained sediments are >20 % likely to remain saturated for at least 2 d following inundation (Fig. 4). Collectively, this area accounts for 9 % of the orchard or 2.85 ha.

Figure 5Vertical cross sections of relative saturation (a–d), coarse fraction (e), pressure head (f–i) and the change in pressure head since day 0 (j–m) for a single simulation for 2 years (730 d) following a single 0.8 m inundation event. Transient increases in saturation within coarse-grained sediments can be observed at 30 d and within fine-grained sediments at 30 and 60 d. The dashed line shows the position of the water table, with obvious water table mounding at 160 and 730 d. Only the upper 50 m of the model domain is shown.

Download

The increase in root zone saturation also drives a minor increase in transpiration within the almond orchard. In the first week after inundation, transpiration is an average of 37 % higher in flood-MAR simulations than in reference simulations (2.16 mm d⁻¹ vs. 1.57 mm d⁻¹), which is reflective of the early growing season in California's Mediterranean climate. After 7 d, the difference decreases, and mean transpiration for the remainder of the growing season is only slightly higher (2.2 % or 0.09 mm d⁻¹). Over the whole simulation period, increased transpiration only accounts for ∼4.8 % (0.038 m) of the total applied water.

4.2.3 Saturated-zone recharge efficiency

Inundation at the surface produces a pressure pulse that cascades through the subsurface, resulting in a transient increase in saturation throughout the vadose zone (Fig. 5). The pulse moves rapidly past the root zone towards the water table. It is important to note, however, that this rapid change in storage corresponds to the transmission of pressure and not to the direct movement of water particles. Water applied at the surface displaces water in the root zone, which displaces water below the root zone and so on. Thus, though saturation may increase at 30 m depth, it does not necessarily indicate that water particles applied at the surface have moved to this depth.

As the pressure pulse moves deeper, it causes the water table to rise (Fig. 5), which represents an increase in saturated-zone storage (storage recharge; Fig. 1). Mounding of the water table induces an increase in discharge out the sides of the domain, which constitutes recharge to the broader aquifer system (flux recharge). Temporal changes in these two components of recharge are plotted for a single simulation in Fig. 6. Storage recharge is the first component to increase following inundation, starting to rise between 34 and 683 d (mean 214 d) after the first inundation event. As storage recharge peaks and slowly decreases, flux recharge decreases as well. However, since flux recharge is cumulative over time (Eq. 3), overall saturated-zone recharge efficiency continues to increase until the end of the simulation (Fig. 6).

Figure 6Temporal changes in the two components of saturated-zone recharge efficiency – storage recharge (blue) and flux recharge (green) – for a single simulation for 10 years following inundation. Flux recharge is cumulative over time, while storage recharge is transient. Thus, even after the water table returns to its pre-MAR level at ∼6 years, water content in the vadose zone remains higher than in simulations without recharge. As this water percolates to the saturated zone, it contributes to flux recharge. The recharge efficiency of this simulation is 46 % after 2 years, which is higher than the median across all simulations (32.7 %).

Download

Among all simulations, deep percolation increases rapidly following inundation. Fourteen days after the final inundation event, the increase in deep percolation accounts for 88.3±17 % of the volume of water applied at the surface (Fig. 7a). By comparison, saturated-zone recharge efficiency is 32.7±26.3 % after 2 years (Fig. 7b). Ten of these simulations were run for an additional 8 years (total of 10 years post-inundation), and results show that saturated-zone recharge efficiency continues to increase beyond 2 years (Fig. 6). However, none of the simulations reach 100 % recharge, and, though some get close, they plateau around 95 % (e.g., Fig. S10). Given this paper's focus on short-term recharge outcomes (<2 years), more long-term simulations are needed to confirm that these results are consistent across the full parameter space.

Figure 7Histogram of deep-percolation efficiency at 14 d (a) and saturated-zone recharge efficiency at 2 years (b) for all 284 simulations. Deep percolation refers to the proportion of water applied at the surface that has moved beyond the root zone within 14 d of the final inundation event. Saturated-zone recharge efficiency ranges from 0.28 % to 87.7 % (median 32.7 %), though the distribution is skewed towards lower values.

Download

Coarse- and fine-grained facies respond differently to recharge, especially between the bottom of the root zone (2 m) and the water table (45 m). Coarse facies in this zone respond rapidly to inundation; on average, water content begins increasing 3.5 d after the first recharge event and peaks 11.2 d after the last event (Fig. 8a). In contrast, water content in the fine facies begins increasing after 22.2 d and peaks 246 d after the final inundation event. Fine-grained sediments are also slow to drain following inundation; 2 years after the first recharge event, storage in fine facies between 2 and 45 m depth is a median of 94 000 m³ higher than in reference simulations, which represents 37 % of the volume of water applied at the surface (Fig. 8b). By comparison, coarse facies hold significantly less water than the fine facies after 2 years (P<0.001, paired t test), which is a median increase of 10 700 m³ (4.2 % of applied water) relative to reference simulations.

Figure 8Storage change in the vadose zone (VZ) below 2 m (a) for a simulation with low saturated-zone recharge efficiency (8.7 % after 2 years). The percent change in vadose-zone storage at 2 years for all 284 simulations is shown as a box plot in panel (b). Results of the paired t test indicate that the difference in the two populations in panel (b) is statistically significant (P<0.001). In both subplots, coarse facies (yellow) are defined as any model cells with coarse fraction ≥ 0.5; fine facies (blue) correspond to coarse fraction < 0.5.

Download