Macro-scale grid-based and subbasin-based hydrologic modeling: joint simulation and cross-calibration

The complexity and heterogeneity of hydrologic interactions over a wide range of spatial and temporal scales (Blöschl & Sivapalan 1995) are one of many challenges in advancing our understanding of the hydrologic cycle. An example of these complexities is how the climate change and land use change affect the hydrologic cycle at various scales (Blöschl et al. 2007). The scales at which the hydrologic processes are investigated depend on the issues to be addressed. Therefore, an increasing need for representing hydrologic processes at different scales has led to three general model classifications: micro-, meso-, and macro-scale hydrologic models. The differences across these three scales demonstrate the trade-offs between spatial capability and complexity of the models.

Micro-scale models take advantage of numerical methods to discretize the space domain into a large number of small cells to provide detailed modeling and in-depth understanding of specific hydrologic processes over relatively small areas. For example, Chu et al. (2013) developed a puddle-to-puddle (P2P) modeling system to simulate micro-scale puddle filling, spilling, merging, and splitting processes dynamics on depression-dominated surfaces. Clearly, intensive high-resolution data are required for such micro-scale modeling. Meso-scale models are developed for simulating dominant hydrologic processes in local watersheds. In a typical meso-scale model, a watershed can be divided into a number of subbasins that are well connected to the channelized drainage system in a semi-distributed fashion.

In contrast, macro-scale models can be associated with regional, continental, or even global scales and they often involve simulations of large-scale climatic and hydrologic processes. The complexity of a macro-scale model varies, mainly depending on data availability, computational resources, and research purposes. Unlike the micro- and meso-scale models, the concept of macro-scale modeling is relatively new (Sood & Smakhtin 2015). Macro-scale models offer a holistic view of hydrologic processes that can be applied for interconnected modeling of Earth systems and climate impact assessments (Gudmundsson et al. 2012a). Increasing availability of remote sensing datasets, high-performance computing resources, and improved representations of hydrologic processes are among the key factors that have contributed to the constant development of macro-scale models (Paiva et al. 2011). The capabilities, limitations, and uncertainties of macro-scale modeling have been investigated before (e.g. Döll et al. 2008; Gudmundsson et al. 2012b; Sood & Smakhtin 2015). For example, Döll et al. (2008) highlighted the advances, difficulties, and research perspectives associated with macro-scale modeling, and indicated that the results from macro-scale models are sensitive to: (1) uncertainties in input data due to measurement, interpolation, and other factors; (2) spatial and temporal discretization and the information density of input data; and (3) model structure.

Two approaches are commonly used for large-scale hydrologic modeling: (1) subbasin-based and (2) grid-based. Concerning the second modeling type, for instance, a grid-based structure has been adopted in macro-scale models such as the variable infiltration capacity model (Liang et al. 1994), the grid-to-grid model (Bell et al. 2007), and the California basin characterization model (Flint et al. 2013), where hydrologic processes are simulated for all grids separately. Tesfa et al. (2014) compared the scalability of surface runoff generation by employing both subbasin-based and grid-based discretization approaches for spatial resolutions of 0.125, 0.25, 0.5, and 1°, indicating that the subbasin-based approach performed better than the grid-based approach, especially for mountainous areas in winters. In addition, calibration and validation of a macro-scale grid-based model by using streamflow data are a challenge because the point-scale streamflow data measured at gauging stations may not be directly used for such purposes. The macro-scale grid-based models require grid-based hydrologic parameters, which are not normally available or measured to be used in the calibration process. Therefore, some parameter regionalization techniques have been developed (e.g. Samaniego et al. 2010) to specify the localized parameters for a group of grids in macro-scale models.

