Abstract

A fully coupled land surface hydrologic model, Flux-PIHM, is developed by incorporating a land surface scheme into the Penn State Integrated Hydrologic Model (PIHM). The land surface scheme is adapted from the Noah land surface model. Because PIHM is capable of simulating lateral water flow and deep groundwater at spatial resolutions sufficient to resolve upland stream networks, Flux-PIHM is able to represent heterogeneities due to topography and soils at high resolution, including spatial structure in the link between groundwater and the surface energy balance (SEB). Flux-PIHM has been implemented at the Shale Hills watershed (0.08 km2) in central Pennsylvania. Multistate observations of discharge, water table depth, soil moisture, soil temperature, and sensible and latent heat fluxes in June and July 2009 are used to manually calibrate Flux-PIHM at hourly temporal resolution. Model predictions from 1 March to 1 December 2009 are evaluated. Both hydrologic predictions and SEB predictions show good agreement with observations. Comparisons of model predictions between Flux-PIHM and the original PIHM show that the inclusion of the complex SEB simulation only brings slight improvement in hourly model discharge predictions. Flux-PIHM adds the ability of simulating SEB to PIHM and does improve the prediction of hourly evapotranspiration, the prediction of total runoff (discharge), and the predictions of some peak discharge events, especially after extended dry periods. Model results reveal that annual average sensible and latent heat fluxes are strongly correlated with water table depth, and the correlation is especially strong for the model grids near the stream.

1. Introduction

The predictability of the atmosphere is limited by the chaotic nature of atmospheric turbulence to a time span of the order of one week or less (Lorenz 1969; Smagorinsky 1969; Lorenz 1982). Earth's surface, however, has “memories” much longer than those of the atmosphere. Significant improvements in short-term climate forecasts as well as weather forecasts can be found by including the modeling of Earth's surface, for example, land surface processes and sea surface temperature (SST), in predictive models (Palmer and Anderson 1994; Beljaars et al. 1996; Koster et al. 2000; Goddard et al. 2001; Ek et al. 2003; Mitchell et al. 2004; Kumar et al. 2008). While covering only 30% of Earth's surface, the land surface plays a distinctive role in weather and climate because of its considerable heterogeneity; its dynamic hydrologic cycle and strong variations of temperature; and its highly variable land use, land cover, and topography (Yang 2004). Land surface processes are critical in the growth of the atmospheric boundary layer (ABL), the formation of clouds and precipitation, and the budgets of heat, momentum, and moisture within the atmosphere. Weather and climate models currently rely on land surface models (LSMs) to represent land surface processes. LSMs provide lower boundary conditions and initialize ground states for numerical weather prediction (Ookouchi et al. 1984; Liang et al. 1994; Betts et al. 1997; F. Chen et al. 1997; Xiu and Pleim 2001; Ek et al. 2003; Kumar et al. 2006; Chen et al. 2007; Niu et al. 2011).

In the past few decades, LSMs have undergone significant development from the “bucket model” (Manabe 1969) to more sophisticated and more physical parameterizations. The evolution of LSMs reflects the community's effort to improve the representation of land surface fluxes, including the benefit of capturing the memory embedded in soil moisture (Ek et al. 2003; Mitchell et al. 2004). Groundwater has a longer memory than soil moisture and has been shown to influence the land surface and the atmosphere (Changnon 1987; Dooge 1992; Skøien et al. 2003; Liang et al. 2003; Maxwell and Miller 2005; Yeh and Eltahir 2005; Kollet and Maxwell 2008). Subsurface waters, however, are not well described in most LSMs. Traditional LSMs are limited to vertical moisture transport in the soil column, and most of them ignore deeper soil moisture processes and lack physical representations of water table or the explicit lateral flow to streams. Therefore, those models have limited ability in representing the contribution of groundwater to the memory of land surface. This groundwater memory effect is particularly challenging to simulate in lower-order watersheds (e.g., first- and second-order streams and upland watersheds), which encompass a large fraction of the land surface (Shreve 1969; Benda et al. 2005), are highly heterogeneous, and thus require hydrologic modeling at spatial resolutions that resolve the complexity of these small-scale systems.

