Spatial Downscaling of Gross Primary Productivity Using Topographic and Vegetation Heterogeneity Information: A Case Study in the Gongga Mountain Region of China

Abstract

Due to the spatial heterogeneity of land surfaces, downscaling is an important issue in the development of carbon cycle models when evaluating the role of ecosystems in the global carbon cycle. In this study, a downscaling algorithm was developed to model gross primary productivity (GPP) at 500 m in a time series over rugged terrain, which considered the effects of spatial heterogeneity on carbon flux simulations. This work was carried out for a mountainous area with an altitude ranging from 2606 to 4744 m over the Gongga Mountain (Sichuan Province, China). In addition, the Moderate Resolution Imaging Spectroradiometer (MODIS) GPP product at 1 km served as the primary dataset for the downscaling algorithm, and the 500 m MODIS GPP product was used as the reference dataset to evaluate the downscaled GPP results. Moreover, in order to illustrate the advantages and benefits of the proposed downscaling method, the downscaled results in this work, along with ordinary kriging downscaled results, spline downscaled results and inverse distance weighted (IDW) downscaled results, were compared to the MODIS GPP at 500 m. The results showed that (1) the GPP difference between the 500 m MODIS GPP and the proposed downscaled GPP results was primarily in the range of [−1, 1], showing that both vegetation heterogeneity factors (i.e., LAI) and topographic factors (i.e., altitude, slope and aspect) were useful for GPP downscaling; (2) the proposed downscaled results (R² = 0.89, RMSE = 1.03) had a stronger consistency with the 500 m MODIS GPP than those of the ordinary kriging downscaled results (R² = 0.43, RMSE = 1.36), the spline downscaled results (R² = 0.40, RMSE = 1.50) and the IDW downscaled results (R² = 0.42, RMSE = 1.10) for all Julian days; and (3) the inconsistency between MODIS GPP at 500 m and 1 km increased with the increase in altitude and slope. The proposed downscaling algorithm could provide a reference when considering the effects of spatial heterogeneity on carbon flux simulations and retrieving other fine resolution ecological-physiology parameters (e.g., net primary productivity and evaporation) over topographically complex terrains.

1. Introduction

The gross primary productivity (GPP), which corresponds to carbon fixation by vegetation at the ecosystem level, is an important indicator to assess the photosynthetic capacity of vegetation and the function of the ecosystem. The accurate estimation of regional, continental and global GPP plays an important role when monitoring vegetation growth conditions [1], the terrestrial carbon budget [2] and the interactions in the soil-vegetation-atmosphere continuum [3].

Carbon cycle models have been developed over many years to predict carbon-climate feedbacks at multiple scale levels [4,5,6], which can be used to estimate spatially and temporally continuous GPP. However, most carbon cycle models are generated from the site scale, without considering the spatial heterogeneity within each modeling grid [7]. Moreover, due to the paucity of the data and computing complexities, the carbon cycle model is usually executed at coarse resolutions, which are based on the simplification of landscape complexities. Such a simplification results in large uncertainties due to spatial heterogeneity, especially in topographically complex terrains. Some researchers have shown that the oversimplification of landscape complexities may inevitably cause model simulations to be considerably biased [8,9,10]. Landscape complexities typically present high spatial heterogeneity, and spatial heterogeneity is scaled by multiple factors, including endogenous and exogenous factors [11,12]. Endogenous heterogeneity refers to the spatial variability of vegetation types and density (i.e., land cover and leaf area index), while exogenous heterogeneity is associated with surface topography (i.e., altitude, slope and aspect). Therefore, it is crucial to consider the effects of spatial heterogeneity on carbon flux simulation by models.