Given the importance of the calibration process of the macro-scale models and to deal with this issue, the main objective of this study is to improve calibration of a macro-scale hydrologic model by proposing a macro-scale modeling framework that links a grid-based hydrologic model (GHM) and a subbasin-based model (i.e. Soil and Water Assessment Tool (SWAT)) for joint simulation and cross-calibration of different models. In this framework, a two-step cross-calibration procedure is also developed. That is, the subbasin-based SWAT model is first calibrated directly by using observed streamflow data. Then, the calibration of the grid-based GHM is indirectly performed by comparing the hydrologic variables simulated by both grid- and subbasin-based models through information transfer between grids and subbasins of the two models. The modeling methodology proposed in this study is tested by an application in the Red River of the North Basin (RRB). Table 1 lists all acronyms used (e.g. GHM, SWAT, and RRB) and their descriptions.

In this study, a GHM was developed for daily macro-scale simulation of hydrologic processes. Figure 1 shows the overall modeling structure, which consists of two main loops: a grid loop and a time loop. Each grid (a basic spatial unit) includes four zones along a profile (i.e. canopy zone, snow zone, surface zone, and vadose zone) (Figure 1). The vadose zone is further divided into three sub-zones (topsoil, root zone, and deep vadose zone). A grid-to-grid (G2G) modeling procedure is implemented at each time step. Within each grid, a zone-to-zone (Z2Z) algorithm (Figure 1) is used to simulate hydrologic processes that vary in a fashion of Z2Z and G2G. The grid-based, bucket-type GHM requires fewer calibration parameters and is efficient in computation. It is particularly suitable for using grid-formatted climatic and hydrologic data. Since most spatially-distributed input data (e.g. rainfall, temperature, land use/land cover (LULC), soil types) may not be available at the same grid scale or resolution, the GHM has a flexible structure that allows users to input gridded data with different resolutions.

The GHM simulates the rainfall-runoff and snowfall-snowmelt processes that are differentiated by a daily mean temperature threshold. In the canopy zone, leaf area index (LAI) is used to partition rainfall into canopy rainfall and direct surface rainfall. In the snow zone, the snowmelt process is simulated by employing the degree-day method in which snowmelt begins when the daily mean temperature is higher than a base temperature. In the surface zone, the curve number (CN) method (USDA 1986) is used to simulate rainfall excess and infiltration, which are linked to the topsoil of the vadose zone. The CN is consistently adjusted based on the water content of the topsoil over the entire simulation period. The CN is also modified by the frozen soil temperature using the approach suggested by Mohammed et al. (2013). In addition, surface depression storage determined by a delineation algorithm is taken into account in the model. To simulate soil moisture content and soil water movement in the vadose zone, the soil moisture accounting (SMA) method (Bennett 1998) is used for the three sub-zones (i.e. topsoil, root zone, and deep vadose zone). The potential evapotranspiration is simulated using the temperature-based method developed by Hargreaves & Samani (1985). Based on the estimated potential evapotranspiration, the actual crop evapotranspiration is then computed using the crop coefficients (Allen et al. 1998). In the modeling of evaporation/evapotranspiration, the water sequentially comes from the canopy zone, surface zone, and vadose zone (including the topsoil and root zone), depending on water availability.

SWAT has been widely used for watershed-scale modeling and assessment of water quantity and quality. In particular, ArcSWAT (Winchell et al. 2013) has been developed and incorporated into ArcGIS (ESRI 2011) as an extension, which provides a user-friendly interface for digital elevation model (DEM)-based watershed delineation, state-of-the-art processing of land use and soil data, HRU (i.e. hydrologic response unit) definition, meteorological data analysis, and preparation of SWAT input files. A variety of modified SWAT models have also been developed to cope with different complex conditions in real applications (e.g. Krysanova & Arnold 2008; Douglas-Mankin et al. 2010). In the current study, ArcSWAT was utilized for model setup.

