Effects of vegetation heterogeneity and surface topography on spatial scaling of net primary productivity

Due to the heterogeneous nature of the land surface, spatial scaling is an inevitable issue in the development of land models coupled with low-resolution Earth system models (ESMs) for predicting land-atmosphere interactions and carbon-climate feedbacks. In this study, a simple spatial scaling algorithm is developed to correct errors in net primary productivity (NPP) estimates made at a coarse spatial resolution based on sub-pixel information of vegetation heterogeneity and surface topography. An eco-hydrological model BEPS-TerrainLab, which considers both vegetation and topographical effects on the vertical and lateral water flows and the carbon cycle, is used to simulate NPP at 30 m and 1 km resolutions for a 5700 km 2 watershed with an elevation range from 518 m to 3767 m in the Qinling Mountain, Shanxi Province, China. Assuming that the NPP simulated at 30 m resolution represents the reality and that at 1 km resolution is subject to errors due to sub-pixel heterogeneity, a spatial scaling index (SSI) is developed to correct the coarse resolution NPP values pixel by pixel. The agreement between the NPP values at these two resolutions is improved considerably fromR2 = 0.782 toR2 = 0.884 after the correction. The mean bias error (MBE) in NPP modelled at the 1 km resolution is reduced from 14.8 g C m −2 yr−1 to 4.8 g C m−2 yr−1 in comparison with NPP modelled at 30 m resolution, where the mean NPP is 668 g C m −2 yr−1. The range of spatial variations of NPP at 30 m resolution is larger than that at 1 km resolution. Land cover fraction is the most important vegetation factor to be considered in NPP spatial scaling, and slope is the most important topographical factor for NPP spatial scaling especially in mountainous areas, because of its influence on the lateral water r distribution, affecting water table, soil moisture and plant growth. Other factors including leaf area index (LAI) and elevation have small and additive effects on improving the spatial scaling between these two resolutions.


Introduction
Remarkable progress has been made in the development of land models as part of Earth system models (ESMs) over the past few decades for predicting land-atmosphere interactions and carbon-climate feedbacks (Bonan et al., 1993;Sellers et al., 1997;Dai et al., 2003;Hong et al., 2009Hong et al., , 2012)).Land models intended for global climate simulations are often executed at coarse resolutions with or without some simple consideration of the spatial variability within each modelling grid, and therefore it is logical to be concerned with errors in the simulated fluxes from ignoring the within-grid variability.The International Land Model Benchmarking (IL-AMB) project (http://www.ilamb.org)identified a set of standard benchmarks to quantitatively evaluate the performance of land models (Luo, et al., 2012).One benchmark for this purpose would be a set of scaling rules to consider the impact of sub-grid heterogeneity on energy and mass fluxes simulated by land models.
Surface heterogeneity stems from both endogenous (biotic) and exogenous (abiotic) sources (Moorcroft et al., 2001).Endogenous heterogeneity is related to the spatial variations in vegetation type and density, while exogenous heterogeneity is often caused by soil and topographical variations.Many studies have examined spatial scaling issues associated with endogenous heterogeneity (Wu and Qi., 2000;Moorcroft et al., 2001;Maayar and Chen, 2006;Gimona et al., 2006;Liu et al., 2010;Hong et al., 2009Hong et al., , 2012)), while there has been little attention given to exogenous heterogeneity for the purpose of spatial scaling over the land surface.
Spatial scaling refers to a process of taking information at one scale and using it to derive information at another scale (Jarvis, 1995).There can be three levels of complexity in spatial scaling (Chen, 1999): (1) obtaining the mean value of a quantity for a land area based on the areal fractions of the cover type and the mean value of each cover type, i.e. a simple weighting process; (2) deriving a static surface parameter (such as the leaf area index) using an algorithm in consideration of the surface heterogeneity, and in this case the scaling depends on the derivation algorithms; and (3) mapping a dynamic variable or a process (such as net primary productivity (NPP) and evapotranspiration (ET)) that depends not only on the surface endogenous and exogenous heterogeneity, but also other conditions, such as meteorology.Chen (1999) addressed the second level of complexity in scaling a surface parameter based on a contexture method.Simic et al. (2004) and Maayar and Chen (2006) tackled the third level of complexity through scaling the NPP and ET, respectively.Simic et al. (2004) developed an algorithm for correcting errors in NPP in coarse resolution pixels by considering endogenous factors such as sub-pixel land cover fractions and the LAI range.Maayar and Chen (2006) made a step forward in ET spatial scaling by considering not only endogenous, but also exogenous factors including sub-pixel soil texture and topographical variations using a distributed hydrological model.In the present study, we attempt to improve the NPP scaling methodology of Simic et al. (2004) by considering both endogenous and exogenous surface heterogeneity, i.e. in addition to land cover and LAI, topographical parameters are also considered in an eco-hydrological model (Chen et al., 2005).This improvement is necessary because topography is often variable over the land surface affecting water redistribution over the landscape, and there has been little quantitative investigation on the influence of topography on the NPP distribution and the mean value at the landscape, regional and global scales (Simic et al., 2004;Sun, et al., 2004;Bertoldi et al., 2010).
Accurate estimation of NPP at regional, continental and global scales is a necessary step in assessing the possible role of these ecosystems in the global carbon cycle.It is, therefore, crucial to investigate and reduce sources of uncertainties associated with large-scale estimates of NPP.It is shown theoretically that vegetation heterogeneity and surface topography are the important factors introducing biases in surface parameter retrieval and ecological modelling (Chen et al., 1999;Grant, 2004;Maayar and Chen, 2006).Considerable errors can result from the overlook of sub-pixel variability of vegetation (land cover, LAI) and surface topography (slope, aspect, elevation) in modelling carbon and water fluxes at regional and global scales (Arora et al., 2001;Kang et al., 2004;Maayar and Chen, 2006).Because topography influences surface and subsurface water redistributions and therefore soil moisture in the root zone (Qiu et al., 2001), it pronouncedly affects the spatial distribution and the magnitude of vegetation productivity and carbon assimilation due to the strong dependence of leaf stomatal conductance on soil water conditions (Jarvis, 1976).Therefore, topography should also be considered as an important factor in NPP mapping and spatial scaling.The main purpose of this study are: (1) to propose a spatial scaling index (SSI) to integrate the influences of vegetation heterogeneity and surface topography, such as land cover type (LC), leaf area index (LAI), elevation, slope, and their interactions, on NPP estimates at coarse resolutions; (2) to investigate the importance of considering the topographical effects on soil water redistribution and NPP distribution; (3) to investigate the relative importance of various vegetation and topographical factors in NPP scaling and the methodology to combine them for large area applications; and (4) to map and compare the spatial variations of NPP at 30 m resolution and 1 km resolution.

