A subbasin-based framework to represent land surface processes in an Earth system model

Realistically representing spatial heterogeneity and lateral land surface processes within and between modeling units in Earth system models is important because of their implications to surface energy and water exchanges. The traditional approach of using regular grids as computational units in land surface models may lead to inadequate representation of subgrid heterogeneity and lateral movements of water, energy and carbon fluxes. Here a subbasin-based framework is introduced in the Community Land Model (CLM), which is the land component of the Community Earth System Model (CESM). Local processes are represented in each subbasin on a pseudo-grid matrix with no significant modifications to the existing CLM modeling structure. Lateral routing of water within and between subbasins is simulated with the subbasin version of a recently developed physically based routing model, Model for Scale Adaptive River Transport (MOSART). The framework is implemented in two topographically and climatically contrasting regions of the US: the Pacific Northwest and the Midwest. The relative merits of this modeling framework, with greater emphasis on scalability (i.e., ability to perform consistently across spatial resolutions) in streamflow simulation compared to the grid-based modeling framework are investigated by performing simulations at 0.125 , 0.25, 0.5, and 1 spatial resolutions. Comparison of the two frameworks at the finest spatial resolution showed that a small difference between the averaged forcing could lead to a larger difference in the simulated runoff and streamflow because of nonlinear processes. More systematic comparisons conducted using statistical metrics calculated between each coarse resolution and the corresponding 0.125 -resolution simulations showed superior scalability in simulating both peak and mean streamflow for the subbasin based over the grid-based modeling framework. Scalability advantages are driven by a combination of improved consistency in runoff generation and the routing processes across spatial resolutions.


Introduction
Land surface processes play an important role in the exchange of energy, moisture, and biogeochemical fluxes in the Earth system.It has been widely recognized that the development of a planetary boundary layer, initiation of shallow and deep convections, and cloud formation and precipitation are sensitive to spatial heterogeneity of hydrologic state variables such as soil water distribution and snow cover (Chen and Avissar, 1994;Quinn et al., 1995;Leung andGhan, 1995, 1998;Pielke, 2001).Hence parameterizations of spatial heterogeneity in land surface models must account for the lateral redistribution of water that subsequently affects the simulated water and energy exchange with the atmosphere (Liang et al., 1996;Niu et al., 2005;Huang et al., 2008;Li et al., 2011).
The most common practice in land surface modeling to resolve spatial heterogeneity is to divide a study domain into regular latitude-longitude or other quasi-uniform rectangular grids for computational convenience.However, Published by Copernicus Publications on behalf of the European Geosciences Union.subbasin-based representation, i.e., dividing the study domain into irregular subbasins, offers distinct advantages over the traditional grid-based representation.First, the subbasinbased representation follows the natural topographic divides and river network structure that strongly govern hydrological processes such as surface runoff and streamflow.Figure 1 shows an example of the Columbia River basin (CRB) in the grid-and subbasin-based representations overlaid by a river network.As highlighted in Fig. 1 by the dashed red line, a single grid cell in the grid-based representation very often crosses over several channel reaches, which leads to great difficulties in parameterizing runoff routing (Guo et al., 2004;Wu et al., 2011Wu et al., , 2012;;Wen et al., 2012).Second, a single grid cell in the grid-based representation also often encompasses areas from several subbasins, which challenges the conceptual basis of runoff schemes such as the TOPMODEL (Topographic model) formulation in which topographic variations at catchment scale are of primary importance to runoff generation (Beven, 1997;Huang and Liang, 2006).For example, the key parameter of such runoff schemes, the topographic index and its statistical distribution (within a spatial unit), essentially measures the accumulated contributing areas from natural divides to the valley and then to the channels, but its physical meaning is confounded in a grid-based representation.Third, at very high resolution, the grid-based approach must be modified to account for lateral redistribution of water from neighboring grid cells, which becomes important in determining the soil moisture states, but in a subbasin-based approach, such requirements are to some extent relaxed because surface water is not redistributed across topographic boundaries of the subbasins.Last but not least, the subbasin-based representation provides a bridge that may enhance co-development of the hydrologic component in land surface models with contributions from the land surface modeling and hydrologic science communities because the latter has mostly focused their theoretical and modeling advances on catchment or subbasin scales.
There have been a few attempts to implement subbasinbased representation in land surface models (Koster et al., 2000;Goteti et al., 2008).Koster et al. (2000) was among the first to adopt this approach to improve parameterizations of spatial variability of soil water in land surface and earth system models.Their study focused more on representing soil moisture; while surface water movements and storages, which are closely related to soil moisture, were not discussed explicitly.Goteti et al. (2008) developed a catchment-based hydrologic and routing modeling system (CHARMS) with explicit treatment of surface water bodies and storages.Despite the important advances, their approach has several limitations in that (1) routing was essentially based on the unithydrograph approach so channel velocity and depth were not directly linked to the discharge; and (2) model inputs including forcing and land surface parameters were remapped, or disaggregated, from the default CLM input data set provided by the National Center for Atmospheric Research (NCAR) at a resolution of 0.5 • (Lawrence and Chase, 2007), which is too coarse compared to the average size of subbasins.Lastly, although catchments were used as the fundamental modeling units in Koster et al. (2000) and Goteti et al. (2008), the subbasin representation has not been systematically compared with the grid-based representation to evaluate its potential advantages in land surface modeling.
The significance of scaling and scale interactions in hydrologic model predictions has been well documented via numerous field and modeling studies over the last several decades (Bloschl and Sivapalan, 1995;Sivapalan and Kalma, 1995;Koren et al., 1999;Schulze, 2000;Kirkby, 2001).However, most modeling studies mainly focused on examining the effects of different spatial resolutions of input data on hydrologic model predictions (Wolock and Price, 1994;Haddeland et al., 2002;Sridhar et al., 2003;Boone et al., 2004;Cerdan et al., 2004), and on identifying a critical spatial resolution for optimal model predictions (Wood et al., 1988;Famiglietti and Wood, 1994;Bruneau et al., 1995;Woods et al., 1995;Liang et al., 2004;Shrestha et al., 2006).While these are important, given the significance of scalable modeling approaches for providing reliable hydrologic predictions under changing climate and environmental conditions, examining modeling approaches from their scalability perspective is very much needed.
Following Koster et al. (2000) and Goteti et al. (2008), we implement another attempt on the use of subbasin-based representation in a land surface model and systematically compared it with the grid-based representation.Comparison on runoff generation revealed that the subbasin-based approach is more consistent across multiple spatial resolutions compared to the standard grid-based land surface modeling approach (Tesfa et al., 2014).The improved scalability of runoff was attributed to scalability of snowfall/rainfall partitioning related to air temperature and surface elevation, and scalability of a topographic parameter that influences rain driven saturated surface runoff.
In this study, we couple a land surface model with a physically based routing model, Model for Scale Adaptive River Transport (MOSART), both using a subbasin-based representation, as a land surface modeling framework.We choose the Community Land Model version 4 (CLM4) (Lawrence et al., 2011) as the basis for our development because it has a large user community and its use of the TOPMODEL approach for parameterizing runoff may allow it to take more advantages of the subbasin-based representation.For brevity, hereafter we denote the subbasin-based representation of CLM as SCLM, while CLM strictly refers to the grid-based representation of CLM4, which, after coupling with the routing model, are denoted as SCLM-MOSART and CLM-MOSART, respectively.Without parameter calibration, we systematically investigated the relative merits of the subbasin-based modeling framework compared to the standard grid-based modeling framework on streamflow simulation, which depends on both runoff generation simulated by the land surface model and river routing represented by MOSART.More specifically we compared the two modeling frameworks (1) to investigate how they simulate streamflow across multiple spatial scales, (2) to determine if the scalability advantage of the subbasin-based approach in simulating runoff is preserved in the coupled SCLM-MOSART framework, and (3) to determine if the understanding of differences in streamflow simulation between the two modeling frameworks is generalizable to other regions.In this paper, we first describe the implementation of the subbasin-based modeling framework, including coupling with the physically based routing model in Sect. 2. Section 3 describes the experimental design, including forcings and land surface parameters for each modeling framework.This is followed by analysis to compare streamflow simulations with observations and between the two modeling approaches in Sect. 4. To the best of our knowledge, this is the first attempt to document comparisons between the two modeling frameworks on their scalability in terms of streamflow simulation.Section 5 closes with summary and conclusions.

