Exploring new topography-based subgrid spatial structures for improving land surface modeling

Topography plays an important role in land surface processes through its influence on atmospheric forcing, soil and vegetation properties, and river network topology and drainage area. Land surface models with a spatial structure that captures spatial heterogeneity, which is directly affected by topography, may improve the representation of land surface processes. Previous studies found that land surface modeling, using subbasins instead of structured grids as computational units, improves the scalability of simulated runoff and streamflow processes. In this study, new land surface spatial structures are explored by further dividing subbasins into subgrid structures based on topographic properties, including surface elevation, slope and aspect. Two methods (local and global) of watershed discretization are applied to derive two types of subgrid structures (geo-located and non-geo-located) over the topographically diverse Columbia River basin in the northwestern United States. In the global method, a fixed elevation classification scheme is used to discretize subbasins. The local method utilizes concepts of hypsometric analysis to discretize each subbasin, using different elevation ranges that also naturally account for slope variations. The relative merits of the two methods and subgrid structures are investigated for their ability to capture topographic heterogeneity and the implications of this on representations of atmospheric forcing and land cover spatial patterns. Results showed that the local method reduces the standard deviation (SD) of subgrid surface elevation in the study domain by 17 to 19 % compared to the global method, highlighting the relative advantages of the local method for capturing subgrid topographic variations. The comparison between the two types of subgrid structures showed that the non-geo-located subgrid structures are more consistent across different area threshold values than the geo-located subgrid structures. Overall the local method and non-geolocated subgrid structures effectively and robustly capture topographic, climatic and vegetation variability, which is important for land surface modeling.


Introduction
Topography plays an important role in land surface processes through its influence on atmospheric forcing, soil and vegetation properties, and river network topology and drainage area.Consequently, accurate climate and land surface simulations in mountainous regions cannot be achieved without considering the effects of topographic heterogeneity (Leung andGhan, 1998, 1995;Ghan et al., 2006).Mountain water resources are particularly sensitive to global warming (e.g., Leung and Ghan, 1999;Ghan and Shippert, 2006;Mote et al., 2007;Kapnick and Hall, 2012).The amplified warming at high elevation due to the lapse-rate effect and snow albedo feedback has a large impact on snowpack accumulation and melt, with consequential effects on runoff and water supply (Leung et al., 2004;McCabe and Clark, 2005;Rasmussen et al., 2011).
Topography has major influence on the spatial pattern of atmospheric forcing, including surface temperature, precipitation, and incoming and reflected solar radiation.Regions characterized by heterogeneous topography generally exhibit diverse hydroclimatic conditions.For example, stable moisture-rich air lifted by the mountains can produce orographic precipitation that dominates the spatial distribution of cold season precipitation in the western United States (Leung et al., 2003).In mid-and high-latitude regions, topography also influences the partitioning of precipitation into snow Published by Copernicus Publications on behalf of the European Geosciences Union.and rainfall.In addition, incoming and reflected solar radiation is highly dependent on the orientation of landscapes, which can also have significant impacts on surface hydrology through the effects of radiation on cloud, precipitation, and snow processes (Lee et al., 2015).
Topography is an important factor in soil formation, exerting dominant control on the spatial patterns of soil properties over watersheds, e.g., soil depth (Tesfa et al., 2009, and references there in).Soils are generally deeper and finer in texture over valleys compared to the shallower and coarser texture over ridges of watersheds.Through its influence on direct and diffuse solar radiation and consequent effects on soil moisture and evapotranspiration, topography affects the spatial pattern of vegetation on a landscape.Different vegetation types grow on different parts of a landscape depending on their water demand and resistance to water stress.In semiarid regions, vegetation types that have high water demand or less resistant to moisture stress grow near streams, while vegetation types resistant to moisture stress can grow further from streams (Tesfa et al., 2011).Topography also determines the topology of river network and drainage areas, which in turn control surface and subsurface flows (Beven, 1997;Chen and Kumar, 2001).Overall, catchment ecohydrology is strongly affected by the topography-mediated interactions among vegetation, soil and river networks (Thompson et al., 2011).Improving representations of land-atmosphere and surface-subsurface interactions affected by fine-scale topography and vegetation have been identified as a grand challenge, motivating the need for hyper-resolution land surface modeling (Wood et al., 2011).While hyper-resolution modeling approaches are being tested at regional (Singh et al., 2015) and continental scales (Maxwell et al., 2015), improving the spatial structures of land surface models to capture the effects of topographic heterogeneity could be crucial to advancing modeling of land-atmosphere interactions in Earth system models.Tesfa et al. (2014a, b) demonstrated improved scalability of simulated runoff and streamflow processes when subbasins instead of structured grids are used as computational units in the Community Land Model.The improvements of the subbasin-based land surface modeling in scalability come from its important conceptual advantages in capturing atmospheric forcing and runoff generation processes, both strongly affected by a topography that defines the boundaries of the subbasins.
Discretization of the subbasins to capture spatial heterogeneity affected by topography may further improve the representation of land surface processes.Ke et al. (2013) evaluated several classification methods to account for subgrid variability of surface elevation and vegetation cover for land surface models with structured grids.To the best of our knowledge, development of subgrid structures for subbasinbased land surface modeling has not been attempted.The purpose of this paper is to explore subgrid structures that capture topographic heterogeneity and its influences on land surface processes for land surface modeling.Such subgrid spatial structures may provide a more realistic spatial distribution of surface properties and their influence on climatic variability, with more reasonable computational requirements compared to discretizing the domain into fineresolution grid-based representations that are reported in the literature (e.g., Singh et al., 2015).
Motivated by the significant influences of topographic heterogeneity on land surface processes, we explore new topography-based spatial structures by further dividing subbasins into subgrid structures or subgrid units (also hereafter denoted as SUs) to take advantage of the emergent patterns and scaling properties of atmospheric, hydrologic and vegetation processes in land surface models.For this purpose, two methods (global and local) of subbasin discretization are applied to derive two types of SUs (geo-located and non-geo-located) over the topographically diverse regions of the northwestern United States.In the global method, the subbasins are discretized into multiple SUs, following the surface elevation classification scheme employed in Leung and Ghan (1998Ghan ( , 1995)), combined with classifications of topographic slope and aspect.The local method utilizes concepts of hypsometric analysis (Willgoose and Hancock, 1998;Sinha-Roy, 2002) combined with the classification of topographic aspect to discretize each subbasin into multiple SUs.We evaluate the two discretization methods and spatial structures for their ability to capture topographic heterogeneity and the implications of this on representations of atmospheric forcing and land cover spatial patterns.For Earth system modeling, the atmospheric model will adopt the subgrid topographic precipitation scheme of Leung and Ghan (1998) and Ghan et al. (2006), which provides subgrid atmospheric forcing for coupling with the topographic SUs of land surface models, representing topographic effects on atmospheric and land surface processes and land-atmosphere interactions.
The remainder of the paper is organized as follows: Sect. 2 describes the study area.Development of new subgrid structures is discussed in Sect.3. The strategy used to evaluate the methods of subbasin discretization and the subgrid structures are discussed in Sect. 4. Section 5 presents the results and discussion, and Sect.6 closes with conclusions and recommendations.

