## 1 Introduction
For example, the environmental scalars of these GPP models were usually generated by the best-fit relationships between environmental indicators and the estimated GPP in a limited number of eddy flux sites [9; 25; 27; 32; 33]. They may not fully capture the actual relationship between a specific factor (e.g., moisture) and the GPP over a large area with substantial spatial heterogeneity [34]. In particular, the environmental stresses on GPP may vary considerably in a given biome under contrasting climate conditions [5]. In addition, a single environmental scalar cannot address the diverse responses of vegetation to environmental stresses, especially water stress, over large areas [3; 17; 33; 34]. Moreover, the satellite-based GPP models do not consider interrelated physiological and ecological processes, and thus tend to neglect the exacerbatedeffects of ecosystem responses during extreme climate events [35; 36]. A comparison analysis of these models is critical for their effective application in estimating ecosystem GPP over large areas. Previous efforts have focused mostly on the comparison of model performances among biomes [9; 17; 18; 37; 38]. The model performances were usually assumed to be the same among vegetation subtypes in a biome [38; 39]. In fact, the carbon fixed by photosynthesis has varied considerably within biomes, and may be larger than that among different biomes [40], suggesting that the model performances may vary substantially inside a biome. This is especially true for the vegetation subtypes limited by contrasting environmental stresses. Although some studies have evaluated satellite-based GPP models at a vegetation subtype level, they have primarily concentrated on individual vegetation subtypes for specific research objectives [4; 33; 36; 41; 42; 43; 33]. In addition, previous studies have typically calibrated and evaluated the models by in situ observations during normal precipitation years [3; 8; 25; 26; 39]. A recent study highlighted that the biases in GPP estimates from the satellite-based models are governed by drought effects [45], indicating that model performances might have been overvalued in previous studies, particularly for areas experiencing frequent drought events. A systematic comparison of model performances among vegetation subtypes in different precipitation years is urgently needed for a better understanding of the satellite-based models, and would contribute to a better application of these models in global carbon cycle research. Grasslands, occupying approximately one-fifth of the global land area [40], play a vital role in providing numerous ecosystem services such as regulating the global carbon cycle and supporting plant and animal biodiversity [46]. For example, grassland ecosystems store about 20% of total global carbon [47]. The status and distribution of grasslands therefore have gained strong interest from scientists and decision-makers and has been the subject of active research [46]. With the rapid advancement of space science and remote sensing techniques, satellite data have been widely used in grassland monitoring/observation, such as the grazing intensity, pasture quality and status, community classification, and carbon sequestration [48; 49; 50; 51; 52]. Among these, the GPP of grasslands has received considerable attention in recent decades due to its significance in understanding the ecosystem carbon cycle and the simplicity of satellite-based GPP models. Using grasslands in northern China as an example, this study evaluated the performances of eight widely used satellite-based GPP models (VPM, MVPM, MODIS-GPP, TG, GR, VI, AVM, and PCM) using in situ GPP estimates at four EC flux towers covering two typical Chinese grassland subtypes (temperate steppe and alpine meadow). We hypothesized that the model performances were largely determined by their powers to explain environmental limitations.