The water cycle, or hydrologic cycle, is central to the Earth system, and water resources are one of the critical environmental and political issues of the twenty-first century (National Research Council 2004). Hydrologic models are important tools to enhance the understanding of hydrological processes and to simulate and predict hydrological events for better decision making. They have evolved from early empirical rainfall–runoff models (e.g., Mulvany 1851; Sherman 1932; Horton 1935) to conceptual models (e.g., Linsley and Crawford 1960; Crawford and Linsley 1966; Dawdy and O'Donnell 1965; Sugawara 1969; Burnash et al. 1973; Burnash 1995; Arnold et al. 1998; Ma and Cheng 1998), to physically based models (e.g., Abbott et al. 1986a,b; Beven and Kirkby 1976, 1979; Beven et al. 1987; Grayson et al. 1992; Fortin et al. 2001a,b; Qu 2004; Qu and Duffy 2007; Kumar 2009). Physically based, spatially distributed hydrologic models have the advantages of taking into account spatial heterogeneity of inputs and characterizing hydrologic variables in space (Beven 1985; Smith et al. 2004; Wood and Lettenmaier 2006; Li et al. 2009). They can also be applied to the watersheds where no adequate data are available for calibrating empirical and conceptual models. Physically based hydrologic models are therefore extremely important for flood/drought forecasting at small-scale, low-order watersheds and for the study and understanding of hydrologic processes.

Evapotranspiration is an important component of the water cycle and an important process in hydrologic models. The accurate forecasting of the timing and magnitude of peak discharge events depends on the forecasting of evapotranspiration, especially after extended dry periods (Kampf 2006). Hydrologic models, however, usually have simplified representations of evapotranspiration. They either take potential evapotranspiration as an external atmospheric forcing (e.g., Beven et al. 1987; Grayson et al. 1992; Burnash 1995; Yeh et al. 2006) or use simplified formulations and parameterizations for calculating evapotranspiration (e.g., Abbott et al. 1986b; Kumar 2009). Studies show that applying different methods for calculating potential evapotranspiration in hydrologic models produces significant differences in evapotranspiration and, consequently, in runoff simulations (Vörösmarty et al. 1998; Kampf 2006; Weiß and Menzel 2008).

Coupled models of land surface and subsurface, which incorporate hydrologic components into LSMs or couple deeper subsurface with the atmosphere, may yield improvements in weather and short-term climate forecasting and flood/drought forecasting. There has been recent interest in incorporating a groundwater component into LSMs or coupling deeper subsurface with the atmosphere to improve the representation of soil moisture at the land surface (e.g., York et al. 2002; Seuffert et al. 2002; Mölders and Rühaak 2002; Liang et al. 2003; Yeh and Eltahir 2005; Maxwell and Miller 2005; Gulden et al. 2007; Maxwell et al. 2007; Kollet and Maxwell 2008; Rosero et al. 2011; Niu et al. 2011). Effects of groundwater and soil moisture on the ABL (Maxwell et al. 2007), as well as improvement in energy fluxes and rainfall predictions (Seuffert et al. 2002; Mölders and Rühaak 2002) were found. Rigon et al. (2006) developed a distributed hydrological model with coupled water and energy budgets (GEOtop) and tested the model performance at a midsized watershed. Bouilloud et al. (2010) coupled a soil–vegetation–atmosphere transfer model [Soil, Biosphere, and Atmosphere (ISBA)] with a conceptual hydrologic model (TOPMODEL) and evaluated the model predictions of discharge at three midsized watersheds. Livneh et al. (2011) developed a unified land model by combining a 1D LSM (Noah LSM) with a spatially lumped conceptual hydrologic model [Sacramento Soil Moisture Accounting model (SAC-SMA)] and tested the coupled model at six midsized watersheds. In the Noah LSM with multiparameterization option (Noah-MP), Niu et al. (2011) implemented a simple groundwater model and a TOPMODEL-based runoff model into the original Noah LSM and evaluated Noah-MP performance at various small watersheds. These coupled modeling studies did not evaluate the subsurface hydrology and land surface fluxes simultaneously. Livneh et al. (2011) and Niu et al. (2011) evaluated subsurface hydrology in some watersheds and land surface fluxes in different watersheds. Rigon et al. (2006) used different simulations to evaluate each model component. A comprehensive evaluation of coupled model performance is needed.

Maxwell and Miller (2005), Maxwell et al. (2007), Kollet and Maxwell (2008), and Rihani et al. (2010) have incorporated a fully coupled, three-dimensional (subsurface and overland) groundwater flow model (ParFlow), with a land surface model [Community Land Model (CLM)] on one hand and a mesoscale atmospheric model [Advanced Regional Prediction System (ARPS)] on the other hand, to study the influences of lateral subsurface water flow and shallow water table depth on the surface energy balance. Results indicate that the coupled land surface groundwater model behaves differently from uncoupled land surface model in synthetic data experiments and reveals strong correlation between land surface and subsurface. Because of the complexity of optimization, the models in those studies were not well calibrated. The improvement in forecasts over uncoupled models, if any, and the interaction between land surface and subsurface, however, still needs further exploration.

Hydrologic model parameters generally need to be calibrated for the model system to reproduce observed hydrologic responses even when a priori parameters can be estimated from existing data (e.g., national soil survey and hydrogeologic databases). Calibration of hydrologic models has been the focus of many studies (e.g., Ibbitt 1970; Johnston and Pilgrim 1976; Pickup 1977; Gupta and Sorooshian 1985; Duan et al. 1992; Sorooshian et al. 1993; Franchini 1996; Vrugt et al. 2003; Wagener et al. 2003; Kollat and Reed 2006; Xie and Zhang 2010). Most of the previous studies focus on the calibration of model discharge, and sometimes water table depth, but tend to neglect other observations. However, “an acceptable model prediction might be achieved in many different ways, i.e., different model structures or parameter sets” (Beven 1993). This nonuniqueness of numerical model parameters and structures is called “equifinality.” Equifinality makes model calibration difficult using only one type of observation. If multiple parameter sets produce equally good discharge predictions, it is very difficult to find the optimal parameter set. One possible solution to equifinality is to use multistate observations to constrain the parameters. Some studies also found that using observations of subsurface conditions (e.g., soil moisture) in addition to discharge observations improves the forecast of streamflow (e.g., Oudin et al. 2003; Aubert et al. 2003; Francois et al. 2003; Camporese et al. 2009; Lee et al. 2011).

Calibrating model parameters using high temporal resolution observations may improve model representation of important hydrologic processes in low-order watersheds. Most previous studies calibrated and evaluated hydrologic models at a daily time scale, ignoring the subdaily changes of hydrologic variables. Because of the impacts of losing streams, infiltration, precipitation, evapotranspiration, and freeze–thaw processes, discharge and groundwater level have diurnal cycles (Lundquist and Cayan 2002; Gribovszki et al. 2010). The diurnal fluctuation of groundwater table is most significant in summer and in highly forested areas and reaches up to 11 cm (Thal-Larsen 1934) in some areas. In low-order watersheds, streamflow and groundwater level exhibit even more temporal variability than in larger basins (Reed et al. 2004) and are influenced by the complication of rapidly varying climatic gradients and topographic effects, as well as complex hydrologic and geomorphic conditions that control basin storage and runoff. Calibration using high temporal resolution observations could help capture the rapidly changing processes and improve the forecasting in low-order watersheds.

In this study, a coupled land surface hydrologic modeling system is developed and evaluated at a small-scale forested watershed. In particular, the Penn State Integrated Hydrologic Model (PIHM; Qu 2004; Qu and Duffy 2007; Kumar 2009) is coupled with the land surface schemes in the Noah LSM (Chen and Dudhia 2001; Ek et al. 2003). PIHM is a fully coupled, spatially distributed, and physically based hydrologic model. This model has advanced model physics, for example, fully coupled surface and subsurface flow, lateral surface and subsurface water flow, and macropore flow, and is capable of resolving small-scale hydrologic modeling at low-order watersheds. The model is implemented at a first-order watershed, the Shale Hills watershed (0.08 km2) in central Pennsylvania. A National Science Foundation (NSF)-sponsored Critical Zone Observatory (CZO), the Susquehanna Shale Hills Critical Zone Observatory (SSHCZO), now exists in this watershed. CZOs are operated at watershed scale to advance the understanding of Earth's surface processes. SSHCZO brings together multiple research disciplines to observe and quantify Earth's surface processes at hillslope to small-watershed scales. The broad array of observations at SSHCZO enables an unprecedented investigation of subsurface–land surface–atmosphere interactions and makes SSHCZO an ideal site for the coupled model test. Multistate observations, including discharge, water table depth, soil moisture, soil temperature, and surface heat fluxes, are used for the comprehensive model calibration and evaluation. The performance of the original PIHM version 2.0 (PIHM V2) is compared with the coupled land surface hydrologic model. Model predictions of hydrologic variables and land surface variables are evaluated at hourly resolution.

2. Description of the coupled land surface hydrologic model

a. The Penn State Integrated Hydrologic Model

PIHM is a multiprocess and multiscale hydrologic model. It simulates evapotranspiration, infiltration, recharge, overland flow, groundwater flow, and channel routing in a fully coupled scheme. The model has been tested at both small-sized (e.g., Qu 2004; Li 2010) and midsized watersheds (e.g., Kumar 2009). Here we summarize the major methodological elements of PIHM. Triangular and rectangular elements are used in PIHM. The land surface is decomposed into unstructured triangular elements while rivers are represented by rectangular elements. The use of an unstructured mesh provides an optimal representation of local heterogeneities in parameters and process dynamics with the least number of elements. It also provides better representation of linear features within the model domain such as river channels and watershed boundaries. Triangular and rectangular elements at the land surface are projected vertically down to bedrock to form prismatic volumes. Surface water is able to flow along topography, infiltrate into soil, or evaporate. Channel flow and overland flow are governed by the 1D (channel flow) and 2D (overland flow) St. Venant equations (Saint-Venant 1871). Infiltration is calculated through the top 10-cm-layer of soil. The subsurface prismatic volume is further subdivided into unsaturated and saturated zones. Soil water is restricted to vertical transport in the unsaturated zone. In the saturated zone, groundwater flow is horizontal with dynamic coupling to the unsaturated zone across the free surface (water table). Unsaturated flow and groundwater flow are described using the Richards equation, and the unsaturated hydraulic conductivities are calculated using the van Genuchten (1980) equations. Macropore flows are simulated in PIHM to account for the rapid water flows through root holes and soil cracks, that is, macropores. Macropores penetrate the soil from land surface to a depth into the soil, defined as the macropore depth. The effective hydraulic conductivity of the subsurface is considered as a weighted average conductivity of the macropores and the soil matrix. The semidiscrete finite volume method (Qu and Duffy 2007) is applied, and the system of governing differential equations is solved simultaneously. PIHM also includes a simple temperature index representation of snowmelt (Kumar 2009). Detailed descriptions and formulations of PIHM are provided by Qu (2004), Qu and Duffy (2007), Kumar (2009), and Li (2010).

b. The land surface scheme

The new land surface hydrologic model Flux-PIHM incorporates a land surface scheme into PIHM. The land surface scheme is mainly adapted from the Noah LSM (Chen and Dudhia 2001; Ek et al. 2003), which has undergone extensive testing (e.g., Chen et al. 1996; T. H. Chen 1997; Liang et al. 1998; Chen and Mitchell 1999; Koren et al. 1999; Schlosser et al. 2000; Ek et al. 2003; Boone et al. 2004) and has been implemented into mesoscale atmospheric models, for example, the fifth-generation Pennsylvania State University–National Center for Atmospheric Research (NCAR) Mesoscale Model (MM5) and the Weather Research and Forecast Model (WRF). The Noah LSM integrates the Penman potential evaporation scheme of Mahrt and Ek (1984), the multiple-layer soil model of Mahrt and Pan (1984), the canopy model of Pan and Mahrt (1987), and the canopy resistance approach of Noilhan and Planton (1989) and Jacquemin and Noilhan (1990) (Chen et al. 2007). The Noah LSM has one canopy layer and four soil layers. It is able to capture diel and seasonal variations in surface heat fluxes and temperature with relatively few parameters. Its relatively simple formulation makes it easy to couple to PIHM.

The Noah LSM is adapted and coupled with PIHM. The soil water retention model, the canopy water capacity formulation, and the canopy drip rate formulation in the Noah LSM are replaced by the PIHM formulations in the coupled land surface scheme. To keep in accordance with PIHM, the Cosby et al. (1984) soil water retention model used in the Noah LSM is replaced by the van Genuchten (1980) soil model in PIHM; in the van Genuchten (1980) equation,

 
formula

where θ is the saturation ratio of soil, α and β are the van Genuchten soil parameters, and h is the hydraulic head. The maximum canopy water capacity Wcmax, which is a constant in the Noah LSM, is formulated as a linear function of the leaf area index (LAI):

 
formula

where S is the reference canopy water capacity, and its default value is set to 0.2 mm. The canopy drip rate is formulated as

 
formula

where the reference drip rate kD = 5.65 × 10−2 m day−1, Wc is the water storage in the canopy, Δt is the time step, and the parameter b has a range from 3.0 to 4.6 as suggested by Rutter and Morton (1977) and is set to be 3.89 here.

The prismatic volume of each grid is divided into several soil layers depending on the depth of the prism. From the ground surface to the bottom, the standard thicknesses of the top four layers are 0.1, 0.3, 0.6, and 1.0 m, respectively, as in the Noah LSM. If the bedrock depth is less than 2 m, the number of soil layers and the thickness of the lowest layer are adapted to match the depth of bedrock. If the bedrock depth is larger than 2 m, additional soil layers are added as needed. Soil moisture contents and soil temperatures of those multiple layers are simulated. The volumetric soil moisture (Θ) equations in the Noah LSM are adapted to allow coupling with PIHM. For the top soil layer,

 
formula

for the bottom soil layer J above the water table,

 
formula

and for the layer j in between,

 
formula

where z represents the vertical direction, is the thickness of the jth soil layer, K is the hydraulic conductivity, D = K(∂h/∂Θ) is the soil water diffusivity, Esoil is the evaporation from soil, and represents the canopy transpiration taken by root in the jth layer. The variable I is the surface water infiltration, and R is the groundwater recharge. When coupled with PIHM, both I and R are provided by the hydrologic model. The Noah LSM uses a parameterized infiltration scheme and a gravitational free-drainage subsurface runoff scheme (Schaake et al. 1996). In comparison, infiltration and groundwater recharge in this model [Eqs. (4a) and (4b)] are calculated using a physical hydrologic model, which should provide more realistic boundary conditions for the LSM. Further, the Noah LSM evapotranspiration formulation, including extraction of transpiration from the rooting zone should provide more realistic upper-boundary conditions for groundwater hydrology. Equation (4) therefore provides the core innovation in this coupled model, Flux-PIHM.

c. Fully coupled land surface hydrologic modeling system

Flux-PIHM is thus a fully coupled land surface hydrologic modeling system. The surface energy balance scheme described in section 2b completely replaces the original evapotranspiration formulation in PIHM. The land surface and hydrologic components are coupled by exchanging water table depth, soil moisture, infiltration rate, recharge rate, net precipitation (i.e., precipitation that is not intercepted by canopy) rate, and evapotranspiration rate between each other, as described in Fig. 1. The 1D land surface scheme is integrated into every model grid, and the passing of information is performed at every time step. The hydrologic component and the land surface component share the same model parameters at subsurface (e.g., hydraulic conductivities). At each time step, the hydrologic component provides the land surface component with water table depth, infiltration rate, recharge rate, and integrated soil water storage over soil column. The soil layers below the water table are set to saturation. Soil moisture equations [Eq. (4)] are applied to the soil layers above the water table. Infiltration rate and recharge rate calculated by the hydrologic component are used as top and bottom boundary conditions for soil moisture transport. PIHM only predicts the total soil water storage within each soil column. The volumetric soil moisture contents simulated by the land surface component are rescaled using the integrated soil water storage provided by PIHM to guarantee mass conservation:

 
formula

and the average volumetric soil moisture content of the unsaturated zone SWC is calculated as

 
formula

where Θj is the volumetric soil moisture at the jth soil layer predicted by the land surface module, is the volumetric soil moisture after rescaling, hunsat is the integrated soil water storage predicted by PIHM, BRD is the bedrock depth, hsat is the groundwater level predicted by PIHM, Θs is the soil porosity, and Θr is the residual porosity. The land surface component then starts the surface energy balance simulation and provides the hydrologic component with the net precipitation rate and evapotranspiration rate.

Fig. 1.

Coupling between the hydrologic model (PIHM) and the land surface energy balance model (adapted from the Noah LSM) yielding the integrated model, Flux-PIHM.

Fig. 1.

Coupling between the hydrologic model (PIHM) and the land surface energy balance model (adapted from the Noah LSM) yielding the integrated model, Flux-PIHM.

Because PIHM is capable of simulating lateral groundwater flow that explicitly resolves low-order channel networks, Flux-PIHM is able to represent the land surface heterogeneities caused by topography and capture the time scales that are generally missing in current land surface models. At the same time, the robust land surface scheme provides accurate sensible heat flux and latent heat flux (evapotranspiration) simulations at all scales in the watershed or river network. The coupled hydrologic and land surface schemes guarantee mass conservation at the subsurface, as well as the land surface, and conserve energy balance at land surface, which provides physical constraints to surface heat fluxes and subsurface water movement.

3. Study site and data

a. The Shale Hills critical zone observatory

The Shale Hills watershed is a 0.08 km2 watershed located in the valley and ridge physiographic province of central Pennsylvania, at 40°39.87′N, 77°54.40′W. The watershed is a small, forested, and temperate-climate catchment carved in shale bedrock. The catchment is V-shaped and is characterized by relatively steep slopes (25%–48%) and narrow ridges. Surface elevation varies from 256 m MSL at the outlet to 310 m MSL at ridge top. The first-order stream within the catchment is a tributary of Shaver's Creek that eventually reaches the Juniata River. The sloping areas and ridges of the watershed are covered by several typical deciduous species. The valley floor and north-facing ridge top are covered by some evergreen species (Wubbels 2010). Five soil series are identified within the watershed (Lin 2006; Lin and Zhou 2008). SSHCZO has been the focus of hydrologic research for several decades (e.g., Lynch 1976; Qu and Duffy 2007; Lin and Zhou 2008; Ma et al. 2010).

Extensive multidisciplinary field surveys have been conducted at SSHCZO, including soil survey, tree survey, surface elevation measurements, and bedrock depth measurements. A real-time hydrologic monitoring network (RTHnet) is operating in SSHCZO. RTHnet provides real-time and high-frequency observations from bedrock to ABL. RTHnet instrument arrays include sensors for soil moisture, soil temperature, soil matric potential, groundwater level and temperature, and snow depth measurement. The arrays also include sap flux measurements, a single above-canopy eddy-covariance flux tower, and a weather station. All available data can be found at http://www.czo.psu.edu/.

b. Model setup and model parameters

The model simulation period is from 0000 UTC 1 January 2009 to 0000 UTC 1 January 2010 with a model time step of 1 min and an output interval of 1 h. The model domain is prepared using the PIHM Geographic Information System (PIHM GIS). The model domain is decomposed into a triangular irregular network of 535 grids with 299 nodes. The average grid size is 157 m2. The stream channel is represented by 20 stream segments. The model grids are presented in Fig. 2. Depending on the bedrock depth, the model grids have either four (if the bedrock depth is equal to or smaller than 2.0 m) or five (if the bedrock depth is larger than 2.0 m) vertical soil layers. A no-flow boundary condition is applied to the watershed boundary (except for the stream outlet) and the aquifer lower boundary. A zero-depth gradient boundary condition (Morris 1979) is applied at the last stream segment, which is the outlet of the watershed.

Fig. 2.

Grid setting for the Shale Hills watershed model domain. The watershed boundary, the stream path, and the locations of RTHnet measurements used in this study are shown. The inset shows the location of SSHCZO within Pennsylvania.

Fig. 2.

Grid setting for the Shale Hills watershed model domain. The watershed boundary, the stream path, and the locations of RTHnet measurements used in this study are shown. The inset shows the location of SSHCZO within Pennsylvania.

Figure 3 shows the vegetation type, soil type, surface elevation, and bedrock depth defined within the model domain. The surface terrain of the Shale Hills watershed is represented by a 1-m resolution digital elevation model (DEM) digitized from airborne light detection and ranging (lidar) data (http://pihm.ics.psu.edu/CZO_NOSL/). The high-resolution DEM is used to develop the numerical mesh necessary to resolve upland surface and subsurface processes.

Fig. 3.

The (a) vegetation type, (b) soil type, (c) surface elevation, and (d) bedrock depth defined in Flux-PIHM for the simulated domain.

Fig. 3.

The (a) vegetation type, (b) soil type, (c) surface elevation, and (d) bedrock depth defined in Flux-PIHM for the simulated domain.

A bedrock depth map was obtained from a field campaign conducted in 2003 (Lin et al. 2006). Bedrock depths range from less than 0.25 m on ridge tops to greater than 3 m along the channel. To account for flow through a deeper weathered shale, an extra 1.5 m is added to measured bedrock depths for all model mesh points. The effective bedrock depths are shown in Fig. 3. The soil map and some soil matrix properties, including vertical saturated hydraulic conductivity, horizontal saturated hydraulic conductivity, and porosity, are also obtained from the field campaign in 2003 (Lin et al. 2006; Lin 2006). Other soil matrix properties used in this study are acquired from the Soil Survey Geographic (SSURGO) national database (http://www.HydroTerre.psu.edu). For soil wilting point (Θw) and field capacity (Θref), the methods in the Noah LSM are adapted. They are calculated as

 
formula

and

 
formula

where θcr is the saturation ratio at which soil matrix vertical hydraulic conductivity equals 0.5 mm day−1. The values of Θw and Θref for different soil types at the Shale Hills watershed can be found in Table 1. Because of the lack of measurements, empirical values are used for soil macropore properties and streambed properties, and those values are subject to calibration. Soil horizontal macropore hydraulic conductivity and soil vertical macropore hydraulic conductivity are assumed to be 100 times their corresponding soil matrix conductivities (Li 2010) and are subject to calibration. The vegetation map is obtained from a tree survey conducted in 2008. Vegetation parameters for different land cover types are obtained from the modified International Geosphere–Biosphere Programme (IGBP) Moderate Resolution Imaging Spectroradiometer (MODIS) 20-category vegetation (land use) data (http://www.ral.ucar.edu/research/land/technology/lsm/parameters/), which are also used in the Noah LSM. Root zone depths are assumed to be 0.6 m for all vegetation types. The model parameters and their a priori values are listed in Tables 1, 2, and 3.

Table 1.

Soil parameters used for the Shale Hills watershed domain. Listed values are the a priori (uncalibrated) parameter values. All parameters in this table are calibrated in the optimization process.

Soil parameters used for the Shale Hills watershed domain. Listed values are the a priori (uncalibrated) parameter values. All parameters in this table are calibrated in the optimization process.
Soil parameters used for the Shale Hills watershed domain. Listed values are the a priori (uncalibrated) parameter values. All parameters in this table are calibrated in the optimization process.
Table 2.

Parameters used for the river segments. Listed values are the uncalibrated parameter values. All parameters in this table are calibrated in the optimization process.

Parameters used for the river segments. Listed values are the uncalibrated parameter values. All parameters in this table are calibrated in the optimization process.
Parameters used for the river segments. Listed values are the uncalibrated parameter values. All parameters in this table are calibrated in the optimization process.
Table 3.

Vegetation parameters used for the Shale Hills watershed domain. Only the parameter Rcmin is calibrated in the optimization process.

Vegetation parameters used for the Shale Hills watershed domain. Only the parameter Rcmin is calibrated in the optimization process.
Vegetation parameters used for the Shale Hills watershed domain. Only the parameter Rcmin is calibrated in the optimization process.

c. Forcing data and evaluation data

Meteorological data from in situ measurements and remote sensing datasets are used to drive and evaluate the model (Table 4). Given the small scale (0.08 km2) of the watershed, spatially uniform forcing is used for this study. To enable the model to respond to realistic temporal variations in meteorology, as many in situ measurements as possible are used to drive the model. Precipitation, air temperature, and relative humidity data are obtained from the RTHnet weather station (http://cataract.cee.psu.edu/czo/rth/). Precipitation is measured with a Hach OTT precipitation gauge, while air temperature and relative humidity are measured with a Campbell Scientific HMP45C temperature and relative humidity probe. Those 10-min interval data are aggregated into hourly data. Downward longwave radiation, downward solar radiation, and surface pressure (prior to April 2009) are not available at SSHCZO and are obtained from the Surface Radiation Budget Network (SURFRAD) Penn State University station. The SURFRAD Penn State University station is located at 40.72°N, 77.93°W, which is 6.48 km away from the Shale Hills watershed. Its surface elevation is 376 m. It is assumed that the radiation and surface pressure at the Penn State SURFRAD station represent the conditions at the Shale Hills watershed. Wind speed data from January 2009 to March 2009 are obtained from the Penn State SURFRAD station. An above-canopy eddy-covariance flux tower was installed in SSHCZO in April 2009. Wind speed and air pressure from the flux tower are used for the period from April to December 2009.

Table 4.

Forcing and evaluation data used for simulation.

Forcing and evaluation data used for simulation.
Forcing and evaluation data used for simulation.

MODIS LAI data are introduced to represent vegetation phenology. MODIS product MOD15A2 provides 8-day composite LAI data at 1-km resolution (Knyazikhin et al. 1999; Myneni et al. 2002). LAI field measurements have also been collected regularly from 25 April to 31 October 2010 at SSHCZO, using a Decagon AccuPAR meter and a LI-COR 2200 leaf area. The MODIS LAI data in 2009 are rescaled based on the comparison between the MODIS product and the CZO field measurements to drive the model:

 
formula

where LAI is the rescaled MODIS LAI, LAIMODIS is the MODIS LAI measurements, LAICZO is the watershed-averaged SSHCZO field-measured LAI in 2010, and LAI0 is a reference value set to be 0.5 m2 m−2.

Evaluation data include outlet discharge, groundwater level, soil moisture, soil temperature, and eddy-covariance surface heat flux measurements (Table 4) of the year 2009. A V-notch weir measures water level at the outlet of catchment. The measured stream stage is converted to discharge rate using a rating curve developed by Nutter (1964) for the V-notch weir at the Shale Hills watershed. Groundwater level and soil moisture data are collected from three RTHnet wells drilled near the stream. Each well is equipped with one water level sensor (Druck pressure transducer CS420-L manufactured by Campbell Scientific) and three soil moisture sensors (Decagon Echo2 probes). The three soil moisture sensors at each well are located at three different levels below ground, with depths of 0.1, 0.3, and 0.5 m below the surface, respectively. Among the three RTHnet wells, Well 1 (Fig. 2) is not deep enough to capture the change of the water table in the summer months, when the water table drops below the screened well depth. For those time periods, only measurements from RTHnet Wells 2 and 3 are used. Two of the nine soil moisture sensors continuously show measurements larger than 1 m3 m−3. Therefore, data from those two sensors are not used in this study. The water table depth (WTD; distance from the land surface to the groundwater table) measurements collected in different wells are averaged to represent the observed WTD at RTHnet wells. The multiple volumetric soil moisture content observations from those three wells are averaged to represent the observed soil water content (SWC) over the soil column. The standard deviations from different water level and soil moisture sensor measurements are calculated as the uncertainty range for observed WTD and SWC. In Flux-PIHM, the model domain is discretized as such that the three RTHnet wells are located at three vertices of one model grid for the convenience of model–data comparison. This model grid surrounded by RTHnet wells is selected for the comparison of WTD and SWC. Soil temperature measurements are collected at different locations and different vertical levels in 2009. Soil temperatures are measured with 229 probes manufactured by Campbell Scientific and 5TE probes manufactured by Decagon. The locations of the five soil temperature sites are shown in Fig. 2. For those five sites, the soil temperature data measured at 5 cm below surface are averaged. Model predictions of soil temperature at the top soil layer (0–10 cm) from the five corresponding model grids are averaged to be compared with the observations.

The above-canopy eddy-covariance flux tower measures wind speed and air temperature with a Campbell Scientific CSAT3 three-dimensional sonic anemometer. Carbon dioxide and water vapor concentrations are measured with a LI-COR LI-7500 CO2/H2O analyzer. Eddy-covariance fluxes are calculated following the quality control methods documented by Vickers and Mahrt (1997). The gaps in processed 30-min sensible and latent heat fluxes are filled using the lookup table method suggested by Falge et al. (2001). The gap-filled sensible and latent heat fluxes are aggregated into hourly data to be compared with hourly predictions of sensible and latent heat fluxes averaged over the model domain.

4. Optimization and evaluation of Flux-PIHM

a. Model optimization and spinup

Flux-PIHM has a large number of tunable model parameters. Some of those parameters are dependent on soil type or land cover type. Thus, calibration could be very complicated if the model domain contains a large number of soil types or land cover types. To simplify the calibration problem, the single global calibration multiplier method (or one-factor method) (Pokhrel and Gupta 2010; Wallner et al. 2012) is adopted. One single global calibration coefficient is used for each model parameter regardless of soil type or land cover type. The global calibration coefficient is a multiplier acting on the corresponding soil- or vegetation-related parameter for all soil or vegetation types. By applying global multiplier calibration coefficients, the dimension of parameter space for calibration is reduced and the ratios between uncalibrated a priori parameters of different soil/vegetation types are preserved. The default value for calibration coefficient is 1.0.

A comprehensive manual calibration is performed using the “trial and error” strategy. In the comprehensive calibration process, outlet discharge, water table depth, soil water content, soil temperature at 5 cm below surface, and eddy-covariance surface heat flux data from June to July 2009 are used to adjust model parameters. This period is chosen because it covers both peak discharge events and a base flow period and exhibits a considerable rate of evapotranspiration. The model parameters that are calibrated and their corresponding calibration coefficients are listed Table 5.

Table 5.

Global calibration coefficient values optimized for Flux-PIHM at the Shale Hills watershed model domain. PIHM V2 does not have parameters Czil, Θref, and Θw. For the other parameters, the same calibration coefficients are used in both Flux-PIHM and PIHM V2.

Global calibration coefficient values optimized for Flux-PIHM at the Shale Hills watershed model domain. PIHM V2 does not have parameters Czil, Θref, and Θw. For the other parameters, the same calibration coefficients are used in both Flux-PIHM and PIHM V2.
Global calibration coefficient values optimized for Flux-PIHM at the Shale Hills watershed model domain. PIHM V2 does not have parameters Czil, Θref, and Θw. For the other parameters, the same calibration coefficients are used in both Flux-PIHM and PIHM V2.

After calibration, Flux-PIHM spinup uses atmospheric forcing from 0000 UTC 20 October 2008 to 0000 UTC 1 January 2009, initialized from the “relaxation mode,” in which the model relaxes from a saturation state. In wet periods such as the winter conditions used for this spinup, it takes approximately three weeks for the models to reach an equilibrium state at the Shale Hills watershed. This 2-month spinup period is used to ensure the model eliminating the effects of initial conditions. Initial surface skin temperature is set to the air temperature and the soil temperatures are obtained from a linear interpolation between the surface skin temperature and the bottom boundary condition.

To evaluate the performance of the coupled land surface hydrologic model, Flux-PIHM is compared with the original PIHM V2 using the same hydrologic and land surface parameters as in Flux-PIHM. PIHM V2 has the identical hydrologic scheme as Flux-PIHM but with a simplified evapotranspiration scheme. The same model parameter sets and initial conditions are used for both Flux-PIHM, and PIHM V2, if applicable, to ensure unbiased comparison.

b. Water budget

Predicted water budgets from Flux-PIHM and PIHM V2 are compared for the period from 1 March to 1 December 2009 and are presented in Fig. 4. Snow-covered conditions are excluded from the model evaluation because Flux-PIHM is using simple snow physics that are not sufficient for accurate hydrologic prediction in winter. Therefore, this study focuses on the warm seasons. Runoff in Fig. 4 is calculated as the total discharge divided by the total area of the watershed. Predicted total runoffs from both models are close to the observation, and Flux-PIHM provides a closer total runoff prediction compared with PIHM V2. Flux-PIHM overestimates the total runoff by only 0.4% (1 mm) while PIHM V2 overestimates the total runoff by 4.1% (11 mm).

Fig. 4.

Comparison of water budget between Flux-PIHM and PIHM V2 from 0000 UTC 1 Mar to 0000 UTC 1 Dec 2009. Observed discharge from RTHnet and the total precipitation during this period is also presented.

Fig. 4.

Comparison of water budget between Flux-PIHM and PIHM V2 from 0000 UTC 1 Mar to 0000 UTC 1 Dec 2009. Observed discharge from RTHnet and the total precipitation during this period is also presented.

The total evapotranspiration predicted by both models are close and only differ by 8 mm. The partitioning of total evapotranspiration, however, is very different in the two models. Flux-PIHM predicts higher canopy evaporation and transpiration compared with PIHM V2, while PIHM V2 predicts higher soil evaporation. Soil evaporation and transpiration extract water from different subsurface layers. Soil evaporation only has access to the top soil layer, while transpiration can transport deep soil water or even groundwater to the atmosphere, depending on the root zone depth and the water table depth. For the six model grids near the stream outlet, transpiration that is extracted directly from the groundwater is calculated. For those six model grids, in Flux-PIHM 34.8% of total evapotranspiration is extracted directly from groundwater, while in PIHM V2, this fraction is only 22.1% because of the lower transpiration prediction in PIHM V2. Such differences will produce different soil moisture profiles and thus affect evapotranspiration and flood/drought forecasting.

c. Evaluation of hydrologic predictions: Discharge, water table depth, and soil water content

To test the model's ability to predict the hydrological states of the watershed, model discharge, WTD, and SWC predictions are compared with in situ measurements. Several criteria are used to quantify model predictions of discharge, including the total bias, the Nash–Sutcliffe model efficiency coefficient (NSE; Nash and Sutcliffe 1970), and the correlation coefficient with observation. The total bias is defined as

 
formula

where Qo is the observed discharge, Qm is the predicted discharge, T is the length of observation, and a superscript t indicates observation or modeled discharge at the time step t. NSE is the most commonly used coefficient to evaluate hydrologic models. It is defined as

 
formula

The correlation coefficient is defined as

 
formula

Figure 5 presents the comparison of hourly outlet discharge among RTHnet measurements and two models from March to December 2009. The statistics of discharge predictions are presented in Table 6.

Fig. 5.

Comparison of hourly outlet discharge among RTHnet measurements, Flux-PIHM, and PIHM V2 predictions from 0000 UTC 1 Mar to 0000 UTC 1 Dec 2009. The inset highlights the peak discharge event from 15 to 19 October, for which Flux-PIHM and PIHM V2 have considerably different predictions.

Fig. 5.

Comparison of hourly outlet discharge among RTHnet measurements, Flux-PIHM, and PIHM V2 predictions from 0000 UTC 1 Mar to 0000 UTC 1 Dec 2009. The inset highlights the peak discharge event from 15 to 19 October, for which Flux-PIHM and PIHM V2 have considerably different predictions.

Table 6.

Bias, NSE, and correlation coefficient of Flux-PIHM and PIHM V2 hourly discharge predictions compared with RTHnet measurements from 1 March to 1 December 2009. NSE is defined as .

Bias, NSE, and correlation coefficient of Flux-PIHM and PIHM V2 hourly discharge predictions compared with RTHnet measurements from 1 March to 1 December 2009. NSE is defined as .
Bias, NSE, and correlation coefficient of Flux-PIHM and PIHM V2 hourly discharge predictions compared with RTHnet measurements from 1 March to 1 December 2009. NSE is defined as .

Flux-PIHM predictions agree well with measurements for both dry period flows and peak flows as illustrated in Fig. 5. Though June and July data are used for calibration, Flux-PIHM still overestimates the peak event in June. In March, April, and May, Flux-PIHM underestimates base flows but overestimates some discharge peaks. The largest errors occur in October, when Flux-PIHM underestimates the highest peak as well as overestimates a series of peak events before that.

Flux-PIHM and PIHM V2 discharge simulations are almost indistinguishable for most of the period. The only notable differences are for the peak events in May, June, and October. The predictions of the two models differ substantially for the 15–19 October flood event, which is highlighted in Fig. 5. The observed peak discharge rate is 277 m3 day−1, the Flux-PIHM predicted peak discharge rate is 778 m3 day−1, and the PIHM V2 predicted peak discharge rate is 1516 m3 day−1. The timings of the predicted peak discharge also differ by 6 h in the two models. As for the average performance, Table 6 shows that Flux-PIHM has improved performance in discharge prediction over PIHM V2, especially at the daily time step.

The comparisons of WTD and SWC among RTHnet well measurements and Flux-PIHM and PIHM V2 predictions are presented in Fig. 6. The biases, correlation coefficients, and root-mean-square errors (RMSEs) of WTD and SWC predictions compared with RTHnet wells are presented in Table 7. Figure 6 shows that Flux-PIHM predictions of WTD and SWC agree well with RTHnet measurements. For most of the evaluation period, the differences between Flux-PIHM predictions and RTHnet measurements are within the observation uncertainty ranges for both WTD and SWC. Like discharge predictions, the largest errors appear in October, when model errors in WTD and SWC fall out of the observation uncertainty ranges. In October, Flux-PIHM continuously underestimates WTD and overestimates SWC, when the error in WTD reaches about −0.5 m and the error in SWC reaches up to 0.1 m3 m−3. In spring (March–May), although model errors of WTD and SWC are always within the observation uncertainty ranges, the model predictions have a consistent wet bias (low bias in WTD and high bias in SWC) compared with observed average WTD and SWC. Figure 6 also shows that the PIHM V2 predictions in WTD and SWC follow the observations better than Flux-PIHM in March, but PIHM V2 is outperformed by Flux-PIHM in almost all following months. The largest differences between Flux-PIHM and PIHM V2 appear in October, when the differences between the two models reach up to 0.25–0.30 m. As indicated in Table 7, Flux-PIHM provides better overall predictions of WTD and SWC than PIHM V2, although PIHM V2 WTD bias is slightly smaller.

Table 7.

Bias, RMSE, and correlation coefficient of Flux-PIHM and PIHM V2 hourly WTD and SWC predictions compared with RTHnet well measurements from 1 March to 1 December 2009.

Bias, RMSE, and correlation coefficient of Flux-PIHM and PIHM V2 hourly WTD and SWC predictions compared with RTHnet well measurements from 1 March to 1 December 2009.
Bias, RMSE, and correlation coefficient of Flux-PIHM and PIHM V2 hourly WTD and SWC predictions compared with RTHnet well measurements from 1 March to 1 December 2009.
Fig. 6.

As in Fig. 5, but for (a) WTD and (b) SWC. The shaded areas indicate the uncertainty range of observations.

Fig. 6.

As in Fig. 5, but for (a) WTD and (b) SWC. The shaded areas indicate the uncertainty range of observations.

d. Evaluation of surface energy balance predictions: Sensible and latent heat fluxes and soil temperature

The aggregated hourly sensible (H) and latent (LE) heat fluxes measured at SSHCZO are compared with model hourly predictions averaged over the model domain. Both H and LE measurements at SSHCZO are available from 1 April 2009. The period from 1 April 2009 to 1 January 2010 is divided into four subperiods to evaluate surface heat flux predictions for different seasons: 1 April to 1 June 2009 as spring, 1 June to 1 September 2009 as summer, 1 September to 1 December 2009 as fall, and 1 December 2009 to 1 January 2010 as winter. Plots of hourly H and LE are not presented to save space: instead, Table 8 lists the statistics for each season based on the comparison between hourly predictions and observations of H and LE. Besides, Fig. 7 presents the comparisons of surface heat fluxes among Flux-PIHM, PIHM V2, and RTHnet eddy-covariance measurements as averaged daily cycles for each season. Flux-PIHM predicts both H and LE, but PIHM V2 only predicts evapotranspiration rate, which is equivalent to LE.

Table 8.

RMSE and correlation coefficient of Flux-PIHM and PIHM V2 hourly H and LE predictions compared with eddy-covariance flux tower measurements for different seasons (from 1 April 2009 to 1 January 2010).

RMSE and correlation coefficient of Flux-PIHM and PIHM V2 hourly H and LE predictions compared with eddy-covariance flux tower measurements for different seasons (from 1 April 2009 to 1 January 2010).
RMSE and correlation coefficient of Flux-PIHM and PIHM V2 hourly H and LE predictions compared with eddy-covariance flux tower measurements for different seasons (from 1 April 2009 to 1 January 2010).
Fig. 7.

Comparison of surface heat fluxes among Flux-PIHM, PIHM V2, and RTHnet measurements as averaged daily cycles at each hour of the day for different seasons.

Fig. 7.

Comparison of surface heat fluxes among Flux-PIHM, PIHM V2, and RTHnet measurements as averaged daily cycles at each hour of the day for different seasons.

Figure 7 shows that Flux-PIHM predictions of midday H and LE have consistent high biases during all seasons. It needs to be noted that Flux-PIHM closes the energy budget at the land surface, while eddy-covariance measurements typically fail to close the energy budget, with the measured H + LE tending to be less than RnG (McNeil and Shuttleworth 1975; Fritschen et al. 1992; Twine et al. 2000; Foken 2008). Therefore, in the calibration process, Flux-PIHM is subjectively calibrated to yield higher H and LE, assuming the eddy-covariance method underestimates surface heat fluxes. Generally, Flux-PIHM prediction of LE has a larger bias in the afternoon than in the morning, while H has a larger bias in the morning than in the afternoon (Fig. 7). Except for spring and winter LE, Flux-PIHM hourly predictions of H and LE have high correlations with RTHnet measurements, with correlation coefficients between 0.91 and 0.94 (Table 8), which shows that Flux-PIHM is capable of capturing the variation of surface heat fluxes well at hourly resolution. Flux-PIHM bias of LE in spring is much higher than in other seasons (Fig. 7). The RMSE for spring LE is larger than the RMSEs in other seasons, and the correlation coefficient for LE in spring is also relatively lower than in summer and fall (Table 8). This substantial high bias in spring agrees with the signs of biases in WTD and SWC predictions in the corresponding months (Fig. 6). Because the calibration of Flux-PIHM only focused on summer months, the seasonality of evapotranspiration is not taken into account in calibration. Using a longer calibration period or including more parameters for calibration would improve the model prediction of the seasonality of evapotranspiration. It is also possible that this high bias in spring is caused by model structural errors or forcing data errors. The correlation coefficient of winter LE is as low as 0.55, but LE measurements of an open-path gas analyzer in winter are less reliable than in other seasons (Burba et al. 2008). The insufficient winter physics in Flux-PIHM may also contribute to the low correlation of LE in winter.

When averaged into daily cycles, Flux-PIHM and PIHM V2 LE are almost indistinguishable in summer and winter (Fig. 7), but predictions of LE by Flux-PIHM are higher in spring and fall. In spring, summer, and fall, Flux-PIHM LE predictions have higher correlation coefficients and lower RMSEs compared to PIHM V2, which indicates that Flux-PIHM is more capable of reproducing the hourly variation in evapotranspiration.

Figure 8 presents the comparison of hourly soil temperature at 5 cm below the surface between Flux-PIHM predictions and RTHnet observations. From March to November, the average bias for soil temperature prediction at 5 cm below the surface is 0.81°C and the RMSE is 1.6°C. The correlation coefficient between model and measurement is 0.97. Figure 8 suggests that the soil temperature predicted by Flux-PIHM has a consistent high bias in summer, but the errors are random in other seasons. At the beginning of March, the error could be as large as ±5°C. This is probably because of the simple snow physics in Flux-PIHM, which may not simulate thaw well.

Fig. 8.

Comparison of hourly soil temperature at 5 cm below surface between Flux-PIHM prediction and RTHnet observations from 0000 UTC 1 Mar to 0000 UTC 1 Dec 2009.

Fig. 8.

Comparison of hourly soil temperature at 5 cm below surface between Flux-PIHM prediction and RTHnet observations from 0000 UTC 1 Mar to 0000 UTC 1 Dec 2009.

e. Correlation between surface heat fluxes and water table depth

The coupled land surface hydrologic model is a good tool for the study of the interaction between the subsurface and the land surface. Annual average sensible and latent heat fluxes as functions of WTD in each model grid from Flux-PIHM simulation for the year 2009 are presented in Fig. 9. Sensible heat fluxes are positively correlated with water table depth and latent heat fluxes are negatively correlated with water table depth. It is also clear that the change of surface heat fluxes with respect to WTD for the Ernest soil type (the diamond-shaped symbols) is much more pronounced than for the other soil types. Because the Ernest soil type is found primarily in the lowest portions of the valley (Fig. 3b), it indicates that impacts of groundwater are stronger on those grids near the stream than the other model grids. For the grids near the stream, the groundwater table is relatively shallow and roots have access to groundwater; thus, groundwater has direct impacts on the surface energy balance. For the grids far from the stream, the groundwater table is deep and roots have no access to groundwater; thus, the impacts of groundwater are weak. For those grids, the effect of groundwater on surface heat fluxes is indirect and is accomplished by affecting the soil moisture profile.

Fig. 9.

Annual average sensible and latent heat fluxes plotted as functions of WTD from Flux-PIHM simulation. Each symbol represents one grid in the Shale Hills watershed model domain. Different soil types are represented by different shapes and different land cover types are plotted using different colors. Blairton soil type is not plotted because only two grids are categorized as Blairton.

Fig. 9.

Annual average sensible and latent heat fluxes plotted as functions of WTD from Flux-PIHM simulation. Each symbol represents one grid in the Shale Hills watershed model domain. Different soil types are represented by different shapes and different land cover types are plotted using different colors. Blairton soil type is not plotted because only two grids are categorized as Blairton.

This correlation between surface heat fluxes and groundwater table has been found in other studies (e.g., Kollet and Maxwell 2008). Kollet and Maxwell (2008) also discovered that when the water table depth is smaller than 1 m or greater than 5 m in the Little Washita watershed, Oklahoma, the land surface and groundwater are decoupled. This decoupling of land surface and groundwater is not found in this study at the Shale Hills watershed.

5. Discussion and conclusions

This paper presents a physically based, land surface hydrologic model, Flux-PIHM, that performs well in a low-order watershed at high temporal resolution when calibrated with multistate, watershed-averaged observations (discharge, WTD, SWC, soil temperature, and surface heat fluxes). Flux-PIHM adds the simulation of the surface energy balance to PIHM while maintaining or improving the predictions of discharge, WTD, SWC, and evapotranspiration. The accurate total discharge prediction is critical for long-term flood/drought forecasting, climate and water resources interactions, and surface energy balance simulation. A high-spatial-resolution (tens of meters), high-temporal-resolution (tens of minutes), multistate, high-fidelity simulation of watershed function in low-order watersheds (a large fraction of the landscape) has been elusive to date. This improved model performance in upland watersheds should also improve our overall understanding of groundwater–land surface–atmosphere interaction over larger basins.

Flux-PIHM's multistate calibration is achieved with high-temporal-resolution data (hourly) with only two months of observations (June and July 2009) and by adjusting a relatively small number of land surface parameters. The model is then evaluated at hourly resolution over the period from March through December of 2009. Compared with other studies, the calibration period for Flux-PIHM is relatively short, the temporal resolution of the calibration and evaluation data is high, and the diversity of data types used for calibration and evaluation is high. Thus, robust model performance is achieved with a relatively small calibration period. Using a longer calibration period and including more parameters for calibration may improve model predictions at other times of the year.

Flux-PIHM simulates interactions between the surface energy balance and hydrologic processes that cannot be achieved without the coupling of these model components. The model results indicate that mean annual surface heat fluxes are highly correlated with water table depth at a significant subset of the grid points in the watershed. The coupling between groundwater and land surface is affected by soil type, land cover type, and topography. The strong coupling between SEB and water table depth implies that LSMs that do not simulate water table are prone (at least at some grid points) to long-term biases as the water table fluctuates.

Flux-PIHM is compared with PIHM V2, which does not have a surface energy balance simulation. Flux-PIHM is shown to have a significantly different water budget for the test watershed, while maintaining similar or slightly improved performance in the components simulated by both models. Flux-PIHM predicts 44% (372 mm) of the total water export from the watershed in March–December 2009 is in the form of transpiration, compared to 37% (311 mm) for PIHM V2, and only 15% (127 mm) in the form of soil evaporation, compared to 23% (192 mm) for PIHM V2. These differences, though not dramatic, have the potential to result in substantially different soil water profiles and responses of the watershed to flood or drought events. Comparison to direct measurements upscaled sap flux and spatially resolved soil moisture measurements would enable future studies to determine which system yields a more accurate simulation of watershed function.

Results indicate that Flux-PIHM does slightly improve the overall hourly predictions of outlet discharge, WTD, and SWC compared with PIHM V2. Although the improvement of model discharge predictions are not clear in Fig. 5, the Flux-PIHM prediction has a smaller total bias, a higher correlation coefficient with observation, and a higher NSE compared with PIHM V2. Differences in discharge predictions are noticeable for some high discharge peaks, especially after the extended dry, low-flow summer period. For a flood event from 15 to 19 October, Flux-PIHM and PIHM V2 provide substantially different predictions in peak discharge rate and timing; the Flux-PIHM simulation is noticeably closer to the observed discharge. Because the hydrologic schemes in those two models are identical, the differences in the predictions of discharge rate and timing are likely caused by the different evapotranspiration calculations. These differences demonstrate the impact of evapotranspiration simulation on flood prediction. We need to point out that the comparisons between Flux-PIHM and PIHM V2 are made by using the same manually calibrated parameter values, and those parameter values may not represent the optimal parameters for both models, although good agreements between model predictions and observations have been shown.

The errors that remain in Flux-PIHM postcalibration illustrate areas for potential future improvements in the modeling system. Notable errors in Flux-PIHM predictions occur in March–May and October. From March to May, the model systematically overestimates SWC and underestimates WTD. Perhaps because of the high bias in soil moisture, Flux-PIHM produces a high bias in LE during this period. Spring simulations are challenging and could be improved by better winter period/snow physics. In October, the model underestimates the highest peak discharge event but overestimates some prior peak discharge events. The largest errors in WTD and SWC during the evaluation period occur in October as well. Improved simulation of the recovery of the watershed from dry, low-flow conditions is another topic of future research.

Acknowledgments

This research was supported by the National Oceanic and Atmospheric Administration through Grant NA10OAR4310166, and the National Science Foundation Susquehanna Shale Hills Critical Zone Observatory project through Grant EAR 0725019. We thank Dr. Henry Lin, Dr. David Eissenstat, and Dr. Kusum Naithani for providing the data and Kelly Cherrey and Colin Duffy for their support of the instrumentation and data collection.

REFERENCES

REFERENCES
Abbott
,
M. B.
,
J. C.
Bathurst
,
J. A.
Cunge
,
P. E.
O'Connell
, and
J.
Rasmussen
,
1986a
:
An introduction to the European Hydrological System—Systeme Hydrologique Europeen, “SHE”, 1: History and philosophy of a physically-based, distributed modeling system
.
J. Hydrol.
,
87
,
45
59
.
Abbott
,
M. B.
,
J. C.
Bathurst
,
J. A.
Cunge
,
P. E.
O'Connell
, and
J.
Rasmussen
,
1986b
:
An introduction to the European Hydrological System—Systeme Hydrologique Europeen, “SHE”, 2: Structure of a physically-based, distributed modeling system
.
J. Hydrol.
,
87
,
61
77
.
Arnold
,
J. G.
,
R.
Srinivasan
,
R. S.
Muttiah
, and
J. R.
Williams
,
1998
:
Large area hydrologic modeling and assessment. Part I: Model development
.
J. Amer. Water Resour. Assoc.
,
34
,
73
89
.
Aubert
,
D.
,
C.
Loumagne
, and
L.
Oudin
,
2003
:
Sequential assimilation of soil moisture and streamflow data in a conceptual rainfall-runoff model
.
J. Hydrol.
,
280
,
145
161
.
Beljaars
,
A. C. M.
,
P.
Viterbo
,
M. J.
Miller
, and
A. K.
Betts
,
1996
:
The anomalous rainfall over the United States during July 1993: Sensitivity to land surface parameterization and soil moisture anomalies
.
Mon. Wea. Rev.
,
124
,
362
383
.
Benda
,
L.
,
M. A.
Hassan
,
M.
Church
, and
C. L.
May
,
2005
:
Geomorphology of steepland headwaters: The transition from hillslopes to channels
.
J. Amer. Water Resour. Assoc.
,
41
,
835
851
.
Betts
,
A. K.
,
F.
Chen
,
K. E.
Mitchell
, and
Z. I.
Janjić
,
1997
:
Assessment of land surface and boundary layer models in two operational versions of the NCEP Eta Model using FIFE data
.
Mon. Wea. Rev.
,
125
,
2896
2916
.
Beven
,
K. J.
,
1985
: Distributed modelling. Hydrological Forecasting, M. G. Anderson and T. P. Burt, Eds., Wiley, 405–435.
Beven
,
K. J.
,
1993
:
Prophecy, reality and uncertainty in distributed hydrological modelling
.
Adv. Water Resour.
,
16
,
41
51
.
Beven
,
K. J.
, and
M. J.
Kirkby
,
1976
: Towards a simple physically based variable contributing model of catchment hydrology. Working Paper 154, School of Geography, University of Leeds, United Kingdom, 11 pp.
Beven
,
K. J.
, and
M. J.
Kirkby
,
1979
:
A physically based, variable contributing area model of basin hydrology
.
Hydrol. Sci. Bull.
,
24
,
43
69
.
Beven
,
K. J.
,
A.
Calver
, and
E. M.
Morris
,
1987
: The Institute of Hydrology distributed model. Tech. Rep. 98, Institute of Hydrology, Wallingford, United Kingdom, 32 pp. [Available online at http://nora.nerc.ac.uk/5977/1/IH_098.pdf.]
Boone
,
A.
, and
Coauthors
,
2004
:
The Rhône-Aggregation land surface scheme intercomparison project: An overview
.
J. Climate
,
17
,
187
208
.
Bouilloud
,
L.
, and
Coauthors
,
2010
:
Coupling the ISBA land surface model and the TOPMODEL hydrological model for Mediterranean flash-flood forecasting: Description, calibration, and validation
.
J. Hydrometeor.
,
11
,
315
333
.
Burba
,
G. G.
,
D. K.
McDermitt
,
A.
Grelle
,
D. J.
Anderson
, and
L.
Xu
,
2008
:
Addressing the influence of instrument surface heat exchange on the measurements of CO2 flux from open-path gas analyzers
.
Global Change Biol.
,
14
,
1854
1876
.
Burnash
,
R. J. C.
,
1995
: The NWS river forecast system: Catchment modeling. Computer Models of Watershed Hydrology, V. P. Singh, Ed., Water Resources Publications, 311–366.
Burnash
,
R. J. C.
,
R. L.
Ferral
, and
R. A.
McGuire
,
1973
: A generalized streamflow simulation system: Conceptual modeling for digital computers. Tech. Rep., U.S. Department of Commerce, National Weather Service, Silver Spring, Maryland, 204 pp.
Camporese
,
M.
,
C.
Paniconi
,
M.
Putti
, and
P.
Salandin
,
2009
: Ensemble Kalman filter data assimilation for a process-based catchment scale model of surface and subsurface flow. Water Resour. Res.,45, W10421, doi:10.1029/2008WR007031.
Changnon
,
S. A.
,
1987
: Detecting drought conditions in Illinois. Illinois State Water Survey Circular 169-87, Department of Energy and Natural Resources, State of Illinois, 36 pp.
Chen
,
F.
, and
K.
Mitchell
,
1999
:
Using the GEWEX/ISLSCP forcing data to simulate global soil moisture fields and hydrological cycle for 1987–1988
.
J. Meteor. Soc. Japan
,
77
,
167
182
.
Chen
,
F.
, and
J.
Dudhia
,
2001
:
Coupling an advanced land surface–hydrology model with the Penn State–NCAR MM5 modeling system. Part I: Model implementation and sensitivity
.
Mon. Wea. Rev.
,
129
,
569
585
.
Chen
,
F.
, and
Coauthors
,
1996
:
Modeling of land surface evaporation by four schemes and comparison with FIFE observations
.
J. Geophys. Res.
,
101
,
7251
7268
.
Chen
,
F.
,
Z.
Janjić
, and
K.
Mitchell
,
1997
:
Impact of atmospheric-surface layer parameterizations in the new land-surface scheme of the NCEP mesoscale Eta numerical model
.
Bound.-Layer Meteor.
,
85
,
391
421
.
Chen
,
F.
, and
Coauthors
,
2007
:
Description and evaluation of the characteristics of the NCAR high-resolution land data assimilation system
.
J. Appl. Meteor. Climatol.
,
46
,
694
713
.
Chen
,
T. H.
, and
Coauthors
,
1997
:
Cabauw experimental results from the project for intercomparison of land-surface parameterization schemes
.
J. Climate
,
10
,
1194
1215
.
Cosby
,
B. J.
,
G. M.
Hornberger
,
R. B.
Clapp
, and
T. R.
Ginn
,
1984
:
A statistical exploration of the relationships of soil moisture characteristics to the physical properties of soils
.
Water Resour. Res.
,
20
,
682
690
.
Crawford
,
N. H.
, and
R. K.
Linsley
,
1966
: Digital simulation in hydrology: Stanford Watershed Model IV. Tech. Rep. 39, Stanford University, Palo Alto, CA, 225 pp. [Available online at http://www.hydrocomp.com/publications/StanfordModelIV.PDF.]
Dawdy
,
D. R.
, and
T.
O'Donnell
,
1965
:
Mathematical models of catchment behavior
.
J. Hydraul. Div.
,
91
,
123
127
.
Dooge
,
J. C. I.
,
1992
:
Hydrologic models and climate change
.
J. Geophys. Res.
,
97
,
2677
2686
.
Duan
,
Q.
,
S.
Sorooshian
, and
V. K.
Gupta
,
1992
:
Effective and efficient global optimization for conceptual rainfall-runoff models
.
Water Resour. Res.
,
28
,
1015
1031
.
Ek
,
M. B.
,
K.
Mitchell
,
Y.
Lin
,
E.
Rogers
,
P.
Grummann
,
V.
Koren
,
G.
Gayno
, and
J.
Tarpley
,
2003
:
Implementation of Noah land surface model advances in the National Centers for Environmental Prediction operational Mesoscale Eta Model
.
J. Geophys. Res.
,
108
,
8851
,
doi:10.1029/2002JD003296
.
Falge
,
E.
, and
Coauthors
,
2001
:
Gap filling strategies for defensible annual sums of net ecosystem exchange
.
Agric. For. Meteor.
,
107
,
43
69
.
Foken
,
T.
,
2008
:
The energy balance closure problem: An overview
.
Ecol. Appl.
,
18
,
1351
1367
.
Fortin
,
J. P.
,
R.
Turcotte
,
S.
Massicotte
,
R.
Moussa
,
J.
Fitzback
, and
J. P.
Villeneuve
,
2001a
:
A distributed watershed model compatible with remote sensing and GIS data. I: Description of model
.
J. Hydrol. Eng.
,
6
,
91
99
.
Fortin
,
J. P.
,
R.
Turcotte
,
S.
Massicotte
,
R.
Moussa
,
J.
Fitzback
, and
J. P.
Villeneuve
,
2001b
:
A distributed watershed model compatible with remote sensing and GIS data. II: Application to Chaudière watershed
.
J. Hydrol. Eng.
,
6
,
100
108
.
Franchini
,
M.
,
1996
:
Use of a genetic algorithm combined with a local search method for the automatic calibration of conceptual rainfall-runoff models
.
Hydrol. Sci. J.
,
41
,
21
39
.
Francois
,
C.
,
A.
Quesney
, and
C.
Ottlé
,
2003
:
Sequential assimilation of ERS-1 SAR data into a coupled land surface-hydrological model using an extended Kalman filter
.
J. Hydrometeor.
,
4
,
473
487
.
Fritschen
,
L. J.
,
P.
Qian
,
E. T.
Kanemasu
,
D.
Nie
,
E. A.
Smith
,
J. B.
Stewart
,
S. B.
Verma
, and
M. L.
Wesely
,
1992
:
Comparisons of surface flux measurement systems used in FIFE 1989
.
J. Geophys. Res.
,
97
(
D17
),
18 697
18 713
.
Goddard
,
L.
,
S. J.
Mason
,
S. E.
Zebiak
,
C. F.
Ropelewski
,
R.
Basher
, and
M. A.
Cane
,
2001
:
Current approaches to seasonal to interannual climate predictions
.
Int. J. Climatol.
,
21
,
1111
1152
.
Grayson
,
R. B.
,
I. D.
Moore
, and
T. A.
McMahon
,
1992
:
Physically based hydrologic modeling. 1: A terrain-based model for investigative purposes
.
Water Resour. Res.
,
28
,
2639
2658
.
Gribovszki
,
Z.
,
J.
Szilágyi
, and
P.
Kalicz
,
2010
:
Diurnal fluctuations in shallow groundwater levels and streamflow rates and their interpretation—A review
.
J. Hydrol.
,
385
,
371
383
.
Gulden
,
L. E.
,
E.
Rosero
,
Z.-L.
Yang
,
M.
Rodell
,
C. S.
Jackson
,
G.-Y.
Niu
,
P. J.-F.
Yeh
, and
J.
Famiglietti
,
2007
: Improving land-surface model hydrology: Is an explicit aquifer model better than a deeper soil profile? Geophys. Res. Lett.,34, L09402, doi:10.1029/2007GL029804.
Gupta
,
V. K.
, and
S.
Sorooshian
,
1985
:
The relationship between data and the precision of parameter estimates of hydrologic models
.
J. Hydrol.
,
81
,
57
77
.
Horton
,
R. E.
,
1935
: Surface Runoff Phenomena. Part 1: Analysis of the Hydrograph. Horton Hydrology Laboratory Publication, No. 101, Edward Bros., 73 pp.
Ibbitt
,
R. P.
,
1970
: Systematic parameter fitting for conceptual models of catchment hydrology. Ph.D. thesis, University of London, 408 pp. [Available online at http://docs.niwa.co.nz/library/public/Ibbrisyst.pdf.]
Jacquemin
,
B.
, and
J.
Noilhan
,
1990
:
Sensitivity study and validation of a land surface parameterization using the HAPEX-MOBILHY data set
.
Bound.-Layer Meteor.
,
52
,
93
134
.
Johnston
,
P. R.
, and
D. H.
Pilgrim
,
1976
:
Parameter optimization for watershed models
.
Water Resour. Res.
,
12
,
477
486
.
Kampf
,
S. K.
,
2006
: Towards improved representations of hydrologic processes: Linking integrated and distributed hydrologic measurements to a physically-based model for a planar hillslope plot. Water Resources Series Tech. Rep. 183, Department of Civil and Environmental Engineering, University of Washington, Seattle, WA, 213 pp. [Available online at http://depts.washington.edu/uwbg/research/Burges/_WRS183-Kampf-2006.pdf.]
Knyazikhin
,
Y.
, and
Coauthors
,
1999
: MODIS leaf area index (LAI) and fraction of photosynthetically active radiation absorbed by vegetation (FPAR) product (MOD15). Algorithm Theoretical Basis Doc., NASA Goddard Space Flight Center, Greenbelt, MD, 126 pp. [Available online at http://modis.gsfc.nasa.gov/data/atbd/atbd_mod15.pdf.]
Kollat
,
J. B.
, and
P. M.
Reed
,
2006
:
Comparing state-of-the-art evolutionary multi-objective algorithms for long-term groundwater monitoring design
.
Adv. Water Resour.
,
29
,
792
807
.
Kollet
,
S. J.
, and
R. M.
Maxwell
,
2008
:
Capturing the influence of groundwater dynamics on land surface processes using an integrated, distributed watershed model
.
Water Resour. Res.
,
44
,
W02402
,
doi:10.1029/2007WR006004
.
Koren
,
V.
,
J.
Schaake
,
K.
Mitchell
,
Q.-Y.
Duan
,
F.
Chen
, and
J. M.
Baker
,
1999
:
A parameterization of snowpack and frozen ground intended for NCEP weather and climate models
.
J. Geophys. Res.
,
104
(
D16
),
19 569
19 585
.
Koster
,
R. D.
,
M. J.
Suarez
, and
M.
Heiser
,
2000
:
Variance and predictability of precipitation at seasonal-to-interannual timescales
.
J. Hydrometeor.
,
1
,
26
46
.
Kumar
,
M.
,
2009
: Toward a hydrologic modeling system. Ph.D. thesis, The Pennsylvania State University, 251 pp.
Kumar
,
S. V.
, and
Coauthors
,
2006
:
Land information system: An interoperable framework for high resolution land surface modeling
.
Environ. Modell. Software
,
21
,
1402
1415
,
doi:10.1016/j.envsoft.2005.07.004
.
Kumar
,
S. V.
,
R. H.
Reichle
,
C. D.
Peters-Lidard
,
R. D.
Koster
,
X.
Zhan
,
W. T.
Crow
,
J. B.
Eylander
, and
P. R.
Houser
,
2008
:
A land surface data assimilation framework using the land information system: Description and applications
.
Adv. Water Resour.
,
31
,
1419
1432
.
Lee
,
H.
,
D. J.
Seo
, and
V.
Koren
,
2011
:
Assimilation of streamflow and in-situ soil moisture data into operational distributed hydrologic models: Effects of uncertainties in the data and initial model soil moisture states
.
Adv. Water Resour.
, 34, 1597–1615.
Li
,
H.
,
L.
Luo
,
E. F.
Wood
, and
J.
Schaake
,
2009
: The role of initial conditions and forcing uncertainties in seasonal hydrologic forecasting. J. Geophys. Res.,114, D04114, doi:10.1029/2008JD010969.
Li
,
W.
,
2010
: Implementing the Shale Hills watershed model in application of PIHM. M.S. thesis, Department of Civil Engineering, The Pennsylvania State University, 97 pp.
Liang
,
X.
,
D. P.
Lettenmaier
,
E. F.
Wood
, and
S. J.
Burges
,
1994
:
A simple hydrologically based model of land surface water and energy fluxes for GCMs
.
J. Geophys. Res.
,
99
,
415
428
.
Liang
,
X.
, and
Coauthors
,
1998
:
The Project for Intercomparison of Land-Surface Parameterization Schemes (PILPS) phase 2 (c) Red-Arkansas River basin experiment: 2. Spatial and temporal analysis of energy fluxes
.
Global Planet. Change
,
19
,
137
159
.
Liang
,
X.
,
Z.
Xie
, and
M.
Huang
,
2003
:
A new parameterization for surface and groundwater interactions and its impact on water budgets with the variable infiltration capacity (VIC) land surface model
.
J. Geophys. Res.
,
108
,
8613
,
doi:10.1029/2002JD003090
.
Lin
,
H.
,
2006
:
Temporal stability of soil moisture spatial pattern and subsurface preferential flow pathways in the Shale Hills catchment
.
Vadose Zone J.
,
5
,
317
340
.
Lin
,
H.
, and
X.
Zhou
,
2008
:
Evidence of subsurface preferential flow using soil hydrologic monitoring in the Shale Hills catchment
.
Eur. J. Soil Sci.
,
59
,
34
49
.
Lin
,
H.
,
W.
Kogelmann
,
C.
Walker
, and
M. A.
Bruns
,
2006
:
Soil moisture patterns in a forested catchment: A hydropedological perspective
.
Geoderma
,
131
,
345
368
.
Linsley
,
R. K.
, and
N. H.
Crawford
,
1960
:
Computation of a synthetic streamflow record on a digital computer
.
Int. Assoc. Sci. Hydrol. Pub.
,
51
,
526
538
.
Livneh
,
B.
,
P. J.
Restrepo
, and
D. P.
Lettenmaier
,
2011
:
Development of a unified land model for prediction of surface hydrology and land-atmosphere interactions
.
J. Hydrometeor.
,
12
,
1299
1320
.
Lorenz
,
E. N.
,
1969
:
Atmospheric predictability as revealed by naturally occurring analogues
.
J. Atmos. Sci.
,
26
,
636
646
.
Lorenz
,
E. N.
,
1982
:
Atmospheric predictability experiments with a large numerical model
.
Tellus
,
34
,
505
513
.
Lundquist
,
J. D.
, and
D. R.
Cayan
,
2002
:
Seasonal and spatial patterns in diurnal cycles in streamflow in the western United States
.
J. Hydrometeor.
,
3
,
591
603
.
Lynch
,
J. A.
,
1976
: Effects of antecedent moisture on storage hydrographs. Ph.D. thesis, The Pennsylvania State University, 192 pp.
Ma
,
L.
,
F.
Chabaux
,
E.
Pelt
,
E.
Blaes
,
L.
Jin
, and
S.
Brantley
,
2010
:
Regolith production rates calculated with uranium-series isotopes at Susquehanna/Shale Hills Critical Zone Observatory
.
Earth Planet. Sci. Lett.
,
297
,
211
225
.
Ma
,
X.
, and
W.
Cheng
,
1998
:
A modeling of hydrological processes in a large low plain area including lakes and ponds
.
J. Japan Soc. Hydrol. Water Resour.
,
9
,
320
329
.
Mahrt
,
L.
, and
M.
Ek
,
1984
:
The influence of atmospheric stability on potential evaporation
.
J. Climate Appl. Meteor.
,
23
,
222
234
.
Mahrt
,
L.
, and
H.
Pan
,
1984
:
A two-layer model of soil hydrology
.
Bound.-Layer Meteor.
,
29
,
1
20
.
Manabe
,
S.
,
1969
:
Climate and the ocean circulation. I. The atmospheric circulation and the hydrology of the Earth's surface
.
Mon. Wea. Rev.
,
97
,
739
774
.
Maxwell
,
R. M.
, and
N. L.
Miller
,
2005
:
Development of a coupled land surface and groundwater model
.
J. Hydrometeor.
,
6
,
233
247
.
Maxwell
,
R. M.
,
F. K.
Chow
, and
S. J.
Kollet
,
2007
:
The groundwater-land-surface-atmosphere connection: Soil moisture effects on the atmospheric boundary layer in fully-coupled simulations
.
Adv. Water Resour.
,
30
,
2447
2466
.
McNeil
,
D. D.
, and
W. J.
Shuttleworth
,
1975
:
Comparative measurements of the energy fluxes over a pine forest
.
Bound.-Layer Meteor.
,
9
,
297
313
.
Mitchell
,
K. E.
, and
Coauthors
,
2004
:
The multi-institution North American Land Data Assimilation System (NLDAS): Utilizing multiple GCIP products and partners in a continental distributed hydrological modeling system
.
J. Geophys. Res.
,
109
,
D07S90
,
doi:10.1029/2003JD003823
.
Mölders
,
N.
, and
W.
Rühaak
,
2002
:
On the impact of explicitly predicted runoff on the simulated atmospheric response to small-scale land-use changes—An integrated modeling approach
.
Atmos. Res.
,
63
,
3
38
.
Morris
,
E. M.
,
1979
:
The effect of the small-slope approximation and lower boundary conditions on solutions of the Saint-Venant equations
.
J. Hydrol.
,
40
,
31
47
.
Mulvany
,
T. J.
,
1851
:
On the use of self-registering rain and flood gauges in making observations of the relations of rainfall and of flood discharges in a catchment
.
Trans. Inst. Civ. Eng. Ireland
,
4
,
1
8
.
Myneni
,
R. B.
, and
Coauthors
,
2002
:
Global products of vegetation leaf area and fraction absorbed PAR from year one of MODIS data
.
Remote Sens. Environ.
,
83
,
214
231
.
Nash
,
J. E.
, and
J. V.
Sutcliffe
,
1970
:
River flow forecasting through conceptual models. Part I: A discussion of principles
.
J. Hydrol.
,
10
,
282
290
.
National Research Council
,
2004
: Groundwater Fluxes Across Interfaces. National Academies Press, 85 pp.
Niu
,
G.-Y.
, and
Coauthors
,
2011
:
The community Noah land surface model with multiparameterization options (Noah-MP): 1. Model description and evaluation with local-scale measurements
.
J. Geophys. Res.
,
116
,
D12109
,
doi:10.1029/2010JD015139
.
Noilhan
,
J.
, and
S.
Planton
,
1989
:
A simple parameterization of land surface processes for meteorological models
.
Mon. Wea. Rev.
,
117
,
536
549
.
Nutter
,
W. L.
,
1964
: Determination of the head-discharge relationship for a sharp-crested compound weir and a sharp-crested parabolic weir. M.S. thesis, Department of Forest Hydrology, The Pennsylvania State University, 87 pp.
Ookouchi
,
Y.
,
M.
Segal
,
R. C.
Kesseler
, and
R.
Pielke
,
1984
:
Evaluation of soil moisture effects of the generation and modification of mesoscale circulation
.
Mon. Wea. Rev.
,
112
,
2281
2292
.
Oudin
,
L.
,
A.
Weisse
,
C.
Loumagne
, and
S.
Le Hégarat-Mascle
,
2003
:
Assimilation of soil moisture into hydrological models for flood forecasting: A variational approach
.
Can. J. Remote Sens.
,
29
,
679
686
.
Palmer
,
T. N.
, and
D. L. T.
Anderson
,
1994
:
The prospects for seasonal forecasting—A review paper
.
Quart. J. Roy. Meteor. Soc.
,
120
,
755
793
.
Pan
,
H. L.
, and
L.
Mahrt
,
1987
:
Interaction between soil hydrology and boundary-layer development
.
Bound.-Layer Meteor.
,
38
,
185
202
.
Pickup
,
G.
,
1977
:
Testing the efficiency of algorithms and strategies for automatic calibration of rainfall-runoff models
.
Hydrol. Sci. Bull.
,
22
,
257
274
.
Pokhrel
,
P.
, and
H. V.
Gupta
,
2010
: On the use of spatial regularization strategies to improve calibration of distributed watershed models. Water Resour. Res.,46, W01505, doi:10.1029/2009WR008066.
Qu
,
Y.
,
2004
: An integrated hydrologic model for multi-process simulation using semi-discrete finite volume approach. Ph.D. thesis, The Pennsylvania State University, 136 pp.
Qu
,
Y.
, and
C. J.
Duffy
,
2007
: A semidiscrete finite volume formulation for multiprocess watershed simulation. Water Resour. Res.,43, W08419, doi:10.1029/2006WR005752.
Reed
,
S.
,
V.
Koren
,
M.
Smith
,
Z.
Zhang
,
F.
Moreda
,
D.-J.
Seo
, and
D. M. I. P.
Participants
,
2004
:
Overall distributed model intercomparison project results
.
J. Hydrol.
,
298
,
27
60
.
Rigon
,
R.
,
G.
Bertoldi
, and
T. M.
Over
,
2006
:
GEOtop: A distributed hydrological model with coupled water and energy budgets
.
J. Hydrometeor.
,
7
,
371
388
.
Rihani
,
J. F.
,
R. M.
Maxwell
, and
F. K.
Chow
,
2010
: Coupling groundwater and land surface processes: Idealized simulations to identify effects of terrain and subsurface heterogeneity on land surface energy fluxes. Water Resour. Res.,46, W12523, doi:10.1029/2010WR009111.
Rosero
,
E.
,
L. E.
Gulden
,
Z.-L.
Yang
,
L. G.
De Goncalves
,
G.-Y.
Niu
, and
Y. H.
Kaheil
,
2011
:
Ensemble evaluation of hydrologically enhanced Noah-LSM: Partitioning of the water balance in high-resolution simulations over the Little Washita River experimental watershed
.
J. Hydrometeor.
,
12
,
45
64
.
Rutter
,
A. J.
, and
A. J.
Morton
,
1977
:
A predictive model of rainfall interception in forests. III. Sensitivity of the model to stand parameters and meteorological variables
.
J. Appl. Ecol.
,
14
,
567
588
.
Saint-Venant
,
B.
,
1871
:
Theory of unsteady water flow with application to floods and to propagation of tides in river channels
.
Proc. French Acad. Sci.
,
73
,
148
154
.
Schaake
,
J. C.
,
V. I.
Koren
,
Q.-Y.
Duan
,
K.
Mitchell
, and
F.
Chen
,
1996
:
Simple water balance model for estimating runoff at different spatial and temporal scales
.
J. Geophys. Res.
,
101
,
7461
7475
.
Schlosser
,
C. A.
,
A. G.
Slater
,
A.
Robock
,
A. J.
Pitman
,
K. Y.
Vinnikov
,
A.
Henderson-Sellers
,
N. A.
Speranskaya
, and
K.
Mitchell
,
2000
:
Simulations of a boreal grassland hydrology at Valdai, Russia: PILPS phase 2(d)
.
Mon. Wea. Rev.
,
128
,
301
321
.
Seuffert
,
G.
,
P.
Gross
,
C.
Simmer
, and
E. F.
Wood
,
2002
:
The influence of hydrologic modeling on the predicted local weather: Two-way coupling of a mesoscale weather prediction model and a land surface hydrologic model
.
J. Hydrometeor.
,
3
,
505
523
.
Sherman
,
L. K.
,
1932
:
Stream flow from rainfall by the unit graph method
.
Eng. News-Rec.
,
108
,
501
505
.
Shreve
,
R. L.
,
1969
:
Stream lengths and basin areas in topologically random channel networks
.
J. Geol.
,
77
,
397
414
.
Skøien
,
J. O.
,
G.
Blöschl
, and
A. W.
Western
,
2003
:
Characteristic space scales and timescales in hydrology
.
Water Resour. Res.
,
39
,
1304
,
doi:10.1029/2002WR001736
.
Smagorinsky
,
J.
,
1969
:
Problems and promises of deterministic extended range forecasting
.
Bull. Amer. Meteor. Soc.
,
50
,
286
311
.
Smith
,
M. B.
,
D.-J.
Seo
,
V. I.
Koren
,
S. M.
Reed
,
Z.
Zhang
,
Q.
Duan
,
F.
Moreda
, and
S.
Cong
,
2004
:
The distributed model intercomparison project (DMIP): Motivation and experiment design
.
J. Hydrol.
,
298
,
4
26
.
Sorooshian
,
S.
,
Q.
Duan
, and
V. K.
Gupta
,
1993
:
Calibration of rainfall-runoff models: application of global optimization to the Sacramento soil moisture accounting model
.
Water Resour. Res.
,
29
,
1185
1194
.
Sugawara
,
M.
,
1969
: The flood forecasting by a series storage type model. Floods and Their Computation, IAHS Publ. 84,
555
559
.
Thal-Larsen
,
J. H.
,
1934
:
Fluctuations in the level of the phreatic surface with an atmospheric deposit in the form of dew
.
Bodenkd. Forsch.
,
4
,
223
233
.
Twine
,
T. E.
, and
Coauthors
,
2000
:
Correcting eddy-covariance flux underestimates over a grassland
.
Agric. For. Meteor.
,
103
,
279
300
.
van Genuchten
,
M. T.
,
1980
:
A closed-form equation for predicting the hydraulic conductivity of unsaturated soils
.
Soil Sci. Soc. Amer. J.
,
44
,
892
898
.
Vickers
,
D.
, and
L.
Mahrt
,
1997
:
Quality control and flux sampling problems for tower and aircraft data
.
J. Atmos. Oceanic Technol.
,
14
,
512
526
.
Vörösmarty
,
C. J.
,
C. A.
Federer
, and
A. L.
Schloss
,
1998
:
Potential evaporation functions compared on US watersheds: Possible implications for global-scale water balance and terrestrial ecosystem modeling
.
J. Hydrol.
,
207
,
147
169
.
Vrugt
,
J. A.
,
H. V.
Gupta
,
L. A.
Bastidas
,
W.
Bouten
, and
S.
Sorooshian
,
2003
:
Effective and efficient algorithm for multiobjective optimization of hydrologic models
.
Water Resour. Res.
,
39
,
1214
, doi:10.1029/2002WR001746.
Wagener
,
T.
,
N.
McIntyre
,
M. J.
Lees
,
H. S.
Wheater
, and
H. V.
Gupta
,
2003
:
Towards reduced uncertainty in conceptual rainfall-runoff modelling: Dynamic identifiability analysis
.
Hydrol. Processes
,
17
,
455
476
.
Wallner
,
M.
,
U.
Haberlandt
, and
J.
Dietrich
,
2012
:
Evaluation of different calibration strategies for large scale continuous hydrological modelling
.
Adv. Geosci.
,
31
,
67
74
.
Weiß
,
M.
, and
L.
Menzel
,
2008
:
A global comparison of four potential evapotranspiration equations and their relevance to stream flow modelling in semi-arid environments
.
Adv. Geosci.
,
18
,
15
23
.
Wood
,
A. W.
, and
D. P.
Lettenmaier
,
2006
:
A test bed for new seasonal hydrologic forecasting approaches in the western United States
.
Bull. Amer. Meteor. Soc.
,
87
,
1699
1712
.
Wubbels
,
J. K.
,
2010
: Tree species distribution in relation to stem hydraulic traits and soil moisture in a mixed hardwood forest in central Pennsylvania. M.S. thesis, Department of Horticulture, The Pennsylvania State University, 39 pp.
Xie
,
X.
, and
D.
Zhang
,
2010
:
Data assimilation for distributed hydrological catchment modeling via ensemble Kalman filter
.
Adv. Water Resour.
,
33
,
678
690
.
Xiu
,
A.
, and
J. E.
Pleim
,
2001
:
Development of a land surface model. Part I: Application in a mesoscale meteorology model
.
J. Appl. Meteor.
,
40
,
192
209
.
Yang
,
Z.-L.
,
2004
: Modeling land surface processes in short-term weather and climate studies. Observation, Theory and Modeling of Atmospheric Variability, X. Zhu et al., Eds., World Scientific Series on Meteorology of East Asia, Vol. 3, World Scientific, 288–313.
Yeh
,
G.-T.
,
G. B.
Huang
,
H.-P.
Cheng
,
F.
Zhang
,
H.-C.
Lin
,
E.
Edris
, and
D.
Richards
,
2006
: A first-principle, physics-based watershed model: WASH123D. Watershed Models, V. P. Singh and D. K. Frevert, Eds., CRC Press, 211–244.
Yeh
,
P. J.-F.
, and
E. A. B.
Eltahir
,
2005
:
Representation of water table dynamics in a land surface scheme. Part I: Model development
.
J. Climate
,
18
,
1861
1880
.
York
,
J. P.
,
M.
Person
,
W. J.
Gutowski
, and
T. C.
Winter
,
2002
:
Putting aquifers into atmospheric simulation models: An example from the Mill Creek Watershed, northeastern Kansas
.
Adv. Water Resour.
,
25
,
221
238
.
Zilitinkevich
,
S. S.
,
1995
: Non-local turbulent transport: Pollution dispersion aspects of coherent structure of convective flows. Air Pollution Theory and Simulation, H. Power, N. Moussiopoulos, and C. A. Brebbia, Eds., Vol. I, Air Pollution III, Computational Mechanics Publications, 53–60.