Study area
To investigate the importance of various watershed discretization methods, the Columbia River basin (located in the US Pacific Northwest) is used as a case study.Figure 1 shows the topographic patterns, subbasins with an average size equivalent of 1/8th degree grid, elevation ranges of the subbasins and two subbasins representing extreme topographic properties.The Columbia River basin encompasses both mountainous and low-lying regions.Climatically, the mountainous regions are characterized by low temperature and higher precipitation dominated by snowfall, whereas the

Development of new SUs for land surface modeling
Two methods of subbasin discretization are implemented to develop land surface subgrid structures that capture the spatial heterogeneity affected by topography.Both methods are applied to derive two types of SUs: geo-located and non-geolocated.The following subsections describe the input data, the two discretization methods and the two types of SUs.

Input data
To derive the subgrid units, the study domain is first delineated into subbasins.We utilize the subbasins equivalent to 1/8th degree grids delineated in Tesfa et al. (2014a, b) using the ArcSWAT (Soil and Water Assessment Tool; Neitsch et al., 2005) with the 90 m digital elevation model (DEM) and the 15 arcsec river network from the Hydrological data and maps based on Shuttle Elevation Derivatives (HydroSHEDS; Lehner et al., 2008).Although DEMs at resolutions of 30 m or finer are available from the United States Geological Survey (USGS), we use the DEM from a global database (i.e., HydroSHEDS) because the main goal of this study is to develop subgrid structures for global land surface models and Earth system models.Along with the delineation of the subbasins, topographic attributes, such as slope and aspect, are derived for the study domain to be used as inputs for the subbasin discretization methods.

Global method
In the global method, the study domain is first discretized into 12 elevation classes based on surface elevation extracted from the 90 m HydroSHEDS' DEM, following the surface elevation classification scheme employed in Leung andGhan (1998, 1995), which uses class intervals of 100 m for surface elevation below 500 m and gradually increases to intervals of 500 and 1000 m for high surface elevations, resulting in 12 elevation classes (see Fig. 2).This method is global because the same elevation classification scheme is used to discretize all subbasins regardless of the elevation spanned by individual subbasins, which can vary substantially.Since topography influences atmospheric and land surface processes through surface elevation, slope and aspect, the global method combines topographic slope and aspect with the elevation classes.For this purpose, the study domain is also partitioned into two classes of topographic slope where slope values less than or equal to 20 Figure 2. The study area classified into elevation bands used in the global method, following the approach described in Leung andGhan (1995, 1998).
sets of SUs derived based on elevation, slope and aspect separately.The SUs derived from elevation, slope and aspect are then intersected to generate SUs based on the combination of topographic elevation, slope and aspect, resulting in a large number of SUs for each subbasin.Since many of the SUs are extremely small in size, and our goal is to capture topographic heterogeneity with only a reasonable number of SUs for computational efficiency, an area threshold value is used to merge SUs with an area smaller than the threshold to their neighboring SUs with a size larger than or equal to the threshold value.The area threshold is defined based on the percentage of the area of each subbasin.Subgrid units that are smaller than the area threshold are merged to their neighboring larger subgrid units.To enable discretization of each subbasin into a reasonable number of subgrid units, a value of normalized area (SU na ) is calculated for each subgrid unit following Eq.( 1) and compared to the area threshold value.
where SU a and Sub a are areas of subgrid unit and subbasin, respectively.The global method has been implemented in Python and utilizes ArcGIS functionalities.In this effort, the global method is applied to derive both geo-located and nongeo-located SUs.