The carbon flux simulation through carbon cycle models usually requires massive forcing data. Atmospheric data are one of the most important inputs for carbon cycle model simulation. Theoretically, atmospheric forcing data, such as air temperature, air humidity and solar radiation, describe the external environment for the carbon cycle process and exhibit strong spatial variation characteristics over complex terrain [13,14,15]. Dense field observations are required to determine the spatial patterns of atmospheric forcing data, especially for rugged regions. However, due to insufficient meteorological stations located in mountainous area, the atmospheric forcing data could not be provided accurately at fine resolution [12]. To overcome the limitation of field meteorological observations for running the carbon cycle model, downscaling is a potential method to obtain high resolution carbon flux simulation in topographically complex terrains. Using the relationships between carbon flux estimation and the surrounding landscape, the disaggregation of low spatial resolution carbon flux estimation and the information from ancillary data could be introduced in the downscaling scheme.

Spatial scaling is a process of using information available at one scale to derive processes at another scale, including upscaling and downscaling. Downscaling refers to an increase in spatial scale following disaggregation of the coarse-resolution dataset and information from ancillary data at a finer resolution, which restores the variations at a finer scale by assuming the values of the coarser are the mean values at the finer scale. In the last decade, the downscaling of the carbon cycle process has been one of the most challenging issues in environmental science, especially over mountainous areas [16,17,18]. In the literature, some studies have been carried out to examine the effects of spatial heterogeneity on the upscaling process of carbon fluxes [8,19], while little attention has been shown towards the issue of downscaling carbon fluxes over complex terrains. Therefore, it is necessary to downscale carbon fluxes using the subpixel information of topography and vegetation heterogeneity, which is an important issue when evaluating the role of ecosystems in the global carbon cycle.

Vegetation heterogeneity and surface topography are important factors that introduce biases into carbon modeling and water fluxes. Vegetation heterogeneity directly influences the photosynthetic capacity per unit of surface area. Topography affects the redistribution of water, radiation and heat, which have significant influences on the process of carbon assimilation. In this paper, the main purposes are (1) to develop a GPP downscaling algorithm using the subpixel information of topography and vegetation heterogeneity and (2) to assess the GPP inconsistency at 1 km and 500 m, to analyze the relationship between spatial heterogeneity and the GPP difference at two scales.

2. Materials and Methods

2.1. Study Area

As shown in Figure 1, the experimental region is northwest of Gongga Mountain on the quaternary sections of the Tibetan plateau. The summit of Gongga Mountain is 7556 m above sea level, which is the highest mountain in the Sichuan Province of China. Mt. Gongga has a mountainous climate that ranges from a cool plateau climate to a subtropical lowland climate [20]. The annual precipitation is approximately 3000 mm at the summit and declines to 1300 mm at lower altitudes [21]. The average annual air temperature decrease with altitude, from approximately 11 °C at lower altitudes to −1 °C at higher altitudes [21]. This experimental region, with an area of 2584 km², has an altitude ranging from 2606 to 4744 m above sea level, and the dominant types of land cover are forestlands (i.e., evergreen forests and mixed forests), grasslands and shrublands.

2.2. Data and Processing

2.2.1. GPP Images at 500 m and 1 km

The Moderate Resolution Imaging Spectroradiometer (MODIS) has produced a regular estimation of GPP since 2000, including the MOD17A2 (8-day, 1 km) and MOD17A2H (8-day, 500 m) products. The algorithm for the MODIS GPP (gC m⁻² day⁻¹) is based on the concept of radiation conversion efficiency and can be described as follows:

$$\mathrm{GPP}={\epsilon }_{\max }\times TMIN\_scalar\times VPD\_scalar\times SWRad\times 0.45\times FPAR$$

(1)

