Decadal change and inter-annual variability of net primary productivity on the Tibetan Plateau

Climate is defined as long-term averages and variations in weather conditions over a period of usually several decades in a region. Climate variations are normally described in two aspects. One aspect describes the short-term (e.g., inter-annual) variability around the mean state and is commonly termed variability. Some of the variability occurs randomly (e.g., extreme events such as heat wave, deluge and drought) and some occurs periodically and relatively regularly (e.g., El Niño and Southern Oscillation or ENSO, Pacific Decadal Oscillation or PDO). Another aspect relates to a shift in the long-term (e.g., decadal) mean, often termed climate change and is represented by a long-term trend, e.g., global warming caused by humans. In this study, we define variation as the combination of long-term change and short-term variability.

Both climate variability and change affect terrestrial ecosystems (Piao et al. 2012; Cuo et al. 2016; De Boeck et al. 2019; Shukla et al. 2019). Extreme temperature and precipitation events can lead to legacy effects in plants years after the original occurrence, indicating the profound influence of climate extremes on terrestrial ecosystems (De Boeck et al. 2019). Browning trends or trends of decreasing photosynthetic activity are projected in regions with increasing drought and heat waves (Shukla et al. 2019). Climate change can lead to plant community range contraction or expansion, biodiversity reduction, and even extinction for those species that are unable to adjust to warming in the alpine region because of reduction in freezing tolerance and restriction in the species-specific phenotypic plasticity (De Boeck et al. 2019). Globally, there is a high confidence that greening trends have increased by 22–32% in the last 2–3 decades in China, India, parts of Europe, central North America, southeast Brazil and southeast Australia due to land use practices, forest conservation and expansive CO₂ fertilization, extended growing season, global warming, nitrogen deposition, and increase of diffuse radiation (Shukla et al. 2019). Among them, China and India stand out with the greatest greening trends mostly due to expanded cropland for both and expanded forest for China (Chen et al. 2019).

Photosynthesis couples the biosphere, atmosphere, hydrosphere and pedosphere together through the exchange of energy, water and geochemicals. Photosynthesis depends on the environmental conditions including air temperature, water, solar radiation and nutrients such as carbon and nitrogen. Primary productivity measures the photosynthesis rate at which a living organism converts light energy or shortwave solar radiation to chemical energy of an organic matter. Gross primary productivity (GPP) is the total light energy assimilated by plants. Net primary productivity (NPP) is GPP minus plant autotrophic respiration of energy used for plant tissue metabolism and is the energy stored in the plant tissue. NPP represents a key integrative process in ecosystem and is an important indicator of plant dynamics and the net carbon exchange between the terrestrial ecosystem and atmosphere. NPP is also linearly correlated with the carrying capacity of an ecosystem. Thus, understanding NPP’s response to climate variability and change will certainly improve our understanding of ecosystem response to the changes that are occurring in the global environment and will also benefit the local livestock and forestry management.

Alpine ecosystems are the communities that exist in the zone above treeline where average growing season temperature is around 6.7 °C or lower (Korner and Paulsen 2004). Alpine ecosystem biomass is low compared to shrubland, woodland and forest; however, it is often located in the headwaters of rivers that sustains life and it resists land surface erosion such as on the Tibetan Plateau (TP) and hence is important for hydrology and water resources. In the alpine region, harsh environment and remoteness often restrict in-situ observations and therefore limited research exists there. Dominant biome of the alpine ecosystems is grassland such as meadow or steppe. As grasslands cover about 30–47% of all terrestrial area (Shukla et al. 2019), they play a critical role in global biogeochemical cycles. To better account for the alpine grassland response to climate variability and change, an integrated approach involving modeling, field observations and remote sensing is necessary.

Relatively abundant studies have focused on the impacts of temperature and precipitation changes on greening and browning using satellite measured normalized difference vegetation index (NDVI) or enhanced vegetation index (EVI) (Sarmah et al. 2018; Du et al. 2019; Feng et al. 2019b; Liu et al. 2019; Qian et al. 2019; Sun et al. 2019; Wang et al. 2019; Li et al. 2020a; Parida et al. 2020; Yuan et al. 2020). Studies using other indicators of vegetation are relatively limited and some of them are mentioned here, e.g., LAI (Piao et al. 2015; Chen et al. 2019), foliar projective coverage (FPC, Cuo et al. 2016), tree ring (Shi et al. 2019), ecosystem indices (Fu et al. 2019), simulated NPP (Zhuang et al. 2010; Piao et al. 2012; Bao et al. 2019; Feng et al. 2019b), modeled carbon and water use efficiency (El Masri et al. 2019), and observed or surveyed phenological states and NPP (Niu et al. 2019; Li et al. 2020b). Few have studied the effects of cloud cover and CO₂ changes (Piao et al. 2012; Cuo et al. 2016) on vegetation. Cloud cover directly affects solar radiation, and it also relates to precipitation, all of which influence photosynthesis.

Observed aboveground grassland NPP at ungrazed sites in the Tianshan mountains in central Asia increased with warming and wetting during 1985–2016 (Li et al. 2020b). By using measured above- and below-ground biomass at alpine grassland sites on the TP during 2013–2015, Niu et al. (2019) found that about 76% of below-ground NPP concentrated in the top 10-cm soil layer and that above-ground (below-ground and total) NPP was mainly affected by temperature (precipitation). Using the Carnegie-Ames-Stanford Approach (CASA) model, Feng et al. (2019b) simulated NPP across China for 1982–2015 and found that spatially averaged annual and seasonal NPP increased due to temperature and solar radiation increase and that there was a shift in the trend around 1998 when global warming hiatus started. By using the ORCHIDEE model, Piao et al. (2012) showed that NPP also increased over 1961–2009 on the TP due primarily to changes in precipitation but secondarily to changes in temperature and CO₂ concentrations. Using the process-based Terrestrial Ecosystem Model with permafrost scheme, Zhuang et al. (2010) revealed that NPP exhibited larger Pearson correlation coefficients with precipitation than those of temperature, soil moisture and soil temperature on the TP during 1901–2002. The responses of various plant functional types (PFTs) to climate variations are rarely investigated on the TP except, e.g., Cuo et al. (2016) and Zhong et al. (2019). Cuo et al. (2016) showed that in the northern TP FPC of temperate needleleaf evergreen trees in the east and alpine meadow/steppe in the northeast and southwest mostly increased, while FPC of barren/sparse grassland in the northwest and perennial temperate summer green scrub/grassland in the east mostly decreased.