Local method
In the local method of subbasin discretization, the subbasins for the study domain are first classified into five groups based on values of elevation range using the natural breaks (Jenks) classification method in ArcGIS.As an example, to derive the hypsometric curves, the two contrasting subbasins shown in Fig. 1 are discretized into 100 elevation contours using elevation data extracted from the 90 m resolution DEM from HydroSHEDS.The relative elevation (RH) and relative area (RA) are calculated for each contour, where relative elevation (RH) is defined as the ratio of the height of the given contour (h) from the base plane of the subbasin to the maximum height of the subbasin (H ), whereas relative area (RA) refers to the ratio of the area above a particular contour (a) to the total area of the subbasin (A).
The hypsometric curves are derived by plotting the relative area (RA) along the abscissa and the relative elevation (RH) on the ordinate axes.In geomorphology, a hypsometric curve is used to characterize the distribution of elevation within  1S and 2S in the Supplement).This method is local as the elevation ranges used to discretize the subbasins vary depending on the topographic variations within each subbasin.Similar to the global method, classes of topographic aspect are extracted for each subbasin and intersected with the corresponding elevation classes to discretize the subbasin into multiple SUs, with some extremely small in size.Since discretizing the subbasins using a hypsometric curve is expected to capture slope variation implicitly, a topographic slope is not used in the local method.With the main goal being to capture topographic heterogeneity with only a reasonable number of SUs for computational efficiency, similar to the global method, area threshold is utilized to merge the SUs with an area smaller than the threshold to their neighboring SUs with a size larger than or equal to the threshold to develop the final SUs.This method has also been implemented in Python and utilizes ArcGIS functionalities.The local method is also applied to derive both geo-located and non-geo-located SUs.In this method, the actual number of SUs of each subbasin depends on the variability of surface elevation and topographic aspect within the subbasin boundary.

Types of subgrid units
Two types of SUs are derived using both the global and local methods: geo-located and non-geo-located.The geo-located SUs are derived by discretizing the subbasins into spatially contiguous structures.They are characterized with explicit geographical location and a single boundary.In this case, SUs with the same topographic characteristics at different locations of the subbasin are treated as separate units.The non-geo-located SUs are developed by discretizing the subbasins into spatially non-contiguous structures.In this case, SUs with the same topographic properties at different locations of the subbasin are treated as a single unit resulting generally in reduced number of SUs compared to the geo-located SUs.
4 Evaluation strategy

Analysis using SUs based only on elevation
Because topographic slope is not explicitly used in the local method, it is logical to ask whether discretizing subbasins using the hypsometric curve is capable of implicitly capturing the variability of topographic slope within the subbasins.To investigate this, geo-located SUs are derived using both global and local methods based on elevation classification only.The number of SUs for each subbasin from both methods is compared to the average values of topographic slope of the subbasins in the study area to determine how topographic slope influences the number of SUs needed to capture subgrid topographic variability in each method.In addition, the spatial pattern of the number of SUs for each subbasin derived using each method is compared to the spatial pattern of topographic slope and elevation range within the subbasins for the study region.An effective subgrid method would allow more SUs in subbasins with complex terrain to capture the subgrid topographic variability and use fewer SUs www.geosci-model-dev.net/10/873/2017/Geosci.Model Dev., 10, 873-888, 2017 in subbasins with small variations of topography.Finally, the relative capability of the two methods in capturing topographic heterogeneity and their sensitivity to the values of area threshold are evaluated, respectively, based on the standard deviation (SD) of the 90 m resolution elevation within the SUs and the variation of statistical metrics (the total number of SUs, mean SU size and SD in SU size) calculated for the study domain across different values of area threshold (1, 2, 3, 4 and 5 %).Methods that are less sensitive to the values of area threshold can provide more robust SUs for representing subgrid topographic heterogeneity.