where ε_max represents the maximum light use efficiency (gC MJ⁻¹) related to the plant functional type (PFT), which can be obtained from the Biome Properties Look-Up Table (BPLUT); TMIN_scalar and VPD_scalar are the scalars for minimum air temperature and vapor pressure deficits (VPD), respectively, ranging from 0 to 1; SWRad represents the incoming short wave radiation (MJ m⁻² day⁻¹), and FPAR represents the fraction of absorbed photosynthetic active radiation, which can be obtained from MODIS product. The MODIS GPP were scaled by multiple environmental factors, and each environmental factor has a comparable level of influence on the final simulated GPP. More details regarding the MODIS GPP algorithm can be found in the literature [22,23]. In addition, some research has been done to evaluate MODIS GPP products across multiple biomes, suggesting that the MODIS GPP works effectively for the majority of PFTs [23,24,25].

The MODIS GPP products MOD17A2 (version 5.0) and MOD17A2H (version 6.0) used in our study were obtained from Julian day 169, 2010 to Julian day 209, 2010, which can be freely acquired at National Aeronautics and Space Administration Distributed Active Archive Center (http://www.modis.ornl.gov/modis/index.cfm). Using spatial heterogeneity information within each 1 km modeling grid, our study was designed to estimate downscaled GPP from a 1 km resolution to a 500 m resolution. Because the 1 km GPP product is executed at a coarse resolution without considering the spatial heterogeneity within each 1 km modeling grid, the GPP product at 1 km served as the primary dataset for the downscaling algorithm in this work. The 500 m GPP product was used as a reference dataset to evaluate the downscaled GPP results at 500 m.

2.2.2. LAI Data

The MODIS LAI products MCD15A2 (version 5.0) at 1 km and MCD15A2H (version 6.0) at 500 m, were selected to describe the vegetation density in this study. The 8-day MODIS LAI product is available at http://wist.echo.nasa.gov in the sinusoidal projection. According to the quality control layer (FparLai-QC) involved in this product, the MODIS LAI consists of a main look-up-table (LUT) algorithm (QC < 64) and a back-up algorithm (64 ≤ QC < 128). The LUT algorithm exploits surface reflectance information in the red and near-infrared bands based on a 3-D radiative transfer equation. When the LUT algorithm fails, the back-up algorithm adopts NDVI-LAI empirical relationships to estimate the LAI [26]. The LUT algorithm usually has a higher reliability than that of the back-up algorithm [26]; the clumping effects at the canopy scale are considered in the LUT algorithm [27].

2.2.3. DEM Data

Some researches have illustrated that the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) global digital elevation model (GDEM) data have a good elevation accuracy over mountainous area [28,29,30], ASTER GDEM were adopted to describe the complex terrain in this work. GDEM data available at a 30 m were averaged to obtain the values at 500 m and 1 km resolutions [7]. The surface slope and aspect were derived from the GDEM data, according to a fast algorithm based on a 3 × 3 pixels window [31].

2.3. Algorithm for Downscaling

As shown in Figure 2, the proposed downscaling algorithm contains two phases: the regression process and area-to-point kriging (ATPK) process. It firstly conducts regression analysis between GPP and heterogeneity factors (i.e., altitude, slope, aspect and LAI), and then ATPK is performed on the downscaling residuals from the regression process.

2.3.1. Problem Formulation

Let GPP_R(x_j) represent the random variable of pixel R centered at x_j (j = 1, 2, …, N, where N is the number of pixels) in 1 km resolution image of GPP, whereas GPP_r(x_i) represents the random variable of pixel r centered at x_i (i = 1, 2, …, M, where M is the number of pixels) in the 500 m resolution image of GPP. The notations R and r represent the 1 km and 500 m pixels, respectively. The purpose of the downscaling method based on topographic factors and the LAI is to predict GPP_r(x) at a 500 m resolution. GPP_r(x) can be described as

$$GP{P}_r\left(x\right)=GP{P}_{r1}\left(x\right)+GP{P}_{r2}\left(x\right)$$

(2)

where GPP_r1(x) and GPP_r2(x) represent the predictions of the regression and ATPK processes, respectively. More descriptions about their calculation can be found in the following sections.

2.3.2. Regression between GPP and Topographic Factors