Worldwide, several studies examined the influence of climate variability on vegetation variability represented by NDVI, EVI and above-ground biomass (Le Houerou et al. 1988; Lotsch et al. 2003; Vicente-Serrano et al. 2013; Thornton et al. 2014; Seddon et al. 2016; Sloat et al. 2018; Bao et al. 2019; Chen et al. 2019; Feng et al. 2019a; Stanimirova et al. 2019), and revealed the dominance of precipitation variability in the arid and semi-arid ecosystems. For example, using the CASA model, Bao et al. (2019) showed that annual NPP was about 265 g C/m² on the Mongolian Plateau and that summer drought played a significant role in determining annual NPP variability. De Boeck et al. (2019) noted that coincident summer time drought and heat waves exerted strong pressure on plant cover and community composition for a long time, indicating the profound effects of extreme climate variability in cold biomes. However, similar studies on the TP are very limited. Only Feng et al. (2019a) found that precipitation could account for the NDVI variability in the permafrost zone of the TP in some periods during 1982–2012.

Being the first of its kind on the TP, this study aims to deconvolve the impacts of annual variability and long-term change in major climate elements on annual NPP by using a dynamic vegetation model with carefully designed scenarios (Table 1). This study will fill the knowledge gap in our understanding of both climate change and variability impacts on NPP.

2.1 Study area

The TP (73–105° E, 26–40° N) is located in the southwest of China. It has an area of 2.5 × 10⁶ km² with elevation ranging in 2000–8800 m, averaging about 4400 m above sea level. The TP is the headwaters of major rivers in Asia including the Dang River, Hei River, Yellow River, Yangtze River, Mekong River, Salween River, Brahmaputra River, Ganges River, Indus River, and Tarim River (Fig. 1a). Altogether, these rivers support about 1.65 billion people (Cuo and Zhang 2017). The TP is influenced primarily by Asia summer monsoons in May–September where 70% of annual precipitation falls and by westerlies in the dry winter season of October–April. Mean annual air temperature was about 1.2 °C and annual precipitation was about 410 mm during 1970–2015, both decreasing from the southeast to the northwest. According to a survey conducted by the Chinese Academy of Sciences in 1981–2001 (CAS 2001), land use/cover of the TP (Fig. 1b) includes perennial alpine meadow (26.9%), perennial alpine steppe (34.8%), Barren/sparse grassland (15.8%), perennial temperate scrub grassland (14.3%), and subtropical needleleaf evergreen trees (4.2%), together occupying about 96% of the TP. From the southeast to the northwest, the ecosystem degrades from forest to desert.

a Geographic location, terrain elevation in 1-km resolution and main river basins and b land cover land use distributions on the Tibetan Plateau. In a, black line outlines the Tibetan Plateau inside China. A The Dang River, B the Hei River, C the Yellow River, D the Yangtze River, E the Mekong River, F the Salween River, G the Brahmaputra River, H the Ganges River, I the Indus River, J the Tarim River. TAR Tibetan Autonomous Region, QH Qinghai Province. Blue lines denote the provincial level political boundaries. Black triangles are NPP observation sites from Zhou et al. (2016). In b Temp temperate, e.g.: evergreen, sg summer green, Subtro subtropical

Full size image

2.2 Model

The Lund–Potsdam–Jena Dynamic Global Vegetation Model (LPJ, Sitch et al. 2003; Gerten et al. 2004) was used in this study to investigate NPP change and variability caused by climate factors. LPJ accounts for vegetation dynamics including establishment, succession and mortality by simulating the process of competition for light, water and nutrients fluxes among PFTs. Establishment is determined by bioclimatic conditions such as lower and upper temperature for CO₂ absorption and photosynthesis and coldest monthly temperature that PFTs can withstand. LPJ simulates photosynthesis, transpiration, soil organic matter and litter dynamics at daily time step. Vegetation production, tissue turnover and population dynamics are simulated at monthly and annual time steps. LPJ has been applied to investigate ecosystem conditions around the world (Sitch et al. 2003, 2005, 2008; Gerten et al. 2004; Steinkamp and Hickler 2015). In this study, LPJ was configured at 0.25° by 0.25° spatial resolution and was driven by monthly gridded air temperature, precipitation and cloud cover at the same resolution and global annual CO₂ data. The model was warmed up by running iteratively for 1000 years during 1957–1986 assuming that a 1000-year spinup could result in a vegetation dynamic equilibrium in the region, which is also a common practice in the LPJ modeling community (Sitch et al. 2003; Cuo et al. 2016, 2020).