Implementation of SCLM-MOSART
CLM applications at regional or global scales involve a large number of computational units.A customized parallel algorithm is embedded in CLM to facilitate such large simulations on clustering computers.To accommodate this parallel algorithm, all CLM units are organized into a twodimensional matrix with each node containing a single grid cell.To take advantage of the parallel algorithm without significantly modifying the original computational structure of CLM, the subbasins of a study domain are also organized into a two-dimensional matrix of subbasins, each treated as a single node, to be consistent with the grid-based representation.Using this subbasin-based representation, all the local land surface processes such as water and energy transfer between the land surface and the atmosphere, as well as subgrid (or within-subbasin) processes such as runoff generation, are represented assuming each subbasin as a pseudo grid cell without significantly modifying the existing CLM modeling structure.Note that the latest public versions of the Community Climate System model and the Community Earth System model (i.e., CCSM4 and CESM1) includes a suite of new coupling capabilities in the CPL7 coupler that allow coupling Earth system components configured on unstructured grids or subbasins (Craig et al., 2012).The SCLM structure has been tested in a preliminary manner at small watersheds (Li et al., 2011;Huang et al., 2013), but without invoking the routing component because the river transport model (RTM) embedded in CLM4 and its supporting parameters are not intended for the subbasin-based representation.The next section introduces the coupling of a new river routing model to SCLM.
The above matrix representation of CLM grids does not distort the real spatial arrangements among the grid cells, i.e., grid cells that are neighbors in a model domain of a region are still neighbors in the matrix because the grids have regular structure (e.g., each rectangular grid has exactly eight neighbors).For SCLM, each subbasin can have different numbers of neighboring subbasins, so the (2-D) matrix structure cannot reflect the real spatial arrangements or linkages among the subbasins.We therefore impose an extra indexing system by assigning a unique index to each individual subbasin.The linkages between the subbasins, i.e., upstream/downstream relationships and other parameters needed for the runoff routing are all preprocessed and identified by their indices defined using the same indexing system and stored in the input data sets as a separate geographical location layer.
MOSART is a large-scale routing model recently developed with explicit treatment of routing processes at hillslope, across tributaries (within a spatial unit) and through the main channel (for more details please refer to Li et al., 2013).MOSART can be applied at both grid-and subbasinbased representations.Li et al. (2013) describes the concept and framework of MOSART and evaluation of its grid-based representation in the US Pacific Northwest region at multiple spatial resolutions.In this work, MOSART is modified to follow the subbasin structure for direct coupling with SCLM.The surface and subsurface runoff produced from the SCLM units is fed into the MOSART units by a one-to-one mapping based on the indexing system described above.MOSART then routes the runoff within and between the subbasins all the way to the ocean or basin outlet.