In SWAT, a watershed is generally divided into a number of subbasins, each of which can be further divided into many HRUs to account for heterogeneity and variability in land covers, soil types, and topography (slopes) in a subbasin. Various hydrologic processes such as surface runoff, infiltration, and groundwater flow are simulated for all HRUs. The HRU-level simulations are then integrated to generate subbasin-level simulations, which are finally routed throughout the drainage network of the entire watershed.

Delineation of a watershed and definition of HRUs require the DEM of the watershed, as well as distributions of LULC and soil types. In addition, meteorological data (e.g. precipitation, temperature, and solar radiation) are needed for model setup. In SWAT, rainfall and snowfall are determined based on a user-defined threshold temperature and utilized to simulate surface runoff, snowmelt, infiltration, and other hydrologic processes. The snowmelt process is simulated using the degree-day method, in which snowmelt is a function of air and snow temperature, melting rate, and areal coverage of snow (Neitsch et al. 2011). Surface runoff and infiltration are modeled by using either the CN method (USDA 1986) or the Green-Ampt (GA) method (Green & Ampt 1911). The former (i.e. the CN method) is used for daily modeling, while the latter (the GA method) is capable of providing sub-daily simulations. Three options are available in SWAT to estimate potential evapotranspiration (Neitsch et al. 2011): Penman–Monteith method (Monteith 1965), Priestley–Taylor method (Priestley & Taylor 1972), and Hargreaves method (Hargreaves & Samani 1985). In SWAT, soil moisture is simulated for all layers and soil water routing is performed, which involves water movement along different pathways such as plant uptake, lateral interflow, and vertical percolation. Simple mass balance equations are used to account for groundwater flow in both shallow and deep aquifers.

The RRB extends from the south of Wahpeton, North Dakota and Breckenridge, Minnesota, crossing the USA–Canada border, to Winnipeg, Manitoba in Canada (Figure 2(a) and 2(b)). The Red River combines with Assiniboine River in Winnipeg and flows into Lake Winnipeg. The RRB has a drainage area of 116,500 km², out of which 103,600 km² is located in the USA (Red River Basin Board 2000), where it drains parts of North Dakota, Minnesota, and South Dakota (Figure 2(b)).

The Red River Valley is the central part of the RRB, which has one of the most fertile soils in North America. According to the National Land Cover Database (NLCD) (Fry et al. 2011), nearly 60% of the RRB within the USA is covered by cultivated croplands, and the developed areas are only 4.6% of the basin area. Lacustrine soil is the dominant soil type (Jin et al. 2008). Also, a major part of the basin belongs to the Prairie Pothole Region (PPR) of North America (Habtezion et al. 2016; Tahmasebi Nasab et al. 2017). The slope of the basin varies from 0.04 to 0.25 m/km (Lin et al. 2015), which makes the RRB very flat and susceptible to frequent flooding. Figure 2(c) depicts the RRB's elevation drop from its headwaters in Wahpeton to a downstream point in Winnipeg. The 65-m elevation difference between these two locations over an 800-km river accentuates the flat topography of the RRB.

The annual mean precipitation and temperature of the basin are 500 mm and 4.3 °C, respectively (Jin et al. 2008; Lin et al. 2015). Precipitation mainly occurs between mid-May and mid-September, partly overlapping the snowmelt period at the beginning and contributing to the spring flooding condition. According to the National Weather Service Flood Stage (Galloway et al. 2011), overtopping of the Red River occurred 47 times during the past 108 years (occurred every year from 1993 to 2011). The key factors that make the RRB prone to flooding are: (1) basin slope, (2) soil type, (3) north-flowing direction, and (4) simultaneous snowmelt and heavy precipitation during the spring (Schwert 2009; Stadnyk et al. 2016).