LPJ was calibrated in our previous studies and the calibrated parameters are shown in Table 2. The calibration was done by examining the impact of each parameter one at a time and then chose the top parameters with the largest impact (Cuo et al. 2011). During calibration, LPJ simulated PFTs for the northern and entire TP were compared to surveyed land cover types of CAS (2001) and Zheng et al. (2008) in Cuo et al. (2016, 2020), respectively. Cuo et al. (2016) found that LPJ simulated alpine grassland matched well to CAS (2001) and Zheng et al. (2008) in the northern TP. Over the entire TP, the similarity in the dominant PFTs between LPJ and CAS (2001) was 70% (Cuo et al. 2020). Cuo et al. (2016, 2020) also evaluated LPJ simulated soil moisture and soil temperature using observations on the TP. Root mean square error (RMSE) for soil moisture ranged from 0.02 to 0.14 m³/m³ at nine observation sites that spread on the TP. RMSE for soil temperature in top 40 cm ranged from 0.36 to 4.51 °C over the northern TP. This study further evaluates LPJ simulated NPP.

2.3 Data

Annual CO₂ and monthly gridded precipitation (P), 2-m air temperature (T) and cloud cover (CC) in 1957–2015 were used as model forcing. Air temperature and precipitation were gridded based on station measurements on the TP and the surrounding areas using procedures describe by Cuo et al. (2013). Cloud cover data were from the Climate Research Unit of the University of East Anglia (Mitchell and Jones 2005) at 0.5° × 0.5° spatial resolution and were regridded to 0.25° × 0.25° resolution assuming uniform distribution within each 0.5° × 0.5° grid cell. Annual CO₂ was from NOAA’s Mauna Loa Observatory. Soil textures were from the Harmonized World Soil Data Base v1.2 produced by the Food and Agricultural Organization of the United Nations and were regridded to 0.25° × 0.25° resolution. Three main soil types on the TP are loamy sand (4%), sandy loam (51%), and loam (42%). Elevation data were from the Shuttle Radar Topography Mission (SRTM) and were interpolated to the 0.25° × 0.25° grids using cubic convolution. Land cover data for 1981–2001 were survey based and were compiled by the Chinese Academy of Sciences (CAS 2001). The Moderate-Resolution Imaging Spectroradiometer (MODIS) NPP product MOD173AH and field derived NPP were used to evaluate the LPJ simulated annual NPP. MOD173AH was regridded to the 0.25° × 0.25° grids for a fair comparison. Field based grassland NPP was obtained from measured fresh grass weight at 11 sites on the TP during May–August of 2003–2013 (Zhou et al. 2016; see Fig. 1a; Table 3 for the sites). The fresh grass weight was measured once every month during the period and the maximum weight during all months was converted to dry grass weight by using a ratio of 0.4–0.6 determined from field experiments. The dry grass weight was further converted to NPP. In order to compare with field derived NPP, the simulated maximum NPP in May–August was chosen for 2003–2013. As the majority of vegetation type on the TP is grassland, the comparison of simulated and observed grass NPP is deemed representative for the TP.

2.4 Analysis

The study period was 1970–2015. The LPJ model has been calibrated and evaluated in terms of PFTs and soil moisture and temperature for the TP in previous studies (Cuo et al. 2016, 2020). In this study, the LPJ simulated monthly NPP was aggregated into annual and growing season NPP for 1970–2015 and was evaluated using MOD173AH and the field data.

To deconvolve the impacts of climate change and variability on annual NPP, 13 scenarios that fell into three categories were designed: variability only (Va), trend only (Tr), and trend and variability combined (Va + Tr) (Table 1). The goal for constructing these scenarios in Table 1 was to represent various aspects of climate variations: (1) short-term variability (Va), (2) long-term changes (Tr), and (3) combined variability with changes (Va + Tr). Within each of the three categories, there was a reference scenario in which consistency in variability, trends, as well as variability plus trends in all four climate factors including air temperature, precipitation, cloud cover and CO₂, were retained. Thus, through comparing each climate factor scenario with the corresponding reference scenario, we would be able to identify to what degree the climate factor contributes to the referenced state. Clearly, the rationale behind this is that the closer a climate factor resolves the referenced state, the more important the climate factor becomes.

For example, in order to analyze the impacts of variability alone, the reference scenario (Reference Va in Table 1) was driven by detrended temperature, precipitation and cloud cover at monthly step and detrended CO₂ at annual step but with their inter-annual variability retained (although CO₂ shows little inter-annual variability); the T Variability scenario was driven by mean precipitation, cloud cover and CO₂ but detrended temperature with its inter-annual variability retained; and similarly for the P Variability and CC Variability scenarios (Table 1). Here, the means of temperature/precipitation/CC/CO₂ were calculated for each month and for the year (CO₂ only) over 1970–2015. This way the seasonal cycle was preserved but the inter-annual variability and long-term changes were removed. Thus, when we used the means of all four climate factors except the one factor whose inter-annual variability or long-term changes were preserved, we could attribute the variations in NPP to the inter-annual variability or long-term changes of that one factor. As a demonstration, Fig. 2 shows the means, variability only, trend only and historical time series of annual temperature, precipitation and cloud cover during 1957–2015 averaged over the TP.

The mean, variability only, trend only and historical time series of annual temperature, precipitation and cloud cover during 1957–2015 averaged over the TP

Full size image

Specifically, taking the T Variability scenario as an example, the trend in each month’s temperature at each grid cell was removed but the inter-annual variability of temperature in each month was preserved, while the mean monthly values of precipitation and cloud cover and mean annual CO₂ were used at the same time. Thus, the comparison in NPP between the Reference Va and T Variability scenarios could reveal the effects of temperature variability only because these two scenarios were driven by the same temperature. Since CO₂ shows a significant annual trend but little annual variability (Fig. 3j), only its trend impacts on NPP were examined in this study.

Means (a, d, g), trends (b, e, h) and standard deviations (c, f, i) of annual air temperature (T, the 1st raw), precipitation (P, the 2nd raw), cloud cover (CC, the 3rd raw) and CO₂ (j) during 1970–2015. Stapled cells represent statistically significant trends at p < 0.05. Annual CO₂ trend is 1.63 ppm/year (p < 0.05)

Full size image