Experimental design
In this study, we applied the two modeling frameworks (CLM-MOSART and SCLM-MOSART) and performed detailed analysis over the topographically diverse region of the Columbia River Basin (CRB), which is located in the US Pacific Northwest (Fig. 1).The CRB region receives the majority of its precipitation during the cold season with its hydrology dominated by snow accumulation and melting.It encompasses both mountainous and low-lying regions (Fig. 1a).The mountains are characterized by low temperature and higher precipitation dominated by snowfall, while the lower elevation regions have higher temperature and lower precipitation mainly as rainfall.To determine if the understanding on differences between the two modeling frameworks is generalizable to other regions, we also apply the models to the US Midwest (MW) study domain, which encompasses the Missouri, Upper Mississippi, and Ohio river basins.Compared to CRB, the MW region is dominated by a less complex topography and the majority of precipitation occurs as rainfall.Located further inland, the MW has a colder winter with less precipitation, but it receives more convective rainfall during summer fed by moisture transported from the Gulf of Mexico.Analysis over the MW region is limited to key aspects identified from analysis over the CRB.
To compare the two modeling frameworks, each is applied at four spatial resolutions (0.125 • , 0.25 • , 0.5 • , and 1 • ) over the CRB and MW.The grid-based framework is applied directly at 0.125 • , 0.25 • , 0.5 • and 1 • grid resolutions.For a fair comparison between the two modeling frameworks, the average sizes of the subbasins delineated in this study are chosen to be equivalent to a 0.125 • , 0.25 • , 0.5 • and 1 • lat/long grid.Both modeling frameworks are driven by the same meteorological forcing and land surface parameters.Details of the subbasin delineation, model inputs, and analysis methods are provided in the subsections below.

Subbasin delineation
To set up SCLM, we utilize the 90 m digital elevation model (DEM) and the 15 arcsec river networks from the Hydrological data and maps based on SHuttle Elevation Derivatives (HydroSHEDS) (Lehner et al., 2008).Although DEMs at resolutions of 30 m or higher over the study area can be obtained from the United States Geological Survey (USGS) database, the goal of our study is to develop a framework suitable for SCLM applications worldwide.Therefore, a global database (i.e., HydroSHEDS) is used in this study.The 15 arcsec river network is reconciled with the 90 m DEM over each study domain for hydrologic conditioning to ensure a consistent delineation of the river network with HydroSHEDS.Using ArcSWAT (soil and water assessment tool; Neitsch et al., 2005), we delineate subbasins within CRB and MW, as well as a river network consistent with the subbasin boundaries at four spatial resolutions, using the hydrologically conditioned DEMs as inputs.For comparison with the grid-based application of CLM4 at 0.125 • , 0.25 • , 0.5 • and 1 • resolutions, the threshold area for the subbasin delineation is adjusted iteratively until the average subbasin size is roughly equivalent to the 0.125 • , 0.25 • , 0.5 • and 1 • grids, respectively.We eventually obtain 5999, 1139, 299 and 75 subbasins for CRB and 18 681, 4019, 1031, and 273 subbasins for MW with average drainage areas equivalent to 0.125 • (140 km 2 ), 0.25 • (770 km 2 ), 0.5 • (3000 km 2 ) and 1 • (1200 km 2 ), respectively.These subbasins are then organized into a 77×78, 38×30, 20×15, and 10×8 matrices for CRB and 137 × 137, 67 × 60, 37 × 28, and 18 × 16 matrices for MW, respectively, where extra grid cells are masked out as nonland cells and therefore excluded from the simulations.Given that subbasins could only be defined over land, any pseudo grid cell in the matrix is either 100 % or 0 % land.This is different from the grid-based applications, in which a single grid can be occupied fractionally by land or ocean.
ArcSWAT also provides the subbasin parameters needed for the routing model (MOSART) including subbasin upstream/downstream dependence information, accumulated contributing area and other channel parameters such as slope and length.The bankfull channel width and depth values are derived based on the empirical hydraulic geometry relationships estimated in Li et al. (2013), consistently for all SCLM-MOSART and CLM-MOSART setups.

USGS stream gauges
The analyses performed in this study require extraction of the USGS stream gauges located within the boundary of each study domain.For the grid-based framework, the stream gauges are co-registered to the dominant river tracing (DRT) flow accumulation grids (Wu et al., 2011(Wu et al., , 2012)).Stream gauges were selected for analyses according to the following criteria: (1) the relative difference between the contributing area determined from the DRT data set and the observed drainage area from USGS is within 20 % across all spatial resolutions and (2) the observed drainage areas are greater than 1500 km 2 (Fig. 1).For the subbasin-based framework, the stream gauges are snapped to the nearest stream network, and in this way the upstream drainage areas are easily preserved across different spatial resolutions.ArcGIS is used to facilitate mapping of the stream gauges to the grids/subbasins across different spatial resolutions, and extraction of input data and simulation results at the stream gauges and the regions draining to the stream gauges.