Analysis using SUs based on elevation, slope and aspect
The two types of SUs are expected to differ in their ability to capture topographic heterogeneity, i.e., the number of SUs, which has important implications to the overall computational burden, and their sensitivity to area threshold values, which is important for defining robust SUs for land surface modeling.Thus, to evaluate the two types of SUs with respect to their applications in land surface modeling, geo-located and non-geo-located SUs for the study area are derived based on elevation, slope and aspect using both global and local discretization methods at different values of area threshold (1, 2, 3, 4 and 5 %).The geo-located and non-geo-located SUs of each method are then compared for their sensitivity across values of area threshold using statistical metrics (total number of SUs, average size of SUs and SD in SU size) calculated over the study domain at different values of area threshold.
The global and local methods are further investigated for their capability in capturing topographic heterogeneity and consistency across different values of area threshold when using the non-geo-located SUs.The relative capability of the non-geo-located SUs from both methods in capturing topographic heterogeneity is evaluated based on the values of SD in surface elevation calculated at each SU across different values of area threshold.In addition, sensitivity of the two methods (global and local) to the values of area threshold when used to derive non-geo-located SUs is evaluated using statistical metrics calculated over the study domain such as total number of SUs, average size of SUs and SD in SU size.

Implications to representation of land surface processes
Since the main goal of this study is to derive land surface structures capable of improving representation of land surface processes in land surface modeling, it is logical to ask how the new structures impact the representation of land surface parameters.For this purpose, the two methods are first evaluated for their relative capability of capturing climatic and land cover variability over the study area using the nongeo-located SUs derived at different values of area threshold.
The capability of capturing climatic variation is investigated by comparing values of SD in precipitation and surface temperature within the SUs derived using the two methods.In this case, the precipitation and surface temperature datasets for the study area are extracted from the 30-year normal annual precipitation and mean annual surface temperature obtained from the PRISM climate datasets (800 m spatial resolution) (http://www.prism.oregonstate.edu/).Similarly, using the normalized difference vegetation index (NDVI) data as a proxy for land cover, the relative capability of the two methods in capturing land cover pattern over the study domain is investigated by comparing values of SD in NDVI calculated within the SUs from the two methods.For this purpose, the NDVI datasets for the study area are obtained from the enhanced Moderate Resolution Imaging Spectroradiometer (eMODIS) data (250 m spatial resolution) portal (http://earthexplorer.usgs.gov/)at the Earth Observation and Modeling Facility (EOMF).Furthermore, to evaluate the relative advantages of the non-geo-located SUs derived using the local method in capturing climatic variability in the study domain as compared to those of the subbasin-based representation, the spatial distributions of precipitation and temperature mapped to the subbasins and the non-geo-located SUs from the local method are compared to the spatial distributions of precipitation and temperature from the original highresolution PRISM grid-based representation.The subbasinbased representation used in this comparison comes from our previous studies (Tesfa et al., 2014a, b), which evaluated the benefits of land surface modeling using the subbasin-based approach against the standard regular grid-based land surface modeling approach, where significant advantages in simulations of hydrologic and streamflow were demonstrated by the subbasin-based approach.
5 Results and discussion

Global versus local methods using elevation-based SUs
Since the main differences between global and local methods are in the way subbasins are discretized into elevation classes and whether topographic slope is included explicitly, the relative capability of the two methods in capturing topographic heterogeneity is investigated using elevation-based SUs. Figure 4 compares how well the global and local subbasin discretization methods capture the topographic slope using elevation-based geo-located SUs derived based on elevation at 1 % area threshold.For this purpose, the numbers of SUs per subbasin resulted from both methods are compared to the average topographic slope calculated over the subbasins.The results show the number of SUs per subbasin from the local method is directly related to the average subbasin slope (r 2 = 0.47); therefore, the steep subbasins are generally discretized into more SUs than the flat subbasins.On the other hand, the number of SUs per subbasin from the global method is not related (r 2 = 0.07) to the average topographic slope of the subbasins.From this comparison, it is clear that discretizing subbasins following the local method (see algorithms in Tables 1S and 2S in the Supplement) using the hypsometric curve characterization within each subbasin is able to capture topographic slope implicitly, making the local method superior to the global method.To investigate how the performance of the global method can be improved when the effect of slope is considered, the number of geo-located SUs per subbasin derived using a combination of topographic elevation and slope is compared to the average topographic slope calculated over the subbasins (Fig. 1s in the Supplement).The results show significant improvement of the capability of the global method in capturing topographic slope (r 2 = 0.44) as compared to the global method using topographic elevation only (r 2 = 0.07); however, the performance is still not as good as that of the local method (r 2 = 0.47).
The Columbia River basin encompasses diverse topography ranging from flat to steep mountainous areas making it an ideal study area for evaluating the relative capability of the two subbasin discretization methods in capturing the spatial pattern of topographic properties.The spatial pattern of the numbers of elevation-based geo-located SUs per subbasin derived using both methods with a 3 % area threshold are compared to the spatial pattern of the average topographic slope and elevation ranges of the subbasins classified based on the natural breaks (Jenks) classification method in Ar-cGIS (Fig. 5).The results suggest that the spatial pattern of the number of SUs per subbasin for the SUs from the local method follows the topographic pattern in the study area better than those of the global method, confirming further the advantages of discretizing the subbasins using the local method.To quantify the correspondence between the pattern of the surface topography and the pattern of the number of SUs, correlation coefficients are calculated between values of surface elevation range within the subbasins and the number of SUs per subbasin.The correlation coefficients are 0.66 and 0.47 for the local and global methods, respectively.Hence, the number of SUs per subbasin from the local method mimics the topographic pattern better than the global method, as more SUs per subbasin are defined over mountainous areas and fewer SUs per subbasin are defined over relatively flat areas of the basin.This enables the model to capture the topographic heterogeneity with a minimum number of SUs over the study domain, which is essential for computational efficiency in land surface modeling.
Figure 6 shows the relative capability of the two methods in capturing subgrid topographic heterogeneity across different values of area threshold using elevation-based geolocated SUs.For this purpose, values of SD in elevation within the geo-located SUs derived using different values of area threshold (1, 2, 3, 4 and 5 %) from both methods are compared.The results again clearly show that the SUs from the local method are able to capture topographic heterogeneity, which is reflected in the smaller SD of topography within each SU across different values of area threshold, better than those of the global method.In addition, the results also show that the local method can capture topographic heterogeneity more consistently across different values of area threshold than the global method, suggesting that the SUs derived using the local method are more robust.
Using the same SUs, the two methods are further investigated for their sensitivity to values of area threshold using the variability of statistical metrics (total number of SUs, mean SU size and SD in SU size) calculated over the whole study domain for different values of area threshold.The results in Fig. 7 show that SUs derived using local method remain more consistent across different values of area threshold than those of the global method, making the local method more robust than the global method for land surface modeling.

Geo-located versus non-geo-located SUs
To evaluate the robustness of the two types of SUs (geolocated and non-geo-located) for land surface modeling, we compare their sensitivity to values of area threshold.For this purpose, geo-located and non-geo-located SUs are derived based on elevation, slope and aspect using both methods at different values of area threshold.The geo-located SUs from each method are then compared to the corresponding nongeo-located SUs derived using the same method based on the statistical metrics calculated over the whole study area.Shown in Fig. 8 are comparisons of the variability of the total number of SUs, average SU size and SD in SU size calculated for the geo-located SUs against those of the non-geolocated SUs for both the global (Fig. 8a, b and c) and the local (Fig. 8d, e and f) methods.In both methods, the results gener-ally suggest that the non-geo-located SUs are more consistent across different values of area threshold than the corresponding geo-located SUs.Thus, in subsequent sections, the two methods of subbasin discretization are evaluated using the non-geo-located SUs only.

Global versus local methods using non-geo-located SUs
Following the evaluation of the two methods using elevationbased geo-located SUs, it is important to investigate whether the advantages of the local method over the global method shown in previous results still apply when the two methods are used to derive non-geo-located SUs based on the combination of multiple topographic properties.Shown in Fig. 9 are comparisons of sensitivity of the global and local meth- ods to values of area threshold when the two methods are applied to derive non-geo-located SUs using the variability of the statistical metrics (total number of SUs, average SU size and SD in SU sizes) calculated over the whole study domain at different values of area threshold.Note that the results in Fig. 9 are intended to evaluate how sensitive the two methods are to the values of area threshold and, unlike the comparison in Fig. 7, the SUs in this comparison are non-geo-located, derived based on a combined classification of elevation and topographic slope and aspect in the global method and elevation and topographic aspect in the local method.Similar to the comparisons in Fig. 7, the results suggest that the SUs from the local method are less sensitive to the values of area threshold, yielding more consistent values of the total number of SUs, average SU size and SD in SU sizes over the study domain than those of the global method.Shown in Fig. 10 are values of SD in elevation within the non-geo-located SUs derived using the global and local methods at different values of area threshold, comparing the capability of the two methods in capturing topographic heterogeneity when used for non-geo-located SUs.Similar to the results shown in Fig. 6, there is a clear difference in the capability of the two methods in capturing topographic heterogeneity across different values of area threshold.The nongeo-located SUs from the local method are able to capture topographic heterogeneity much better than those of the global method across different values of area threshold.The results in Table 1 show that the local method reduces the SD of sub-  and  3).The local method has been designed to minimize computational demand by discretizing mountainous areas into more subgrid units and flat areas into fewer subgrid units; therefore, its advantages are expected to be more pronounced when it is applied over topographically heterogeneous regions as opposed to regions characterized by homogenous topography.Furthermore, the results from Fig. 9a show that the area threshold value resulting in the same number of nongeo-located subgrid units from the two methods lies between 1 and 2 %.Hence, the results from Fig. 10 implicitly suggest that the local method still performs better even when we compare the two methods that produce the same number of SUs.The results also suggest that the capability of the non-geo-located SUs from the local method in capturing topographic heterogeneity remains more consistent at different values of area threshold than those of the global method, confirming the superior advantages of the local method.
From the results shown so far, relative to the global method, the SUs from the local method are superior in capturing topographic heterogeneity yielding more SUs per subbasin over mountainous areas and fewer SUs per subbasin over flat areas, which is essential for more realistic representations of the spatial distributions of precipitation and snow cover in mountainous areas and computational efficiency in land surface modeling.Furthermore, the SUs from the local method are more consistent across different values of area threshold than those of the global method.Subsequently, it is important to examine whether similar advantages exist for the local method in capturing climatic and land cover variability as compared to the global method.The following section focuses on the implications of the non-geo-located SUs in the representations of climatological and land cover variability in the study area.