To examine the impacts of climate trend, five scenarios, i.e., Reference Tr, T trend, P trend, CC trend and CO₂ trend, were designed (Table 1). Reference Tr was the scenario where all forcings’ variability was removed but their long-term trends retained. This was done by first obtaining the long-term means in each month for temperature, precipitation and cloud cover and then adding the trends in each month (regardless of significance) back to the monthly mean time series of the corresponding climate element. This process effectively removes annual variability but preserves trend. As an example, in order to examine temperature trend impacts alone (T trend in Table 1), the scenario was driven by the monthly mean time series of all climate elements but with the monthly trends added for temperature only. Then comparing the simulated NPP from the Reference Tr and T trend scenarios could reveal the impacts of temperature trend. This was done similarly for precipitation, cloud cover and CO₂ although in the case of CO₂ annual time series were used.

Next, four scenarios, i.e., historical, T trend and variability, P trend and variability, and CC trend and variability (Table 1, no CO₂ scenario due to its little variability) were designed to investigate the impacts of historical climate variations that include both interannual variability and long-term changes. The historical or reference scenario was driven by historical temperature, precipitation, cloud cover and CO₂ that include both long-term trends and short-term variability. As an example, to investigate temperature variation (T trend and variability) impacts, historical temperature together with the monthly mean time series of precipitation and cloud cover and mean annual CO₂ were used to drive the model. By contrasting the simulated annual NPP from the historical and T trend and variability scenarios, the impacts of temperature variation on annual NPP could be revealed. This was done similarly for precipitation and cloud cover.

For temperature, precipitation and cloud cover, the variation at monthly time step was used to enable seasonal vegetation growth. For CO₂, the input time step in the model is annual and its mean and trends were all at annual step. The Mann-Kendall analysis and Sen’s slope (Sen 1968) were employed to test and obtain the annual and monthly trends of the simulated NPP, temperature, precipitation, cloud cover and CO₂ for each cell on the TP in 1970–2015. Conventional statistical matrixes including mean, correlation coefficient (R), coefficient of determination (R²), coefficient of variation (CV), standard deviation (Std) and root mean square error (RMSE), were also examined. For fair comparisons, the standardized anomalies of NPP of all scenarios were used.

The analysis of model uncertainty including forcing, parameter and model structure uncertainties could enhance the reliability of the study. Cuo et al. (2020) evaluated model forcings’ uncertainty by comparing them to station observations and NOAA CMORPH (Climate Prediction Center Morphing Technique) precipitation and multimodel GLDAS (Global Land Data Assimilation System) products, and showed that the forcings exhibited fairly good quality given the currently available station distributions on the TP. LPJ parameter uncertainty was examined in this study by changing the values of the sensitive vegetation and soil parameters by ± 10% relative to the calibrated values for all the scenarios. The changed vegetation parameters were gmin, GDDs, GDD5min, Int and Ws-m (see Table 2). The other vegetation parameters in Table 2 were not changed to avoid alteration in the climate zones and disruptions in vegetation spatial patterns. The changed soil parameters included percolation, water holding capacity, heat diffusivity, soil bubbling pressure, exponential decrease rate of soil moisture at freezing temperature, soil bulk density, soil density, and quartz content. The simulations with ± 10% changes in vegetation and soil parameter values were compared with the simulations with the calibrated vegetation and soil parameters. The parameter change induced uncertainty in annual NPP variability and change was also examined. However, LPJ model structure uncertainty has not been examined to date and this could be conducted in the future by utilizing multiple dynamic vegetation models.

3.1 Climate on the TP

Spatially averaged temperature, precipitation and cloud cover are about 1.2°C, 410 mm, and 53%, respectively, over 1970–2015. Spatially averaged air temperature over the TP increases at 0.044 °C/year (p < 0.05), while precipitation increases at 0.50 mm/year (p < 0.1) and cloud cover decreases insignificantly at 0.02%/year during 1970–2015. Spatial patterns of long-term annual means, trends and variability of air temperature, precipitation, cloud cover and CO₂ on the TP over 1970–2015 are shown in Fig. 3. Annual mean air temperature ranges from − 4 to 18 °C, low in the center and high in the southeast (Fig. 3a). Air temperature increases significantly at most cells on the TP, with a few cells even showing increase rates in excess of 0.10 °C/year (Fig. 3b). Annual temperature has high standard deviation in the central TP (Fig. 3c). Annual precipitation varies from less than 100 mm in the west to 1200 mm in the southeast (Fig. 3d). Worthy of mentioning here is the isolated 600-mm area along the Himalayas probably the artifact of the extrapolation using limited nearby stations. This can be fixed only when more stations become available. Annual precipitation displays increase trends over the most part of the TP, with significantly increasing trends occurring mainly in the center and the northeast (Fig. 3e). Significantly decreasing trends in annual precipitation are found only in isolated spots in the east and southeast. The standard deviation of annual precipitation decreases from the southeast to the northwest, following the spatial pattern of mean annual precipitation (Fig. 3f).

Annual mean cloud cover ranges from 40 to 77%, high in the east and west and low in the north and southwest where desertification occurs (Fig. 3g, Cuo et al. 2020). Cloud cover shows decrease trends over the most part of the TP, with significantly decreasing trends located in the central to southeast and a few significantly increasing trends in the south and west (Fig. 3h). Annual cloud cover has high standard deviation in the northeast and southeast (Fig. 3i). Annual mean CO₂ is 359 ppm during the period and annual CO₂ increases steadily and significantly at 1.63 ppm/year (Fig. 3j). Annual CO₂ exhibits little interannual variability and this is the reason that the effects of CO₂ variability on annual NPP is not examined.

3.2 Evaluation of simulated NPP