According to the USGS hydrologic unit code (HUC), river basins in the USA are delineated and mapped by using a hierarchical system of nested hydrologic units at different scales. In the SWAT model in this study, the RRB was divided into 178 HUC-10 subbasins ranging from 159.5 to 1,038.7 km². In the GHM, the study area was divided into 6,351 4 × 4 km grids, based on the availability of main gridded input data (e.g. precipitation and temperature). Table 2 lists the major input data and their sources for setting up the SWAT model and the GHM. The daily gridded (4 × 4 km²) precipitation and minimum and maximum temperature data used in the GHM were obtained from the Parameter-elevation Regressions on Independent Slopes Model (PRISM) (Daly et al. 2008, 2015). In addition, daily precipitation and minimum and maximum temperature data from 12 NOAA weather stations were also obtained and used in the SWAT model. These datasets were compared to ensure that they were consistent in both models.

Soil layers and their physical characteristics play a significant role in simulating vadose zone processes. Soil type distribution was retrieved from the Digital General Soil Map of the United States or STATSGO2 (NRCS 2017). Based on the STATSGO2 data, 25 different soil types were identified in the study area. However, SWAT reclassifies the original soil data and merges similar soil types together to create a modified soil map. Similarly, to use the STATSGO2 data in the GHM, original soil data were reclassified into nine different soil types, indicating that 44% of the study area is covered by loams and only 1.7% of the area is covered by sandy clay loams. Table 3 lists the soil hydraulic property parameters for the nine soil types. These parameters include saturated hydraulic conductivity (K_s), capillary head (h_s), saturated water content (θ_s), field capacity (θ_FC), wilting point (θ_p), residual water content (θ_r), and soil water retention parameters (n and α).

LULC data were obtained from the NLCD (Fry et al. 2011). The original LULC data were reclassified into 15 categories (including open water, open space developed, low intensity developed, medium intensity developed, high intensity developed, barren land, deciduous forest, evergreen forest, mixed forest, shrub/scrub, herbaceous, hay/pasture, cultivated crops, woody wetlands, emergent herbaceous wetlands) that were incorporated into the GHM. Similarly, SWAT reclassified the original LULC data into 15 categories. The LAI and canopy interception storage capacity values for the study area were obtained from Hasan et al. (2014) and Suárez (2005), respectively, and used for different LULC classes. In addition, crop coefficient (Kc) estimates for different LULC classes were extracted from Allen et al. (1998) for three different growth stages, namely, Kc_ini for the initial stage in spring, Kc_mid for the middle stage in summer, and Kc_end for the end stage in fall and winter.

Using the 30-m DEM of the RRB, surface depressions were identified by employing a simple topography-based delineation method suggested by Planchon & Darboux (2002). As a result, 39 depressional zones were included in the GHM. In the delineation process, depressions and their storage values were determined for each zone separately. Considering the existence of spurious pits or artifacts of the DEM, a special effort was made to remove those spurious pits and the adjusted depressional data were generated. In the GHM, the maximum depth of depression storage over a grid was used for each depressional zone.

Figure 3 depicts the required steps for the joint simulation and cross-calibration of the GHM. After all input data are prepared (as described above), the basin is delineated into subbasins and HRUs for the SWAT modeling and grids for the GHM modeling. Then, both models are applied to the RRB and auto-calibration is performed for the SWAT model (Figure 3). Eventually, using the calibrated SWAT results, the GHM is cross-calibrated (Figure 3). The cross-calibration process and the selected calibration parameters for both models are detailed as follows.

Note that the RRB region has suffered from frequent flooding and drought events in the past decade. Over the study period, historical floods occurred in 2008, 2009, and 2011, followed by an extreme drought event in 2012 (Lin et al. 2015). However, this study focused on demonstrating the proposed modeling framework and strategy, instead of examining and showing the capabilities of the GHM to simulate hydrologic extremes. Thus, a 5-year period from 2003 to 2007 (1,826 daily time steps) was selected for the GHM modeling.