Implications to the representation of land surface processes
Topography can influence land surface processes through its impact on atmospheric forcing and vegetation variability.Consequently, it is essential to examine the implications of the new SUs on representations of climatic and vegetation variability.This is particularly important as our goal is to couple land surface and atmosphere models both with topographic subgrid units to provide the largest improvement for capturing subgrid variability of land surface processes.As an example, the subgrid orographic precipitation scheme of Leung and Ghan (1995Ghan ( , 1998) ) has been shown to improve simulations of surface temperature, precipitation and snowpack in mountainous areas by representing the impact of subgrid topography on airflow and cloud processes.Shown in Fig.  are values of SD of the 30-year normal annual precipitation (Fig. 11a) and mean annual surface temperature (Fig. 11b) obtained from the high-resolution PRISM dataset calculated within the non-geo-located SUs derived using the global and local methods at different values of area threshold.This comparison is intended to evaluate the capability of the two methods to reduce climatic variability within the SU boundaries.
The results show generally lower values of SD in both precipitation and temperature for the SUs derived using the local method than those of global method across all values of area threshold.Consistent with the comparison to the capability of capturing topographic heterogeneity shown in Figs. 6 and 10, these differences reflect the dominant control of topography and the impact of spatial structure on precipitation and surface temperature, suggesting improved capability of capturing climatic variability for the local method.Furthermore, shown in Fig. 12 are values of SD of NDVI calculated at the non-geo-located SUs from both global and local methods at different values of area threshold.In this comparison, the NDVI is used as a proxy for land cover information during spring (Fig. 12a) and summer (Fig. 12b) extracted from the eMODIS dataset, showing the relative capability of the two methods in capturing land cover variability in the study area.The results generally show lower values of SD for the SUs derived using the local method than those of the global method across all values of area threshold, suggesting that the SUs from the local method have a better capability of capturing land cover variation in the study domain, which is essential to representation of land cover in land surface modeling.In all the results shown so far, the SUs from the local method have demonstrated clear advantages in capturing topographic heterogeneity and climatic and land cover variation compared to those of the global method over the study domain.Therefore, we further examined how representation of climatological forcing improves when using the nongeo-located SUs derived using the local method at 3 % area threshold value compared to the subbasin-based representation.This comparison is intended to evaluate the improvement in representing the spatial pattern of precipitation and temperature from the high-resolution PRISM datasets when using the non-geo-located SUs from the local method as compared to those of the subbasins without a subgrid clas- sification.Figure 13 compares the spatial pattern of the 30year mean precipitation represented based on subbasins at roughly 1/8th degree resolution (Fig. 13a), non-geo-located SUs within the subbasins (Fig. 13b) and the original grid representation from the PRISM dataset at 800 m resolution (Fig. 13c).Note that the Canadian part of the study domain is missing from the map because the PRISM data are only available for the United States.The results show that the SUbased representation yields a similar spatial pattern of precipitation to that of the original PRISM grids with no visually discernible difference.The spatial pattern of precipitation for the subbasin-based representation has noticeable differences from those of the SUs and original PRISM grid representa-tions, especially over the mountainous areas.With the ability to better capture the spatial heterogeneity of precipitation, land surface models that use the SU-based representation are expected to produce more realistic distribution of snow cover over the mountains compared to the subbasin-based representation, and it is considerably more efficient computationally compared to modeling land surface processes using hyper-resolution such as the PRISM grids.Further comparison of the three representations of precipitation using statistical metrics (average and SD values) reveals that both the mean and SD values from the SU-based representation are much closer to those of the original PRISM grids as compared to the subbasin-based representation (Table 2).
Similar comparisons are also shown in Fig. 14 for surface temperature.Similar to the results for precipitation, there is no visually noticeable difference in the spatial pattern of temperature between the SU-based and original PRISM gridbased representations, whereas the subbasin-based representation misses important variability indicated in the PRISM data.The advantage of the SU-based representation in capturing the spatial pattern of temperature is more pronounced over the mountainous areas.Comparison using statistical metrics of temperature (Table 2) confirms the advantages of the new SUs representation.Further comparison of statistical metrics of NDVI for spring and summer are shown in Table 3, comparing NDVI representations using the nongeo-located SUs from the local method generated at 3 % area threshold and the subbasins against those of the original high-resolution NDVI grid representation.The results show some improvement for the SU-based representation in the SD values, whereas the mean values of the SUs are generally higher than those of the subbasin and original grid representations.The latter is caused by the fact that mountainous areas are generally characterized by higher values of NDVI and they are discretized into more numbers of SUs compared to flat areas that generally have lower values of NDVI, so the mean value calculated over the study domain with the SUs is generally higher.Similar behavior is also observed in Table 2 where the average value of precipitation from the SUs calculated over the study domain is higher than that of the PRISM grid representation.
The advantages demonstrated for the SUs derived using the local method in representing topographic features are expected to be significant for land surface modeling in mountainous areas such as the Columbia River basin, where topography has dominant control on precipitation and temperature characteristics that translate to differences in runoff and streamflow characteristics (Tesfa et al., 2014a, b).
Also shown in the Supplement are similar comparisons of the non-geo-located SU-based representation from the global method against the subbasin-based representation of meteorological forcing and NDVI (Tables 3s and 4s and Figs.2s  and 3s).The results show general improvements from the use of non-geo-located SUs as compared to those of the subbasin-based representation.

