In August 2016, the National Weather Service Office of Water Prediction (NWS/OWP) of the National Oceanic and Atmospheric Administration (NOAA) implemented the operational National Water Model (NWM) to simulate and forecast streamflow, soil moisture, and other model states throughout the contiguous United States. Based on the architecture of the WRF-Hydro hydrologic model, the NWM does not currently resolve channel infiltration, an important component of the water balance of the semiarid western United States. Here, we demonstrate the benefit of implementing a conceptual channel infiltration function (from the KINEROS2 semidistributed hydrologic model) into the WRF-Hydro model architecture, configured as NWM v1.1. After calibration, the updated WRF-Hydro model exhibits reduced streamflow errors for the Walnut Gulch Experimental Watershed (WGEW) and the Babocomari River in southeast Arizona. Model calibration was performed using NLDAS-2 atmospheric forcing, available from the NOAA National Centers for Environmental Prediction (NCEP), paired with precipitation forcing from NLDAS-2, NCEP Stage IV, or local gauge precipitation. Including channel infiltration within WRF-Hydro results in a physically realistic hydrologic response in the WGEW, when the model is forced with high-resolution, gauge-based precipitation in lieu of a national product. The value of accounting for channel loss is also demonstrated in the Babocomari basin, where the drainage area is greater and the cumulative effect of channel infiltration is more important. Accounting for channel infiltration loss thus improves the streamflow behavior simulated by the calibrated model and reduces evapotranspiration bias when gauge precipitation is used as forcing. However, calibration also results in increased high soil moisture bias, which is likely due to underlying limitations of the NWM structure and calibration methodology.
From 1984 to 2013, flood events in the United States cost the nation $7.95 billion per year and resulted in an average of 85 fatalities per year (National Weather Service 2014), making flooding one of the deadliest and costliest hazards for the United States. Flash floods occur when atmospheric conditions are favorable for the occurrence of sustained (often convective) heavy precipitation, where sufficient ascent and moisture are both present (e.g., Doswell et al. 1996). Flash flood events are difficult to predict in semiarid environments, such as the southwestern United States, due to the difficulty of having good estimates of antecedent soil moisture, and due to the relatively short durations and highly localized character of extreme precipitation events. Initial land surface soil moisture conditions can be highly uncertain due to limited soil moisture observations, while elevated terrain hampers accurate estimation of precipitation using weather radar (e.g., Zamora et al. 2014). Therefore, flash flood prediction in these semiarid domains is particularly challenging and is further complicated by the dry antecedent conditions of ephemeral channel beds, resulting in large transmission losses (Goodrich et al. 1997). Previous efforts to improve the physical process representation in hydrologic forecasts within these environments have focused on the implementation of spatially distributed hydrologic models, making use of gridded high-resolution precipitation data (e.g., Blöschl et al. 2008; Looper and Vieux 2012; Broxton et al. 2014; Hardy et al. 2016; Reed et al. 2007; Gourley et al. 2017). In the present study, the added value of accounting for channel infiltration in a spatially distributed hydrologic model is evaluated for two basins in the southwest United States.
To ameliorate the uncertainties associated with hydrologic forecasts, and to provide longer forecast lead times for flash flood events, the National Weather Service (NWS) Office of Water Prediction (OWP), of the National Oceanic and Atmospheric Administration (NOAA), recently implemented the operational National Water Model (NWM), based on architecture of the Weather Research and Forecasting (WRF) WRF-Hydro hydrologic model (Gochis et al. 2015). The NWM is a continental-scale distributed hydrologic model that produces streamflow forecasts for 2.7 million stream reaches across the contiguous United States (CONUS) based on observed and forecasted precipitation. For the semiarid southwestern United States, the model is subject to significant errors (Dugger et al. 2017) and, accordingly, the present study focuses on improving the performance of the WRF-Hydro model (configured as NWM v1.1) for this region.
In the southwest United States, much of the streamflow that occurs in ephemeral channels is due to surface runoff associated with convective rainfall occurring during the North American monsoon season (NAM) (e.g., Maddox et al. 1995; McCollum et al. 1995; Adams and Comrie 1997). During the NAM, precipitation can either be phase-locked to the high terrain or propagate into the low deserts as squall lines or mesoscale convective systems (MCSs; e.g., Pytlak et al. 2005; Bieda et al. 2009; Finch and Johnson 2010; Newman and Johnson 2012; Seastrand et al. 2015; Lahmers et al. 2016; Luong et al. 2017), resulting in significant amounts of precipitation and surface runoff.
Throughout southeastern Arizona, the groundwater levels tend to be deep, so that the surface runoff that flows through the channel network infiltrates into the soil, where it becomes a source for groundwater recharge (e.g., Blasch et al. 2004). Goodrich et al. (2004) note that the depth to groundwater in the U.S. Department of Agriculture (USDA) Agricultural Research Service (USDA-ARS) Walnut Gulch Experimental Watershed (WGEW; Moran et al. 2008), a tributary to the San Pedro River in southeast Arizona with no perennial channels, ranges from approximately 50 to 145 m in the lower and middle basin, respectively. However, local recharge from channel infiltration can make up a significant portion of the local water balance, as recharge estimates (based on groundwater modeling, surface water balance measurements, and evaluation of deep groundwater storage using microgravity measurements) suggest that channel infiltration accounted for 15%–40% of groundwater recharge in the San Pedro basin during wetter than average monsoon seasons (Goodrich et al. 2003, 2004). Virtually all of the runoff generated in the WGEW is the result of monsoon convective precipitation (Goodrich et al. 1997; Stone et al. 2008).
As currently implemented, WRF-Hydro does not permit water to exit the channel network once it enters a channel from the terrain routing grid (Gochis et al. 2015). Consequently, the operational version of the NWM does not account for channel infiltration losses. Because of the importance of channel infiltration processes in semiarid regions like the southwest, we modified the WRF-Hydro structure to include channel infiltration (e.g., Goodrich et al. 2004). Specifically, we used a scheme similar to the one implemented in the KINEROS2 distributed watershed model (Goodrich et al. 2012), which was originally developed for the southwestern United States and has been widely used for hydrologic forecasting in semiarid domains. KINEROS2 has been tested and used extensively (e.g., Goodrich et al. 2012; Yatheendradas et al. 2008) in the WGEW and other basins in southern Arizona.
Previous work has attempted to quantify and estimate the nature and magnitude of channel infiltration at various scales. Lane (1983) used an ordinary differential equation to approximate the rate of change in runoff volume with channel flow distance; observed inflow–outflow discharge from ephemeral reaches were then used to compute equation coefficients via regression analysis. Noorduijn et al. (2014) modeled infiltration from an artificial channel using Philip’s equation, accounting for both the gravity and pressure terms in the sorption equation, while Callegary et al. (2007) quantified the recharge potential of ephemeral channel reaches in the southwest CONUS using field data. Redistribution of water across the land surface caused by infiltration may also affect the local water balance; Zampieri et al. (2012) showed, using a modified version of the Community Land Model (CLM), that parameterization of the redistribution of water through infiltration out of the channel network was needed for the model to reproduce Oklahoma Mesonet soil moisture observations.
In addition to incorporating a representation of channel infiltration, the present study investigates the model calibration problem; this is important due to the inherent uncertainty in the model parameter estimates and their impact on the simulated hydrologic response (e.g., Yatheendradas et al. 2008). For spatially distributed hydrologic models (potentially having tunable parameters at each grid point), calibration is complicated by the high dimensionality of the model (e.g., Smith and Gupta 2012). This problem can be ameliorated using spatial regularization (Gupta et al. 2008, 2009; Samaniego et al. 2010), wherein model parameters are associated with a priori nonlinear functions of observed surface data. Doing so results in a smaller number of hyperparameters (i.e., parameters controlling the degree of modification to the actual model parameters through transfer functions) that must be calibrated to ensure that the original model parameters are physically consistent with catchment properties (e.g., Gupta et al. 2008, 2009). For example, Pokhrel et al. (2012) used a spatial regularization approach to calibrate the NWS Hydrology Laboratory Research Distributed Hydrological Model (HL-RDHM), and Vergara et al. (2016) used spatial climate and land data to derive a priori estimates of routing parameters for the kinematic wave routing model across the CONUS to be used for flash flood forecasts. Advanced spatial regularization techniques that have been proposed include the multiscale parameter regionalization (MPR) method, which uses spatial transfer functions to compute model parameters based on high-resolution spatial data that account for spatial heterogeneity across grid cells (e.g., Samaniego et al. 2010).
In this study, we followed the approach presented by Pokhrel et al. (2012) and calibrated a select set of WRF-Hydro parameters (determined by sensitivity analysis). This methodology is consistent with the calibration methodology developed for the NWM v1.1 by the National Center for Atmospheric Research (NCAR).
We evaluate the combined impacts of accounting for channel infiltration and parameter calibration on the simulated hydrologic response of WRF-Hydro for the WGEW and the Babocomari River basin, which are both tributaries to the San Pedro River. This work was focused primarily on the WGEW (149 km2), where relatively accurate precipitation forcing data are available from a network of 88 weighing rain gauges. The same calibration approach was extended to the larger Babocomari River basin, where the impacts of channel infiltration are greater due to its large area, despite precipitation forcing uncertainties. The NWM model structure and datasets used are described in section 2 of this paper, and the calibration approach is outlined in section 3. Calibration results are presented in section 4. Results from the model calibration and evaluation and possible future changes to the NWM are further discussed in section 5, and conclusions are included in section 6.
2. WRF-Hydro model and user modifications
a. NWM WRF-Hydro model structure and configuration
WRF-Hydro (Gochis et al. 2015) is a parallelized distributed hydrologic model that can either be forced offline using prescribed atmospheric forcing variables, or coupled to the Advanced Research version of the WRF (WRF-ARW) atmospheric model (Skamarock et al. 2008). A simplified schematic of WRF-Hydro is shown in Fig. 1. Atmospheric forcing data needed to execute WRF-Hydro offline include incoming shortwave radiation, incoming longwave radiation, specific humidity, air temperature, surface pressure, and near surface wind (both u and υ components). The NWM is a particular configuration of WRF-Hydro, which uses the Noah-MP land surface model (LSM; Niu et al. 2011) to resolve vertical fluxes within the soil column and exchanges with the atmosphere. Noah-MP is configured using gridded NWM soil parameters, which govern the drainage of water through the soil column, with 1-km grid resolution. NWM WRF-Hydro resolves horizontal surface and subsurface fluxes on a 250-m grid resolution routing grid (with 10-s and 60-min time steps, respectively). Since the routing grid cells are 4 times smaller than the Noah-MP grid cells, spatially varying quantities on the 250-m routing grid are aggregated back to the 1-km grid during model time steps when Noah-MP is called (every 60 min) and disaggregated back to the 250-m routing grid.
The subsurface flow module on the WRF-Hydro grid computes changes to the water table in the 250-m soil column (which is assumed to be 2 m deep everywhere), using Dupuit–Forcheimer assumptions. This assumes that the hydraulic gradient is based on differences in the groundwater table depth along the steepest gradient in eight possible directions around a routing grid point (Gochis et al. 2015). If subsurface flow causes a model grid point to become saturated, exfiltration is computed and this resultant water ponding is combined with infiltration excess and routed as surface runoff. Surface flow is computed using diffusive wave routing based on the steepest gradient around each grid point (Julien et al. 1995; Ogden 1997). Details of the surface and subsurface routing schemes of WRF-Hydro are discussed in detail in Gochis et al. (2015).
When surface flow reaches a grid cell that is designated as a channel, it is mapped to the vector channel network and routed downstream (with a 5-min time step). The NWM channel network is based on the National Hydrography Dataset (NHD) Plus Version 2 (NHDPlusV2) (McKay et al. 2012). In the present study, we added a channel infiltration parameter that is physically representative of the channel bed conductivity (m s−1) (ChannK; Table 1). Flow in the channels is computed using an iterative Muskingum–Cunge function for each reach. This vector routing scheme is more computationally efficient than other WRF-Hydro channel routing configurations (e.g., 1D diffusive wave routing) and has the capability to be mapped to specific rivers and reaches that may be of interest to emergency managers and stakeholders; however, it cannot resolve backwater flow (Gochis et al. 2015).
Deep baseflow (below the soil column) in WRF-Hydro is computed using a conceptual exponential bucket model. All water that drains out of the Noah-MP LSM 2-m soil column is mapped to a groundwater catchment, which corresponds to the NHDPlusV2 channel reach/catchment topology. There are a few relatively small reaches in the NWM domain that are not mapped directly to an NHD bucket. Water from this conceptual groundwater bucket is gradually returned to the channel reach that directly corresponds to its underlying catchment, where the rate of discharge is controlled by three tuning parameters in an exponential equation.
WRF-Hydro has parameters that may be specified through input tables and grids, and these parameters may be adjusted or calibrated depending on the study region and hydrologic behavior of interest. Priority Noah-MP and WRF-Hydro parameters, that were selected by the NCAR to regionally calibrate NWM version 1.1 (Dugger et al. 2017), are shown in Table 1. For all modeling described herein, model code and a priori parameters are based on NWM v1.1, except where modifications are noted. A priori Noah-MP soil parameter values in WRF-Hydro are based on Rawls et al. (1982), from the STATSGO2 soil texture dataset (available at http://websoilsurvey.nrcs.usda.gov/). Channel parameters, which correspond to individual reaches based on stream order and other channel characteristics (with the exception of the channel infiltration parameter; see next section) and that were evaluated after the first round of calibration in the WGEW (see section 3a), are shown in Table 1.
b. Model domains
In the present study, the effects of modeled channel infiltration are evaluated for the WGEW and the Babocomari River basin. Figure 2 shows the elevation grid for the 250-m NWM routing grid for the entire San Pedro basin; however, the simulations performed in this study use cutout grids for the Babocomari basin and WGEW.
For all calibration simulations, atmospheric forcing is derived from NLDAS-2 atmospheric variables that are regridded to the model domains; however, the model is calibrated with both regridded 1/8° NLDAS-2 precipitation and NCEP 4-km Stage IV precipitation, which is based on WSR-88D radar and gauge precipitation (Lin and Mitchell 2005). To reduce the influence of possible forcing uncertainty and bias, WRF-Hydro in WGEW was forced with gauge-based precipitation (88 gauges in the 149 km2 watershed area) interpolated to the WGEW-domain model grid. WGEW is useful for analyzing the performance of WRF-Hydro because it includes 20 soil moisture measurement sites, two AmeriFlux (Department of Energy 2018) flux towers, and subwatersheds with supercritical flumes for measuring runoff (Fig. 3).
Figure 4 demonstrates the profound role of channel infiltration in the lowermost reach of the WGEW, where streamflow during two runoff events from the two upstream flumes that come together to form a single channel (runoff gauges 2 and 7), is highly attenuated by the time it reaches runoff flume 1, situated 7 km downstream, defining the basin outlet. Model topography and two initial soil parameters in WGEW, which are a function of soil type, are shown in Fig. 3. This figure shows that soil saturated hydraulic conductivity based on NWM v1.1 grids in WGEW is relatively homogeneous, as it is derived from STATSGO2 soil texture class which is homogenous in this domain. Higher-resolution SSURGO soil texture data (available at https://websoilsurvey.nrcs.usda.gov/) in WGEW shows considerably more variability, which is considered in our discussion (section 5); however, this study uses the NWM STATSGO2 grids to stay consistent with the NWM v1.1 configuration. DKSAT, according to the NWM a priori parameters, is constant at 3.37 × 10−6 m s−1 everywhere in WGEW, except in the extreme upper basin. Vegetation across WGEW includes primarily desert shrubs and grasslands; however, the town of Tombstone, Arizona, is partially urbanized. WGEW has oak forests in its extreme upper reaches.
The Babocomari River basin characteristics are shown in Fig. 5. This basin has a larger drainage area, and much of the headwaters are at high elevations. Like WGEW, NWM STATSGO soil characteristics likely underestimate the variability of the soils, per the SSURGO dataset. The Babocomari basin has a flux tower and three soil moisture observation sites with reliable data through the calibration period (Fig. 5), so model states at these sites are also evaluated.
c. Baseflow bucket parameter modifications
The current WRF-Hydro conceptual bucket model, which assumes one-way direct connection between attenuated baseflow from a groundwater basin and the overlaying channel, provides a poor representation of baseflow in event-driven, ephemeral channels in semiarid environments. Depth to groundwater is often substantial in the southwest CONUS, such that surface water–groundwater processes can become decoupled. In addition, water from channels often infiltrates to recharge the local aquifer (e.g., Blasch et al. 2004), and this source of recharge is currently not represented in the NWM.
To prevent unrealistic simulations of baseflow in the channel network, we disabled the baseflow bucket model everywhere in the model domain by setting the Noah-MP bottom drainage scaling (SLOPE) parameter to zero, imposing a no-flow condition at the bottom of the Noah-MP LSM in both the WGEW (Fig. 3b) and the Babocomari basin (Fig. 5b). This assumption is consistent with NWM calibrated parameters for this region, where a small SLOPE parameter effectively eliminates deep baseflow to the channel. Neither WGEW nor the Babocomari basins contain perennial channels fed from deep groundwater, and there are no perennial channels shown in either basin per the NHDPlusV2 dataset attributes, so they are not good candidates for applying the conceptual baseflow bucket module in its current form.
d. Channel infiltration loss function
The WRF-Hydro channel routing scheme assumes a trapezoidal channel geometry, and the length and slope of specific reaches is specified in the NHDPlusV2 dataset. A cross section of a WRF-Hydro channel is shown in Fig. 1. If the volume of water in a reach is known, the height of the water h (m) and the wetted perimeter p (m) can be calculated by the following argument. First, the cross-sectional area a (m2) of the water may be computed by dividing the channel volume by the length of the reach, assuming a constant height along the entire length of the reach. Based on the trapezoidal shape assumption, cross-sectional area is equivalent to
We can compute ws (m) as a function of the riverbank slope s (h/ws in Fig. 1) and water height. Note that s is equivalent to the WRF-Hydro ChSlp parameter. This may be written as
Since area can be derived for a channel volume, for a given bottom width w and bank slope s, Eq. (2) can be solved for h. It is possible to derive the wetted perimeter p from h using
Due to the effects of uneven channel bed topography, wetted perimeter will be reduced during periods of low flow, as water will only flow within the lowest portions of an irregular channel. To account for this tendency in trapezoidal channels, p is reduced during periods of low flow using a conceptual model similar to what is used in KINEROS2 (Woolhiser et al. 1990). KINEROS2 computes a corrected wetted perimeter pe for a channel using the function
where b (unitless) is set to a constant value of 0.15. Channel infiltration I (m3 s−1) may then be derived as
In this function, k (ms−1) is the saturated conductivity of the channel bed (the ChannK parameter), and l (m) is the length of a channel reach. Channel infiltration is accounted for in the Muskingum–Cunge routing scheme of WRF-Hydro, as an added sink in the iterative calculation, by deriving the effective flow height and wetted perimeter for a channel reach volume and assigning a saturated conductivity for the channel bed below. We assume b to be constant, while the saturated conductivity of the channel bed (ChannK) is initially set equivalent to the saturated soil conductivity of Noah-MP (DKSAT). Water that infiltrates out of the channels is assumed to contribute to deep groundwater recharge and is therefore removed from the model. During calibration, initial ChannK is adjusted by a scalar multiplier, a simple form of spatial regularization (e.g., Pokhrel et al. 2012).
3. WRF-Hydro calibration intervals
We calibrated the NWM with three different configurations of increasing complexity to evaluate: 1) the importance of channel parameters following the addition of the channel infiltration function, 2) the added value of the channel infiltration function in an ephemeral catchment with a relatively simple hydrologic response (i.e., WGEW), and 3) the calibrated model performance with two nationally available forcing products (Table 2).
Optimization of daily streamflow was performed using the dynamically dimensioned search (DDS) algorithm (Tolson and Shoemaker 2007), which is capable of converging to near optimal parameter sets with fewer (approximately 100–500) iterations (e.g., Lespinas et al. 2017) than the widely used Shuffled Complex Evolution function (e.g., Duan et al. 1992) that can require ~10 000 iterations to converge to an optimal solution. This makes DDS better suited for calibration of computationally expensive, distributed, physically based, models like the WRF-Hydro model. For all calibration exercises, except the first (section 3a), we use 500 iterations of DDS with the updated NWM v1.1 utilizing channel infiltration. Lespinas et al. (2017) show that improvements to model skill are greatest between 100 and 500 DDS iterations and that less value is added beyond 500. As the first exercise (section 3a) is simply intended to bring the NWM parameters to an acceptable state so that channel parameters can be evaluated, we only performed 250 iterations for this step. For all calibration simulations, NWM v1.1 with channel infiltration is calibrated to daily streamflow to stay consistent with NWM v1.1 calibration methods, and we consider hourly streamflow when evaluating the NWM in WGEW since most streamflow events occur at the hourly time scale.
All of the calibrations reported here are based on optimization of the Kling–Gupta efficiency (KGE) performance metric (Gupta et al. 2009), which equally weights correlation, water balance, and variance errors. Like the Nash–Sutcliffe efficiency (NSE), KGE is optimal when equal to 1, and negative values of KGE are considered to have low skill. KGE is optimal when the Euclidian distance from an ideal point for the ratio of modeled to observed standard deviation α, ratio of modeled to observed mean β, and correlation coefficient r are minimized (Gupta et al. 2009):
For all basins, initial model states were derived by executing WRF-Hydro with default parameters with their calibration precipitation forcing (WGEW-gauge, Stage IV, or NLDAS-2) for an 8-yr spinup period during water years (WY) 2007–15. The model state at the end of this period was used as the initial “warm” state for calibration, consistent with the practices of the NCAR WRF-Hydro NWM development team. This ensures a long-term model spinup of multiple years to allow state variables to reach equilibrium. Then for each calibration iteration, model parameters were given one year of additional spinup to equilibrate before the calibration and evaluation periods (described in detail below).
a. Stage 1: Channel parameter sensitivity
To evaluate the sensitivity of the channel routing parameters, WRF-Hydro was first calibrated to eliminate water balance errors using WGEW gauge precipitation forcing. This step was only performed for the WGEW basin. As bias is a component of KGE, this objective function permitted us to reduce water balance errors and simulate an otherwise realistic hydrologic response (since correlation and variance errors are also accounted for). To select parameters for this calibration, sensitivity analysis on daily KGE was first performed on the updated NWM, using linear adjustments applied to default prior estimates (e.g., DKSAT shown in Fig. 3c) of the model parameters (Table 1). Selected calibration parameters included ChannK (channel conductivity), DKSAT (soil conductivity), REFKDT (infiltration scaling), and SMCMAX (soil porosity). These parameters were chosen because sensitivity analysis revealed that they had the greatest impact on model KGE. ChannK, SMCMAX, and DKSAT were computed by multiplying initial NWM parameters by a constant, and REFKDT was assumed constant in the whole model domain. DKSAT was adjusted by an addition constant (before being adjusted by a multiplicative constant).
For this first step of calibration, the Noah-MP time step was reduced from 60 min (the standard for the operational NWM) to 15 min, so that the model produced outputs, including streamflow, at 15-min temporal resolution. This permitted us to critically analyze the sensitivity of the channel infiltration function at fine resolution (not shown). This calibration improved WGEW model KGE (with daily temporal resolution) from −1.57 to 0.80 for WY2009–11 (the calibration period; spinup period, WY2008) by reducing water balance errors that caused excessive runoff and low ET biases (not shown).
Following initial calibration, a simple one-at-a-time parameter sensitivity analysis was performed on the WRF-Hydro channel parameters, including BtmWdth (channel width), ChSlp (outer channel slope), and N (Manning’s N; see Table 1). For each parameter, initial parameter values in the WGEW domain were adjusted by either addition or multiplication to the initial parameter with a constant. A priori parameter values were assigned by the NCAR NWM WRF-Hydro development team. This exercise permitted evaluation of the model hydrologic response to perturbations of the channel routing parameters. The main evaluation metric of interest for this analysis was the hourly correlation coefficient because its accuracy is a proxy for errors in the timing of streamflow events (e.g., Gupta et al. 2009). This analysis showed that ChSlp and N, which affect the transit time of streamflow, impact the model correlation coefficient (Table 3), with correlation coefficients increasing when these parameters were multiplied by 2.5 and 1.5, respectively. However, in further analyses, we chose to only vary ChSlp, while keeping Manning’s channel roughness N constant. We could have also varied this latter parameter, as both parameters (ChSlp and N) affect the model correlation coefficient (Table 3). However, as both parameters have a similar impact (though, not necessarily identical) on NWM streamflow transit time, it was decided that including an additional parameter to calibrate would have been somewhat redundant.
b. Stage 2: Evaluation of channel infiltration in WGEW
To demonstrate the added value of channel infiltration, NWM WRF-Hydro was recalibrated for WGEW with the same parameters used for Stage 1 of calibration (see last section) but with the addition of the ChSlp parameter, for WY2011–13 (WY2010 for spinup) with channel infiltration either active or disabled. The change to the calibration period from Stage 1 was intended to include both wet and dry years in the calibration. The ChannK parameter multiplication constant was constrained between 0.0 and 1.0, based on the assumption that saturated conductivity (DKSAT) sets an upper limit for channel infiltration. This implicitly assumes that the channel bed properties for a region will tend to be consistent with local soil characteristics, with soil conductivity being maximized when saturated. The Noah-MP time step for this calibration was set to 60 min, consistent with the operational NWM v1.1.
All adjustment constants for selected parameters and their ranges are shown in Table 2. For parameters where a priori values are adjusted by multiplication or addition constants (e.g., DKSAT), the ranges of the adjustment constants are shown. Note that these adjustment factors are not the parameters themselves. For example, if SMCMAX has an adjustment factor of 0.8, the a priori values of SMCMAX will be multiplied by 0.8 everywhere to compute the new parameters, thus preserving the spatial patterns of the parameters as they are calibrated. REFKDT is left constant throughout the model domain, so the possible ranges of the actual values for this parameter are shown in Table 2. The final parameter adjustment factors and calibration configuration for these simulations are shown in Table 4. As noted above, 500 iterations of DDS were used, and daily streamflow KGE was the objective function. To evaluate whether the calibration of the NWM improved the entire physical model state and not just simulated streamflow, soil moisture and evapotranspiration (ET) fluxes within WGEW were evaluated after calibration, with and without channel infiltration.
c. Stage 3: Calibration with Stage IV and NLDAS-2 precipitation
To demonstrate the potential application of the aforementioned channel infiltration function and calibration methods in the operational NWM, WRF-Hydro was also calibrated using NLDAS-2 and NCEP Stage IV precipitation as forcing instead. Calibration periods were 3 years long, with one extra prior year for spinup (i.e., WY2011–13 with spinup WY2010 and WY2009–11 with spinup WY2008 for the WGEW and Babocomari basin, respectively) and selected to include both wet and dry years for each basin (Table 2). To demonstrate the effects of channel infiltration, the model was calibrated with and without channel loss.
Following calibration, WRF-Hydro, configured as NWM v1.1 with channel loss, was run with both default uncalibrated parameters (as a benchmark), calibrated parameters without channel infiltration, and calibrated parameters with channel infiltration for WY2009–16. WY2008 was used as spinup, and 3 years of this 8-yr period were calibration years. The results from this evaluation are presented herein.
a. Calibrated model performance and added value of channel infiltration in WGEW
In WGEW, using rain gauge observations, calibration eliminated water balance errors, largely by eliminating spurious flashy peaks in the uncalibrated simulation, and performed better when channel infiltration was included, yielding slightly higher correlation coefficients (Fig. 6). The resulting cumulative simulated streamflow for WY2009–16 is less biased, demonstrating the added value of calibration with channel loss for reducing water balance errors. Adding channel loss reduced the negative percent bias of the modeled versus observed coefficient of variation (labeled as “cvdiff” in Fig. 6), suggesting that the model with channel loss had less tendency to underestimate the variance of the streamflow. To quantify the amount of channel infiltration in the calibrated model, we executed the NWM with the same parameters as the calibrated solution with channel infiltration, but with channel infiltration disabled. Adding channel infiltration reduces the output at the outlet of WGEW by 21.73% over the model with channel infiltration disabled (see figure in the online supplemental material), indicating the fraction of water that becomes channel infiltration. Table 5 shows that when WRF-Hydro is calibrated with channel infiltration, the model KGE is 0.9431, and the correlation coefficient is 0.9442 for daily streamflow during the calibration period. KGE was reduced to 0.8679 when the model was calibrated without channel loss. These values show slightly less skill when zero values of modeled and observed streamflow are omitted from the analysis (see supplemental material). As would be expected, model skill decreases in the evaluation period (WY2009–10 and WY2014–16), with a notably lower KGE (0.4794 versus 0.6596) and higher bias (40.85% versus 7.69%) when the model was calibrated without channel loss, versus with channel loss (Table 5). For hourly resolutions, the calibrated model results with channel loss had a KGE of 0.83 (0.55) and a correlation coefficient of 0.83 (0.80) for the calibration (evaluation) periods (not shown here). For two sample hydrographs, the bottom panels of Fig. 6 show that calibration improves modeled hydrograph shape even though the model struggles to represent specific hourly events. The calibrated model without channel loss produces slightly more intense spurious runoff events than the model with loss.
b. Soil moisture and flux tower ET evaluation in WGEW
To go beyond runoff only and consider the physical process representation of WRF-Hydro, model soil moisture and ET were evaluated against observations in WGEW starting in WY2009 and ending in April 2015. WGEW 5-cm soil moisture observations were compared to the area average of Noah-MP 0–10-cm soil moisture, based on 20 soil moisture sites in WGEW including the two AmeriFlux flux tower sites (i.e., Kendall Grassland and Lucky Hills). At these flux tower sites, observed hourly soil moisture is compared to model soil moisture averages from the Noah-MP grid point closest to the site and the two points surrounding it in all directions (25 points total), to minimize the impact of spurious model simulations for a single grid point. Calibration increased the positive bias of near-surface soil moisture, averaged throughout the basin; however, including channel infiltration in the calibration did not increase the soil moisture bias by as much (Table 6). Calibrating without channel loss nominally increased correlation coefficients for the basin average and two sites, and adding channel infiltration had little effect on the correlation coefficients after calibration. Soil moisture from WY2014, which had a wetter than average NAM season, is shown in Fig. 7. This figure also indicates that the modeled soil moisture solutions have a greater bias after calibration. We computed KGE for model soil moisture, which also reflects the increased bias after calibration more than the nominal gains in correlation coefficient. These results show that the calibrated WRF-Hydro model (with and without channel infiltration) increases the total amount of water inside the upper soil layer by 4.8% (computed basin average soil moisture in control vs calibration with infiltration, with WGEW gauge forcing). This increase in soil moisture in desert catchments is not observed in the observation datasets. All modeled solutions have a positive bias and slow dry-down during dry periods, which might be from limiting the bottom drainage scaling factor to eliminate deep groundwater baseflow. Increased values of REFKDT from calibration in all the basins likely also increased soil moisture, by increasing infiltration. Analysis of 15-cm soil moisture data (not shown) at Kendall Grassland and Lucky Hills reveals that soil moisture is less biased compared to observations, which may suggest that the Noah-MP 0–10-cm layer might be better representative of deeper soils. We acknowledge that representativeness errors in the soil moisture observation sites may also be possible.
Benefits of calibration for the water balance can be identified when considering observed daily ET fluxes from the Kendall Grassland and Lucky Hills AmeriFlux flux tower sites from WY2009 to December 2015 (Table 6), compared to averages of corresponding model grid points (same as for soil moisture). WRF-Hydro, without calibration has a negative ET bias at Lucky Hills and Kendall Grassland. This negative ET bias is consistent with the uncalibrated model’s tendency to overestimate streamflow, indicating that too little precipitation infiltrates at the land surface (due to low REFKDT and DKSAT values), forcing too much water to flow into the channel network as overland flow. Calibration of WRF-Hydro with channel infiltration, which eliminates the positive streamflow bias, further reduces the magnitude of the ET bias at both Kendall Grassland and Lucky Hills (compared to the calibrated simulation without channel loss). This implies that when WRF-Hydro is executed with accurate forcing precipitation, calibration can improve the accuracy of the simulated streamflow and reduce the ET bias, despite increasing soil moisture bias, to improve the water balance. Figure 7 shows that the high soil moisture bias might be associated with too much ET during the dry season, while ET is less during the wet season. This tendency may be related to the soil moisture bias discussed above.
c. Extension of calibration methods with NLDAS-2 and Stage IV precipitation
WRF-Hydro was calibrated with NLDAS-2 and Stage IV precipitation in WGEW, and streamflow bias errors were consistently reduced with calibration. However, correlation coefficients demonstrate little improvement outside the calibration period regardless of whether channel loss is added (Fig. 8, Table 5). This is true with both Stage IV and NLDAS-2 forcing. NLDAS-2 forced results are shown in the supplemental material. Bias also increases outside of the calibration period. Adding channel infiltration does reduce coefficient of variation errors and produces a realistic hydrograph, despite timing errors. The fact that correlation coefficients and bias for the same study area are degraded when Stage IV or NLDAS-2 precipitation is used suggests that precipitation forcing may be a source of uncertainty for hydrologic models in this region. Table 7 implies that soil moisture bias errors persist regardless of precipitation forcing, but correlation coefficients remain high regardless of calibration. As with the gauge forcing, KGE values seem to reflect bias errors. Model ET (at Kendall Grassland and Lucky Hills) has a positive bias, which increases with calibration with both datasets (Table 7), implying that calibration might be compensating for precipitation biases and the distribution of Stage IV precipitation events (see discussion). Given the low runoff/rainfall ratios common in the semiarid Southwest, small errors in rainfall forcing typically result in large errors in runoff predictions (Goodrich et al. 2012).
Calibrating WRF-Hydro without channel loss did improve the KGE and bias metrics in the Babocomari basin (regardless of forcing), but it did not yield a realistic hydrologic response (see Fig. 9 for Stage IV results and supplemental material for NLDAS-2 results). The same figure shows that calibrating with channel loss did produce a realistic hydrologic response. Correlation coefficients were low during the evaluation period (WY2012–16), outside the calibration period (WY2009–11), consistent with WGEW. Cumulative streamflow in Fig. 9 shows that the later years (WY2014–16) were wetter and had higher bias, suggesting that the hydrologic response may not be stationary (i.e., having constant statistical properties over time; Wilks 2006) over the calibration and evaluation periods. The poor hydrologic response in these later years may also be due to errors in the Stage IV precipitation. Despite the limitations of the NLDAS-2 and Stage IV precipitation products in the Babocomari basin, the added value of channel infiltration is clear, as calibrating the model without channel loss cannot completely eliminate the positive bias (Fig. 9, Table 5). The calibrated model without channel loss produces excess flow during dry periods, as the flashy runoff peaks are reduced by increasing infiltration into the soil column, leading to excessive baseflow entering the stream network, which may be associated with a low coefficient of variation bias (possibly implying more frequent flow with less variance).
Calibrating WRF-Hydro with Stage IV forcing with channel infiltration increases Noah-MP level 1 soil moisture bias at three soil moisture sites operated by the NOAA Hydrometeorology Testbed (HMT; Table 8) and has little effect on model soil moisture correlation coefficient (Fig. 10, Table 8). When channel infiltration is disabled, calibration improves soil moisture bias, but this is clearly at the expense of hydrologic response (see above). Observed ET at the Audubon Research Ranch AmeriFlux Site from WY2009 to December 2011 is also compared to model ET, which generally reflects the positive ET bias, also shown in WGEW with Stage IV precipitation (Fig. 10; Table 8). For deeper soil levels (20 and 50 cm), the model underestimates soil moisture variability and has a consistent high bias (see Fig. 11 and supplemental material) at all three HMT sites. This high bias decreases the soil moisture KGE. This is consistent with the soil moisture bias in the upper layers, and it suggests that too much water is entering the soil moisture column, to compensate for the runoff biases in the NWM. The filling of the soil columns is also likely due to the imposing of no-flow conditions at the base of Noah-MP to shut off the baseflow bucket model or due to other limits of Noah-MP parameters.
a. Calibration and precipitation forcing errors
Calibration with channel loss improves streamflow skill in the Babocomari basin with NCEP Stage IV and NLDAS-2 precipitation forcing, while calibrating the NWM without channel loss cannot produce a realistic hydrologic response in the same basin, unlike in WGEW. To compensate for the lack of channel infiltration, the optimal value for the REFKDT parameter increased in the Babocomari basin, to enable increased soil infiltration before the water reached the channels (Table 4). This relationship was less consistent in WGEW. Despite this compensating effect, the model without channel loss still had higher bias (Table 5), particularly in the evaluation period. Thus, it appears that parameter compensation can partially mask out the biases from channel infiltration in WGEW; however, this is not possible in a larger basin with more drainage area (and more area for water to enter the channels that must eventually infiltrate) like the Babocomari basin.
In the San Pedro basin, where streamflow in excess of small base flows is primarily driven by warm season convection, the calibrated model with NCEP Stage IV and NLDAS-2 forcing is associated with lower correlation coefficients, particularly during the evaluation period. As Stage IV remains the more reliable product (in the present study), this discussion centers on assessing the errors associated with only Stage IV precipitation. Figure 12 demonstrates the likely cause of the poor streamflow correlation, as observed precipitation from a NOAA HMT gauge at Elgin, Arizona, sometimes poorly matches Stage IV precipitation, despite the climatology of the two datasets being similar. Figure 12 shows that the evaluation period was wetter than the calibration period, which may partly explain the high model bias after 2013 in the WGEW and the Babocomari basins. Similar results from two other Babocomari basin rain gauges are shown in supplemental material.
One possible reason for this correlation error is that precipitation estimation by radar in the southwest United States is affected by elevated terrain (causing beam blockage and thus leading to the need for high measurement elevations; e.g., Zamora et al. 2014). Radar and available gauge observations are used to derive NCEP Stage IV precipitation (Lin and Mitchell 2005). Limitations in the quality of Stage IV data may cause the reduced streamflow correlation coefficients in the Babocomari River and WGEW. A detailed analysis of differences in WGEW gauge and WSR-88D rainfall data can be found in Morin et al. (2003, 2005).
Regridded 4-km Stage IV and 1/8° NLDAS2 forcing may also spatially spread precipitation over a larger area, therefore buffering precipitation over the landscape and reducing localized high-intensity events. As surface runoff only occurs when there is sufficient precipitation intensity to exceed the infiltration capacity of the soil, spreading of precipitation could reduce surface runoff that might otherwise occur over a small area associated with locally heavier precipitation intensity. This might explain why Stage IV forcing produces more precipitation and less streamflow over the WGEW domain (Fig. 12). Spatial spreading of precipitation, and subsequent reduction of surface runoff, may also be part of the cause of the high soil moisture and ET bias.
b. Suggestions for future work
The addition of channel infiltration and subsequent calibration of WRF-Hydro enables it to produce a more realistic hydrologic response with reduced water balance errors in two basins in the southwest United States. Calibration of WRF-Hydro in WGEW with gauge precipitation produces hourly streamflow with KGEs of 0.83. While this calibration does increase the positive bias of near surface soil moisture in WGEW, it does maintain the soil moisture correlation coefficients and reduces ET bias. Stage IV or NLDAS-2 forcing, which is derived by WSR-88D radar data and gauge precipitation or model and satellite remote sensing products, respectively, is subject to greater error and uncertainty than gauge precipitation, which reduces the skill of the model, irrespective of the channel infiltration function. Given the limitations of the NWM with channel infiltration and its calibration routines described in this work, we make several suggestions for future development of the NWM and other distributed hydrologic models used for forecasting.
Short-term solutions that may ameliorate the calibration challenges for the NWM presented in our results include:
New precipitation forcing datasets: NOAA is now developing the high-resolution (~1 km) Analysis of Record for Calibration (AORC) precipitation and temperature dataset. While testing is still ongoing, this dataset may add some value over the NLDAS-2 and NCEP Stage IV forcing products, and partially ameliorate some soil moisture and ET bias, described above.
Additional calibration parameters and metrics: The increased soil moisture bias, despite reduction in streamflow errors from calibration, may indicate potential deficiencies in the Noah-MP structure, including some hard-coding of parameters (e.g., Mendoza et al. 2015). Cuntz et al. (2016) showed that the water partitioning of Noah-MP was dependent upon several hard-coded parameters. Efforts are now ongoing to update the Noah-MP code to allow for more multicomponent calibration and adjustment of these hard-coded parameters, to improve its physical process representation.
The limitations of the NWM and its calibration approach may also be further addressed with more permanent based solutions in the long term, including:
Multivariable calibration: Addressing the soil moisture biases may additionally require calibrating the model to quantities other than streamflow, including soil moisture. Basins around the United States with soil moisture observations, including those maintained by NOAA HMT in Arizona and California (Zamora et al. 2011), would be good candidates for this. Calibrating WRF-Hydro to soil moisture may permit researchers to reduce soil moisture errors, thus making the model’s representation of streamflow and ET more physically consistent.
Land and soil datasets: More detailed soil datasets (e.g., SSURGO) could also be used to better constrain soil parameters, as Figs. 3 and 5 indicate that a priori NWM parameters likely underestimate soil heterogeneity compared to the SSURGO data shown in the same figures. SSURGO data, which are based on local soil surveys, were used to analyze output from the Noah LSM in Zamora et al. (2014). Use of a SSURGO dataset may also reduce the soil moisture errors and biases shown earlier in this study.
Alternatives to the baseflow bucket model: Soil moisture errors in the Babocomari basin in deeper layers may also have been caused by modifying the free drainage scaling parameter to essentially a no-flow boundary on the bottom of the Noah-MP LSM. The current baseflow bucket model’s use of NHD catchments keeps baseflow relatively local by returning it to the same channel. This is not realistic in basins with groundwater decoupled from surface water, so it was disabled for this work. A WRF-Hydro configuration that is coupled to a physically based groundwater model, such as ParFlow (e.g., Maxwell et al. 2015) could be a possible solution to this.
Further enhancements to channel infiltration: Another area for future analysis is to make the channel infiltration function more representative of physical processes, like suction from dry soil (Smith and Goodrich 2000; Smith et al. 2002). The infiltration scheme implemented in the NWM for the present study only accounts for the impact of gravity. Soil suction results in an increase of infiltration into the dry soil at the start of a runoff event (Smith and Parlange 1978; Parlange et al. 1982); however, its impact has not yet been considered. Accounting for this process in a future version of the model, may yield further improvements to WRF-Hydro channel infiltration. Another potential improvement would be if the channel infiltration function were coupled back into the Noah-MP LSM, permitting the model to resolve the effects of channel infiltration on the surface fluxes in riparian areas, even when water is not flowing. A saturated soil column could also limit channel infiltration during wetter periods.
The operational NWM is subject to systematic errors in the southwest CONUS, in part due to lack of a representation of ephemeral channel infiltration. These errors reduce the skill of the NWM as an operational tool in the southwest CONUS by the NWS. Addition of a channel infiltration function and subsequent calibration permits the model to produce a more realistic hydrologic response in low-elevation regions that are fed by convective rainfall. To fully realize the benefits of this infiltration scheme, future work is needed to eliminate the soil moisture biases of the Noah-MP LSM in WRF-Hydro. This might be accomplished through direct calibration to soil moisture, but may also require changes to the Noah-MP structure, and modification or elimination of the WRF-Hydro baseflow bucket scheme to permit deep groundwater recharge in arid regions. Spatial spreading of precipitation by the Stage IV precipitation product likely also contributes to soil moisture and ET biases.
Analysis of modeled streamflow and associated forcing precipitation data suggests that uncertainties associated with national precipitation datasets will need to be considered in future model calibration efforts. A logical step for evaluating model forcing uncertainty is to use the NOAA AORC forcing dataset for calibration of the NWM. Heterogeneity of the soil column depth is also not yet resolved in WRF-Hydro. These limitations of WRF-Hydro could be addressed in subsequent releases of the operational NWM. The observed improvements to modeled ET fluxes, when WRF-Hydro was forced with gauge data and calibrated, suggests that the model could potentially resolve surface fluxes in a realistic manner. Calibration in the larger Babocomari basin, clearly demonstrates the added value of channel infiltration, as calibration without it yields an unrealistic hydrologic response.
We acknowledge Michael Smith (NOAA/NWS Office of Water Prediction) for providing us feedback as we configured the NWM WRF-Hydro Domains in the southwest United States. This work was supported by a University Corporation for Atmospheric Science (UCAR) COMET Cooperative Project and NOAA Joint Technology Transfer Initiative (JTTI) Federal Grant NA17OAR4590183. We acknowledge the USDA-ARS Southwest Watershed Research Center for development and operation of the WGEW research facilities and diligent long-term collection of high-quality hydrometeorological data. We additionally acknowledge our reviewers for their constructive feedback on this manuscript.
Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/JHM-D-18-0064.s1.