Model inputs
The meteorological forcing in this study was extracted from the phase-two North America Land Data Assimilation System (NLDAS-2) at an hourly time step from 1979 to 2008 (Xia et al., 2012), including precipitation, shortwave and long-wave radiation, air temperature, humidity, surface pressure and wind speed at 0.125 • resolution derived from the 32 km resolution 3-hourly North American Regional Reanalysis (NARR).For grid-based applications, the NLDAS-2 forcing data are either applied directly at 0.125 • resolution or spatially aggregated to the corresponding coarser resolutions.For SCLM, an area-average algorithm is applied to remap the NLDAS-2 forcing to the subbasins defined by their boundaries at each spatial resolution.That is, the algorithm computes the value of each meteorological variable in a subbasin as the average of the corresponding variable from all the 0.125 • grid cells that intersect with the subbasin weighted by the overlapping areas.ArcGIS is used to link the subbasins to the intersecting grids and to compute the fractional areas of the intersecting grids.Figure 2 shows the spatial distributions of annual mean precipitation and surface temperature of CLM and SCLM in the CRB domain at the finest (0.125 • ) resolution.It can be seen from the figure that the differences between the two representations are very small.
Land surface parameters and leaf area index (LAI) are derived from the global land parameter data set developed by Ke et al. (2012) at 0.05 • resolution based on the most recent MODIS (Moderate Resolution Imaging Spectroradiometer) land cover and improved MODIS LAI products.Soil texture is generated based on a hybrid of 30 arcsec State Soil Geographic Database (STATSGO, now referred to as the US General Soil Map) (for CONUS) and 5 min Food and Agriculture Organization (outside CONUS) 16-category two-layer soil-type data (Chen et al., 2007;Miller and White, 1998).The two-layer soil-type data is then converted to composition of clay and sand (Cosby et al., 1984;Dai et al., 2003) within each 30 arcsec grid cells and interpolated to 10 vertical layers down to 3.8 m depth.Soil depth is a notoriously difficult input parameter for hydrology modeling.It is thus greatly simplified in CLM4 by assuming that soil depth has a globally uniform value of 3.8 m (Oleson et al., 2010) in both modeling frameworks.Even though it is a bold assumption, it is typical in the field of land surface modeling due to lack of global soil-depth data.Other land surface parameters such as soil color and soil organic matter are derived from the default 0.5 • CLM4 global input data set provided by NCAR.For the grid-based CLM simulations, the CLM4 preprocessing package (Oleson et al., 2010) is used at each spatial resolution in this study.
The SCLM soil, land cover, and vegetation parameters are derived by overlaying the subbasin boundaries of each spatial  resolution over the grids as described above.Similar to the forcing parameters, ArcGIS is used to link the subbasins to the grids and calculate the area weights of the intersecting grids.Consistent with the CLM4 preprocessing package (Oleson et al., 2010), soil properties such as percent clay, percent sand etc. are calculated using an area-dominant algorithm, where each parameter value for the subbasins is assigned to the value covering the largest fraction of the subbasin.Land cover characteristics and plant functional types (PFTs) for each subbasin are determined using the areaaverage algorithm described above for the atmospheric forcing.The LAI and stem area index (SAI) parameters for each subbasin are calculated using a PFT-weighted area-average algorithm.Between the two representations, the differences of most land surface parameters are small due to the high resolution of source data.But this is not the case for soil color since the resolution of source data is 0.5 • , which is much coarser (figure not shown).
Topography exerts an important control on the lateral redistribution of soil moisture and therefore affects the generation of saturation excess runoff and baseflow.In CLM4, two important parameters are used to represent the effects of topography on runoff generation, the maximum saturated area fraction within a spatial unit, f max , and an empirical coefficient to describe the variation of actual saturated area fraction with the groundwater table, C s (Hou et al., 2012;Huang et al., 2013;Li et al., 2011;Niu et al., 2005).In this work f max is derived following the algorithm described in Niu et al. (2005) and Niu and Yang (2006) for both SCLM and CLM simulations.Based on the DEMs, compound topographic indices (CTIs) are first derived following the definition used in TOP-MODEL (Beven;1997;Quinn et al., 1995) using ArcGIS.The CTIs are then fitted to a distribution, within each subbasin or grid, to estimate the topography-relevant hydrologic parameters (i.e., f max and C s ).In CLM4 and its previous versions (e.g., CLM3, CLM3.5), these parameters are derived from coarse resolution (e.g., 1 km) DEMs (Niu et al., 2005) due to the lack of higher-resolution DEMs with a global domain.However, as discussed in our previous study (Li et al., 2011), the estimation of these parameters using 1 km DEMs is problematic due to its inconsistency with hydrology theory.Interested readers are referred to Li et al. (2011) for details.With the newly available HydroSHEDS (Lehner et al., 2008) database, 90 m DEM data are now available globally.We therefore estimated f max values using the 90 m DEM from HydroSHEDS.