Calibration of a macro-scale grid-based model is a challenge. Satellite data can be utilized to calibrate different hydrologic processes (e.g. Lohmann et al. 1996; Güntner 2008); however, the temporal and spatial resolutions of such data are often limited. The current study involved two distinct modeling approaches: grid-based GHM and subbasin-based SWAT. It should be noted that this study was not intended to show the superiority of these two methods over each other. Instead, a cross-calibration strategy was implemented in this study to take full advantage of both modeling methods. Specifically, auto calibration of the subbasin-based SWAT model was performed by using SWAT-CUP (Abbaspour 2013) and comparing the SWAT simulations against the observed discharges at the 16 gauging stations. The information from the well-calibrated SWAT was then transferred to the GHM for its calibration. To do so, a cross-model comparison scheme (Figure 4) was designed to relate the different hydrologic processes simulated by the subbasin-based SWAT model to those by the grid-based GHM. As shown in Figure 4, the 178 HUC-10 subbasins in the SWAT model were overlaid with the 6,351 GHM grids, from which the hydrologic quantities (e.g. snowmelt, snowpack water equivalent (SWE), surface runoff, and infiltration) simulated by SWAT were extracted and used to calibrate the GHM. Typically, each subbasin included a number of GHM grids (Figure 4). The major calibration parameters for the GHM included CN, snowmelt threshold temperature, and initial abstraction coefficient.

In the cross-calibration process, the grid-averaged CN values in the GHM were adjusted to be consistent with the corresponding subbasin CN value in the calibrated SWAT model. This process was repeated for all subbasins and their associated grids in the two models. When simulating snowpack and snowmelt, SWAT used several parameters to account for various factors such as surface topography, snow temperature, and drifting, while GHM adjusted the threshold temperature of snowmelt for each grid to account for heterogeneity. In contrast, SWAT used one value for the snowmelt threshold temperature for the entire basin (see Table 4). The initial abstraction (I_a) in the CN method is generally estimated by I_a = αS, in which S is the potential retention and the default initial abstraction coefficient α = 0.2. Since the GHM had separate simulations for interception of the canopy zone and the impact of surface depressions, it was logical to calibrate α, instead of using a fixed value of 0.2.

Using the streamflow data observed at the 16 USGS gauging stations in the RRB, 14 SWAT parameters were calibrated (Table 4). Figure 5 shows the calibrated monthly average flows at the Doran (a few miles south to Wahpeton, Figure 2), Fargo, Grand Forks, and Drayton gauging stations (Figure 2) along the main channel of the Red River. Overall, SWAT provided reasonable simulations for low flow conditions at the four stations (Figure 5(a)–5(d)). For high flow conditions, the simulated discharges also matched the observed data very well in terms of both distribution and magnitude for the years before 2009, except for over-predicting three snowmelt events at the upstream Doran station in 2005–2007 (Figure 5(a)). For 2009–2011, however, SWAT failed to capture the peaks caused by the historical extreme floods and underestimated the high flows at the four stations during the flooding periods (Figure 5(a)–5(d)). In addition, a severe drought occurred in the RRB in 2012, resulting in very low streamflow levels. As shown in Figure 5, the simulated monthly average discharges for 2012 were consistently higher than the observed data.

To evaluate the performance of the SWAT model, NSE and PBIAS were also calculated for both monthly average and daily discharges at the four gauging stations along the main channel of the Red River (Table 5). According to Moriasi et al. (2007, 2015), NSE values greater than 0.50 and PBIAS values between −25 and +25% indicate an acceptable level of model performance. As shown in Table 5, both NSE and PBIAS demonstrated the good and acceptable performance of the SWAT model for the RRB. It is worth noting that the SWAT model tended to overestimate the streamflow at Doran whereas it slightly underestimated that at the other three stations (Table 5 and Figure 5). Nevertheless, it can be concluded that the calibrated SWAT model was able to provide acceptable simulations of streamflow in the Red River. In addition to these four stations, NSE and PBIAS were calculated for another 12 stations within the RRB. The results for those 12 stations suggested that the PBIAS values for all stations fell into the acceptable range from −25 to +25%. In addition, eight out of 12 gauging stations had an NSE of greater than 0.6, with the highest NSE of 0.83 for Sand Hill River at Climax, MN.

