Increasing gross primary productivity under soil warming and wetting on the Tibetan Plateau

The soil freeze-thaw process has undergone significant changes on the Tibetan Plateau (TP) in the context of global change, resulting in the changes of soil physical and chemical properties, thereby affecting the vegetation phenology and photosynthesis through affecting the utilization capacity of CO2 and light by vegetation. However, little is known about how soil temperature (ST) and soil moisture (SM) affect the gross primary productivity (GPP) on the TP at different seasons and elevations. In this study, the spatiotemporal variation patterns of GPP, ST, and SM were analyzed based on the Community Land Model version 5.0 (CLM5.0) simulations in order to illustrate the impacts of ST and SM in surface (0–10 cm) and root zone soil (0–100 cm) on GPP between 1979 and 2020. The results showed that the CLM5.0-based GPP and ST were in good agreement with in situ observations. ST, SM and GPP increased at the rates of 0.04 °C a−1, 2.4 × 10−4mm3 mm−3 a−1, and 5.36 g C m−2 a−2, respectively. SM dominated the variations of GPP in winter (64.3%), while ST almost was the dominant factor in other periods, especially spring (99.9%) and autumn (94.7%). The explanatory power of ST and SM for GPP increased with elevation, especially for ST. The relative contributions of ST and SM to GPP at different time scales in root zone soil were similar to those in surface soil. This study provided a new understanding of how soil freeze-thaw affected GPP changes on the TP in the context of the intensification of warming and humidification.


Introduction
The intensification of fossil fuel combustion, human activities, and land use changes exacerbate global warming, with a significant impact on the soil freezethaw process and the productivity of terrestrial ecosystems.In particular, the Tibetan Plateau (TP), which is extremely sensitive to climate change, has experienced almost twice the degree of global warming during the same period (Yan et al 2016, 2020, Palazzi et al 2017, You et al 2017).Permafrost is widely distributed on the TP, and the soil freeze-thaw process is one of the most significant land features.The carbon storage in permafrost regions is the largest component of the terrestrial carbon pool, almost twice that of the atmosphere (Schuur et al 2015, Wu et al 2017, Natali et al 2021, Wu et al 2022, Wang et al 2023a).The soil freeze-thaw process has changed significantly under the background of climate change, characterized by delay of the freeze start-date (FSD), advance of freeze end-date (FED), and thickening of the active layer (Wang et al 2015, Luo et al 2020, Liu et al 2022, Peng et al 2023).Furthermore, the changes in hydrothermal states and biochemical characteristics of soil improve the physiological state of vegetation and increase the photosynthesis rate, thus increasing the gross primary productivity (GPP) of the terrestrial ecosystem (Schimel et al 2014, Shen et al 2015, Wang et al 2019), which plays a significant role in regulating the ecosystem function of the TP (Wang et al 2012, Chen et al 2013).Therefore, studying the impacts of the soil freeze-thaw process on GPP will provide theoretical support and a scientific basis for predicting the changes in carbon sources and sinks in terrestrial ecosystems in the future.
The GPP is defined as the gross carbon fixed by green plants through photosynthesis per unit time and area (Wu et al 2014, Zhang et al 2018), which can not only reflect ecosystem production capacity (Beer et al 2010, Li et al 2016b), but is also an important indicator of the terrestrial carbon cycle.Previous studies used in situ observations and remote sensing datasets to explore the effects of soil temperature (ST) on GPP in arctic ecosystems.Mauritz et al (2017) used CO 2 flux data to discuss the effects of soil warming on GPP in arctic permafrost, and showed that soil thawing stimulated vegetation growth and increased GPP.In the arctic, the beginning of the phenological period in many ecosystems was determined by snowmelt date (Semenchuk et al 2016).Leblans et al (2017) and Descals et al (2020) discussed the effect of soil thawing on length of the growing season in the alpine biomes and subarctic grasslands of the northern hemisphere, respectively.The results indicated that soil warming led to advances in soil thawing, and vegetation took the phenomenon of soil thawing as a key signal to restore physiological activities (Descals et al 2020).
Different soil moisture (SM) patterns and vegetation responses might be related to soil freeze-thaw dynamics (Jiang et al 2018).The soil freeze-thaw process redistributes SM by affecting the movement of soil water, resulting in changes in soil physical properties, which affect vegetation phenology (Jiang et al 2018).There are extremely complicated energy and hydrologic exchange processes in soil, vegetation, and atmosphere during the alternating process of the soil freeze-thaw cycle.However, previous studies have mainly focused on arctic ecosystems, and relatively little attention has been paid to the effects of ST and SM on the GPP over the TP in the context of global warming.Therefore, investigating the dominant factors affecting GPP at different periods on the TP will help to further understand the underlying mechanisms of atmosphere-soil-vegetation interaction (Wang et al 2019, You et al 2020).
In this study, high resolution simulations based on the Community Land Model version 5.0 (CLM5.0)were used to investigate the spatiotemporal variation characteristics of GPP, ST, and SM in surface (0-10 cm) and root zone soil (0-100 cm) of the TP.Then, the effect of ST and SM on GPP at different time scales was quantitatively analyzed.Moreover, the impact of elevation on the explanatory power of ST and SM for GPP was also examined.