Methods of analyses
For a fair comparison of the two modeling frameworks, all simulations are driven by the same meteorological forcing (NLDAS-2 1979(NLDAS-2 -2008) ) and land surface parameters, generated using the same methods described above, and spun up until the state variables such as soil moisture and temperature reached equilibrium states.The relative merits of CLM-MOSART and SCLM-MOSART land surface modeling frameworks in streamflow simulations are then investigated using streamflow simulated after the spin up.The two modeling frameworks are compared using two approaches.
In the first comparison, we investigate how the two modeling frameworks simulate streamflow across multiple spatial scales (defined as the size of upstream drainage area of gauge station) compared to the naturalized streamflow data from the Surface Hydrology Group, University of Washington (http://www.hydro.washington.edu/2860).The sources of differences in simulated streamflow at the highest resolution (0.125 • ) between the two modeling frameworks are explored in different climatic regions of the CRB and classified based on rainfall/snowfall fractions.The USGS stream gauges are also classified based on the dominant climate regime in their catchment areas for insight on the role of climate in the scalability differences of the modeling frameworks.In the second comparison, the two modeling frameworks are evaluated for their scalability across multiple spatial resolutions in both CRB and MW.Scalability is measured by the ability of the models to produce simulations that asymptotically approach the high-resolution simulations as model resolution increases.Model scalability is thus demonstrated by using the simulations performed at 0.125 • resolutions (CLM0125-MOSART and SCLM0125-MOSART) as the "reference" solution for comparison with simulations performed at increasingly coarser spatial resolutions.
In both comparisons, we compare the streamflow simulated by the two modeling frameworks across multiple spatial scales at the USGS stream gauges selected to represent major rivers of the CRB.Absolute differences in specific peak flows (ADP) (long-term average peak streamflow normalized by the corresponding contributing area) are calculated at all the stream gauges of the CRB and selected using the procedure described in Sect.3.2 between each coarse-resolution simulation and the corresponding finest-resolution (0.125 • ) simulation.The ADP values are used to evaluate the scalability of the two modeling frameworks in simulating peak flow at multiple temporal scales.Furthermore, Nash-Sutcliffe efficiency (NSE) values are calculated at the same stream gauges between the streamflows simulated by each coarse resolution and the corresponding reference (finest resolution) simulation as shown below: where, F r , F c , and F rave are streamflow simulated by the reference simulation, coarse-resolution simulation, and average streamflow from the reference solution, respectively.To understand the role of contributing area on model scalability, the NSE values are also calculated using streamflow values normalized by the respective contributing area.A nonparametric statistical significance test is employed using the CRB stream gauges to evaluate the significance of the scalability differences between the two modeling frameworks.Finally, NSE values are calculated at the USGS stream gauges of the MW domain that are selected based on the procedure described in Sect.3.2 (Fig. 1) to investigate scalability of the two modeling frameworks at the spatial resolutions that showed statistically significant differences in scalability in the CRB study domain.The comparison at the MW domain is further extended to the MW domain basins (MSRB, UMRB and OHRB, respectively, Missouri, Upper Mississippi and Ohio river basins) to gain more insight on any differences at a regional scale.Results from these analyses are discussed in the following section.

Runoff and streamflow at the finest spatial resolution
Figure 3 compares how well SCLM-MOSART and CLM-MOSART applied at the finest resolution of 0.125 • simulate streamflow at a number of major USGS stream gauges representative of the CRB study domain (see Fig. 1).There is a systematic phase shift between the streamflow simulated by the two frameworks, with SCLM-MOSART producing higher streamflow in January-March and lower streamflow in April-June.Also shown in Fig. 3 is the naturalized streamflow with human influences such as reservoir operation removed.Due to the availability of the naturalized streamflow data, the analysis period is chosen as October 1979-September 1989 for Fig. 3. Compared to the Since the northern part is wetter and provides more runoff than the southern part (e.g., streamflow at CHIEF is much higher than HCANY although they have comparable drainage area), both SCLM-MOSART and CLM-MOSART slightly underestimate streamflow at the DALLE stream gauge, which is downstream of the confluence where the two parts join each other.Also, at most stream gauges, both SCLM-MOSART and CLM-MOSART produce peak flow 1 month earlier than the naturalized streamflow.This can be attributed to the parameterization of snow processes in CLM that very often leads to earlier snowmelt (Wang et al., 2008;Li et al., 2011).Due to the earlier snowmelt, one may conclude that SCLM-MOSART does not necessarily perform better than CLM-MOSART in simulating streamflow.One reason why SCLM-MOSART is not performing better than CLM-MOSART is that the runoff simulations in SCLM and CLM are both governed by the same set of hydrological formulations and parameters originally calibrated for CLM based on the grid-based configuration at the global scale.If reproducing the observed streamflow is the target, a fair comparison between SCLM-MOSART and CLM-MOSART should be conducted with separate parameter calibration for each.The effective and meaningful parameter calibration of SCLM and CLM, however, is itself challenging particularly over large regions.This has been a topic of research in separate studies (Huang et al., 2013;Sun et al., 2013), and is beyond the scope of this study.Here, our main objective is to investigate the differences between the two modeling frameworks in streamflow simulation caused purely by different approaches to delineating the fundamental spatial units without changes in model parameters or adjustments of model parameterizations to take advantage of one representation over the other.
Streamflow is a direct product of runoff routing processes that are fed by, and therefore directly controlled by, runoff generation in terms of both magnitude and timing.Runoff generation itself is controlled by the interactions between climate and landscape properties and the latter two are very often closely interrelated to each other.Thus, to explain the differences in streamflow simulation between the two modeling frameworks, we first explore their differences in simulating runoff generation in different climate regimes.For this purpose, the subbasins/grids of the finest (0.125 • ) resolution in the CRB domain are grouped into different regimes by rainfall fraction (ratio of rainfall to the total precipitation) as snow dominated (areas with rainfall fraction ranging between 0.1 and 0.5), intermediate (areas with rainfall fraction ranging between 0.5 and 0.75), and rain dominated (areas with rainfall fraction ranging between 0.75 and 1.0) regimes (Fig. 4).The grids and subbasins are classified based on the same criteria, which result in spatial distributions of rainfall fraction largely consistent with the spatial distribution of elevation in the basin in both representations.This is not surprising since rainfall/snowfall partitioning of precipitation is dominated by near-surface air temperature, which is closely related to elevation variation (Tesfa et al., 2014).The total area of each regime in the CRB is listed in Table 1.Using different thresholds of 0.1, 0.4, 0.7 and 1.0 does not change the conclusions except the snow-dominated area is rather small as the two middle thresholds decrease.Hence subsequent analysis is based on the classification with thresholds of 0.1, 0.5, 0.75, and 1.0.In the subsequent sections, the climate regimes are used to investigate the differences of the two modeling frameworks in runoff generation.
Shown in Fig. 5 are the seasonal variation of runoff of the two modeling frameworks averaged over the whole CRB domain and the three different climate regimes.The phase shift between the two frameworks is consistent among different runoff components, i.e., surface runoff, subsurface runoff, and the total.This phase shift is more evident in the snowdominated areas than in the rain-dominated areas and more evident in the surface runoff than the subsurface runoff.In   the snow-dominated areas, SCLM-MOSART produces less subsurface runoff due to drier soil (Fig. 5).The latter is because the evaporation from bare soil and canopy simulated by SCLM-MOSART is overall slightly higher than that simulated by CLM-MOSART, which affect the soil moisture simulations.Compared to the phase shift in the simulated streamflow shown in Fig. 3, the phase shift in the simulated runoff is less significant.However, the transformation from runoff to streamflow is captured by the routing process, which is nonlinear in nature.Therefore one could infer that it is this transformation that has amplified the phase difference between the runoff simulated by the two modeling frameworks.It is then logical to ask whether and how this phase difference exists in the climatic forcings that are major drivers of runoff generation processes.
Figure 6 shows the seasonal variation of precipitation, temperature, and the partitioning of precipitation into rainfall and snowfall in the two modeling frameworks.From the plots for the whole CRB, one can see that there is no difference between the mean precipitation and temperature averaged over all subbasins and grids, which is expected because the remapping from grids to subbasins conserves the area-averaged forcings used as inputs to the models.However, the total precipitation is noticeably larger in CLM than SCLM in the snow-dominated areas, which is compensated by slightly smaller precipitation in CLM than SCLM in the rain-dominated areas, as the latter occupies a much larger fraction of the total area of the CRB.Similarly, the differences in rainfall and snowfall are more noticeable in the snow-dominated areas, with smaller differences in rainfall also noted in the rain-dominated areas.The models partition the total precipitation into rainfall/snowfall depending on air temperature.From the plots of air temperature for the different regimes, the difference between the two modeling frameworks is barely discernible in any climate regime.However, even very small differences in air temperature can lead to noticeable differences in the partitioning of rainfall/snowfall in areas with very high total precipitation.Hence larger total precipitation in CLM in the snow-dominated areas translates to larger snowfall in the cold season and larger rainfall in the warm season compared to SCLM, with opposite compensating effects in the rain-dominated areas.These differences reflect the dominant control of topography, hence air temperature, on the precipitation regimes; therefore, the model's spatial structure has an impact on precipitation regimes that translate to differences in runoff (Fig. 5) and streamflow (Fig. 3) due to the runoff generation and river routing processes.

Scalability of streamflow simulations
In this analysis, we explore how the two modeling frameworks simulate streamflow at different spatial resolutions.Figure 7 compares the streamflow simulated by SCLM-MOSART and CLM-MOSART at all spatial resolutions (0.125 • , 0.25 • , 0.5 • , and 1 • ) at the USGS stream gauges selected to represent major river basins in the CRB domain (Fig. 1).Results show that streamflow simulated by both frameworks has considerable variations across spatial resolutions.However, at these stream gauges, SCLM-MOSART is more consistent across spatial resolutions compared to that of CLM-MOSART in that SCLM-MOSART-simulated streamflow is consistently lower as spatial resolution coarsens, while the streamflow from CLM-MOSART tends to be less consistent across spatial resolutions and produces outliers at 1 • resolution.This is particularly true over smaller drainage areas such as manifested at the HCANY, BROWN, CJSTR, MILNE and ARROW stream gauges.This issue is thus investigated further at the stream gauges selected in Sect.3.2 in the subsequent sections.
Shown in Fig. 8 are scatterplots and statistics boxplots of the absolute differences in specific peak flow (ADP) calculated, as described in Sect.3.4, between each coarseresolution simulation and the corresponding reference simulation at 0.125 • resolution in each modeling framework at 3-hourly (Fig. 8a, d), daily (Fig. 8b, e) and monthly (Fig. 8c, f) timescales.In the figure, symbols represent spatial resolutions, while colors are used to identify the dominant climate regimes in the catchment area of the stream gauges.From the scatterplots, it is obvious that the differences between the coarse simulations and the reference simulation from SCLM are generally smaller than that from CLM, especially in snow-dominated areas (e.g., the blue symbols are more often below the 1 : 1 line than the red symbols), consistent with the finding of Tesfa et al. (2014) for runoff.Hence in both types of plots (scatterplots and boxplots), SCLM-MOSART tends to show some scalability advantages compared to CLM-MOSART at all temporal scales, which becomes more evident as one goes from 3-hourly to monthly temporal scales, particularly at the 0.5 • and 1 • resolutions.In general, these results suggest that improved scalability in runoff generation combined with the routing processes resulted in better scalability in peak flow simulation for SCLM-MOSART compared to CLM-MOSART at the stream gauges within the CRB domain.
We also calculated NSE values between each coarseresolution simulation and the reference simulation to compare scalability of the two approaches.Figure 9  for each coarse resolution (0.25 • , 0.5 • , and 1 • ) at the USGS stream gauges of CRB at 3-hourly (Fig. 9a), daily (Fig. 9b) and monthly (Fig. 9c) temporal scales.Similar to Fig. 8, spatial resolutions are identified by the symbols, while the colors represent the climate regimes dominating the catchment area of the stream gauges.Results show that SCLM-MOSART has superior scalability compared to CLM-MOSART at all the coarse spatial resolutions (0.25 • , 0.5 • , and 1 • ) and all temporal scales.The scalability advantage of SCLM-MOSART in streamflow simulation becomes more pronounced at the coarser spatial resolutions and monthly temporal scales.Since the routing process is less important at a monthly scale, these results suggest the importance of the scalability advantage of SCLM in runoff generation, particularly in snow-dominated and intermediate areas, at monthly scale.Note that the snow-dominated regions in this study domain largely overlap with the mountainous regions defined based on topographic steepness in Tesfa et al. (2014).This generally suggests that the scalability advantages of SCLM in runoff generation, discussed in Tesfa et al. (2014), are preserved even after coupling with the routing model (MOSART).As streamflow simulated at the stream gauges depends on the contributing areas, it is important to determine if contributing area differences between the two approaches may play a role in the scalability differences and how scalability differences may vary across spatial scales (i.e., upstream drainage areas of the stream gauges).We calculated another set of NSE values by normalizing the simulated streamflow by the corresponding contributing area and compared against the NSE values calculated using nonnormalized streamflow at the same stream gauges at daily (Fig. 10a, d) and monthly temporal scales (Fig. 10b, e).Results show similar contrast in scalability between the two modeling frameworks using both sets of NSEs, suggesting a minor role of contributing area in the scalability differences of the two modeling frameworks in streamflow simulation.Results also show more clearly that the scalability advantages of SCLM-MOSART are more significant at the 0.5 • and 1 • spatial resolutions.Comparison of the two sets of NSEs across spatial scales shows that the slight differences between the area normalized (Fig. 10c) and nonnormalized (Fig. 10f) streamflow occur mostly at stream gauges with smaller drainage areas (less than ∼ 10 5 km 2 ), suggesting that the role of contributing area on the scalability differences of the two modeling frameworks diminishes as the spatial scale increases.The results also show that both modeling frameworks have a threshold behavior with increasing spatial scale (drainage area), that is, in both sets of NSE values, the ability to reproduce the finestresolution (0.125 • ) simulations generally improve with increasing catchment size at all spatial resolutions (0.25 • , 0.5 • , and 1 • ).But, in both sets of NSE values, SCLM-MOSART results converge to the reference simulation faster with increasing drainage areas than that of CLM-MOSART at all coarse spatial resolutions.Comparisons at 3-hourly and daily temporal scales show the same pattern (figures not shown).
Following the results discussed so far, it is logical to ask whether the scalability advantages of SCLM-MOSART have any statistical significance.Table 2 shows the p value results from a nonparametric statistical significance t test (Bauer, 1972) on the NSEs calculated from the nonnormalized streamflow at each coarse (0.25 • , 0.5 • , and 1 • ) spatial resolution at 3-hourly, daily and monthly temporal scales.Using a confidence level of 95 %, the results show (1) significant differences in NSEs between SCLM-MOSART and CLM-MOSART at the 0.5 • and 1 • spatial resolutions at all (3-hourly, daily and monthly) temporal scales, (2) insignificant difference in NSEs at 0.25 • resolution at all temporal scales, and (3) the significance of the differences in NSEs between the two modeling frameworks increases from 3-hourly to monthly temporal scales.
To determine if the scalability difference in the CRB, which is dominated by topographic control on precipitation and runoff, may be generalized in other regions, similar analysis is performed for the MW study domain for streamflow at stream gauges selected using the procedure described in Sect.3.2.Analysis is only performed at the spatial resolutions that showed statistically significant scalability differences in the CRB domain.For this purpose, NSE values are calculated using the nonnormalized streamflow simulated at the coarse (0.5 • and 1 • ) resolutions and the corresponding reference (0.125 • resolution).Figure 11a shows minimal scalability differences between the two frameworks in MW.Since the MW domain is large and more heterogeneous in climate and topographic regimes than the CRB, we further compared the NSE values at the three major river basins: OHRB, MSRB, and UMRB of the Midwest (Fig. 11b-d  and the CRB domains, the results at the whole MW domain, MSRB, as well as UMRB are not surprising.The MW domain is dominated by flat topography and precipitation occurs mainly as rain, while the CRB domain is dominated by mountainous topography and precipitation occurs mainly as snow.Thus, these results are generally consistent with the findings in Tesfa et al. (2014), which showed the scalability advantages of the subbasin-based land surface modeling in runoff generation to be dominated by its superior scalability in mountainous and snow-dominated regions due to better consistency in representing mountainous topography and snow over complex topographic regions.However, the results over OHRB deserve further investigation.
To get a better understanding of the causes of the scalability differences in OHRB, we further explore how the topographic slope and rainfall fractions compare over the three major river basins of the MW domain (OHRB, MSRB and UMRB).Figure 12 shows the statistics of the topographic slope (Fig. 12a) and rainfall fractions (Fig. 12b) of the three river basins.The results reveal that both topographic slope and rainfall fractions of OHRB are quite different from those of MSRB and UMRB.It is interesting that, compared to MSRB and UMRB, OHRB is dominated by a much steeper topographic slope and dominantly higher rainfall fraction, suggesting that the scalability advantage of SCLM-MOSART in streamflow simulated in OHRB comes from a combination of improved scalability in the saturated fraction of surface runoff driven by rain over steep topography and the routing processes.As discussed in Tesfa et al. (2014), the subbasin-based approach (SCLM) is more consistent than the grid-based approach (CLM) in representing mountainous topography, which plays important role in the scalability of the rain driven saturated component of surface runoff.These results thus generally suggest that the scalability differences in streamflow simulation between the two modeling frameworks are generalizable to other regions.

Summary and conclusions
In this study, we have implemented a subbasin-based representation of CLM called SCLM, coupled with a physically based river routing model (MOSART).The relative merits of the subbasin-based modeling framework (SCLM-MOSART) in streamflow simulation are compared to the grid-based modeling framework (CLM-MOSART) over topographically and climatologically contrasting regions: the Columbia River Basin (CRB) and US Midwest region (MW).For this purpose, the two modeling frameworks are applied at four spatial resolutions (0.125 • , 0.25 • , 0.5 • , and 1 • ) in both the CRB and MW, and streamflow simulated by SCLM-MOSART and CLM-MOSART are compared with each other and with naturalized streamflow data.
We found that in the CRB where topography dominantly controls the precipitation regimes, small differences between the averaged atmospheric forcing for the two modeling frameworks could lead to larger differences in simulated runoff and streamflow at the finest (0.125 • ) resolution because of the nonlinear runoff generation and streamflow routing processes.Our results showed that simply by using a spatial structure that follows subbasin boundaries defined by topography without any change in model parameterizations, SCLM-MOSART exhibits improved scalability in simulating both peak and mean streamflow compared to CLM-MOSART.The scalability advantages of SCLM-MOSART are more apparent in snow-dominated and intermediate climate regimes, and in areas with steeper topography.This suggests that the scalability advantages of SCLM in runoff generation, discussed in Tesfa et al. (2014), are preserved even after coupling with the routing model (MOSART).Both modeling frameworks showed a threshold behavior with spatial scales (drainage area of the stream gauges); i.e., for drainage area larger than a threshold, spatial resolution becomes less important so the coarse-resolution simulations resemble the fine-resolution simulations, but SCLM-MOSART converges to the high-resolution simulations faster with increasing drainage area than CLM-MOSART.Lastly, we found that the understanding on scalability differences in streamflow simulation between the two modeling frameworks is generalizable to other regions.
The scalability results presented in this study suggest that the subbasin-based representation is more robust than the grid-based representation across spatial scales.This reduced sensitivity to model resolution for both peak and mean flow is an important advantage for reliable hydrologic predictions.Given that the scalability advantages have been identified for both runoff generation and streamflow simulations, it would be valuable to further examine how the scalability advantages partition between the two nonlinear processes and further contrast scalability differences in different topographic and climate regions.This would include (1) analyses of runoff generation over the catchment area of each stream gauge; and (2) a detailed investigation of river routing parameters such as drainage density, topographic slope and channel geometry, which are potential sources of differences between the streamflow simulated by the two modeling frameworks.Furthermore, given that the topography-relevant runoff generation parameter, f max , is derived from the HydroSHEDS (Lehner et al., 2008) 90 m DEM database, which is at a considerably finer resolution compared to the 1 km data provided with CLM4, it would be valuable to examine its relative merits on runoff generation/streamflow simulation in the two modeling frameworks across different spatial resolutions.
Since soil depth is a notoriously difficult parameter in hydrologic modeling (Tesfa et al., 2009), it is greatly simplified in CLM4, assuming a globally uniform soil-depth value of 3.8 m (Oleson et al., 2010).Even though it is a bold assumption, it is typical in the field of land surface modeling due to lack of global soil-depth data.Effective incorporation of spatially distributed soil depth in CLM requires two major steps: (1) deriving a large-scale soil-depth map; and (2) significant modifications of CLM source code to accommodate new soil-layer delineation and associated hydrology, thermodynamics and biogeochemical processes.Both of these steps demand substantial efforts and cannot be achieved in a short period.We therefore leave the soil-depth issue for future study.
The analyses presented in this study, covering a wide range of model resolutions, are particularly useful as the land surface modeling community is exploring the feasibility and advantages of ultra-high resolution (e.g., 1 km resolution globally) (Wood et al., 2011).Conceptually, the differences between the two approaches may be larger as lateral water redistribution becomes increasingly important in determining soil moisture states at smaller spatial scales.The results in this study also suggest that without proper calibration, SCLM-MOSART may not necessarily perform better than CLM-MOSART.A fair comparison between the two modeling frameworks in reproducing the observed streamflow requires a separate proper parameter calibration for each; thus, future research to include parameter calibration of SCLM and CLM on smaller basins with good forcing and evaluation data would be useful to test this point.Given that CLM is the land component of an Earth system model and can interact with the atmosphere component and ocean component, it would be valuable to examine how the subbasinbased representation of terrestrial processes may affect the global cycling of water and energy between land, atmosphere and ocean.Such a scientific pursuit is now supported by recent progresses made in software engineering.With the latest public versions of the Community Climate System model and the Community Earth System model (i.e., CCSM4 and CESM1) there is a suite of new coupling capabilities in the CPL7 coupler that allow more flexibility and extensibility for very-high-resolution modeling and coupling Earth system components configured on unstructured grids (Craig et al., 2012).

1Figure 1 .
Figure 1.(a) The Columbia River basin (CRB) and US Midwest (MW), elevation, stream network, USGS stations with drainage area greater than 15 000 km 2 and Midwest regions; (b) subbasin delineation overlaid with regular grids at 0.125 • resolution (highlighted with red dashed line are example grids containing multiple river channels). 1

Figure 2 .
Figure 2. Comparison of precipitation and air temperature representation in the two modeling frameworks.

Figure 4 .
Figure 4. Climate regions based on model simulated rainfall/snowfall partitioning.

Figure 5 .
Figure 5. Seasonality of runoff and soil water over the climate regions.Note that soil water is included here because it is closely related to the runoff generation.

1Figure 6 .
Figure 6.Seasonality of forcing over the climate regions.

Figure 8 .
Figure 8. Specific peak flow comparison at 3-hourly, daily and monthly temporal scales at the USGS stations with contributing area larger than 15 000 km 2 : (a), (b) and (c) comparison over the climate regions and (d), (e) and (f) comparing their statistics.

1Figure 9 .
Figure 9. Scalability comparison using a NSE of streamflow at 3-hourly (a), daily (b) and monthly (c) calculated between each coarse scale and the corresponding fine scale in each modeling framework at the USGS stations with contributing area larger than 15 000 km 2 .

Figure 10 .
Figure 10.Statistics of (a and b) area-normalized and (d and e) nonnormalized NSEs of streamflow at different temporal scales and (c and f) NSEs of monthly streamflow across the spatial scale over USGS stations with drainage area greater than 15 000 km 2 in the CRB.

Figure 11 .
Figure 11.NSE values from monthly streamflow compared at the USGS stations with contributing area larger than 15 000 km 2 located in the whole Midwest (a), and Missouri (b), Upper Mississippi (c) and Ohio (d) regions.

Figure 12 .
Figure 12.Topographic slope and rainfall fraction over the Missouri, Upper Mississippi and Ohio regions.

Table 1 .
The portion of areas within three climate regimes (km 2 ).

Table 2 .
Nonparametric t-test p values on NSE.