Summary and conclusions
Topography exerts a major control on land surface processes through its influence on atmospheric forcing, soil and vegetation properties, network topology and drainage area.Thus, the spatial structure of land surface models that captures spatial heterogeneity affected by topography may improve modeling of terrestrial water cycle and land-atmosphere interac-tions.In both land surface and atmospheric modeling, such spatial structures are very much needed for accurate simulations of land surface processes in Earth system models.
While there are similar efforts to improve representation of the impacts of topography in atmospheric models, this study focuses exclusively on the development of new land surface spatial structures to improve representation of land surface processes in land surface models by further discretizing sub- basins into subgrid units (SUs) based on their topographic attributes (e.g., elevation, topographic slope and aspect).Two methods of watershed discretization (local and global) have been developed and applied over the Columbia River basin in northwestern United States to derive two types of topography-based subgrid structures (geo-located and nongeo-located).In addition, the two methods have been evaluated for their consistency, capability of capturing topographic heterogeneity and climatic and land cover variability of the study domain using both types of subgrid structures.
In the global method, the study domain is initially discretized into 12 elevation classes following the surface elevation classification scheme employed in Leung andGhan (1998, 1995).Then, following the subbasin boundary, the elevation classes are intersected with classes of topographic slope and aspect, discretizing each subbasin into multiple subgrid units.The local method utilizes concepts of hypsometric analysis to first discretize each subbasin into elevation classes using algorithms developed in this study, which are then merged with classes of topographic aspect to divide the subbasin into multiple subgrid units.In both methods, values of area threshold are used to merge small subgrid units into the neighboring large subgrid units, yielding a reasonable number of subgrid units per subbasin.Both methods are applied to derive two types of subgrid structures: geolocated (spatially contiguous) and non-geo-located (spatially non-contiguous).Furthermore, using both types of SUs, the two methods of subbasin discretization are investigated for their capability of capturing topographic heterogeneity, their implications on representations of climatic and vegetation variability in the study area, as well as their sensitivity to the area threshold values.
Using elevation-based geo-located subgrid units, comparison of the two methods showed that the local method is able to capture the topographic variability better than the global method.Taking advantage of hypsometric analysis, the local method can capture slope variability implicitly, and therefore it generally requires fewer SUs to represent subgrid topographic variability.The local method more effectively captures the topographic pattern across the region than the global method by discretizing steep subbasins into more subgrid units and flat subbasins into fewer subgrid units.Using the local method, the SD of surface elevation within the subgrid units is noticeably smaller and less sensitive to the values of area threshold than the global method.Hence, the local method is clearly more effective and robust for representing subgrid elevation variability for land surface modeling.
Comparing the two types of subgrid structures derived using the global and local methods revealed that the nongeo-located SUs are more consistent than the geo-located SUs across different area threshold values.Further investigation of the relative capability of the two methods, with non-geo-located subgrid units representing multiple topographic features (elevation, slope and aspect) based on the SD in surface elevation within the subgrid units and statistical metrics calculated over the whole study domain, further demonstrated superior capability and consistency for the local method compared to the global method.Similarly, investigation of the relative capability of the two methods in capturing climatic and land cover variability based on the highresolution PRISM precipitation and surface temperature and NDVI data, respectively, reveals that the local method is generally better than the global method.Finally, comparing the precipitation and surface temperature over the study area when represented using non-geo-located SUs from the local method against those of the subbasin-based and original PRISM grid-based generally showed the spatial pattern and statistical values of the subgrid units are much closer to those of the original PRISM grids than those of the subbasins.
In summary, this study demonstrated that adopting the hypsometric curve characterization for discretizing subbasins yields improved capability of capturing topographic heterogeneity and consistency across different values of area threshold.This resulted in improved representation of climatic and land cover variability in land surface modeling.The improved capability of capturing subgrid variability of atmospheric forcing, surface topography and vegetation cover with a nominal increase in computational requirement is essential for improving simulations of land surface modeling in mountainous regions.The focus in this paper is the development and evaluation of the methods and new spatial structures.Future efforts will implement the non-geo-located SUs from the local method in a land surface model based on the Community Land Model subgrid structure to investigate how the addition of topographic subgrid units to the subgrid hierarchical structure may translate to improved simulations of evapotranspiration, soil moisture, snowpack, and runoff and streamflow.Coupling land surface models with the non-geo-located SUs to an atmosphere model with a subgrid parameterization of orographic precipitation may further improve modeling of land-atmosphere interactions in topographically diverse regions.

