AbstractLand surface models intended for large‐scale applications are often executed at coarse resolutions, and the sub‐grid heterogeneity is usually ignored. Here, a spatial scaling algorithm that integrates the information of vegetation heterogeneity (land cover type and leaf area index) and surface topography (elevation, slope, relative azimuth (Raz) between the sun and the slope background, sky‐view factor, and topographic wetness index), was proposed to correct errors in gross primary productivity (GPP) estimates at a coarse spatial resolution. An eco‐hydrological model named BEPS‐TerrainLab was used to simulate GPP at 30 and 480 m resolutions for 16 mountainous watersheds selected globally. Results showed that an obvious improvement on GPP estimates at 480 m resolution was achieved after the correction in comparison with GPP modeled at 30 m resolution, with the determination coefficient increased by 0.38 and mean bias error reduced by 203gCm−2 yr−1. The combination of all the seven factors made the largest improvement for GPP estimation at 480 m resolution, suggesting that a larger improvement would be achieved when more factors of surface heterogeneity are considered. More specifically, our results indicated that five factors, including land cover type and leaf area index regarded as integrated outcomes of all the environmental conditions, Raz and sky‐view factor associated with radiation redistribution, and slope related to soil water redistribution, were especially important in the spatial scaling procedure. This study suggests that incorporating the information of surface heterogeneity into the spatial scaling algorithm is useful for improving coarse resolution GPP estimates over mountainous areas.