The regression takes advantage of valuable spatial heterogeneity information from the ancillary data at 500 m spatial resolution. The relationship between topographic factors and MODIS GPP at 1 km is modeled by multiple linear regression:

$$GP{P}_{R\_obs}\left(x\right)=b_0+b_1{H}_R+b_{2}COS(S_R)+b_{3}COS(A_R)+b_{4}LA{I}_R$$

(3)

where GPPR_obs(x), H_R, S_R, A_R and LAI_R represent the MODIS GPP(gC d⁻¹), altitude (m), slope (deg), aspect (deg) and LAI (m² m⁻²) at a 1 km resolution, respectively. Since GPPR_obs(x), H_R, S_R, A_R and LAI_R in Equation (3) are known, b₀, b₁, b₂, b₃ and b₄ can be estimated. Currently, a wide range of fitting methods has been applied in the theoretical framework of regression, such as generalized least squares (GLS) and ordinary least squares (OLS). In this paper, the OLS algorithm was adopted to estimate the regression coefficients.

The relationship in Equation (3) is assumed to be universal at different resolutions [32,33], and thus the relationship modeled at a 1 km resolution can be implemented at a 500 m resolution. Based on this assumption, the regression prediction GPP_r1(x) is calculated as

$$GP{P}_{r1}\left(x\right)=b_0+b_1{H}_r+b_{2}COS(S_r)+b_{3}COS(A_r)+b_{4}LA{I}_r$$

(4)

where H_r, S_r, A_r and LAI_r represent the altitude (m), slope (deg), aspect (deg) and LAI (m² m⁻²) at a 500 m resolution, respectively. Variables b₀, b₁, b₂, b₃ and b₄ are regression coefficients.

2.3.3. ATPK for Downscaling Residuals

If the result of the first phase is perfect, there should be no bias between it and the original GPP at a 1 km resolution. However, it is unrealistic to obtain such an ideal result, and there are inevitable residuals from the first process. The residuals at a 1 km resolution are calculated as:

$$GP{P}_{R\_residual}\left(x\right)=GP{P}_{R\_obs}\left(x\right)-(b_0+b_1{H}_R+b_{2}COS(S_R)+b_{3}COS(A_R)+b_{4}LA{I}_R)$$

(5)

In this work, the residuals at a 1 km resolution were considered for downscaling, and the ATPK process served as the second phase when downscaling the residuals, GPP_{R_residual}(x), at 1 km to those at 500 m, GPP_{r_residual}(x). In the ATPK approach, GPP_r2(x) represents a linear combination of residuals from N 1 km pixel, which can be described as:

$$GP{P}_{\mathrm{r}2}\left(x\right)={\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{N}}{\lambda }_{\mathrm{i}}}GP{P}_{\mathrm{R}\text{\_}residual}\left(x_{\mathrm{i}}\right)\begin{array}{ccc},& s.t.& \end{array}{\displaystyle \sum _{\mathrm{i}}^{\mathrm{N}}{\lambda }_{\mathrm{i}}}=1$$

(6)

where λ_i represents the weight for the residual of the ith 1 km pixel, which is centered at x_i. As shown in Equation (6), the ATPK phase takes the spatial correlation among pixels into account, which is not considered in the above equations.

The main issue of the ATPK phase is obtaining the weights, λ, which are estimated by minimizing prediction error spatial variations at a 1 km resolution. The kriging system can be expressed as:

$$\left[\begin{array}{cccc}{\mathrm{C}}_{\mathrm{RR}}\left(x_1,x_1\right)& \cdots & {\mathrm{C}}_{\mathrm{RR}}\left(x_1,x_N\right)& 1\\ \cdots & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots \\ {\mathrm{C}}_{\mathrm{RR}}\left(x_N,x_1\right)& \cdots & {\mathrm{C}}_{\mathrm{RR}}\left(x_N,x_N\right)& 1\\ 1& \cdots & 1& 0\end{array}\right]\left[\begin{array}{c}{\lambda }_1\\ \cdots \\ \cdots \\ \cdots \\ {\lambda }_N\\ \mu \end{array}\right]=\left[\begin{array}{c}{\mathrm{C}}_{rR}\left(x,x_1\right)\\ \cdots \\ \cdots \\ \cdots \\ {\mathrm{C}}_{rR}\left(x,x_N\right)\\ 1\end{array}\right]$$