The grid-based GHM was cross-calibrated based on the well-calibrated, subbasin-based SWAT simulations. Figure 6 shows the comparisons of four monthly hydrologic variables, including snowmelt, SWE, infiltration, and surface runoff, simulated by the subbasin-based SWAT and the grid-based GHM for the RRB. A visual comparison between the SWAT and GHM results in Figure 6(a)–6(d) indicates that both magnitude and timing of the peaks from the two models are in good agreement. For example, the snowmelt peaks simulated by the SWAT and the GHM for March 2007 were 44.54 and 46.91 mm, respectively (Figure 6(a)) and the snowmelt peaks for other years were also similar. For SWE, however, the SWAT peaks occurred earlier than the GHM peaks (Figure 6(b)). This can be attributed to the underlying methodologies used for computing SWE in the two models. Specifically, SWAT utilized an exponential equation to account for the impacts of drifting, shading, and topography on the snow cover within each subbasin, while GHM did not consider these factors and only employed a simple threshold temperature-based approach. The infiltration and surface runoff simulated by SWAT and GHM were also comparable (Figure 6(c) and 6(d)).

Table 6 provides a summary of the four main hydrologic variables (i.e. snowmelt, SWE, infiltration, and surface runoff) simulated by the SWAT and the GHM for the RRB. For the 5-year simulation period, the differences of the yearly values were less than 17%, except for the SWE in 2003–2004 and the snowmelt in 2003. The differences between the annual averages of the four hydrologic variables ranged from 3.87 to 16.43%.

Figure 7 depicts the spatial distributions of rainfall, surface runoff, and infiltration over the RRB on July 7, 2004. The average rainfall over the entire RRB was 13 mm and the southern and eastern parts of the basin received more rainfall than the northwestern basin (Figure 7(a)). The average surface runoff over the RRB was 1.62 mm on July 7, 2004, and more surface runoff was generated in the middle and southern parts of the basin (Figure 7(b)). This trend can be directly attributed to the spatial distribution of CN values controlled by soil types and LULC in the RRB. For example, the dominant silty clay soil and the high percentage of impervious land surfaces in the area close to Fargo led to a much higher amount of surface runoff (24.6 mm). Figure 7(c) illustrates the distribution of infiltration over the RRB and the average infiltration over the basin on July 7, 2004 was 10.84 mm. A visual comparison between Figure 7(a) and 7(c) suggests that rainfall and infiltration followed a similar spatial distribution pattern where both decreased from southeast to northwest.

Snow accumulation and snowmelt are prominent hydrologic processes in the RRB in winters. Figure 8(a)–8(f) exemplify the spatial distributions of snowpack SWE and snowmelt over the RRB for three selected dates: January 22, March 22, and March 29, 2005. Figure 8(a)–8(c) display a progressively decreasing trend of the snowpack SWE from southwest to northeast on the three selected days in 2005. The average snowpack SWE values were 39.08, 39.14, and 22.74 mm on January 22, March 22, and March 29, respectively. Understandably, these variations were mainly affected by the variations in temperature. The average temperatures on March 22 varied from +2.86 °C in Wahpeton (south) to −3.12 °C in Pembina (north), resulting in a reduction in snowpack depth in the southern part of the basin. In addition, Figure 8(d)–8(f) reveal the snowmelt trend in the RRB. The snow-melting process started from the southwestern side of the basin and moved towards the northeastern side (Figure 8(d)–8(f)). The average snowmelt over the RRB on January 22, March 22, and March 29 was 0, 0.19, and 7.89 mm, respectively.

