Examining subsurface response to an extreme precipitation event using HYDRUS‐1D

North‐central Colorado experienced an extreme precipitation event (EPE) in September 2013, during which the equivalent of 80% of the region's annual average precipitation fell in a few days. Widespread flooding occurred above ground, but the short‐ and long‐term subsurface response remains unclear. The objective of the study is to better understand the dynamic subsurface response, namely how the water table and soil water storage responded to a large amount of infiltration in a short period of time and how the hydrologic properties of the subsurface influence the response. Better understanding of subsurface response to EPEs is expected to increase with the advent of more intense and frequent EPEs in the coming decades. A one‐dimensional subsurface flow model using HYDRUS‐1D, was built to simulate and examine infiltration of an EPE at a site in the Boulder Creek Watershed, Colorado. Model calibration was conducted with local field data to fit site observations over a 6‐yr period. A rapid water table depth response in field observations was observed, with the modeled subsurface storing water for 18 mo acting as a hydro‐buffer during recovery. To examine influence on model results, a sensitivity study of soil hydraulic parameters was conducted. The sensitivity study found that changes in n, an empirical parameter related to pore‐size distribution, most significantly affects water table depth. The implications are that one‐dimensional models may provide useful estimates of water table fluctuations and subsurface hydro‐buffer capacities in response to EPEs, which could be of use to regions preparing for EPE effect on water resources.

dollars of property and infrastructure damage and the tragic loss of eight lives (Coffman, 2013). The above-surface flooding response to similar extreme events has been documented and photographed in Colorado for over 125 yr (BASIN, 2005). However, the subsurface physical response to EPEs is not as easily observed or measured in real-time, making it one of the more poorly understood hydrogeologic topics of the 21st century (Vereecken et al., 2015).
Subsurface response to EPEs may involve rapid fluctuations of soil water storage and abrupt water table depth changes (Freeze & Witherspoon, 1967;Jasechko & Taylor, 2015;Tashie et al., 2016). French et al. (1996) examined subsurface response to regular and intense precipitation at a highelevation study site in the southwest region of the United States. The shallow soils at this high-elevation site extended 1 m below ground level to fractured bedrock. The hypothesis stated that if infiltrating water could penetrate this 1-m physical transition, then it could likely result in groundwater recharge (French et al., 1996). Examination of soil water data found that fall and winter events (October-April) more often infiltrated below a 1-m depth. It was suggested that this was due to (a) the longer duration of the precipitation events and snowmelt and (b) lower evapotranspiration rates. In contrast, summer events (May-August) were observed to be of short duration and affected by high evapotranspiration, diminishing infiltration past the 1-m depth. The study concluded that it was unclear how soil profiles deeper than 1 m may respond to varying precipitation events, either normal or intense, and more research was suggested. Ng et al. (2010) studied the effects of different climate predictions on diffuse episodic recharge for a study site in the southern High Plains of the United States. They found that high-rainfall periods, equivalent to EPEs, were more likely to result in recharge during the winter months (December-March), when evapotranspiration is lower and plant roots are dormant. At the same time, the study acknowledged that EPEs and interannual variability were not represented, which may have underestimated a significant fraction of the total recharge predicted. They called for future studies to use field measurements of interannual variability, including EPEs, where possible, especially for predicting recharge in arid environments, as is the case for the southwestern United States, the Middle East, most of Australia, and northern African. Shao et al. (2018) found that consecutive wet years promoted groundwater recharge more significantly than years with average precipitation. This is an important finding because the numbers of wetter years and drier years are expected to increase in the future, while years of average precipitation are expected to decrease (Lehmann et al., 2015;Trenberth, 2011;Wasko et al., 2016). In particular, precipitation events are expected to shorten in duration and increase in intensity (Pendergrass & Knutti, 2018;Pfahl et al., 2017;Prein et al., 2017). This highlights an urgent need to move