CLM5.0 simulations
The CLM released by the National Center for Atmospheric Research is the land component of the Community Earth System Model (Dickinson et al 2006, Swenson et al 2012).In this paper, the latest version of CLM5.0 was selected to simulate GPP, ST, and SM during 1979-2020 over the TP.High spatial land surface heterogeneity in CLM is represented as a nested subgrid hierarchy.Each grid cell is composed of multiple land units, each land unit can have a different number of columns, and each column has multiple plant functional types (PFTs).We chose the PFT data from Ran and Li (2019) and soil properties data from Shangguan et al (2013) for CLM5.0 simulations instead of the default datasets.The elevation data (figure 1(a)) were collected from the Resource and Environment Science and Data Center (www.resdc.cn).The CO 2 concentration data for 1979-2020 were collected from the National Oceanic and Atmospheric Administration (NOAA) Mauna Loa Observatory (https://gml.noaa.gov/ccgg/trends/data.html).A high resolution meteorological forcing dataset with 1/30 • × 1/30 • was used to force CLM5.0.This included precipitation, temperature, downward longwave radiation, incident shortwave radiation, pressure, wind speed, and specific humidity.The forcing data were generated by combining a short-term high-resolution Weather Research and Forecasting Model simulation with ERA5 reanalysis and then merging observations using the Climatology Aided Interpolation and Random Forest methods (Jiang et al 2023).Before the simulation, a long-term cyclically spin-up run (2500 years) was prepared to allow the simulated data to reach their equilibrium state.After that, a normal transient simulation was initiated and run with a spatial resolution of 0.1 • × 0.1 • and temporal resolution of 1800 s.

In situ observations
This study used in situ ST observations at a depth of 0.1 m at three stations of the National TP Data Center (figure 1): BJ (91.9  2015, 2017-2018, 2010, and 2014-2015, respectively.

Satellite-derived and upscaled GPP products
In this study, two other GPP products were used to compare with CLM5.0 simulations.The first GPP

Identification of soil freeze-thaw status
The analysis of soil freeze-thaw process mainly included the determination of FSD, FED, and thaw duration (TD).The definitions of FSD and FED were adopted from Guo et al (2014,2018).When the ST changed from above 0 • C to below 0 • C for three consecutive days, the first day was defined as the FSD.
When the ST changed from below 0 • C to above 0 • C for three consecutive days, the first day was defined as the FED.The 3-day threshold was chosen to avoid the potential impact of random weather process on ST.The TD was the duration from FED to FSD in the current year.

Analysis methods
Because ST and SM are not independent, a partial regression method was used to analyze the contributions of ST and SM to GPP (Heikkinen et al 2005, Tang et al 2018).The specific formulas follow: where, a1, b1, c1, d1, a2, b2, c2, a3, and b3 are regression parameters.It is assumed that a, b, and c represent the individual impact of ST, SM, and the synergistic impact of ST and SM on GPP, respectively.The coefficients of determination of equations ( 1)-(3) were a + b + c, a + c, and b + c, respectively.Consequently, a is the difference between the coefficient of determination of equation ( 1) and that of equation ( 3), and b is the difference between the coefficient of determination of equation ( 1) and that of equation ( 2).After obtaining the values of a and b, c could be calculated by the difference between the coefficient of determination of equation ( 1) and the sum of those of a and b.In this way, the individual contribution and synergistic impact of ST and SM on GPP were obtained.

Validation of CLM5.0 simulations
To validate and evaluate the simulation results, the monthly mean ST (figure 1(b)) at depth of 0-10 cm and GPP (figure 1(c)) were compared with in situ observations.The simulated ST (R 2 = 0.96, p < 0.01) and GPP (R 2 = 0.63, p < 0.01) were in good agreement with the observations, and the scatter points were practically distributed around the line y = x, especially at MAWORS, QOMS, and Zoige stations (R 2 = 0.98, p < 0.01).Then the spatial distribution patterns of annual mean GPP based on satellite and upscaled products were compared with the simulation results for 2001-2018 (figure S1).The GPP simulated by CLM5.0 had similar spatial distribution characteristics to the other datasets, showing a gradually increasing trend from northwest to southeast, but the values in some central regions were larger than that of other products.The values of GPP might differ due to different PFTs (figure S2).

Spatiotemporal variation characteristics of ST, SM, and GPP
Figure 2 shows the spatiotemporal variation characteristics and trends of ST, SM, and GPP.The spatial patterns of ST at depth of 0-10 cm generally increased gradually from north to south (figure 2(a)) with SM and GPP increasing from northwest to southeast (figures 2(c) and (e)).The spatial patterns of ST and SM at depth of 0-100 cm were similar to those of surface soil (figure S3).The surface ST and SM warmed and wetted the soil on the TP at average area rates of 0.04 • C a −1 and 2.4 × 10 −4 mm 3 mm −3 a −1 , respectively (figure 3), and the corresponding trends in the root zone soil were 0.04 • C a −1 and 2.7 × 10 −4 mm 3 mm −3 a −1 (figure S4).Driven by the warm and wet environment, the FED was advanced with an area-mean value of 0.28 d a −1 , and the FSD was delayed 0.16 d a −1 , resulting in an extension of the TD at a rate of 0.44 d a −1 , thus prolonging the TD by 18.48 d during the study period (figure S5).Due to the extension of TD, the amount of organic carbon fixed by vegetation through photosynthesis increased.The spatial distribution of the trend of GPP was similar to that of ST and SM (figures 2(b), (d), and (f), and increased at an average area rate of 5.36 g C m −2 a −2 (figure 3).The ST, SM, and GPP all had large declining trends in 1997, so 1997 was used as the boundary to separate the time series of ST, SM, and GPP into two different periods.During 1979-1997, ST, SM, and GPP increased at rates of 0.01 • C a −1 , 1.4 × 10 −4 mm 3 mm −3 a −1 , and 3.31 g C m −2 a −2 , respectively.Corresponding increase rates during 1997-2020 were 0.03 • C a −1 , 3.1 × 10 −4 mm 3 mm −3 a −1 , and 3.46 g C m −2 a −2 (figure 3).Compared with the previous period, there were significant increases in GPP, ST, and SM in both surface and root zone soil for 1997-2020, so it was necessary to further explore the relationships between them.

Relationships of GPP with ST and SM
The relationships of GPP with ST and SM were considered for different time scales.The GPP was significantly positively related with ST (R 2 = 0.57, p < 0.01) and SM (R 2 = 0.42, p < 0.01) at annual scale (figures 4(a) and (b)), and the relationships of GPP with ST (R 2 = 0.60, p < 0.01) and SM (R 2 = 0.54, p < 0.01) in the root zone soil were similar to those for 0-10 cm (figures S6(a) and (b)).The warm and wet conditions provided by soil are beneficial for vegetation growth and promote the increase of GPP.In this paper, these relationships were considered on the seasonal scale by constructing exponential equations of GPP and ST and regression equations of GPP and SM,  respectively.The results showed that GPP was positively correlated with ST and SM in the four seasons, especially in spring (R 2 (ST) = 0.93, R 2 (SM) = 0.46, p < 0.01) and autumn (R 2 (ST) = 0.97, R 2 (SM) = 0.91, p < 0.01) (figures 4(c) and (d)).The ST explained more of the GPP than SM in spring and autumn, but SM explained more than ST in winter, and the explanatory power of both was comparable in summer.

Effects of ST and SM on GPP
The relative contributions of ST and SM to GPP were calculated (figure 5).The SM dominated the variation in GPP for winter, accounting for 64.3% of the relative influence, with comparable effects in summer, while ST made the main contribution in   other periods: spring (99.9%), and autumn (94.7%).
In spring, the increase of ST was beneficial to the advance of soil thawing of large areas on the TP (Li  2023b).In contrast, the soil was already completely frozen in most areas in winter, and the increase of SM could provide the SM conditions required for vegetation growth (Wang et al 2018), so SM accounted for the main contribution in winter.In summer, the continuous increase of ST resulted in the soil completely thawing in large areas, and the increase of precipitation led to increased SM and soil water availability (Ji and Fan 2020, Weijers 2022), which promoted vegetation growth.The ST in the root zone soil remained dominant for most periods, except winter, consistent with the characteristics of surface soil (figure S7).

Spatial patterns of dominant factors in controlling GPP
After analyzing their relative contributions, the spatial patterns of dominant factors controlling GPP were discussed (figure 6).1(a)).The results demonstrated that the explanatory power of ST (R 2 = 0.21-0.51)and SM (R 2 = 0.12-0.48)for GPP increased with greater elevation (figure 7), with similar characteristics for root zone soil (figure S9).To verify the result, the elevation interval was converted from 1000 to 50 m, and the relationships between the explanatory power of GPP by ST and SM and the elevation were examined (figure 8).The relationships of GPP with ST and SM were significantly dependent on elevation, especially ST (R 2 = 0.71, p < 0.01).This may be due to the significant decrease in ST with increased elevation, resulting in a higher demand and stronger dependence on ST for vegetation growth (Shi et al 2020, Cui et al 2022).

Conclusions
This study quantified the relative contributions of ST and SM on GPP, based on CLM5.0 simulations.Model simulations were first evaluated using in situ observations and then used to discuss the relationships of GPP with ST and SM on the TP.During the study period, soil of the TP showed characteristics of warming and wetting, with a significant increase in GPP.The SM dominated the variations of GPP in winter, with equivalent effects in summer, while ST was the dominant factor in other periods, especially spring and autumn.Spatially, SM dominated the western TP, while ST dominated the central and eastern regions.The explanatory power of ST and SM for GPP increased with elevation, especially ST.The impacts of ST and SM in root zone soil on GPP at different time scales were similar to those for surface soil.
The variations in ST and SM had a significant effect on the GPP, especially ST, providing new insights into how GPP changes in the context of climate change.Our future work will consider the trends in autotrophic and heterotrophic respiration, which will provide a scientific basis for predicting the changes in carbon sources and sinks in ecosystems on the TP.

Figure 1 .
Figure 1.(a) Elevation distribution of the Tibetan Plateau and locations of gross primary productivity (GPP, purple circles) and soil temperature stations (ST, red triangles).Comparison of simulated and observed monthly mean (b) soil temperature (ST) at three stations and (c) gross primary productivity (GPP) at four stations.Different symbols represent different stations.The black dotted line is y = x.** Significance at α = 0.01.

Figure 2 .
Figure 2. Spatial patterns of surface (0-10 cm) (a) soil temperature (ST) and (c) soil moisture (SM), and (e) gross primary productivity (GPP) during 1979-2020.Spatial patterns of the trends of (b) ST, (d), SM, and (f) GPP, where areas with significance exceeding α = 0.05 are denoted with dots.The units are in the upper right corner of each subgraph.

Figure 3 .
Figure 3.Time series of annual mean gross primary productivity (GPP), surface (0-10 cm) soil temperature (ST), and soil moisture (SM) during 1979-1997 and 1997-2020, where k represents the trends for 1979-2020, k1 represents the trends for 1979-1997, and k2 represents the trends for 1997-2020.Green indicates GPP, red indicates ST, blue indicates SM, and gray vertical line in the figure shows the dividing line with 1997 as the boundary.Dashed lines represent the linear trend.* Significance at α = 0.05, and * * significance at α = 0.01.

Figure 8 .
Figure 8. Relationships between the coefficients of determination of surface (0-10 cm) (a) soil temperature (ST) and (b) soil moisture (SM) on gross primary productivity and elevation.The value of elevation is the mean value of each elevation level.** Significance at α = 0.01.

Table 1
These results indicate that CLM5.0-basedGPP simulations in this study are relatively reliable.

Table 1 .
Statistics of different gross primary productivity products at four stations.MB is the mean bias, and RMSE is the root mean square error, the unit is g C m −2 month −1 .
al 2016a), which promoted microbial decomposition and released nutrients conducive to vegetation growth, resulting in the advance of spring phenology (Descals et al 2020).Similarly, the increase of ST in autumn delayed the FSD, thus delaying the autumn phenology of vegetation (Menzel et al 2006, Yu et al 2010, Zhang et al 2013).The changes in phenology further affected vegetation productivity.

. Dependence of relationships of GPP with ST and SM on elevation Elevation
The ST was the dominant factor in most periods, especially spring and autumn, and the areas dominated by ST accounted for 88.3% and 82.7% of the TP, respectively.This may be because surface soil freezes in autumn (approximately 91.9%) and thaws in spring (approximately 73.4%) in most regions of TP, which is consistent with previous studies showing that the area-averaged value of FSD is the 290th day and the FED is the 107th day on the TP (Guoand Wang 2014, Luo et al 2020).Therefore, ST is the most important factor controlling the soil freeze-thaw state (Jia et al 2021, Li et al 2021, Luo et al 2021, Zhao et al 2021, Fu et al 2022), which controls vegetation growth.In winter, approximately 79.0% and 90.0%(figureS8) of the regions on the TP were dominated by SM in surface and root zone soil, respectively, likely because the extensive soil freezing limits the liquid water content of soil and affects the carbon sequestration capacity of vegetation.The SM dominated the variations of GPP mainly in the western region, especially for summer.This was mainly because there was less precipitation in the TP western region, i.e. the water-limited area, while the central and eastern parts were relatively temperature limited.Therefore, the limiting factor for vegetation growth in the western region was SM, but was ST in the central and eastern regions (Wang et al 2023b).hasasignificant effect on the soil freezethaw process(Peng et al 2016, Cao et al 2017, Luo et al  2020, Liu et al 2022, Xue et al 2023).