LPJ simulated NPP displays high values on the order of 800 g/cm²/year in the southeast where tropical and subtropical forests exist and low values of 0–100 g/cm²/year in the northwest, similar to what MODIS NPP shows (Fig. 4a, b). However, LPJ simulated NPP is higher over the central TAR (Tibetan Autonomous Region) but much lower along the eastern periphery of QH (Qinghai Province) compared to MODIS NPP, the reason of which may be due to higher precipitation over the central TAR and lower precipitation in the eastern QH in the forcing as indicated in Fig. 3d. LPJ simulated NPP also exhibits an isolated area of high value along the Himalayas that is not supported by MODIS NPP due to overestimated precipitation (Fig. 3d). The regionally averaged NPP during 2000–2014, 120 g/cm²/year in the model vs. 126 g/cm²/year in MODIS product, is comparable. Spatial pattern correlation between Fig. 4a and b is 0.58 (p < 0.05). Scatter plot between field derived and LPJ simulated NPP reveals a generally linear correspondence with a statistically significant coefficient of determination of 0.41 (Fig. 4c), but there is clear underestimation in LPJ simulated NPP which is perhaps due to the coarse resolution of the LPJ model compared to the point measurements. Our previous studies on the TP showed a satisfactory performance of LPJ in simulating PFTs for the region and soil temperature and moisture at multiple sites across the TP (Cuo et al. 2016, 2020). Thus, despite some deficiencies that are likely related to the forcings due to limited observations, LPJ simulations are generally reasonable on the TP. The model simulations are subsequently used for climate variation analyses.

Spatial distributions of the LPJ simulated (a) and MODIS derived (b) average annual NPP during 2000–2014; scatter plot of LPJ simulated annual NPP vs. field derived annual growing season NPP from Zhou et al. (2016) for 2003–2013 (c). *Statistical significance at p < 0.05

Full size image

3.3 Correlation of NPP with forcings at monthly time step

Concurrent and lagged correlations between LPJ simulated NPP from the historical scenario and observation-based climate factors are calculated. Here, monthly temperature/precipitation/cloud cover were paired with 0-, 1-, 2-, 3-, and 4-month lagged NPP and the correlations were obtained for the entire TP over 1970–2015. In general, only the concurrent and 1-month lagged correlations exhibit significantly large R² (Fig. 5), ranging from 0.52 to 0.77 for the former (largest for temperature) and 0.72–0.74 for the latter (largest for cloud cover). As the lagged time progresses from one month to four months, the R² values become gradually smaller and smaller, e.g., 0.35 (0.54) for 2-month lag for temperature and precipitation (cloud cover).

Scatter plots between the anomalies of gridded monthly air temperature (precipitation, cloud cover) and LPJ simulated concurrent NPP (a, c, e) and 1-month lagged NPP (b, d, f) from the historical scenario during 1970–2015

Full size image

3.4 Annual NPP variability

Regionally averaged time series of the standardized anomalies of annual NPP simulated by the four variability scenarios (Table 1, with the calibrated values and ± 10% changed values in vegetation and soil parameters) are shown in Fig. 6a–c. NPPs using the calibrated values and ± 10% changed values completely overlay on top of each other, indicating low uncertainty in the simulated NPP due to uncertainties in vegetation and soil parameters for all scenarios. The standard deviations of the normalized regionally averaged annual NPP for Reference Va, T/P/CC variability are the same, i.e., 1.14 g/cm²/year. RMSE of simulated annual NPP between the Reference Va scenario and the T/P/CC variability scenarios, reflecting the closeness between them, are 0.53, 0.99 and 1.56 g/cm²/year, respectively. T variability simulated NPP most closely follows the reference NPP throughout the period, with P variability coming in second and CC variability showing the least match.

Regionally averaged standardized anomalies of annual NPP simulated from the reference Va, T variability, P variability and CC variability (a–c); from reference Tr, T trend, P trend, CC trend and CO₂ trend (d–g); and from the historical, T trend and variability, P trend and variability, CC trend and variability (h–j). Simulated NPP using the calibrated values and ± 10% changes in vegetation and soil parameter values are shown for all scenarios. Thick black lines represent NPP from reference simulations and are the same in each column. Gray lines represent NPP from the designed scenarios and are different between the panels. The first (second, third and forth) raw shows the impacts of temperature (precipitation, cloud cover and CO₂, respectively). Designed scenarios are listed in Table 1

Full size image

Spatial distributions of coefficients of variation (CV) of annual NPP from the four variability scenarios with the calibrated values and ± 10% changed values in vegetation and soil parameters over the TP are shown in Fig. 7. For Reference Va, high CV is found in the west, central-eastern and northern TP, and low CV in the northwest, north, southeast and east (Fig. 7a–c), indicating that (1) low (high) NPP tends to have high (low) variability (see Fig. 4) and (2) modeled NPP shows low uncertainty with changes in vegetation and soil parameters. Over the areas with little vegetation growth, e.g., the northwestern TP where the Chaidamu Desert, Qilian Mountain range and East Kunlun Mountain range are located, NPP displays zero variability.

Coefficients of variation (CV) of annual NPP using the calibrated values and ± 10% changes in vegetation and soil parameter values from reference Va (a–c), T variability (d–f), P variability (g–i), and CC variability (j–l). The first column is from calibrated vegetation and soil parameters, the second (third) column is from − 10% (+ 10%) perturbation in vegetation and soil parameter values. Scenario definitions are listed in Table 1

Full size image