Code availability
The code developed to generate the subgrid structures is available upon request.Please contact Teklu K. Tesfa at teklu.tesfa@pnnl.gov.
The Supplement related to this article is available online at doi:10.5194/gmd-10-873-2017-supplement.
Competing interests.The authors declare that they have no conflict of interest.

Figure 1 .
Figure 1.The topographic distribution (left) and subbasin delineation (right) of the study area (Columbia River basin).Two subbasins selected to represent the extreme classes of elevation ranges are shown on the far right.

Figure 3 .
Figure 3. Hypsometric curves of two subbasins with extreme contrast of elevation variability discretized into three parts following Willgoose and Hancock (1998) and Sinha-Roy (2002): the head, body and toe, as used in the local method.

Figure 4 .
Figure 4.The number of elevation-based geo-located SUs (LUs) plotted against the average topographic slope for each subbasin derived using the global (a) and local (b) methods with 1 % area threshold.

Figure 5 .
Figure 5. Spatial patterns of the number of elevation-based geo-located SUs per subbasin derived using the global (c) and local (d) methods compared to the spatial pattern of the topographic slope (a) and elevation ranges of the subbasins in the study area.

Figure 6 .
Figure 6.The standard deviation (SD) in elevation within the elevation-based geo-located SUs derived using different values of area threshold.On each box, the central mark (notch) is the median (q2), the edges of the box plot are the 25th (q1) and 75th (q3) percentiles, and the whiskers extend to the most extreme data points (q3 + 1.5x interquartile range (q3 − q1) and q1 − 1.5x interquartile range (q3 − q1); outliers are not considered).

