Dynamic characteristics and synergistic effects of ecosystem services under climate change scenarios on the Qinghai–Tibet Plateau

Study area

The QTP is located in southwestern China (25° ~ 40°N, 75° ~ 103°E), with a total area of 2.5 million km² and an average elevation above 4000 m (Fig. 7). The QTP is mainly covered with permafrost and grassland, with areas of glacier and desert⁴⁸. The QTP, also known as the “Asian Water Tower”⁴⁹, is the source of 13 major Asian rivers (e.g., the Indus, Ganges, Brahmaputra, Yangtze, and Yellow Rivers). The QTP has a clod, arid climate, with an annual average temperature below 0 °C and an annual mean precipitation of 400 mm. The seasonal distribution of precipitation is uneven, with most precipitation concentrated in the period June to September. There is a decreasing trend in precipitation from the southeast to the northwest of the plateau⁵⁰. Known as the “Roof of the World” and “Third Pole”, the QTP is also an area that is sensitive to global climate change, showing increasing warming and humidification in recent decades⁵¹. In addition, the QTP contains a diversity of ecosystems and fosters a historic ecological security barrier, which nurtures the development of animal husbandry and diverse cultures.

Data sources

RCP scenarios and climate change dataset

The RCP scenarios released by the IPCC 5th Assessment Report⁵² supply a forecasting standard for climate change research. RCP values ranging from 2.6 to 8.5 reflect radiation forcing values in 2100 relative to the beginning of the Industrial Revolution in 1750⁵³. Different radiative forcing scenarios represent different future climate scenarios. RCPs consist of one high-emission scenario (8.5 ({text{W}} cdot {text{m}}^{ – 2}), RCP8.5), two medium-emission scenarios (6.0 ({text{W}} cdot {text{m}}^{ – 2}), RCP6.0; 4.5 ({text{W}} cdot {text{m}}^{ – 2}), RCP4.5), and one low-emission scenario (2.6 ({text{W}} cdot {text{m}}^{ – 2}), RCP2.6)⁵⁴. In this study, we adopted the RCP2.6, RCP4.5 and RCP8.5 climate change scenarios choosing RCP4.5 to represent the medium emission scenario in consideration of increasing activity through global initiatives in response to climate change. Specific descriptions of each scenario are shown in Table 1.

We adopted the climate change dataset outputs from five global circulation models(GCMs) (namely GFDL-ESM2M, HadGEM2-ES, IPSLCM5ALR, MIROC-ESM-CHEM, and NorESM1-M) within the fifth phase of the Coupled Model Intercomparison Project (CMIP5)⁵⁵. The dataset outputs from GCMs were downscaled to a resolution of 0.5° and bias-corrected with Water and Global Change (WATCH) data (Integrated Project Water and Global Change, http:/www.eu-watch.org/data_availability)56. The baseline period of the dataset is 1950–2005 and the forecast period is 2006–2099.The climate change dataset included daily precipitation, air pressure, solar radiation, air temperature, maximum air temperature, minimum air temperature, wind speed, and relative humidity.

Auxiliary data