To examine the impacts of temperature, precipitation and cloud cover variability on annual NPP variability spatially, the spatial distributions of CV from the T/P/CC variability scenarios are shown in Fig. 7d–l. Compared with Reference Va, it is apparent that T variability produces the closest match which is also consistent with what is revealed in Fig. 6a–c. P variability and CC variability exhibit very small CV of NPP across the region and with little spatial variation. However, P (CC) variability produces an area of high CV in the northern sparsely vegetated region (along the west edge of the TP) similar to what is shown in Fig. 7a–c, indicating strong influence of P (CC) variability there. R² of the spatial patterns of annual NPP between Reference Va and T/P/CC variability are 0.74, 0.0 and 0.04, respectively. This analysis indicates that the temporal and spatial annual NPP variability is dictated primarily by temperature variability.

3.5 Annual NPP long-term change

Averaged over the TP, the standardized anomalies of annual NPP simulated by the trend scenario Reference Tr increases significantly at 0.086 g/cm²/year² (p < 0.01, black lines in Fig. 6d–g). Comparisons between Reference Tr and the T/P/CC/CO₂ trend scenarios clearly show a match in annual NPP between the CO₂ trend scenario and Reference Tr, indicating that the influence from CO₂ increases dominates the annual NPP trend across the region. Annual NPP from the T/CC trend scenarios also follow the reference annual NPP although for T/CC trend there are small but clear variations due to seasonal fluctuations superposed on the increasing trends. For the P trend scenario, annual NPP exhibits even larger variations due to seasonal fluctuations obscuring long-term trend seen from other scenarios.

For the spatial patterns of long-term change, annual NPP simulated from Reference Tr with the calibrated values and ± 10% changed values is identical and all exhibits predominantly increase trends across the TP (76% positive trends vs. 8% negative trends) (Fig. 8a–c). The T/P/CC trend scenarios all display mixed trends spatially, and among them, the T trend scenario shows the most negative trends (47%) that spread across the region. The P and CC trend scenarios correspond respectively to 23% and 25% negative trends that are mostly distributed in the south. On the other hand, CO₂ increases result in predominantly positive trends totaling 83% of the entire region and the spatial patterns between the CO₂ trend scenario and Reference Tr are quite similar. The correlation coefficient (R) between spatial pattern of reference Tr and T/P/CC/CO₂ were 0.50, 0.53, 0.26 and 0.60, respectively. Clearly, CO₂ is the most important factor in long-term annual NPP change in the region, followed by precipitation and temperature.

Simulated annual NPP trends from reference Tr (a), T trend (d), P trend (g), CC trend (j), and CO₂ trend (m). The first raw is from reference simulation, the second (third, forth and fifth) raw is from T trend (P trend, CC trend and CO₂ trend, respectively). The uncertainty caused by − 10% (+ 10%) perturbations in vegetation and soil parameter values for each scenario are shown in the second (third) columns. Stapled cells denote statistically significant trends at p < 0.05. Scenario definitions are listed in Table 1

Full size image

3.6 Annual NPP variation

Spatially averaged standardized anomalies of annual NPP from the Historical and T/P/CC trend and variability scenarios (i.e., historical forcings with both variability and long-term change retained) are shown in Fig. 6h–j. In general, annual NPP from the T/P/CC trend and variability scenarios largely follows that from the Historical scenario. However, differences are noted in correlation and RMSE in annual NPP between the Historical and T/P/CC trend and variability scenarios in that correlation coefficient ranges from 0.67, 0.57 to 0.25, and RMSE from 0.93, 1.06 to 1.40 for temperature, precipitation and cloud cover, respectively. Hence, in terms of the combined trend and variability of the regionally averaged annual NPP, temperature influence is stronger than precipitation while cloud cover exerts the least impact among these three climate elements.

The spatial patterns of annual mean NPP, annual NPP trend and annual NPP CV during 1970–2015 from the Historical scenario are shown in Fig. 9. Annual mean NPP is in the range of 0–1000 g/cm²/year and decreases from the southeast to the northwest on the TP (Fig. 9a). The long-term trends of annual NPP vary spatially, with statistically significant decrease trends located over the arid west and northwest and mixed trends over the eastern TP (Fig. 9b). Annual NPP CV is large in the west where vegetation is sparse and low in the east where vegetation is relatively dense (Fig. 9c). Among the dominant PFTs on the TP (Table 4), subtropical needleleaf evergreen trees (barren/sparse grassland) shows the highest (lowest) mean annual NPP while perennial alpine steppe (subtropical needleleaf evergreen trees) displays the greatest (least) CV of 0.46 (0.24). All dominant PFTs exhibit increase trends in annual NPP, with barren/sparse grassland and subtropical needleleaf evergreen trees having significant trends (p < 0.05).

Annual mean NPP (a), NPP trend (b), and coefficient of variation (c) simulated by historical temperature, precipitation, cloud cover and CO₂ during 1970–2015. Staples in b represent statistically significant trends at p < 0.05

Full size image

We next use correlation coefficients (Rs) of annual NPP between the Historical and T/P/CC trend and variability scenarios to examine the influence of temperature, precipitation and cloud cover on the trend and variability of NPP across the TP together with vegetation and soil parameters uncertainty (Fig. 10). Same as previous analysis, Fig. 10 reveals low parameter uncertainty in the model. CO₂ is not examined here due to little annual variability (Fig. 3j). The results reveal that the entire region (99% area) is dominated by positive R in excess of 0.6 for temperature scenario and the reference (Fig. 10a–c), indicating the dominant role of temperature. Precipitation (Fig. 10d–f) and cloud cover (Fig. 10g–i) present mixed and spatially varying R patterns, with 77% (67%) of the region showing positive R and statistically significant positive R located mainly in the south and central TP (central and western TP) for precipitation (cloud cover).

Correlation coefficients of simulated NPP anomalies between the historical and T trend and variability (a, b, c), P trend and variability (d, e, f), and CC trend and variability (g, h, i). The uncertainties caused by − 10% (+ 10%) perturbations in vegetation and soil parameter values for each scenario are shown in in the second (third) column. Gray area represents no R values due to zero NPP where desert or high elevation exists throughout the study period