Figure 7 .
Figure 7. Sensitivity of the global (gray) and local (black) methods to different values of area threshold for the total number of SUs (a), average SUs size (b) and standard deviation in SU size (c) of the elevation-based geo-located SUs derived using different values of area threshold.

Figure 8 .
Figure 8.Comparison of the geo-located (black) versus non-geo-located (gray) SUs derived based on elevation, slope and aspect using the global (a-c) and local (d-f) methods, in terms of their sensitivity to different values of area threshold for the total number of SUs, average SU size and standard deviation in SU size over the study domain.

Figure 9 .
Figure 9.Comparison of the two methods (global and local) using non-geo-located SUs in terms of their sensitivity to different values of area threshold for the total number of SUs (a), average SU size (b) and standard deviation in SU size (c) over the study area.SUs are constructed based on elevation, slope and aspect. 11

Figure 10 .
Figure 10.Similar to Fig. 6, but for the capability of the global and local methods to capture topographic heterogeneity based on the standard deviation in elevation, within the non-geo-located SUs derived based on elevation, slope and aspect using different values of area threshold.

Figure 11 .
Figure 11.Similar to Fig. 5, but for the capability of the global and local methods to capture climatic variability based on SD of the PRISM 30-year normal precipitation (a) and surface temperature (b), within the non-geo-located SUs derived based on elevation, slope and aspect across different values of area threshold.

Figure 12 .
Figure 12.Similar to Fig. 6, but for the capability of the global and local methods to capture land cover variation based on standard deviation values of eMODIS NDVI during spring (a) and summer (b), within the non-geo-located SUs based on elevation, slope and aspect across different values of area threshold.

Figure 13 .
Figure 13.PRISM 30-year normal precipitation represented using the subbasins (a) and non-geo-located SUs from the local method using 3 % area threshold (b) compared to those of the original PRISM grids (c).The Canadian territory of the study area is not represented in the PRISM dataset.

Table 1 .
Comparing the local and global methods in capturing topographic heterogeneity when using non-geo-located SUs.
grid surface elevation over the study domain by 17 to 19 % across different values of area threshold compared to those of the global method, highlighting the relative advantages of the local method over the global method in capturing topographic heterogeneity.The improved capability of the local method shown in this comparison comes from the advantage of performing elevation discretization based on hypsometric curve characterization in the local method (see Figs. 6

Table 2 .
Comparing the SUs generated using 3 % area threshold from the local method and subbasin representations against the original PRISM grid representation using statistical summary of mean annual precipitation and surface temperature calculated over the study domain.

Table 3 .
Comparing the SUs generated using 3 % area threshold from the local method and Subbasin representations against the original NDVI grid representations using statistical summary of spring and summer NDVI values calculated over the study domain.