In addition to the surface hydrologic processes, the GHM also provided simulation results for subsurface processes. Figure 9(a) and 9(b) show the spatial distribution of topsoil water storage (TWS) in the RRB on January 1, 2006 and April 2, 2006, respectively. The average TWS for the entire basin on January 1, 2006 was 67.44 mm, whereas it increased to 78.89 mm on April 2, 2006. This trend also occurred for all other years, which can be mainly attributed to the spring snowmelt as a direct result of higher temperature in March and April. Figure 9(a) and 9(b) also reveal a spatial pattern which affected surface runoff and infiltration. A closer look at the red box in Figure 9(b) indicates that the middle part of the basin along the Red River Valley accounted for the higher TWS values due to the silty clay and silty clay loam soils in the Red River valley. This unique soil distribution also impacted infiltration and surface runoff. Figure 9(c) and 9(d) highlight the spatially distributed surface runoff and infiltration on April 2, 2006. It can be observed that more runoff was generated, and less water infiltrated into soils along the Red River Valley. The average surface runoff for the entire basin on April 2, 2006 was 0.48 mm; however, the surface runoff values in the areas of Wahpeton and Grand Forks, dominated by silty clay and silty clay loam soils, were as high as 4.23 and 2.67 mm, respectively.

The GHM can also be used to simulate the temporal distributions of hydrologic processes for the entire basin to reveal basin-wide trends. Figure 10(a) shows the basin-wide average snowmelt and SWE, indicating that most snowmelt peaks occurred right after the corresponding SWE peaks. For example, the maximum average SWE occurred on March 25, 2006, which was followed by the maximum average snowmelt on March 30, 2006. Moreover, Figure 10(b) shows the variations in the average infiltration and the soil water content (SWC) of the root zone. The SWC peaks coincided with the infiltration peaks. SWC varied within a range of 0.36–0.42 cm³/cm³.

In this study, a macro-scale modeling framework was proposed for joint simulation and cross-calibration of a new GHM and the widely-used SWAT model. This modeling strategy was applied to the RRB. A two-step cross-calibration procedure was developed, which involved (1) direct calibration of the subbasin-based SWAT model using the observed streamflow data and (2) indirect calibration of the GHM by transferring the well-calibrated, subbasin-based SWAT simulations to the GHM grids. The simulation results from the two models were analyzed and compared. The performance of the SWAT model was quantitatively assessed by the NSE and PBIAS for the simulated and observed hydrographs. The NSE values were greater than 0.5 and the PBIAS values were within a range of ±25%, showing acceptable model performance.

This study presented a new way of calibrating a macro-scale GHM by taking full advantage of two different modeling structures and approaches. The simulation results demonstrated the feasibility of the two-step cross-calibration procedure that combined the subbasin-based and grid-based modeling methodologies, as well as the capability of the GHM in grid-based simulations of surface runoff, infiltration, snowmelt, and other hydrologic processes. Some discrepancies in the simulations from the subbasin-based SWAT and grid-based GHM can be mainly attributed to the differences in the modeling methods used in the two models for specific hydrologic processes.

The gridded simulation results from the calibrated macro-scale GHM can be used to better understand the localized (up to a grid) hydrologic processes (e.g. overland flow generation, infiltration, depression storage, snowpack, snowmelt, soil moisture condition, and subsurface flow), and further analyze their variations across space and time, as well as their spatio-temporal trends at a macro scale. The joint modeling framework, grid-based hydrologic modeling method, and the two-step cross-calibration procedure developed in this study can be readily applied to other macro-scale basins. Although this cross-calibration approach makes the calibration of the macro-scale, grid-based model much easier, there are some limitations. Caution should be taken if this method is used for very large grids that might not fit into subbasin boundaries. In addition, the uncertainties related to the subbasin-based SWAT model can be transferred to the GHM.

This material is based upon work supported by the National Science Foundation under Grant No. NSF EPSCoR Award IIA-1355466. The North Dakota Water Resources Research Institute also provided partial financial support in the form of graduate fellowships.