Full size image

The following equation is often used for studying attribution (Nayak et al. 2013; Feng et al. 2019a, b). This equation implies that the examined system is linear. Here we applied this equation to examine long term change and variations of NPP using LPJ simulations and gridded temperature, precipitation and cloud cover and observed CO₂, and compared with our modeling results in which system nonlinearity is preserved.

where X_i is the ith climate factor, n total number of climate factors, R residual, \(\frac{\partial NPP}{\partial {X}_{i}}\) the sensitivity term, and \(\frac{{\partial X}_{i}}{\partial t}\) the trend of individual climate factor.

Figure 11 shows long-term change in NPP with variability removed derived from Eq. (1). The sensitivity of NPP changes to temperature, precipitation and cloud cover, and the trends of temperature, precipitation and cloud cover in Eq. (1) are obtained from the least-square fitting lines. NPP long-term trends (Fig. 11a) from Tr reference simulation are calculated from Sen’s slope. Temperature change induced NPP change (Fig. 11b) is the product of NPP from the Reference Tr scenario (Table 1) and gridded temperature trend. The same procedure applies for precipitation, cloud cover and CO₂ induced NPP change shown in Fig. 11c–e. Figure 11, similar to Fig. 8, shows the dominant role played by CO₂ change in long term NPP change as revealed over the majority of the TP by (1) consistent positive NPP change between historical and CO₂ induced NPP change (Figs. 8a, m, 11a, e) and (2) more than 80% CO₂ induced NPP change contribution to historical NPP change (Fig. 11i). Temperature also plays a secondary role. There are clearly differences for precipitation and cloud cover induced NPP change between Figs. 8 and 11, perhaps due to the linearity assumption implicit in Eq. (1).

Long-term change in NPP simulated from Tr reference scenario (a), and changes in temperature (b), precipitation (c), cloud cover (d) and CO₂ (e) induced NPP changes, and contributions from changes in temperature (f), precipitation (g), cloud cover (h) and CO₂ (i) to historical NPP change computed using Eq. (1). Here, contribution is computed as the ratio of climate element induced NPP changes to Tr trend scenario simulated NPP changes

Full size image

Figure 12 shows historical NPP variations based on Eq. (1) and Sen’s slope. Temperature induced NPP variation is the product of the sensitivity of NPP from Historical scenario (see Table 1) to temperature and the linear trend of temperature. The same procedure is applied to precipitation and cloud cover. Clearly, temperature induced variation in NPP displays rather similar spatial patterns to historical NPP variation calculated from Sen’s slope (Fig. 12a, b), and accounts for more than 40% of historical NPP variation across the TP, larger than precipitation and cloud cover contributions in general (Fig. 12e–g). Hence, this analysis corroborates our finding shown in Fig. 10a in that temperature is the major factor for NPP variation over the TP. The spatial patterns of contributions from precipitation and cloud cover (Fig. 12f, g) are not consistent with those of correlation coefficients shown in Fig. 10d, g, likely related to the shortcoming of linearity assumption in Eq. (1). We noticed that the residual term in Eq. (1) can sometimes be several times greater than actual NPP variation.

Historical trend and variability of NPP (a), trend and variability in temperature (b), precipitation (c), cloud cover (d) induced NPP trend and variability. Contributions from temperature (e), precipitation (f) and cloud cover (g) to NPP trend and variability computed using Eq. (1). Here, contribution is computed as the ratio of climate element induced NPP trend and variability to historical NPP trend and variability

Full size image

This study revealed the effects of inter-annual variability and long-term change in temperature, precipitation, cloud cover and CO₂ (long-term change only) on annual NPP variability and trend on the TP. Temperature variability exerts dominant influence on annual NPP variability both temporally and spatially, while precipitation and cloud cover display much smaller impacts in general and their influences are limited to certain areas. In the past six decades, the TP has experienced increase trends in extreme temperatures (Ding et al. 2018). Also, in the past 50 years, precipitation changes showed heterogeneous spatial patterns, with more frequent extreme precipitation occurring in the northeast but fewer precipitation days in the northwest of the northern TP (Cuo et al. 2013). As global warming continuous, it is likely that heat waves, drought and deluge will occur more frequently around the world (IPCC 2013). The extremes in temperature and precipitation variability are likely to occur more frequently on the TP as well. Given the strong influence of temperature variability on NPP variability across the region in general, and specifically in the Pamir Plateau, central and northern TP, and strong influence of precipitation variability in northern TP, frequent occurrences of extreme temperature and precipitation could pose a great danger to the local ecosystems. As was shown in the Mongolian plateau, more drought occurrences during 1982–2011 had detrimental effects on annual NPP (Bao et al. 2019).

Some studies that examined the regional NDVI/EVI/ANPP (above-ground NPP) trends on the TP showed the dominance of precipitation change effects on NPP up to 2009 (e.g., Zhuang et al. 2010; Piao et al. 2012); however, this is not the case during 1970–2015 as shown in this study. It is possible that the discrepancies in the findings may be due to the different trends for different time periods and also the different approaches used. On average, the TP experienced a precipitation trend of 0.61 mm/year and a warming trend of 0.041 °C/year during 1970–2009, whereas the precipitation trend was 0.50 mm/year and the warming trend was 0.044 °C/year during 1970–2015. Furthermore, Piao et al. (2012) calculated the climate and CO₂ change impacts on NPP by comparing individual simulations forced by detrended climate elements (as targeted scenarios) and historical climate (as reference scenario). Here we used variability removed climate forcing to drive the LPJ model and found that CO₂ is the dominant factor in long-term annual NPP change. Temperature and cloud cover display secondary effects on the regionally averaged NPP long-term change which is in line with the study from Feng et al. (2019b) for the entire China. Precipitation due to its great seasonal variability and less significant trends exerts small effects on the long-term NPP change. The modeling approach applied in this study effectively removed the contamination of the effects of annual variability in the forcing.