(7)

where C_RR(x_i, x_j) is the residual covariance between a 1 km pixel centered at x_i and a 1 km pixel centered at x_j, and C_rR(x, x_j) is the residual covariance between a 500 m pixel centered at x and a 1 km pixel centered at x_j. Variable μ represents the range multiplier. To estimate weight λ according to Equation(7), these two types of residual covariance need to be calculated in advance.

Assuming each pixel as a point, s is intended to represent the Euclidean distance between two pixels. The 1 km to 1 km residual covariance, C_RR(s), and the 500 m to 1 km residual covariance, C_rR(s), can be estimated as:

$${\mathrm{C}}_{RR}\left(s\right)={\mathrm{C}}_{rr}\left(s\right)\ast h_R(s)\ast h_R(-s)$$

(8)

$${\mathrm{C}}_{rR}\left(s\right)={\mathrm{C}}_{rr}\left(s\right)\ast h_R(s)$$

(9)

where C_rr(s) represents the 500 m to 500 m residual covariance, and h_R(∗) represents the point spread function (PSF). Variable −s indicates that the distance from a 500 m pixel (A) to another 500 m pixel (B) within a 1 km pixel is opposite to the distance from a 500 m pixel (B) to a 500 m pixel (A) (i.e., s).

Supposing the value of the coarse pixel is the mean of the 500 m pixels within it, and F is the ratio between the resolutions of 1 km pixel and 500 m pixel, then, PSF can be described as:

$$h_R(x)=\{\begin{array}{cccc}& \frac{1}{F^2}& & x\in R(x)\\ & 0& & otherwise\end{array}$$

(10)

Given Equations (8)–(10), the calculation of C_RR(x_i, x_j) and C_rR(x, x_j) can be simplified as:

$${\mathrm{C}}_{RR}\left(x_i,x_j\right)=\frac{1}{F^4}{\displaystyle \sum _{m=1}^{F^2}{\displaystyle \sum _{m\prime =1}^{F^2}{C}_{rr}(s_{mm\prime })}}$$

(11)

$${\mathrm{C}}_{rR}\left(x,x_j\right)=\frac{1}{F^2}{\displaystyle \sum _{m=1}^{F^2}{C}_{rr}(s_m)}$$

(12)

where S_mm’ is the distance between each pair of 500 m pixels from within two 1 km pixels centered at x_i and at x_j, respectively, and s_m is the distance between 500 m pixels centered at x and any 500 m pixels within the 1 km pixel centered at x_j. Therefore, the critical issue when estimating the weight is the estimation of the 500 m to 500 m residual covariance, C_rr(s), which is essentially semivariogram modeling [32,34].

2.4. Result Validation and Method Evaluation

In this study, the performance of the proposed downscaling method was evaluated by MODIS GPP at a 500 m resolution. The coefficient of determination (R²) and the root mean square error (RMSE) between the downscaled results and the 500 m MODIS GPP were adopted as indicators of validation. In addition to the proposed downscaling approach, other three downscaling methods, including ordinary kriging (OK) [35], spline [36] and inverse distance weighting (IDW) [37], were used to provide a systematic comparison and illustrate the advantages of the proposed downscaling method.

3. Results and Analyses

3.1. Topographic and Vegetation Heterogeneities at Two Scales

Figure 3 shows the spatial distributions of topographic factors and the LAI on Julian day 193 at 500 m and 1 km resolutions. Although the values of these variables at two resolutions presented similarity on the whole distributions, the spatial variations at a 500 m resolution are more pronounced than those at a 1 km resolution.

The statistics of the topographic factors and the LAI on Julian day 193 at 500 m and 1 km are listed in

Table 1. Comparisons of topographic factors and at 500 m and 1 km.

Pixel Resolution (m)	Min		Max		Mean		STD
	500	1000	500	1000	500	1000	500	1000
Altitude (m)	2606	2647	4744	4671	3784	3787	501	489
Slope (deg)	0	0	38	28	18	12	8	6
Aspect (deg)	0	0	360	360	178	177	101	100
LAI (m2 m−2)	0	0	7.0	5.7	1.88	1.58	1.93	1.61

. The altitude at a 500 m resolution varied between a minimum of 2606 m and a maximum of 4744 m; at a 1 km resolution, the altitude varied between a minimum of 2647 m and a maximum of 4671 m. The mean altitudes at the two resolutions are similar, while the range and deviation in altitude at 500 m is larger than that at a 1 km resolution. The slope had a maximum of 38° at a 500 m resolution, which varied from the 28° maximum at a 1 km resolution. The mean slope value at 500 m (18°) was greater than that at a 1 km resolution (12°), showing that there was a considerable loss of topographic information at a coarse resolution. The statistical results also show that the LAI on Julian day 193 at 1 km had a smaller deviation than that at 500 m, and the mean value of the LAI at 500 m (7.0 m² m⁻²) was greater than that at a 1 km resolution (5.7 m² m⁻²). In addition, the aspect at 500 m had a similar range, mean value and deviation to those at 1 km.

3.2. GPP Difference of the Two MODIS Products

To assess the GPP inconsistency at 1 km and 500 m, a pixel-by-pixel comparison between the MODIS GPP at two resolutions was implemented with a time interval of 8 days during Julian days 169–209. Figure 4a presents the relationships between MODIS products at 500 m and 1 km during the study period. From Julian day 169 to Julian day 209, the 1 km MODIS GPP explained between 24% and 47% of MODIS GPP at 500 m, with RMSEs varying between 1.34 and 2.14 gC m⁻²d⁻¹. In general, the MODIS GPP at 1 km had a significant inconsistency with MODIS GPP at 500 m (R² = 0.36, RMSE = 1.62) during the whole study period, showing that the simplification of complex terrain within each modeling pixel might cause the model results to be considerably biased.

To investigate the effects of topography on GPP, a mismatched MODIS GPP between a 1 km resolution and a 500 m resolution was analyzed. Figure 5 illustrates the relationships between topographic factors (i.e., altitude, slope) and the inconsistency in the MODIS GPP at 500 m and 1 km during Julian days 169–209. Altitude and slope were divided by 500 m and 10° intervals, respectively, and the relative error (RE) was adopted to evaluate the inconsistency in the MODIS GPP at 500 m and 1 km. It was found that the inconsistency in MODIS GPP at two resolutions increased with the increase of slope and altitude. In addition, there was not an obvious relationship between the aspect and the inconsistency in MODIS GPP at two resolutions.

3.3. Downscaled Results

To evaluate the quality of the downscaled GPP, a pixel-by-pixel comparison between the MODIS GPP at 500 m and the downscaled GPP results was implemented during Julian days 169–209 (Figure 4b). From Julian day 169 to 209, the downscaled GPP for all pixels could be explained between 85% and 91% of 500 m MODIS GPP, with RMSEs varying between 0.71 and 1.30 gC m⁻²d⁻¹. The points in the scatterplots clustered around the 1:1 lines (R² = 0.89, RMSE = 1.03), and the downscaled GPP were slightly higher than MODIS GPP at high GPP values. In general, the downscaled GPP showed a strong consistency with the MODIS GPP at 500 m.

Figure 6 illustrates the absolute differences of GPP values between the 500 m MODIS GPP and the downscaled GPP results. During Julian days 169–209, the percentages of absolute differences in ranges of [0, 1], [1, 2], [2, 3], [3, 4], [4, 5] and [5, 6] were 70.78%, 23.02%, 3.47%, 1.61%, 0.89% and 0.23%, respectively. The pixels whose absolute GPP differences between the 500 m MODIS GPP and the downscaled GPP results were in the range of [0, 2], accounting for 93.80% of the total pixels in the study area, showing that the downscaling method proposed in this study could be effectively used to obtain the GPP at 500 m resolution. Figure 7 shows the density distributions concerning the GPP differences between the 500 m MODIS GPP and the downscaled GPP results. The mean values of differences between downscaled GPP and 500 m MODIS GPP varied between 0.60 and 1.03 gC m⁻²d⁻¹, with standard deviations changing between 0.47 and 1.10 gC m⁻²d⁻¹, which were smaller than those between MODIS GPP at 500 m and 1 km (mean = 1.36, std. = 1.34). It was worthwhile to note that, the proposed downscaling method could effectively solve the GPP underestimation of 1 km MODIS GPP.

To illustrate the advantages and benefits of the proposed downscaling method, the downscaled results in this work, along with the kriging downscaled results, spline downscaled results and IDW downscaled results, were compared against the MODIS GPP at 500 m (Figure 4). In general, the proposed downscaled results (R² = 0.89, RMSE = 1.03) had a stronger consistency with the 500 m MODIS GPP than the kriging downscaled results (R² = 0.43, RMSE = 1.36), the spline downscaled results (R² = 0.40, RMSE = 1.50) and the IDW downscaled results (R² = 0.42, RMSE = 1.10) for all Julian days. Hence, it can be concluded that the proposed downscaling method has a potential application for GPP downscaling to obtain 500 m GPP distributions, especially in mountainous areas. The results also illustrated that vegetation heterogeneity factors and topographic factors were effective for downscaling land surface fluxes, showing good matches with other related studies [7,8,38].

4. Discussions

4.1. The Topographic Effects on GPP Inconsistency at Two Resolutions

The time series GPP estimation is very important for monitoring vegetation growth conditions. Mountainous regions typically present high spatial heterogeneity, which is characteristic of steep slopes, altitude variations and notably biological diversity [11,39,40]. The results in this study showed the topography had a significant relationship with MODIS GPP inconsistency at two resolutions, as the topography influences the surrounding environment of the carbon assimilation procedure, such as air humidity, air temperature, water redistribution and solar radiation. Generally, air temperature increases linearly with the decrease of altitude. However, in most carbon cycle models, air temperature is assumed to be homogeneous within each modeling grid. This assumption neglects the impact of inner altitude variations in each grid, bringing uncertainties to the simulation over complex terrain. Besides air temperature, solar radiation, air humidity and precipitation also face the same problems. Related studies have illustrated that the spatial variations of altitude and slope greatly influenced hydrologic procedure [13,14,41], and hence affect the carbon assimilation procedure [42,43,44]. In general, the patchy ecosystem structures in mountain areas usually follow topographic and climatic gradients, and the remarkably spatial heterogeneity has considerable effects on the spatial distributions of GPP estimation.

4.2. Evaluations of the Proposed Downscaling Method

Serving as an effective tool in the estimation of GPP, the carbon cycle model is typically executed at coarse resolutions due to the paucity of the data and computing complexities. Landscape complexities are assumed to be simplified in a carbon cycle model by ignoring spatial heterogeneities within each modeling grid. It has been illustrated that the oversimplification of landscape complexities may inevitably cause model simulations to become considerably biased [8,9,10]. In this approach, massive forcing data were required to run the carbon cycle model, including the spatial distributions of temperature, radiation, vapor pressure deficits and so on. In mountainous areas, the regional estimation of the meteorological element was a challenge, especially at fine resolution [12,13,45].

In this study, a downscaling algorithm of GPP was developed to obtain the GPP distribution at 500 m resolution. The proposed downscaling method contains two phases: a regression process and an ATPK process. The first phase was used to reduce uncertainties caused by topographic effects and vegetation heterogeneities. A regression analysis was first conducted between the GPP estimation and spatial heterogeneity factors; then, the ATPK was performed to downscale residuals from the regression process.

This study tested the thesis that the 500 m GPP could be obtained from 1 km GPP and spatial heterogeneity information, which avoided complex calculation and the limit of forcing data through the carbon cycle models. It was found that the proposed method could be effectively used to obtain the GPP at 500 m resolution, showing significant advantages over the kriging method, spline method and IDW method. Our results also illustrated that both endogenous factors (i.e., LAI) and exogenous factors (i.e., altitude, slope and aspect) were useful for GPP downscaling. Exogenous factors have significant effects on the various instantaneous land surface fluxes due to the redistribution of water and their interactions with atmospheric elements [43,46,47]. The endogenous factors can be regarded as the accumulated outcomes of the surrounding environmental conditions [48], such as topography, climate and soil property. Therefore, both exogenous and endogenous factors should be considered in the downscaling scheme to remove biases within the coarse modeling grid.

4.3. Limitations of the Current Work and Prospects for Future Studies

The proposed downscaling method in this work fully utilized valuable topographic and vegetation information from ancillary data at a 500 m spatial resolution to reduce uncertainties in spatial heterogeneities over rugged terrain. However, due to practical constraints, there were also several limitations in this study. First, the proposed downscaling algorithm was based on the MODIS GPP, DEM and LAI data, and any errors in the source data would propagate into the final downscaled results [49,50]. Second, the resampling method for 30 m GDEM data influenced the quality of DEM data at 500 m and 1 km resolutions, and hence affected the accuracy of the topographic information [51,52]. In addition, the choice of algorithm for slope and aspect also resulted in uncertainty [53,54,55], besides slightly affecting the downscaling procedure. Third, although our algorithm considered the main factors (i.e., vegetation density, altitude, slope and aspect) for GPP downscaling, several other factors (e.g., land cover, soil texture, and soil depth) were not included in the algorithm, which might result in slight uncertainties. Finally, due to the lack of reference data for validation at smaller resolution, MODIS GPP at 1 km resolution were scaled only down to a 500 m resolution. This study tested the thesis that the 500 m GPP could be obtained from 1 km GPP and spatial heterogeneity information. Based on this thesis, the relevant work of obtaining the GPP at a finer resolution (e.g., 30 m) will be presented in an additional study.

To improve the performance of the proposed downscaling method, future studies should consider more factors in the downscaling algorithm, and attention should be given to the effects of accurate ancillary data (i.e., DEM and LAI data) on the final GPP downscaled results. Moreover, the proposed downscaling method should be applied in more heterogeneous areas to test the feasibility and stability of the algorithm.

5. Conclusions

In this work, a downscaling algorithm was developed to retrieve the GPP at 500 m resolution in a time series over rugged terrain, which considered the effects of spatial heterogeneity on carbon flux simulations. The GPP product at 1 km served as the primary dataset of the downscaling algorithm, and the 500 m GPP product was used as the reference dataset to evaluate the downscaled GPP results. It was found that the proposed downscaling method had advantages over the kriging method, spline method and IDW method. Our results also showed that vegetation heterogeneity factors (i.e., LAI) and topographic factors (i.e., altitude, slope and aspect) were useful for GPP downscaling.

Our study developed a reliable downscaling method to derive a 500 m resolution GPP for time series over mountainous areas, which is an important issue when evaluating the role of ecosystems in the global carbon cycle. Furthermore, our proposed downscaling algorithm could provide a reference when considering the effects of spatial heterogeneity on carbon flux simulations and retrieving other fine resolution parameters (e.g., net primary productivity and evaporation) in topographically complex terrains.