Model used
In this research, the BEPS (Boreal Ecosystem Productivity Simulator) model (Liu et al., 1997(Liu et al., , 1999(Liu et al., , 2002(Liu et al., , 2003) ) is coupled with a hydrological model named TerrainLab (Chen et al., 2005).The coupled model is capable of simulating the spatial distributions of carbon and water fluxes between terrestrial ecosystems and the atmosphere under the influence of topography.The model is driven by remotely sensed vegetation parameters (LAI, land cover type), climate, soil texture, and digital elevation model (DEM).Climatic inputs required to run the model include daily maximum, minimum, mean and dew point temperatures, daily precipitation, and incoming solar radiation.Daily solar radiation was calculated from daily amplitude of air temperature, precipitation, and water vapour pressure.BEPS was developed based on the Forest BioGeochemical Cycles (Forest-BGC) model (Running and Coughlan, 1988) with the following improvements.The daily canopy photosynthesis is simulated by up-scaling instantaneous photosynthesis at the leaf level (Farquhar et al., 1980) through a temporal and spatial integration scheme after sunlit and shaded leaf separation (Chen et al., 1999).This scaling approach was developed to consider the nonlinear effect of the diurnal variation in incoming solar radiation on daily photosynthesis.Radiation intercepted and absorbed by the canopy is simulated as a function of LAI, clumping index and solar zenith angle (Chen and Cihlar, 1995;Liu et al., 1997).Compared with daily gross photosynthesis derived from tower flux data, this method is a large improvement over big-leaf models (Chen et al. 1999;Sprintsin et al., 2012).Through the temporal and spatial integration, this model greatly enhances the computation efficiency and therefore suitable for the applications to large areas.BEPS has been used in different ecosystems with a proven ability to provide reliable estimates of NPP (Amiro et al., 2000;Matsushita and Tamura, 2002;Sun et al., 2004).TerrainLab is a spatially distributed hydrological model, which was designed to provide realistic simulations of daily soil moisture content and water table for modelling carbon balance of terrestrial ecosystems.This model carefully describes topographic effects on the spatial variations of climatic variables and the movement of soil water via a subsurface saturated flow mechanism.The extraction of transpirational water from two soil layers (saturated and unsaturated) is partitioned in terms of the relative abundance of active roots using the methodology of Jackson et al. (1996).Each pixel is linked with its surrounding 8 pixels by saturated subsurface baseflow, which is computed according to local slope and water table (Wigmosta et al., 1994).Terrainlab has been modified and validated in a forested watershed in Canada (Chen et al., 2005).The formulations of BEPS and TerrainLab were introduced in detail in previous publications (Liu et al., 1997(Liu et al., , 1999(Liu et al., , 2002(Liu et al., , 2003;;Chen et al., 1999Chen et al., , 2005;;Sonnentag et al., 2008;Govind et al., 2010).

Simulation procedures
Two sets of simulations were conducted to analyse the effects of vegetation heterogeneity and topography on NPP calculation.In the first simulation, the model was run at the 30 m resolution, whereas the second run was made at the 1 km resolution after all inputs at the 30 m resolution were aggregated to the 1 km resolution.To generate the 1 km simulation datasets, input data, except for land cover types, available at the 30 m resolution were averaged to obtain input datasets at the 1 km resolution (Fig. 1).As usual, the land cover type of each 1 km pixel was assumed to be the dominant type, i.e. having the largest percentage of coverage based on the land cover image at the 30 m resolution.For the NPP run at the 1 km resolution, the vegetation was assumed to be uniform with a unique cover type within each 1 km pixel.NPP distributions at the 1 km resolution are therefore obtained in two ways: (1) by averaging NPP modelled at the 30 m resolution to the 1 km resolution, hereinafter referred to as distributed NPP (NPP d ), and (2) modelling NPP at the 1 km resolution, referred to as lumped NPP (NPP l ).

Algorithms for spatial scaling
The theory used in this research is the extension of contexture-based scaling methodology presented by Chen (1999).The algorithm assumes that NPP d represents the reality, while NPP l is an approximation and can be improved through a scaling procedure.The reason is that the effect of surface heterogeneity within a coarse pixel is mostly preserved in the distributed calculation, while it is not the case in the lumped calculations (Fig. 1).We, therefore, attempt to develop a spatial scaling index (SSI) to integrate the effects of vegetation heterogeneity and topography on NPP estimation, in order to reduce errors in lumped NPP estimation.The relationship between NPP d and NPP l is expressed as follows: SSI is needed for the adjustment of NPP l to the correct distributed value.It can be a function of two scaling indices where SSI s and SSI T represent the scaling index related to vegetation heterogeneity and topography, respectively.For the first index, we further express it as a function of land cover (LC) and leaf area index (LAI), A mathematical function is to be developed to quantify the influence of subpixel land cover and LAI variations on the mean NPP value of the pixel simulated without considering these variations.For the second index, we consider it as a function of elevation (ELE) and slope (SL):  NPP simulations using BEPS-TerrainLab, the mean elevation and slope of each coarse pixel are considered, but these simulations are only approximations of those made at 30 m resolution and need to be corrected with the function.NPP l is biased because input data used in the coarse resolution simulation are the averages in each 1 km 2 pixel and only the ecosystem model parameters associated with the dominant type are used.It is known that the responses of NPP to topographical and LAI variations are nonlinear (Campbell and Norman, 1998;Kenward et al., 2000), and hence NPP simulated using the mean input values would differ from that simulated with ranges of values about the mean values.These differences are to be removed using these scaling indices.For the convenience of mathematical operation, we expressed SSI as follows: where SSI LC , SSI LAI , SSI ELE , and SSI SL , are functions that correct the biases in the lumped calculation of NPP for the individual factors of land cover, LAI, elevation, and slope, respectively.Our scaling algorithm follows the methodology of Maayar and Chen (2006) who assume that the influence of these factors on lumped ET calculations are independent and multiplicative, and therefore can be treated separately in the algorithm development.
3 Site description and data processing

Site description
The experimental site is located in the southwest slope of the Taibaishan Natural Reserve, in the middle of the Qinling Mountain.This mountain significantly alters the climate of interior China and its nearby areas in the southwest of Shanxi Province (31 2).This site has elevations ranging from 518 m to 3767 m above the sea level and a transitional mountainous

Data processing
Inputs data to BEPS-TerrainLab include LAI and land cover maps derived from remote sensing, climate interpolated from measurements at Taibai county weather stations, DEM, and soil texture.All datasets are firstly prepared at the 30 m resolution on a UTM projection for the distributed run and then resampled to the 1 km resolution for the lumped run.LAI was measured at 9 sampling plots (yellow squares in Fig. 2) using the TRAC instrument (Chen and Cihlar, 1995) in June 2004, and tree rings were taken at the same time.The area of each sampling plot is about 30 m × 30 m, which is equal to the grid size of the Landsat TM data.The geographical locations of these sites were obtained using a Magellan global positioning system (GPS) with an accuracy of 5 m.LAI was calculated according to the measured canopy gap fraction after taking into account the element clumping index derived from the measured gap size distribution (Chen and Cihlar, 1995).
The following variables were measured: bark thickness and diameter at 1.3 m breast height (DBH), tree height, crown base height (using a Vertex hypsometer), and crown radius (estimated from beneath the canopy).The measurements of DBH and bark thicknesses had a precision of 1 mm.The height and width of tree crowns were measured with a precision of 0.1 m.Data for validating the simulated annual NPP are obtained from two sources: tree-ring data (9 plots) and inventory data (black dots in Fig. 2).Tree-ring samples were collected with a 5-mm incremental borer applied at DBH.The cores were immediately glued onto wooden mounts and subsequently polished in laboratory with several successively finer sandpapers.Tree-ring series were dated following a standard procedure (Stokes and Smiley, 1968).At least 9 dominant trees were cored at each site in plots with no evidence of tree mortality.Two cores (one from the north and one from the south) were taken at DBH for each tree.In total, about 162 tree cores were collected.Tree age ranged between 24 and 58 yr.
We adopted the relative growth method to measure the biomass (Feng et al., 1999).The biomass of each component of a tree (stem, branch, leaf, root) was estimated from tree height and DBH using a formula derived from field measurements.The increments of DBH were obtained from the tree ring samples, which were then used to estimate the increment of each biomass component of each tree from 2002 to 2003.The annual NPP of a sampling plot was calculated according to the number of trees in each plot and the biomass increment of each tree between these two years, which was estimated from the increment of DBH.In the study area, there was a good relationship between measured LAI and NPP (Fig. 3).Overall, annual NPP increased by 206 g C m −2 yr −1 per unit LAI, and the model was able to capture 84 % of variability in NPP measured using tree rings and inventory data (Fig. 4).
The remote-sensing, climate, soil and DEM data and their processing methods were described in detail in Chen et al. (2007).
The study area was divided into 20 patches to measure soil moisture with Time Domain Reflectometers (TDR) during 1 July to 31 August in 2004.In each patch, one plot containing four sampling points is set up.At each sampling point, soil moisture was measured at three depths: 5, 10 and 20 cm, and a time series of soil moisture values is recorded at every depth at every sampling point.The sampling plots are square or rectangular in shape, ranging from 100 m 2 (10 m × 10 m) to 600 m 2 (20 m × 30 m) in size, and covering different topographic settings, from flat to sloping terrain.In Fig. 5b, soil moisture in situ values at 5 cm depth is the average daily value of each plot.Soil moisture in the unsaturated zone simulated by BEPS-Terrainlab for the sampling plots correlate well with measured moisture at 5 cm depth, with R 2 = 0.84, RMSE =2.5 % (Fig. 5a).The model systematically underestimates the soil moisture possibly because errors in estimating the field capacity based on soil texture.However, the spatial variability of soil moisture is well captured by the model, judging from the close correspondence between modelled and measured values at these plots distributed widely in the watershed (Fig. 5b).The ability of the model in capturing the soil moisture spatial variability gives us confidence in our NPP scaling.

Comparison of inputs and NPP at two spatial resolutions
The statistics of input data (LAI, LC, slope, elevation, aspect) and NPP calculated at the fine and coarse resolutions are shown in Table 1.At the 30 m resolution, three land cover types are distinguished, including coniferous forest (CF), mixed (conifer and broadleaf) forest (MF), and small patches of open land (OL).MF is the dominant land cover type accounting for 60 % of the total watershed area, followed by CF (25 %) and OL (15 %).After aggregating the land cover types to 1 km resolution, the coverage of OL marginally reduced to 12 %, and the proportions of MF and CF changed to 70 % and 18 %, respectively.The value of open land LAI is zero.In general, the average of mixed forest LAI was 3.1 with a range from 1.5 to 3.4, and the average of conifer forest was 3.7 with a range from 3.2 to 5.9.The LAI of most pixels are about 4.0.The slope at the 30 m resolution varied with a mean of 28.7 • and a maximum of 85.4 • .However, the slope at the 1 km resolution varied with a mean of 26.6 • and a maximum of 43.5 • , suggesting a considerable loss of slope information at the coarse resolution.The elevation varied between a minimum of 518 m and a maximum of 3767 m at the 30 m resolution, and varied between a minimum of 836 m and a maximum of 3159 m at the 1 km resolution.The elevation range is larger at 30 m (3249 m) than at 1 km (2323 m).In general, the mean values of the corresponding datasets shown in Table 1 are very similar, but the ranges are reduced for all variables at the coarse resolution.
Figure 6 shows a comparison of spatial distributions of NPP at the 30 m and 1 km resolution.Although NPP values at 30 m and 1 km resolution show similar spatial distributions in general, the spatial variations of NPP at the 30 m resolution were more pronounced than that at the 1 km resolution.Moreover, NPP values simulated at the 1 km resolution were smaller than those at the 30 m resolution in the southern part of the study area.

Relationship between NPP d and NPP l before correction
Lumped NPP values compare relatively well with distributed NPP values (R 2 = 0.78) for all pixels (Fig. 7).The strongest relationship is observed within conifer-labelled pixels (R 2 = 0.81), and mixed-labelled pixels also have

Scaling based on land cover type
The algorithm is presented as follows: where f LC ij is the area fraction of a non-dominant cover type i in a pixel labelled as land cover type j ; C LC ij is the correction coefficient for a non-dominant cover type i present in a coarse pixel (1 km in this study) labelled as (or dominated by) land cover type j , and n is the number of non-dominant land cover types within that pixel.Because a coarse resolution pixel is treated with only one dominant cover type in the lumped calculation, the fractions of non-dominant cover types that exist in the pixel can cause errors in the simulated mean NPP for the pixel because of the differences in NPP between the non-dominant cover types and the dominant cover type.These differences are corrected using C LC ij (Table 2), which is obtained through the regression of NPP d /NPP l against the fraction of the cover type i for pixels dominated by cover typej (Fig. 8).C LC ij is taken as the slope of the regression as it signifies the deviation of NPP d /NPP l from 1 as the non-dominant cover type fraction increases.Figure 8a-b are for the 1 km pixels dominated by CF, Fig. 8c-d for the 1 km pixels dominated by MF, and Fig. 8e-f for the 1 km pixels dominated by OL.It can be seen that the ratio between NPP d and NPP l become less variable with increasing areal fraction of a cover type, in particular conifer and mixed forests.In CF dominated pixels, the mean value of this ratio does not change much with MF fraction (Fig. 8a) because the NPP values of MF and CF are similar.The same reason can explain the small variation of the ratio with CF fraction in MF dominated pixels (Fig. 8c).For CF and MF dominated pixels, this ratio decreases with the OL fraction in these pixels (Fig. 8b and d) because OL has lower NPP than CF and MF.In OL dominated pixels (Fig. 8e and f) NPP d /NPP l increases with increasing MF and CF fractions, because CF and MF have much larger NPP values than OL.
Figure 9 shows the outcome of the land cover type correction for NPP l .The correction significantly improved the agreement between NPP l and NPP d , particularly for open land (Fig. 9c) with R 2 improved from 0.26 to 0.52.However, a number of pixels containing large fractions of open land are widely scattered.As discussed in Sect.4.2, this suggests that it is necessary to further consider LAI variability (see Sect. 4.4).

Scaling based on both cover type and LAI
To correct the errors caused by aggregating LAI, an additional LAI correction (SSI LAI ) was developed based on the relative deviation of LAI values of the dominant land cover type in a lumped pixel from the mean LAI value of the pixel.
where L i is the relative deviation of LAI values for a nondominant cover type i within a lumped pixel; L j is the relative deviation of LAI values for the dominant cover j within a lumped pixel; and α is a nonlinearity factor to describe the NPP-LAI relationship for the dominant land cover type    becomes more nonlinear at higher LAI values.As shown in Eq. ( 9), the LAI correction depends on both the variability (L i ) and the mean value (L j ) of LAI in the pixel.Similar to land cover correction, the correction for the LAI variability for a non-dominant cover type is proportional to the areal fraction of the non-dominant cover type and the difference in NPP between the non-dominant and the dominant cover type.
If the LAI-NPP relationship is linear, α is zero.It should generally be between 0 and 1.The importance of LAI variations within the same cover type in the scaling process is demonstrated in Fig. 10.Although there is only a slight improvement in the correlation between NPP l and NPP d , pixels containing large fractions of open land fall closely to the regression line after the addi-tional LAI correction, suggesting that LAI plays a significant role in capturing the interaction between open land and other cover types.

Scaling based on slope and elevation
According to the pattern of NPP of forest (excluding OL) changing with slope at the 1 km resolution (Fig. 11), slope was classified into three ranges, that is 0-20  The algorithm to correct the effects of slope on NPP is similar to that used in the LC correction, i.e.
where f SL ij is the areal fraction of a non-dominant slope range i in a pixel labelled as slope range j ; C SL ij (Table 2) is a regression coefficient for a particular non-dominant slope range i present in a coarse pixel labelled as (or dominated by) slope range j ; and n is the number of non-dominant slope ranges within that pixel.The logic of Eq. ( 11) is similar to that of Eq. ( 8), because the impact of the distorted slope range on the lumped NPP calculation depends on the areal fraction of the non-dominant range and the difference in NPP between the non-dominant and dominant slope range.
The correlation between NPP d and NPP l is improved after SSI SL correction (Fig. 11), especially after combining it with LC and LAI corrections for open land pixels (Table 3) in comparison with the case before any correction, indicating that slope is a relative dominant topographic factor influencing NPP spatial scaling.
Because of NPP peaks at locations at about 1350 m above sea level and decreases with further increase of elevation from this height (Chen et al., 2007), two elevation sections were analysed as follows: where f ELE ij is the areal fraction of a non-dominant elevation range i in a pixel labelled as elevation range j at the 1 km resolution; C ELE ij (Table 2) is a regression coefficient for a particular non-dominant elevation range i present in the coarse pixel labelled as (or dominated by) elevation range j ; and n is the number of non-dominant elevation range within that pixel.The SSI ELE correction based on Eq. ( 13) improved by a small extent the correlation between NPP d and NPP l (Fig. 12), indicating that elevation is also a significant topographic factor to consider in NPP spatial scaling.The effectiveness of individual corrections based on slope, elevation and their combinations with LAI + LC is given in Table 3.In general, the independent elevation correction is smaller than the slope correction.However, the elevation correction combined with LAI + LC is slightly larger than that of slope correction combined with LAI + LC except for open land pixels, suggesting a low level interaction of these factors in influencing the NPP spatial pattern.

Scaling based on all factors
Large improvements of NPP l are achieved for OL, CF and DF pixels after applying corrections for all factors (land cover, LAI, slope, and elevation) (Fig. 13, and Tables 3  and 4).The improvements are particularly significant for OL pixels.The mean bias error (MBE) in NPP modelled at the 1 km resolution is reduced from 14. 4.8 g C m −2 yr −1 in comparison with NPP modelled at 30 m resolution, where the mean NPP is 668 g C m −2 yr −1 (Table 4), and about 68 % mean bias error is reduced for the entire watershed with these corrections.In term of percentage, the NPP l MBE was reduced from 2 % to 0.7 %, the total correction resulted in a decrease of average mean bias error (MBE) in NPP l from 1.9 % to 0.6 % for the coniferous forest and from 2 % to 0.1 % for mixed forest (Table 4).Table 4 also indicates that these corrections improved minimum and maximum of NPP l values.Vegetation heterogeneity and surface topography were of almost equal importance in NPP spatial scaling (Table 3) 5 Discussion Ambroise (1995) pointed out that large topography-driven lateral redistributions of water and considerable topographyrelated heterogeneities on all scales are the two main characteristics of mountainous regions.Topography alters the water cycle through its direct effects on the spatial patterns of precipitation and potential ET and its indirect effect on the spatial patterns of soil water content and streamflow generation.All these spatial patterns exert influences on the spatial pattern of vegetation productivity.On a slope, the amount of rainfall reaching the ground varies according to the slope and aspect (Fig. 14), and the general tendency of increasing soil moisture downslope is modulated by variations in the depth of saturated or less-permeable layers.This tendency can be amplified by a convergence of flows in laterally   concave zones or attenuated by their divergence in laterally convex zones (Dunne and Black, 1970;Anderson and Burt, 1978).The wettest conditions thus usually occur on shallow soils and in convergence zones at the bottom of concave slops (Fig. 14).The spatial pattern of soil moisture, which is controlled by the topography, determines the NPP spatial pattern to a large extent.One of the main objectives of coupling BEPS with Terrainlab is to estimate soil moisture precisely (Chen et al., 2005), because a reliable simulation of soil moisture is a must for reliable simulations of the net carbon exchange between terrestrial ecosystems and the atmosphere.The water redistribution in the soil is controlled by topography (Moore et al., 1988;Grayson et al., 1997;Western et al., 1999;Qiu et al., 2001), and BEPS-Terrainlab can therefore simulate the spatial variations of carbon and water fluxes under the influence of topography.NPP is sensitive to soil moisture in the rooting zone, and therefore the NPP spatial variation with slope is mostly induced by the influence of slope on soil moisture in the rooting zone.This may be the reason that slope is an important factor to consider in NPP spatial scaling.The resulting variations in soil moisture in the rooting zone could be used to evaluate variations in NPP, and the soil moisture scaling will be valuable to interpret NPP spatial scaling (Figs. 15 and 16).
As shown in Fig. 15, NPP at 30 m resolution increases with slope from 6 • to 25 • , varies slightly from 26 • to 60 • , and then decreases considerably from 61 • to larger slopes, while NPP at 1 km resolution increases with slope from 6 • to 20 • , decreases slightly from 21 • to 40 • , and then decreases sharply from 41 • to larger slopes.Corresponding to these NPP variation patterns, soil moisture also shows the similar patterns at both 30 m and 1 km resolutions, respectively (Fig. 16).Specifically, at 30 m resolution, the average soil moisture in the rooting zone during the growing season shows the following spatial patterns: (i) it has the largest values at locations with small slopes (generally at the bottom of hills) due to convergence of water from upper slopes; (ii) it has medium  values at locations with medium slopes at which the gain of water from upper slopes may be balanced by the transfer of water to lower slopes; and (iii) it decreases with increasing slope after a threshold value of about 63 • at which ground water is drained to lower positions rapidly after rainfall.At Fig. 15.NPP variation with slope at 30 m and 1 km resolutions.The range of slope at 1 km resolution is plotted at the same scale as that at 30 m, but is much smaller because of the averaging process. 1 km resolution, the largest slopes disappear due to the averaging process, and this distribution pattern is distorted.The correction of NPP l based on slope is to remove the difference in NPP caused by the difference in this soil moisture spatial pattern between these two resolutions.In Fig. 15, the NPP values at small slopes are smaller than those at medium range slopes because the soil was usually saturated in flat areas (Fig. 16) and the plant growth is impeded under the partially anaerobic conditions (Chen et al., 2005).This suggests that, on the one hand, a reliable NPP simulation depends closely on reliable soil moisture simulation, and on the other hand, slope is the most important topographic factor affecting water redistribution and soil moisture, and hence affecting the NPP spatial distribution.Slope is, therefore, an important topographic factor to be included in NPP spatial scaling in mountainous areas.
In general, the changes in NPP and soil moisture with slope are larger at finer resolutions (at the 30 m resolution, NPP ranges from 443 to 642; soil moisture ranges from 14.2 % to 24.0 %) than at coarser resolutions (at the 1 km resolution, NPP ranges from 573.4 to 652.4; soil moisture ranges from 17.1 % to 24.0 %).In addition to the topographical effect on water redistribution, forest was also disturbed by human activities, such as timber harvest and clear-cutting for farming.Therefore, NPP values are low in the locations with the lowest slopes.At the highest slopes, soil moisture is the lowest (Fig. 16), and the effect of limited available water on forest growth causes NPP to be the lowest.In addition, largest slopes generally occur at high elevations, and low temperature becomes a negative factor in tree growth.Therefore, NPP values are often the lowest at locations with the highest slopes.NPP values are high at locations with medium slopes where the soil moisture is optimum for tree growth.
Our results show that both exogenous factors, such as slope and elevation, and endogenous factors, such as land cover and LAI, are useful for NPP spatial scaling, i.e. to apply information at one scale to another.The spatial distributions of endogenous factors may be regarded as integrated outcomes of all environmental conditions including climate, soil, topography and disturbance and therefore are theoretically most important factors to consider in spatial scaling as they intrinsically include the accumulated influences of endogenous factors over a large time scale (a few years to a few decades).In this regard, the methodology adopted in some land models (Bonan et al., 1993;Chen et al., 2012) to compute fluxes for various cover type fractions having different mean LAI values within a grid is an effective way to capture the first-order effects of within grid heterogeneity on the grid-level computation.However, exogenous factors have strong influences on the various instantaneous fluxes between land and atmosphere due to the redistribution of water over the landscape and their interactions with meteorological variables (at least radiation, but also often precipitation, temperature and humidity).Moreover, eco-hydrological models attempted for hydrological and ecological simulations at coarse resolutions would suffer from the inevitable distortion to the topographical parameter fields, and this would also be an issue for spatial scaling.Therefore, both endogenous and exogenous factors should be integrated within a spatial scaling scheme for the purpose of removing biases in the coarse resolution results, and this scheme may be included as a benchmark of land models.The methodology and theories proposed in this study will be helpful for developing this scheme.
In this study, soil texture and depth are not included in the scaling algorithm, although they are considered in the hydro-ecological model.Soil texture has strong influence on the land cover type (e.g.pine forests usually grow on sandy soils), and therefore its influence might have been mostly included in the scaling using land cover information.Soil depth is not considered in this study because it is over 0.5 m at most locations in the watershed and therefore not a strong limiting factor for plant growth.In different landscape settings, it may be worthwhile to consider these two additional parameters.
The scaling algorithm developed in this study is transferrable to other landscapes.However, some of the coefficients in the algorithm need to be adjusted for local values, such as the mean differences in the flux of concern (NPP in this study) among different cover types and topographical classes.These differences can be estimated using a processbased model or a simple model depending on the required www.biogeosciences.net/10/4879/2013/Biogeosciences, 10, 4879-4896, 2013 accuracy for the flux.The application of our spatial scaling methodology requires high-resolution maps of land cover, LAI, slope and elevation.Attention should be given to the accuracy of these high-resolution maps as their errors would propagate to the final scaling results.

Conclusions
The BEPS-TerrainLab model, which is an eco-hydrological model coupling an ecosystem model BEPS with a distributed hydrological model TerrainLab, was applied to a mountainous forested watershed, in Baohe River Basin in China.
Through studying the NPP results simulated at fine and coarse resolutions, we draw the following conclusions regarding to NPP spatial scaling methodology for the purpose of correcting errors in NPP modelling at coarse resolutions: 1.A spatial scaling index (SSI), which integrates the effects of vegetation heterogeneity (land cover, leaf area index) and topographical variations (elevation, slope) on NPP, is effective for correcting the errors in coarse resolution NPP estimation.
2. Vegetation heterogeneity and topography are of almost equal importance in determining the NPP spatial distribution pattern and therefore its spatial scaling.The effect of vegetation heterogeneity on NPP spatial scaling is generally larger than that of topography based on the correlations between the distributed NPP and lumped NPP.
3. The accuracy in NPP spatial scaling is incremental as more factors are considered, suggesting that the effects of these factors can be considered individually.
4. LAI and slope are especially significant factors for removing errors in coarse pixel NPP estimation for open land.Land cover is generally more important than LAI in NPP spatial scaling.In scaling with topographical parameters, slope is generally more effective than elevation.
5. NPP spatial scaling in complex terrains depends much on the amount of the distortion of the soil moisture field at the coarse resolution.Topographical variations are often greatly suppressed as the resolution decreases, causing various biases in the estimation of radiation and water flows.The spatial redistribution and movement of ground water in complex terrains tightly control the NPP distribution.Results of this study suggest that it is indeed necessary to consider topography in NPP modelling and spatial scaling in addition to considering vegetation heterogeneity, although such effort would require the use of a distributed hydrological model.
6.In general, both endogenous and exogenous factors are useful for spatial scaling of land surface fluxes and they can be integrated into a comprehensive scaling scheme for considering the influences on sub-grid heterogeneity on the grid level computation.A scaling scheme of this type may be used to produce a new benchmark for evaluating land models.

Fig. 2 .
Fig. 2. Study area with the locations of ground sites in Qinling Mountain in southern Shanxi Province in China.
climate of a warm temperate region in the northern sub-tropic region.The annual precipitation is about 630 mm at low elevations and increases to about 1000 mm at the top of the Taibaishan Mountain.The mean annual air temperature decreases from +9.5 • C at the bottom to −4.5 • C at the top of the mountain.Influenced by topography, the Qinling Mountain has different vegetation covers at different elevations and terrain aspects, including pine (Pinus armandi Franch) and broadleaf mixed forest at elevations from 1330 to 3420 m on the northwest slope, conifer (Abies fabric (Mast.)Graib., Pinus tabulaeformis Carr.) and broadleaf mixed forest at elevations from 518 m to 1300 m on the north slope, conifer forest in the alpine zone at elevations from 1250 m to 3400 m on the west slope, and pine (Pinus armandi Franch) and broadleaf mixed forest (Quercus spp.) at elevations from 912 m to 2428 m on the south slope.The pine and broadleaf mixed forests are the dominant vegetation types.The vegetation density varies with topography, providing large dynamic ranges for analysing the effect of vegetation heterogeneity and topography on the spatial scaling of NPP.The forest ecosystems are disturbed in various ways because the local people cut trees and shrubs at low altitudes to feed live stocks and to make furniture and building materials.Areas with low slopes up to 20 • near rivers are mostly converted to irrigated farmland. 1

Figure 3 .
Figure 3. Relationship between measured annual NPP and measured LAI.NPP was 2

Fig. 3 .Figure 4 .Fig. 4 .
Fig. 3. Relationship between measured annual NPP and measured LAI.NPP was calculated based on tree ring increment and biomass measured in the field, and LAI was measured using TRAC in 2004.

Fig. 5 .
Fig. 5. Simulated soil moisture (SM) at 30 m resolution using BEPS-TerrainLab in comparison with measured SM using TDR during 1 July to 31 August in 2004 (a); comparison of BEPSsimulated SM with measured SM at various sites (b).The solid line is the 1 : 1 line.

Figure 6 .Fig. 6 .
Figure 6.Comparison of s patial distribution of NPP simulated at 30 m resolution and 1 km resolution

Fig. 7 .
Fig. 7. Regression between distributed (NPP d ) and lumped (NPP l ) annual values before correction: (a) pixels dominated by the coniferous forest (CF-pixels) at the 1 km resolution; (b) pixels dominated by the mixed forest (MF-pixels) at the 1 km resolution; (c) pixels dominated by open land (OL-pixels) at the 1 km resolution; and (d) all pixels for all land cover types (all-pixels).

Figure 8 .Fig. 8 .
Figure 8. Variation of average NPP d /NPP l as a function of the fraction of

Fig. 10 .
Fig. 10.Regression between distributed (NPP d ) and lumped (NPP l ) annual values after land cover and LAI (SSI LC+LAI ) correction: (a) pixels dominated by the coniferous forest (CF-pixels); (b) pixels dominated by the mixed forest (MF-pixels); (c) pixels dominated by open land (OL-pixels); and (d) all pixels for all land cover types (all-pixels).

Figure 11 .Fig. 11 .
Figure 11.Regression between distributed (NPP d ) and lumped (NPP l ) annual values , similar to the finding of Maayer and Chen (2006) for ET scaling.The combination of three factors of land cover, LAI and slope (SSI LC + LAI + SL ) made the largest improvement to NPP l for open land (R 2 = 0.7118), suggesting that slope is especially a significant topographic factor influencing NPP spatial scaling in a region with mixed open land and other cover types.The result shows that different factors of vegetation heterogeneity and topography have different contribu-tions to NPP l correction.These contributions are generally additive, i.e. larger improvements are achieved with more factors considered.

Fig. 12 .
Fig. 12. Regression between distributed (NPP d ) and lumped (NPP l ) annual values after elevation (SSI ELE ) correction: (a) pixels dominated by the coniferous forest (CF-pixels); (b) pixels dominated by the mixed forest (MF-pixels); (c) pixels dominated by open land (OL-pixels); and (d) all pixels for all land cover types (all-pixels).

Figure 14 .
Figure 14.Schematic diagram illustrating the processes of water redistribution at 5

Fig. 14 .
Fig. 14.Schematic diagram illustrating the processes of water redistribution at different positions of a slope.

Fig. 16 .
Fig. 16.Variations of the average soil moisture in the rooting zone with slope during the growing season at 30 m and 1 km resolutions.

Table 1 .
Comparison of inputs and NPP results at 30 m and 1 km resolutions (CF: coniferous forest; MF: Mixed forest; OL: open land).NPP d denotes distributed NPP at 1 km resolution based on averages of 30 m resolution NPP values.NPP l denotes lumped NPP.
* Open land fraction excluded.

Table 2 .
C LC The subscripts con, mix, and ol refer to coniferous forest, mixed forest, open land, respectively.The ELE low , ELE high , SL low , SL mid and SL high refer to elevation less than and equal to 1350 m, greater than 1350 m , slope 6-20 • , slope 21-40 • and slope 41-45 • , respectively.
ij , C SL ij , C ELE ij coefficients used in the land cover correction (SSI LC ), elevation correction (SSI ELE ), and slope correction (SSI SL ).

Table 3 .
Comparison of the correlations between NPP d and NPP l before and after various corrections.

Table 4 .
Statistics of the distributed and lumped NPP.Shown are averaged values of coniferous and mixed forest excluding open land in 2003.MBE is the mean bias error.All values are expressed in g C m −2 yr −1 (CF: coniferous forest; MF: Mixed forest; OL: open land).