The differences in NPP response to temperature/precipitation between the TP and other regions are likely due to the differences in environmental conditions of the ecosystems which in turn lead to differences in energy/water limitation factors. Some researchers use Budyko’s dryness index to assess whether the environment is water or energy limited, e.g., Donohue et al. (2009). The Budyko’s theory states that when mean precipitation is greater than mean potential evapotranspiration during a given period, the environment is energy limited; otherwise, it is water limited. Zhang et al. (2009) showed that in the eastern and southern TP mean annual potential evapotranspiration is generally greater than mean annual precipitation during 1971–2004 at meteorological stations located mostly in valleys where human inhabit, implying water limitation in these areas. However, in the vast western TP and other un-inhabited areas there are essentially no meteorological stations available for assessing the Budyko’s theory. Our study using the interpolated forcing data and the LPJ model suggested that temperature is more important than precipitation when averaged over the entire TP for NPP annual variability and variability and trend (e.g., Figs. 6, 7, 10, 12). In contrast, the tropical region and Australia are blessed with abundant solar radiation year around, and these regions are likely water limited (Nayak et al. 2015; Donohue et al. 2009). However, with increasing CO₂, the water or energy limiting factors might not be the sole natural reasons for vegetation change, as demonstrated by Donohue et al. (2009) in Australia.

The advantage of using the modeling approach is that any exaggerated influences that are resulted from higher inter-annual variability in precipitation and cloud cover could be removed. This is because on the TP historical precipitation and cloud cover displayed higher inter-annual variability (standard deviation of 21.6 mm for precipitation and 1.40% for cloud cover) but much smaller trends (0.50 mm/year for precipitation and 0.02%/year for cloud cover) than temperature (standard deviation of 0.64 °C and a trend of 0.044 °C/year). The higher inter-annual variability could result in greater impacts on NPP with the consequence of inflated responses for long-term trend impacts even though the long-term trends are small. By using scenarios of including variability only, trend only, and combined trend and variability for individual climate elements, we can isolate and identify the impacts of the inter-annual variability and long-term trend of individual climate elements on annual NPP. This approach can also be used in the other impact investigations through separating variability and long-term trend and their effects.

This study identified that CO₂ and temperature are the most important factors for long-term change and inter-annual variability of annual NPP on the TP, respectively. With continued warming and increasing CO₂, the positive correlation between the regionally averaged annual temperature/CO₂ and NPP may hold an optimistic future for ecosystem development on the TP in general. However, across the region, temperature trend impacts on NPP vary spatially and negative impacts are found in the sparsely vegetated western region due to intensified evapotranspiration and enhanced drought related to significant temperature increases (Figs. 7d, 8d). Thus, the benefit of the dominant effect of temperature change also depends on the baseline conditions that the region has. Then a logical question to ask is: how much of the detrimental effects of the enhanced precipitation and temperature variability (i.e., extreme events) on NPP will be canceled out or strengthened by the effects of the long-term temperature and CO₂ changes? This needs to be addressed in future studies.

NPP is the product of integrative processes involving temperature, precipitation, cloud cover and CO₂, and it is difficult to isolate impacts from each climate factor in the field measurements or remotely sensed products. A modeling study has the advantage of directly revealing the dominant influencing factors or causality between impacts and climate factors by using carefully designed scenarios. However, this study only uses the LPJ model and the analysis results will need to be verified using other models as well. Nevertheless, unlike other studies that focused on either climate change impacts or climate variability impacts, this study considered both aspects of climate variations and examined the effects of shifts in long-term mean conditions and short-term variability around the mean conditions. This kind of design should be useful for producing synergetic knowledge in the relationship between climate and NPP and for revealing the underlying mechanisms, which could provide better assistance for local Tibetan nomads and many local communities in designing short-term and long-term mitigation strategies and policies. On the TP, most of these communities are entirely dependent on rain-fed cropland and grassland and the natural grassland is also managed only by using fences. Clearly, both short-term variability and long-term change in climate have direct impacts on the livelihood of the local communities and a better understanding of the impacts and informed mitigation approaches are important.

MODIS and field observed NPP was reasonably well simulated by the LPJ model as both LPJ and MODIS had similar regionally averaged mean annual NPP (modeled 120 vs. observed 126 g/cm²/year), significant spatial correlation of 0.58 (p < 0.05) in mean annual NPP between LPJ and MODIS. Significant R² of 0.41 in growing season maximum NPP between LPJ and field observations was also obtained. Regionally averaged inter-annual variability of NPP represented by standard deviation was 1.14 g/cm²/year, and long-term trend of annual NPP was 0.086 g/cm²/year².

At monthly time scale, NPP was significantly correlated concurrently (R² in 0.52–0.77) and at 1-month lag (R² in 0.72–0.74) with temperature, precipitation and cloud cover. At annual time scale, relative to mean, annual NPP variability was high (low) in the areas where mean annual NPP was low (high). The influences of the variability of temperature, precipitation and cloud cover on annual NPP variability were spatially heterogeneous, and temperature variability dominated NPP inter-annual variability over the TP as reflected by high R² of 0.74 between annual mean NPP and temperature. On the other hand, long-term NPP changes were strongly correlated with CO₂ increase, with a spatial correlation of 0.60. For the impacts of the trend & variability in climate elements, temperature, with the highest spatial correlation coefficient of 0.67, was still the major factor, followed by precipitation and cloud cover on the TP. Overall, parameter uncertainty had very limited effects on the LPJ model simulated NPP on the TP. These findings will benefit local communities and stakeholders by assisting in designing short-term and long-term polices and strategies for mitigating the potential detrimental impacts of climate variations.