Core Ideas
• Effects of extreme precipitation events on subsurface processes were explored through modeling. • A model was built to fit field data and simulate varying water table depth over a 6-yr period. • The modeled subsurface stores water for 18 mo, acting as a hydro-buffer during recovery. • Model results are most sensitive to n, an empirical value related to pore-size distribution.
beyond annual precipitation and use comprehensive interannual variability, including EPEs, in modeling efforts to better understand subsurface response. A better understanding of subsurface response to EPEs can improve future planning of groundwater resource allocations (Gurdak et al., 2009;Kløve et al., 2013). Using the HYDRUS-1D subsurface flow model with average soil hydraulic parameters estimated by Schaap et al. (2001), Corona et al. (2018) found that the prescribed flux (precipitation) and period (30, 180, 365, and 730 d) were the most statistically significant predictors of whether an infiltration flux became steady or transient recharge. The study examined the combinations of daily precipitation rates and soil types that could lead to recharge, finding that daily precipitation of lower intensity and finer-grained soils resulted in little to no recharge, whereas daily precipitation of greater intensity and coarser-grained soils, like sand, resulted in greater recharge. A sensitivity study of parameter influence on infiltration fluxes in the vadose zone was also conducted but did not consider precipitation variability within a period and the subsequent soil response. Where available, field observations of water table depths and soil water content changes are a useful guide to better understand the EPE-subsurface connection (Jasechko & Taylor, 2015;Thomas et al., 2016). Where data are limited, numerical models coupled with available field data can provide some understanding of the EPEsubsurface connection (Dadgar et al., 2020;Mo'allim et al., 2018).
The objective of this study was to build a HYDRUS-1D model with local field data and examine how the subsurface responded to the EPE that affected north-central Colorado in September 2013. This study addressed the following questions: How does the water table fluctuate in response to the EPE? What can a sensitivity study show about parameter uncertainty? How does soil water storage respond to the EPE?
Exploring these questions can shed new light on infiltration flux through the subsurface, dynamic changes in subsurface water storage, and the temporal extent of subsurface system response to EPEs.

Study area and field data
The Gordon Gulch drainage basin is located 30 km west of Boulder, CO (2,400-2,800 m asl) ( Figure 1). The basin has a total area of ∼3.6 km 2 ; the upper basin has an area of ∼1.0 km 2 , and the lower basin has an area of ∼2.6 km 2 . Gordon Gulch lies in a montane climate zone, with an average annual precipitation of 520 mm yr -1 (BcCZO, 2020). The basin, hereinafter "Gordon Gulch," was chosen due to the extensive data available. In 2011, the Boulder Creek Critical Zone Observatory (BcCZO) installed six wells in the upper basin of Gordon Gulch. The wells have been monitored and maintained by the BcCZO since December 2011. Wells 1, 2, and 6 have working pressure transducers that record water table depth variations at 10-min intervals. Well 1 was chosen because it has the largest vadose zone extent (∼10 m) of the three wells measuring water table depth, making it the ideal candidate for examining pressure head and soil moisture response to EPEs above the water table. Well 1 is at a horizontal distance of 150 m away from the small ephemeral stream and 12 m higher in elevation above the streambed.
The flow record of the nearest stream gauge has not shown evidence of stream influence on Well 1 (Anderson & Ragar, 2021a;Henning, 2016;Salberg, 2021). In contrast to Well 1, Well 2 is influenced by nearby streamflow, whereas Well 6 is affected by lateral flow and upslope infiltration (Anderson & Ragar, 2021b;Henning, 2016;Salberg, 2021). Wells 2 and 6 are henceforth omitted. Well 1 is screened from a depth of 9.4 m to the bottom of the well at 18.55 m, with an average water table depth of ∼9.6 m. The soil lithology of Well 1 is considered representative of the subsurface of the study site based on the geophysical surveys conducted by Befus et al. (2011). For this study, the depth of the well penetrating 10 m of the unsaturated zone and 10 m of saturated zone makes it a good candidate for studying the dynamics of water table fluctuation.
The daily precipitation data are derived from a meteorological station ∼3 km south of the well site, at the Sugarloaf Station CO94, managed by the National Atmospheric Deposition Program (NADP, 2020). Although not co-located, the station experienced the same amount of precipitation as the well site (Uccellini, 2014). It has a record of daily precipitation from 1986 to 2017 (NADP, 2020). To align with the available daily precipitation record, only December 2011-December 2017 water table depth data are used in this study.

Unsaturated flow in the vadose zone
Subsurface processes are difficult to observe and quantify in real-time. Numerical models, such as the public domain HYDRUS source code (Šimůnek et al., 2008), solve Richards' equation to examine one-dimensional (1D) water flow in an unsaturated-saturated porous medium and calculate the overall water mass balance. Ignoring the air-phase flow and thermal effects, Richards' equation has the following form (Richards, 1931): where θ is the water content, t is time (T), ψ is the pressure head (L), K = K(ψ) is the unsaturated hydraulic conductivity dependent on the pressure head (LT −1 ), and z is the downward distance from the ground surface (L). HYDRUS-1D implements the van Genuchten (1980) equations that use Mualem's (1976) pore-size distribution model. The van Genuchten (1980) equations are a set of closed-form analytical expressions that provide continuous functional relationships for the soil water retention, θ(ψ), and the unsaturated hydraulic conductivity, K(ψ), of a soil: where S e is the effective saturation: where θ r and θ s denote the residual and saturated water content, respectively; K s is the saturated hydraulic conductivity; α is a parameter inversely related to the air-entry pressure; n is the pore-size distribution; and l is a pore-connectivity parameter (van Genuchten, 1980;Mualem, 1976). Initially, HYDRUS-1D solves equation 1 for ψ(z). The unsaturated hydraulic conductivity as a function of pressure head, K(ψ) in Equation 1, is obtained from Equations 2-5. During each time step, K(ψ) and θ(ψ) values are obtained iteratively, where ψ(z) values from the prior time step and specified soil parameters (K s , θ r , θ s , α, n, m, l) in Equations 2-5 are used to compute K(ψ) at every depth and then used to solve Equation 1 for ψ(z). Once the ψ(z) values between iterations converge, HYDRUS-1D proceeds to the next time step.

Model setup
The HYDRUS-1D model can be used to analyze water movement in partially saturated and fully saturated porous media (Šimůnek et al., 2013). The model domain is a 1D vertical column extending downward from the land surface to a depth of 50 m. In the model domain, the pressure heads change from negative in the unsaturated zone to positive in the saturated zone. The water table position is found where the pressure head is zero. To reflect the depth of the average water table at the site (Salberg, 2021), the modeled water table was initialized at a 10-m depth. An initial sensitivity study (not shown) was conducted to identify whether a varying soil column length (z = 20, 50, or 100 m) affected water table fluctuations. Model runs showed that the water table fluctuations (where ψ = 0 at z = 10 m) were very similar when comparing the 50-and 100-m lengths; as a result, a 50-m length was chosen. The soil column is discretized into 241 nodes. A sensitivity study of refining soil profile discretization (101, 201, 241, 301, and 501 nodes) found no significant differences in model results with profile discretization of finer than 241 nodes. The area of interest in this study is the unsaturated zone (z = 0-10 m), with denser node spacing in the first 10 m of the profile (z = 0-10 m), with a spacing of 0.072 m between nodes. From z = 10 to 50 m, node spacing is less dense and increases linearly from 0.072 to 0.72 m.
A prescribed flux is applied as the top boundary condition ( Figure 2). A deep drainage flux is applied at the bottom of the soil column. The drainage flux out of the column, q(ψ), is approximated by the following expression (Hopmans & Stricker, 1989): The variable q(ψ) (LT −1 ) is the flux across the bottom boundary, and A (LT −1 ) and B (L −1 ) are adjustable empirical parameters. The ψ bottom (L) is the pressure head at the bottom boundary, and GWL (L) is a long-term equilibrium water table position relative to the bottom boundary, where GWL = 50 m for this study. We calibrated the A and B parameters iteratively to fit the available water table data following the methodology of Neto et al. (2016). In this study, the unit of length (L) is meters (m), and the unit of time (T) is days (d). Model setup allows for a conceptualization of the system such that model results can explain the field observations. The lithology is characterized by four soil types: soil, saprolite, weathered bedrock, and unweathered bedrock ( Figure 2).

Vadose Zone Journal
F I G U R E 2 One-dimensional model setup. The top boundary is defined by a prescribed flux condition. The bottom boundary is defined by a deep drainage condition. The soil lithology is characterized by the following thickness: 1 m soil, 13 m saprolite, 3 m weathered bedrock, and 33 m unweathered bedrock This composition has been identified by geophysical surveys, soil pit hand-dug records, and well lithology records (Anderson et al., 2013;Befus et al., 2011;Shea, 2013). The initial parameters, K s , θ r , θ s , α, n, (Table 1), are obtained from a previous calibrated model and field study of the upper Gordon Gulch drainage basin (Henning, 2016). Figure 2 shows the soil stratigraphy of the well, where a portion of the saprolite layer (10-14 m), the entire weathered bedrock layer (14-7 m), and the entire unweathered bedrock layer (17-50 m) are below the 10-m water table and considered fully saturated. Because there is no pore space for air to enter in the saturated zone, changes to the α value for the unweathered and weathered bedrock layers have little influence on model results. Similarly, adjusting the n value for the fully saturated layers would not influence model results.
During initial calibration, the tortuosity parameter l was calibrated in conjunction with the other parameters. A literature review found that a tortuosity value of l = 0.5 led to poor predictions of unsaturated hydraulic conductivity in 235 soil samples of varying textures (Schaap & Leij, 2000). Schaap and Leij (2000) suggested that the tortuosity be optimized at values of −1 or lower. Yates et al. (1992) suggested that optimal values for l can range from −3 to >100. Thus, simulations were run with the initial value of 0.5 and in increments/decrements of 0.5 from −10 to 10. A tortuosity of l = −2 was determined to be the optimal value for this study.
The model is initialized with a prescribed pressure head distribution that linearly decreases from ψ = −10 m at the surface (z = 0 m) to ψ = 40 m at the bottom of column (z = 50 m). For model spin-up we use the daily average precipitation minus evapotranspiration as the recharge boundary condition on the model top. To account for evapotranspiration, we examined previous studies that estimated potential evapotranspiration for the Gordon Gulch basin and found that potential evapotranspiration values may range from 31 to >100% of the annual average precipitation (Hale, 2022;Langston et al., 2015;Salberg, 2021). Most recently, Salberg (2021) calculated monthly total evapotranspiration loss for a catchmentscale water budget of Gordon Gulch. They suggested that ∼435 mm of the average 580 mm annual precipitation of the Gordon Gulch drainage basin was lost to evapotranspiration, a ∼75% average. Thus, an evapotranspiration rate of 75% is used, which is also backed by regional climate model estimates of evapotranspiration (55-85%) (Sanford & Selnick, 2013;Reitz et al., 2017). Given the 75% loss to evapotranspiration, the model recharge is 25% of the precipitation amount.
The model was spun-up for 400 d to allow the model to equilibrate to a steady state, considering the initial condition. The initial condition is the state from which the model's transient simulations can initiate. For the transient simulations, we apply the 2011-2017 precipitation record minus evapotranspiration. The transient model is run and calibrated by iteratively adjusting the soil hydraulic parameters K s , θ r , θ s , α, and n. Calibration results in a parameter scenario that allows the modeled water table to best match the field observations. We note that though these parameters may not be the only ground truth parameter scenario in the field, they represent the best T A B L E 1 van Genuchten parameters used to define the four soil layers for the base case

Statistical indicators
The R 2 and the RMSE are used in this study to provide a firstorder assessment comparing the modeled and observed water table fluctuations. The R 2 describes the proportion of the variance of the field observation data that can be explained by the model. The R 2 can range from 0 to 1, with higher values indicating less probability of error variance; R 2 values >.50 are acceptable (Moriasi et al., 2007). The RMSE index quantifies the error of a model in predicting observations by measuring the residual spread from the observations. In Equation 7, P i denotes the model predicted values, and O i denotes the field observations for a sample n.
The RMSE is the square root of the average of the squared errors. Thus, a lower RMSE typically suggests a lower chance of error, with an RMSE of zero suggesting a perfect fit between the predicted and observed. The RMSE is commonly used because it calculates the error of a comparison in the units of the constituent of interest (Moriasi et al., 2007;Reusser et al., 2009). It has been proposed that RMSE values less than half of the standard deviation of the observations may indicate a low probability of error (Moriasi et al., 2007;Singh et al., 2005).

Sensitivity analysis
Once the base case is constructed, we conduct (a) a sensitivity study of parameters to understand parameter influence on results and (b) post-processing analysis of simulated soil water storage. The sensitivity study focuses on the local unsaturated zone, where we consider only the soil (0-1 m) and the saprolite (1-14 m) layers. Six sensitivity simulations are conducted for three parameters (two per parameter): the residual water content (θ r ), the empirical parameter inversely related to the air-entry pressure (α), and the pore-size distribution (n). For each sensitivity simulation, the respective parameter is increased or decreased for the soil and saprolite layers of the base case model ( Table 2). The chosen range of values reflects low and high averages across the 12 soil textural classes predicted by the ROSETTA Soil Catalog (Schaap et al., 2001). For θ r , a decrease from the base case (θ r = 0.01) suggests less water remaining in a soil pore at high tension. An increase (θ r = 0.15) suggests more water remaining in a soil pore at high tension, which may be indicative of a clay-rich soil (van Genuchten , 1980). For α, a decrease (α = 0.10) suggests a higher minimum matric suction required for air to enter pore spaces, whereas an increase (α = 2.00) suggests a lower minimum matric suction required. For n, a decrease (n = 1.25) represents a wider pore-size distribution (larger variation in pore sizes in the soil), whereas an increase (n = 5.50) represents a narrower pore-size distribution.
A preliminary sensitivity analysis of changes in the saturated water content (θ s ) showed similar trends in water table fluctuations to changes in residual water content (θ r ). The residual water content is often overlooked and difficult to assess, in contrast to θ s , which is often easier to characterize (van Genuchten, 1980;Vanapalli et al., 1998). In a landmark paper, van Genuchten (1980) suggested that the poor matching between the predictive soil water retention curve and the observed curve could be attributed to the θ r value, which was estimated to be zero. van Genuchten (1980) suggested that future studies consider the importance of having an independent procedure for estimating θ r . Despite decades of progress, correctly assessing the θ r for a soil remains a challenge (Nimmo, 2006;Vanapalli et al., 1998;Vogel et al., 2001). Including it in the sensitivity analysis as in the present study could help us better understand the consequences of changing the residual water content in subsurface flow modeling. Changes to the saturated hydraulic conductivity, K s , and tortuosity, l, had little effect on model results and are henceforth omitted.

Base case
Water  (Figure 3a, red arrow) is an exception, where the base case model predicts a shallower water table depth than the field observations. A linear regression between the model predictions and the field observations suggest that the model predicts shallower water table depths relative to the field observations ( Figure 3b). For example, during the EPE, an ∼8.1-m field measurement was predicted to be ∼7.6 m by the base case. Although the model overestimates water table depths at certain times during the EPE, the changes remain within ∼6% of the field observations, suggesting an overall good fit.
The green squares indicate a rainy May in 2015, and the base case model predicts deeper water table depths compared with the field observations. Around May 2015 (Figure 3a), the base case predicts a water table depth of ∼8.3 m, compared with the shallower depth of ∼8.0 m of the field observations. This underestimation is within 4% of the field observation. The linear regression gives a R 2 value of .56 for the 6-yr time series. Singh et al. (2005) published guidelines stating that RMSE values less than half the standard deviation of the field observation data could be interpreted as indicating a good fit of the model to the field observations. The standard deviation of the field observations for the 6-yr period is 0.40 m. Following the guidelines of Singh et al. (2005)

Sensitivity analysis of hydraulic parameters
With the base case established, we examine parameter uncertainty in model results by conducting a sensitivity analysis of

(a) (b) (c) (d)
F I G U R E 4 (a) Time series of field observations and model predictions of water table depths for varying residual water content, θ r . The "decrease in θ r " represents a 0.09 decrease from the base case to 0.01 for both Layers 1 and 2. The "increase in θ r " represents a 0.05 increase to 0.15 for both Layers 1 and 2. (b-d) Linear regression of the model predictions versus the field observations for the decrease (b) and increase (d) in θ r three van Genuchten parameters: θ r , α, and n for the soil layer and saprolite layer, respectively. For each model run, one of the three parameters is increased or decreased (Table 2). Figure 4a shows model sensitivity to changes in residual water content, θ r . When θ r is decreased for both Layers 1 and 2, the modeled water table depth consistently predicts decreases in water table depths (overestimates). The decrease in θ r negatively affects the correlation between the field observed and base case, lowering the R 2 to .45. The respective RMSE (0.24 m) indicates an acceptable model fit. Increasing θ r results in a dampened response where water table fluctuations are slightly subdued (Figure 4a,d). An increase in θ r results in a slightly exaggerated water table response. The R 2 remains at .56, and the RMSE (0.30 m) indicates an acceptable model fit at the cusp of being a poor model fit. Figure 5a shows model sensitivity to changes in α (m −1 ), an empirical parameter that is inversely related to the airentry pressure value. The α value chosen for the base case is α = 0.18 (m −1 ). For the sensitivity study, α is initially decreased for the soil and saprolite layers (Layers 1 and 2) to 0.10 (m −1 ). A decrease in α allows the model to conform to field observations but consistently predicts slightly deeper water table depths (Figure 5a). Figure 5b shows that a decrease in α slightly improves the R 2 (.57). The respective RMSE (0.28 m) indicates an acceptable model fit. When α values for are increased to 2.0, water table fluctuations are subdued (Figure 5a). The R 2 value decreases (.21 m) with an increase in α, suggesting a poor correlation (Figure 5d) between the model results and field observations. The respective RMSE (0.28 m) indicates an acceptable model fit. Figure 6a shows the model sensitivity to changes in n (1), an empirical parameter that characterizes pore-size distribution. In HYDRUS-1D, n must be >1 (Equation 4). The default values for n are set by HYDRUS-1D and vary by soil type. For the base case, n = 1.50. For the sensitivity study, n is decreased to 1.25. A decrease in n allows the model to better conform to field observations but also consistently predicts slightly deeper water table depths (Figure 6a). The decrease

Soil water storage of the base case
In HYDRUS-1D, the soil water storage, V (m), is defined as the volume of water per unit area at a point in time. The V is calculated as: where θ i and θ i+1 are water contents evaluated at elements i and i + 1, and Δz i is the size of the element (Šimůnek et al., 2008, 2013). The summation in Equation 8 is taken over the 241 elements in the flow domain. Figure 7 shows the base case variability in soil water storage, V, from December 2011 to December 2017 by water year (WY). For example, WY2013 denotes the water year from 1 Oct. 2012 to 30 Sept. 2013. The dashed lines represent the average soil water storage for the respective water year. The soil water storage, V, for the profile ranges from 6.50 to 6.62 m (Figure 7). The V values (e.g., 6.50 m) are the product of the average water content across the entire domain (e.g., 0.13) and the total column length (50 m; 0.13 × 50 m = 6.50 m). As such, a higher V suggests that a greater portion of the available pore space is saturated, indicating that the subsurface is wetter than average (Figure 7). A lower V suggests that less of the available pore space is saturated, indicating that the subsurface is of average wet conditions or drier. A wetter subsurface may result in recharge, whereas a drier subsurface may result in little or no recharge. 6-yr record. The following WY2014 and WY2015 (dotted lines) exhibit the highest average V per water year, indicating that the EPE influenced subsurface processes for two water years after its occurrence.

How does the subsurface respond to the 2013 EPE in terms of water table fluctuations?
The September 2013 EPE had a long-term consequence in the subsurface, as shown by the Gordon Gulch field observations of water table fluctuations. The 1D, four-layered, homogeneous base case model represented the best scenario based on available data and provided model output that could most closely match the field observations. Visually, a comparison of the field observations and base case model showed good compatibility. Figure 3a shows a consistent downward trend in the water table depth that reaches its deepest point (∼9.6 m) every March of every year except 2014. In March 2014 (after the EPE), the deepest point is ∼9.25 m, a shallower depth than every other year. It is not until March 2015 that the water table deepens to ∼9.6 m again. We speculate that the contrast between the consistent water table depths of ∼9.6 m and the shallower water table depth of March 2014 (∼9.2 m) provide evidence that the water table remained shallower in large part due to the EPE footprint, which remained until at least March 2015. The field observations and model results agree that the subsurface continued to respond to the EPE infiltration flux for at least 18 mo after the event, which is longer than previously suggested (Henning, 2016).
For further comparison between the field observations and model results, we calculated the R 2 to measure the goodness of fit between the field observations and base case model. In addition, the RMSE was used to calculate the square root of the variance of the residuals to indicate how close the observed data are to the model results. The base case model is considered an acceptable model fit, as indicated by the R 2 (.56) coefficient and the RMSE (0.23 m) value (Table 3). The results present opportunities for improvements while highlighting the limitations of the 1D modeling approach. For example, the 1D modeling approach does not simulate lateral flow process at the hillslope scale or regional scale, which could affect the goodness of fit. The greatest deviation in correlation occurred during the EPE, where the base case predicted a higher water table (∼7.5 m) than the field observations (∼8.2 m), though this value was still within a 10% margin of the field data. Although out of the scope of this study, a two-dimensional or three-dimensional model accounting for lateral flow may improve the goodness of fit between the field observations and modeled results.

Sensitivity analysis of hydraulic parameters
A sensitivity analysis examined how model response may be affected by parameter change. The R 2 value calculated for each sensitivity analysis ranged from .21 (poor correlation) to .58 (acceptable correlation). The RMSE value calculated for each sensitivity analysis ranged from 0.24 to 0.30 m (Table 3).
All RMSE values were indicative of acceptable model fit, with one model (increase in θ r ) at the cusp between acceptable and poor, despite an acceptable R 2 value. The lack of RMSE values indicating a good fit model could be attributable to the high sensitivity that the RMSE has to outliers (i.e., the largest differences between field observations and model results). Two EPEs of different temporal extents-the September 2013 EPE and the May 2015 month-long rain event-resulted in large differences between the field observations of water table depths and the model results. These EPE-derived outliers skewed the RMSE away from indicating a good model fit. Future studies examining EPEs may benefit from statistical methods that are not strongly biased toward outliers.
The visual outcome of the sensitivity runs can be described by two general responses: dampened or exaggerated (Table 3). An exaggerated response indicates shallower and deeper water table depths relative to the base case. A dampened response would indicate the opposite of exaggeration (i.e., more smoothed, tempered variation).

4.2.1
Residual water content From the sensitivity analysis, decreasing the residual water content, θ r (Figure 4d), yields a dampened water table response. Increasing the θ r (Figure 4b) yields an exaggerated water table response. These responses are likely due to the local effective porosity of the material (Horton et al., 1988). The effective porosity, also thought of as the "drainable porosity," is defined as the percentage of interconnected void space with respect to the bulk volume (Brooks & Corey, 1964).
where φ is the effective porosity (1), V p is the total volume of interconnected voids (m 3 ), and V b is the bulk volume (m 3 ). A soil with a higher φ has a larger total volume of V p relative to V b . Decreasing θ r (with θ s held constant) increases the total volume of V p relative to V b , indicating a higher φ.
Decreasing θ r can result in water being more readily held in pore spaces (higher φ), slowing the rate of flow. Water held in pore spaces may result in a slower drainage out of the pore spaces, which can dampen fluctuations of the water table. In contrast, increasing θ r (with θ s held constant) decreases the total volume of V p relative to V b , indicating a lower φ. A lower φ can be indicative of a faster draining soil with minimal available pore space for water to fill. A low φ may allow faster flow through the subsurface, thus resulting in more rapid water table response (Figure 4b).

Air-entry pressure
The alpha parameter, α (m −1 ), is inversely related to the air-entry pressure, the matric suction value required to fill (empty) pore spaces (Nimmo, 2004). Its purpose is to serve as an approximation of the steepest section of the soil water retention curve (van Genuchten, 1980). Where the soil water retention curve is steepest, the water content (typically plotted on the y-axis) is most sensitive to changes in the pressure head (typically plotted on the x-axis). This is also mathematically evident in the van Genuchten equations (Equation 2), where α as part of the denominator influences the fraction from which the quotient determines water content as a function of pressure head, θ(ψ) given ψ < 0 (Brooks & Corey, 1964;Mualem, 1976;van Genuchten, 1980). The base case α = 0.18 m −1 is considered the best fit for the field data.
Results from the sensitivity study show that decreasing α (relative to the base case) causes the water table to fluctuate more rapidly (Figure 5a). Conceptually, a decrease in α translates to an increase (due to the inverse relation) in the minimum matric suction value that air must attain to enter a pore space. While pressure builds so that air can attain the higher matric suction to enter a pore space, water can enter (drain from) these same pores with greater ease. The ease at which water can fill (drain) pores in the unsaturated zone is known as the unsaturated hydraulic conductivity, K (m d −1 ). A decrease in α can allow water to fill (empty) pores with greater ease, resulting in a larger K value, which consequently results in more rapid downward flow, and thus more dramatic water table fluctuations, as seen in Figure 5a.
An increase in α translates to a decrease (due to the inverse relation) in the minimum matric suction value that air must attain to enter a pore space. A lower minimum matric suction value means that air can more easily enter (exit) pores. Water has difficulty entering pore spaces now relative to air, and as a result the unsaturated hydraulic conductivity, K of the soil decreases. A smaller K value implies slower and delayed downward flow to the water table. The downward flow is dampened with time as it slowly moves downward, resulting in smoother water table fluctuations (Figure 5a).

4.2.3
Pore-size distribution Changes to the pore-size distribution parameter, n, result in higher correlation between the field observations and model results. In the van Genuchten soil hydraulic functions used to determine θ (Equation 2), both α and ψ are raised to the power of n. The n parameter is further used to determine the empirical parameter, m (1) (Equation 4), where m is an exponent used to solve for the unsaturated hydraulic conductivity and the soil water retention (Equations 2 and 3). Figure 6 shows that a decrease in n causes an exaggerated response, whereas an increase in n causes a dampened response. Physically, n represents the allowed abundance of varying pore sizes in a volume of soil (Nimmo, 2004). When water infiltrates a soil with a narrow and uniform distribution of pore sizes, the water flux can more easily fill (or empty) pores at the same matric suction. In the subsurface, matric suction is defined as the difference between pore air pressure and pore water pressure. Conventionally, pore air pressure is equal to atmospheric pressure and is ignored (Chiorean, 2017). Expanding the allowed distribution of pore sizes (higher n) increases the possible variation of pore sizes. A soil with more highly varying pore sizes requires highly varying matric suction for water to fill (empty) the varying size pores, generally retarding downward flow (Nimmo, 2004;Zhang et al., 2019).

Vadose Zone Journal
In response, the water flux becomes dampened in the subsurface, visually translating to a smoother water table response ( Figure 6). In contrast, narrowing the allowed distribution of pore sizes (smaller n) reduces the possible variation of pore sizes. Reducing the possible variations in pore sizes allows water to fill (empty) pores more easily with the same matric suction. Such ease allows water to flow downward at a faster rate, resulting in a more exaggerated water table response (Figure 6a).

4.3
How does soil water storage respond to the EPE?
The modeled changes in soil water storage suggest new developments that affect our understanding of how the subsurface responded to the 2013 EPE. First, there was a rapid increase in soil water storage in late WY2013. After the EPE, soil water storage remained elevated through WY2014 and into WY2015. The early part of WY2014 had comparatively higher soil water storage during the winter months (December 2013-March 2014) relative to all other water years for the same time frame. The heightened soil water storage during this time frame may be a strong indicator that recharge occurred for several months after the EPE, especially because minimal evapotranspiration occurs in the area during the winter season.
The increase in V (Figure 7, dotted lines) post-EPE is sustained through WY2015. During WY2016 (February 2016), the soil water storage once again reaches a winter low as seen during WY2012 (February 2012) and WY2013 (February 2013), pre-EPE. The modeled changes in soil water storage suggest that a 2-yr recovery occurred in response to the EPE-induced infiltration flux.

CONCLUSIONS
The 2013 Colorado EPE not only flooded the surface and rivers downstream but resulted in rapid infiltration and heightened water table response. Here are the conclusions we draw from this study: 1. Both the field observations and model results showed a water table rise after the EPE, which persisted for ∼18 mo before the water table recovered to pre-EPE levels. 2. Average annual soil water storage post-EPE for WY2014 and WY2015 was higher than all other water years in the record, indicating a wetter subsurface post-EPE. 3. The post-EPE could serve as a hydrologic buffer that stores a portion of extreme precipitation for various seasons. 4. A sensitivity study of model parameters showed that the modeled water table was most sensitive to changes in the empirical parameter that represents the pore size distribution value, n. Pore-size distribution cannot be measured in the field, it is essential to scrutinize the values to which empirical parameters are set in simulations.
Given the characteristics in geology, hydrology, and geographic considered herein, the model setup and results could be applicable to regions of similar characteristics. By assessing the potential for unsaturated zone profiles to serve as natural storage space for EPE-induced infiltration, this study could provide a scientific basis for water managers to timely use the stored water that may be released to streams over time. More research regarding local subsurface response to EPEs is needed because EPEs are predicted to occur more frequently worldwide (Lehmann et al., 2015;Trenberth, 2011;Wasko et al., 2016). Understanding the effects of individual EPEs on the subsurface could also provide the basis for predicting aggregated effects over longer time scales. From another viewpoint, in headwater regions, snowmelt is the primary source of groundwater recharge. Although the rate of snowmelt is typically not as dramatic as EPEs, snow could occur at an accelerated rate under warming (Pepin et al., 2015). This study could be informative for projecting the potential hydrologic consequences of accelerated snowmelt. More broadly, the results of this study contribute to a better understanding of how the subsurface can act as a long-term hydrologic buffer for infiltration from an EPE before recharge occurs.

A C K N O W L E D G M E N T S
This research was supported by the National Science Foundation, Division of Earth Sciences EAR (Grant 1834290). The authors appreciate the constructive comments from three anonymous reviewers and the associate editors.

C O N F L I C T O F I N T E R E S T
The authors declare no conflict of interest.