The auxiliary data for our research include the following. (1) The land use and land cover (LULC) map was obtained from the Resource and Environment Science and Data Center (RESDC), Chinese Academy of Sciences (https://www.resdc.cn) for 1980, 1990, 1995, 2000, 2005, 2010, 2015 and 2020 at a 1 km resolution. The LULC data have six major classes: cropland, grassland, forestland, water, built-up land and barren land. (2) The spatial distribution of soil type data, digital elevation model (DEM), watershed boundaries and normalized difference vegetation index (NDVI) data with a resolution of 1 km were obtained from the RESDC. (3) Soil physical and chemical property data (available soil water capacity, absolute depth to bedrock, silt content, clay content, sand content and soil organic carbon content) were obtained from the International Soil Reference and Information Centre (ISRIC Data Hub) (https://data.isric.org) with a 1 km spatial resolution. (4) During 1986–2005 and 1986–2098 (RCP2.6; RCP4.5; RCP8.5), the permafrost datasets in the Northern Hemisphere (https://doi.org/10.12072/ncdc.CCI.db0032.2020) and the response of the alpine grassland ecosystem to climate change (RCP2.6, RCP4.5, and RCP8.5) in the permafrost region of the Qinghai-Tibet Plateau from 1981 to 2099 (https://doi.org/10.12072/ncdc.CCI.db0006.2020) were provided by the National Cryosphere Desert Data Center (https://www.ncdc.ac.cn).

Future land use simulation and validation

In this study, we used the Future Land Use Simulation model (FLUS) to simulate the LULC in 2030, 2050 and 2100 under the three RCP scenarios. This model was developed by⁵⁷ and is available for download at (www.geosimulation.cn/flus.html). The FLUS model is an efficient land use simulation tool and has been widely used^58,59. We selected the DEM, slope, precipitation, temperature, soil type, and permafrost distribution to calculate the suitability probability. Based on the land use transfer from 2010 to 2015, we calculated the total land use in 2030, 2050 and 2100 under three RCP scenarios by the Markov model. To validate the FLUS model, we set 2010 as the starting year and simulated the land use in 2015. The output results were compared with the real 2015 land use data, and we calculated the Kappa coefficient as follows:

where (P_{0}) is the number of pixels converted correctly,(P_{C}) is the correct number of pixels to be converted in the random case, and (P_{P}) is the correct number of pixels to convert under ideal conditions.

Assessment of ecosystem services under different RCP scenarios

This study assessed four ESs namely WY, SR, CS, and RMP, under climate change scenarios in 1980, 1990, 1995, 2000, 2005, 2010, 2015, and 2030 (short-term); 2050 (medium-term); and 2100 (long-term). We adopted the Integrated Valuation of Environmental Service and Tradeoffs (InVEST)⁶⁰ model to assess the WY, SR, and CS ecosystem services. The InVEST model developed by the Natural Capital Project(www.naturalcapitalproject.org) is an effective model to evaluate ESs⁶¹ and is widely used in ES research on the QTP^22,23,24,25. All spatial data were processed into a 1 km resolution and Albers projection by ArcGIS 10.2 before input into the InVEST model. The data requirements of the InVEST model and its processing are shown in Table S1. We use net primary productivity (NPP) to evaluate the RMP, and NPP can be used to represent the richness of biomass and the supply of organic materials. We adopted the Carnegie-Ames-Stanford Approach (CASA)⁶² model to estimate NPP.

Water yield

Water yield is a key ecosystem service. It refers to the annual quantity of water available for human use, as measured by the supply of surface water per unit area⁶³. We adopted the InVEST 3.9.0 water yield model to estimate WY services in the QTP region. The water yield model is based on the water balance principle⁶⁴. The biophysical parameter table required by the model is shown in Table S2. The parameters in the biophysical table come from the published literature^26,63,65. The annual WY is calculated as follows:

where (Y_{xj}) is the annual WY of land cover type j in pixel x; (P_{x}) is the annual average precipitation of pixel x; and (AET_{xj}) is the actual evapotranspiration of land cover type j in pixel x.

where (omega_{x}) is a dimensionless nonphysical parameter representing soil properties under natural climate conditions. The calculation method is as follows:

where Z is a seasonal rainfall factor representing the regional precipitation distribution and other hydrogeological characteristics. The higher the Z value is, the less the seasonal constant Z affects the model results⁶⁶. Since the QTP region belongs to the arid and cold climate zone in China, the Z value is set as 9. (AWC_{x}) is the soil effective water content of pixel X, which is determined by the soil depth and physical and chemical properties. (R_{xj}) is the Budyko dryness index, which is calculated as follows:

where, (K_{xj}) is the reference crop evapotranspiration and (ET_{0}) is the reference evapotranspiration in pixel x. We adopted the modified Hargreaves method to calculate (ET_{0}).

In the above formula, (T_{av}) represents the average daily maximum temperature and minimum temperature, (TD) represents the difference between the daily maximum temperature and minimum temperature, (RA) represents astronomical radiation (MJm^–2d^-1) and P represents precipitation (mm/month).

Soil retention

Soil retention refers to the ability of various land cover types to prevent soil erosion. The InVEST 3.9.0 sediment delivery ratio (SDR) was employed to estimate SR services in the QTP region. The SDR model is based on the Revised Universal Soil Loss Equation (RUSLE)⁶⁷, and the model is calculated as follows:

where SR is the total amount of soil retention (tons ha^-1 a^-1), LS is the topographic factor, and LS is calculated from the slope length factor (L) and slope steepness factor (S). C is the vegetation and management factor. P is the support practice factor. C and P are shown in Table S2. R is the rainfall erosivity index(MJ mm ha^-1 h^-1 a^-1), which was calculated via monthly precipitation²⁸. K is the soil erodibility, which was calculated from the sand, silt, clay and organic soil moisture contents⁶⁸. R and K are calculated as follows:

where P_i is the precipitation in month i. SAN, SIL, CLA, and SOM are the contents of sand, silt, clay and organic moisture, respectively. Other parameters are shown in Table S1.

Carbon storage

Carbon storage services refer to the carbon that ecosystems store in vegetation, soil and debris. The InVEST 3.9.0 carbon model uses a simple method to estimate CS based on land use data. The carbon pools in this model include four categories: aboveground carbon, belowground carbon, soil organic carbon and dead organic matter. This model simplifies the carbon cycle, and the change in carbon storage is mainly caused by change in land use69. The carbon pools for land use types were set according to published literature^70,71,72. The carbon storage is calculated as follows:

where ({text{C}}_{{{text{total}}}}), (C_{above}), (C_{below}), (C_{soil}) and (C_{dead}) are the total carbon storage, aboveground carbon, belowground carbon, soil organic carbon and dead organic matter, respectively.

Raw material provision

Raw material supply refers to the organic matter provided by the ecosystem for human production and life, such as pasture and wood. In this study, RMP was quantified by the annual NPP. The NPP in the QTP region is calculated by the CASA model, which is a light use efficiency model driven by climate and remote sensing data^73,74. The CASA model has been widely used to estimate NPP in terrestrial ecosystems^75,76. In the CASA model, NPP is calculated as follows:

where, (APARleft( {x,t} right)) is the photosynthetically active radiation(MJ m^-2) absorbed by pixel x in month t, (varepsilon left( {x,t} right)) is the actual light energy utilization rate(gC MJ^-1), and the (APARleft( {x,t} right)) calculation method is as follows:

In the formula, (SOLleft( {x,t} right)) is the total solar radiation in pixel x in month t(MJ M^-2); (FPARleft( {x,t} right)) is the absorption ratio of photosynthetically active radiation by vegetation, which is determined by the normalized difference vegetation index (NDVI); and the constant 0.5 is the proportion of photosynthetically active radiation to the total radiation. (SOLleft( {x,t} right)) can be calculated by the solar shortwave radiation as follows:

where, (R_{s}) is the solar shortwave radiation(MJ M^-2 d^-1), n is the actual sunshine time(hours), N is the time of day(hours), and (frac{n}{N}) is the relative sunshine time; The constants (a_{s} = 0.25) and (b_{s} = 0.5).

And the (varepsilon left( {x,t} right)) is calculated as follows:

where, (T_{varepsilon 1}) and (T_{varepsilon 2}) are the stress factors of cold and heat, respectively; (W_{varepsilon }) is the water stress factor, reflecting the influence of water conditions; (varepsilon_{max}) is the maximum light use efficiency(gC MJ^-1) under the optimal conditions, in this study, (varepsilon_{max}) is 0.389.

Trend analysis

The Mann–Kendall nonparametric test and Sen’s slope estimator were used to analyze the trend of ESs in the QTP region. The Mann–Kendall method is widely used to analyze climatic and hydrological time series variation trends⁷⁷. The advantage of the Mann–Kendall test is that it does not require the sample to follow a certain distribution, allows the existence of missing values, is not affected by a small number of outliers, and has strong quantitative ability⁷⁸. The Mann–Kendall test is as follows:

For time series data, i.e., {x₁, x₂, …, x_n}, n is the length of the data, and (sgnleft( {x_{j} – x_{i} } right)) is derived as:

In this study, we set the significance level of (alpha = 0.05), when (left| Z right| le Z_{1 – alpha /2}) accepts the null hypothesis. Otherwise, the null hypothesis is rejected, and the trend is statistically significant.

where p is the number of nodes in the dataset and (t_{j}) is the length of the nodes.

Sen’s slope estimator is an estimation method based on the median and its insensitivity to outliers⁷⁸.

Trade-offs and synergy analysis

Synergies and trade-offs were used to describe the relationships among the ESs. A trade-off analysis was conducted to reflect the difference in ESs and their responses to climate change. Trade-offs are when ESs change in the opposite direction. Synergies are when ESs change in the same direction⁷⁹. Correlation analysis is often used to evaluate trade-offs and synergies between ESs². To analyze the trade-offs and synergies of ESs at different administrative and natural scales, we allocated the ES values at the 10 km (pixel), county and watershed scales by the “zonal statistic” module of ArcGIS 10.2, and conducted minimum–maximum normalization in R4.0.3 (www.R-project.com). To analyze the relationship between any two of the four ES types, the R package PerformanceAnalytics was adopted to measure the Spearman correlation matrix at different scales.

Source: Ecology - nature.com