Cooling effects of increased green fodder area on native grassland in the northeastern Tibetan Plateau

With increasing livestock production due to high demand for consumption, the planted area of green fodder, an essential livestock supplement, has grown rapidly and will continue to grow in China. However, the climate feedback of this rapid land cover conversion is still unclear. Using multisource data (e.g. remote sensing observation and meteorological data), we compared the land surface temperature of green fodder plantation areas and native grassland in the northeastern Tibetan Plateau. The green fodder area was detected to be cooler than the native grassland by −0.54 ± 0.98 °C in the daytime throughout the growing season. The highest magnitude (−1.20 ± 1.68 °C) of cooling was observed in August. A nonradiative process, indicated by the energy redistribution factor, dominated the cooling effects compared to the radiative process altered by albedo variation. The results indicate the potential cooling effects of increasing green fodder area on native grassland, highlighting the necessity of investigating climate feedback from anthropogenic land use change, including green fodder expansion.


Introduction
Over the past few decades, livestock consumption in China has grown rapidly by 93% per capita, as aggregated by official statistical institutions (www.stats.gov. cn/), causing an increase in demand for grazing supplies. Forage production that merely relies on native grassland is becoming unsustainable for growing forage demands, especially in the current state of practices regarding native grassland degradation [1,2]. These current practices require supplementation of native grasslands to complement forage and protein requirements, among which green fodder stands out, with richer nutrients and higher water content compared to native grass [3].
The northeastern Tibetan Plateau, as an important livestock origin in China, has shown strong momentum toward increasing green fodder area, driven by both socioeconomic demand and policy support [4]. From 2010 to 2019, the green fodder area rapidly expanded by 119.8 km 2 [5], which drove a significant land cover transition in this region. Along with the rapidly increasing green fodder area, most research related to green fodder has focused on forage yield [6,7], leaving the climate feedback of increased green fodder area on native grassland relatively uninvestigated. Additionally, most studies on anthropogenic climate feedback mainly focus on more intuitive land cover conversion processes, such as the decrease of land surface temperature (LST) along with afforestation processes [8]. However, understanding of the feedback from land cover conversion between agricultural fields and climate, especially climatic effects via land surface biophysical processes that are dominated by the surface energy budget, is still in the early stages [9]. The Tibetan Plateau is widely known as a fragile environment that is sensitive to climate change [10,11]. Thus, investigations into the climate feedback of anthropogenic land use and land cover change in this area will be necessary.
Green fodder species, such as oat (Avena sativa), alfalfa (Medicago sativa), and rapeseed (Brassica napus), have different phenologies during the growing season when compared to native grassland [12]. Specifically, green fodder is planted later when the native grassland has already sprouted and greened up, and agricultural management, such as harvests, is implemented to promote productivity [13]. Because of this, native grassland has an earlier peak growth stage, as indicated by vegetation indices (e.g. the normalized difference vegetation index, the enhanced vegetation index, and the normalized difference water index) in the growing season [5]. The green fodder planted in the northeastern Tibetan Plateau is mainly rainfed, and the main forage types are fodder oat and fodder rapeseed, which are planted and harvested during similar time windows. Overall, the green fodder shows higher water content and brighter green reflectance values from stems and ears in the late growth stage (from July to September). Changes in the phenology of vegetation alter land surface characteristics, which further change the surface energy budget via radiative (namely the changes in albedo) and nonradiative (namely the changes in the energy redistribution process) phenomena, and ultimately provide feedback to climate factors like LST [14].
Statistical results can help us to understand the differences in LST between different land cover types [8]. However, the underlying mechanism of change in LST can be better understood by modeling. Variation of LST between different land covers is affected by a combination of factors, including the differences in input radiative energy and the redistributed energy via latent heat and sensible heat fluxes. The temperature response model, initially developed by Lee et al [15] and further revised by Bright et al [16], has been demonstrated to be a dependable method of measuring temperature variation caused by land cover conversion. Based on the temperature response model, a previous study found nonradiative dominated cooling effects on the land surface along with paddy rice (Oryza sativa) expansion [17]. The temperature response model could help to further investigate the functional mechanisms of LST variation induced by different energy distribution processes.
Although green fodder expansion is still at an early stage in China, there is great growth potential for its cultivation driven by the quickly growing forage demands from the livestock industry [5]. Thus, the magnitude of climate feedback via biophysical processes will become more and more significant with the increasing scale of green fodder plantations. A better understanding of these biophysical processes is needed to assess the impacts of the climate feedback of green fodder expansion at a large scale.
To quantify the climate impact of increased green fodder area on native grassland, we compared differences in LST between pure green fodder and native grassland pixels, which were derived from green fodder tracing production and land use classification products at 30 m resolution [5]. The LST in green fodder fields and native grassland, both in the daytime and nighttime, were compared based on a pairwise comparison analysis [8,15,18]. We also explored the mechanism of LST variation along with the conversion from native grassland to green fodder by utilizing a temperature response model [16], which illustrated the potential climate feedback of increasing green fodder area on native grassland.

Study area
The study area is located in the northeastern Tibetan Plateau, which covers Guinan, Zeku, and Tongde counties in Qinghai Province of China (figures 1(a) and (c)). Belonging to the plateau continental climate zone, the climatic condition in this region is relatively cold and dry. The mean annual air temperature is around 2.0 • C, and average annual precipitation is around 450 mm. According to the Yearbook data of statistical departments, the overall monthly air temperature from May to September is parabolic (figure 1(d)) (http://tjj.qinghai.gov.cn/nj/2020/ indexch.htm). Alpine meadow is the dominant vegetation type in the region. As one of the main livestockproducing regions in China, expanding native grassland degradation was detected in the northeastern Tibetan Plateau due to management issues, such as overgrazing. Green fodder planting, as an efficient supplementary method for the shortage of forage required by livestock production, has been widely popularized and implemented by the local government ( figure 1(b)). The area of green fodder increased dramatically by 737% from 16.26 km 2 in 2010 to 136.09 km 2 in 2019 [5]. Undoubtedly, the increase of green fodder area in this region is in its primary stage and will quickly expand over large areas [19]. The green fodder planted in the study area is mainly rainfed, and the main forage types are fodder oat and fodder rapeseed, which are planted and harvested during similar time windows [4,20].

Satellite data and land cover map
Optical and thermal data from the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor in 2019 were used, specifically the average 8 day 1 km LSTs (MYD11A2) from the Aqua satellite. A noteworthy feature of the MYD11A2 data product is that the observation is applied twice a day, representing both the daytime (local solar time ∼01:30 PM) and nighttime (∼01:30 AM) LSTs [21,22]. The daily 500 m albedo (MCD43A3) is generated from a combination of the Terra and Aqua satellite observations. The emissivity of the land surface is estimated from the average of emission band 31 and emission band 32 data [23,24]. The albedo data were further utilized as the essential factor for the net radiation energy calculation in this study. We took the average of the shortwave band (0.3-0.5 µm) reflective of the black sky and the white sky in each pixel as the albedo value [25,26].
The green fodder pixels were derived from the 30 m green fodder detection dataset in the Tibetan Plateau [5], while the native grassland pixels were extracted from the 30 m global land cover data product (GLC_FCS30-2020: Global Land Cover with Fine Classification System at 30 m in 2020), which were derived from the 2019 to 2020 time-series Landsat surface reflectance data [27]. The accuracy of these two datasets was 96.6% and 82.5%, respectively. To match the spatial resolution of the MODIS-derived LST data, all the pixels were aggregated to 1 km spatial resolution grids. The internally homogeneous grids (more than 50% of the area of a certain land cover type) were defined as pure pixels, and a pure pixel map was further generated based on the pure pixels of the two land covers ( figure 1(b)).

Meteorological data
Meteorological data were derived from the Famine Early Warning Systems Network's Land Data Assimilation System (FLDAS) dataset, the TerraClimate monthly reanalysis data, and a 1 km resolution monthly air temperature dataset for China in 2019 (GPRChinaTemp 1 km: 1 km monthly minimum air temperature for China from January 1951 to December 2020). Specifically, the FLDAS dataset used here was generated based on a Noah's version 3.6.1 surface model that provides daily downward long-wave radiation flux values [28]. One specific band in the TerraClimate monthly reanalysis dataset provides information about downward shortwave radiation [29], which can be further used in temperature response model simulations. The air temperature data were derived from the GPRChinaTemp 1 km reanalysis dataset, which is generated using the Gaussian process regression model based on machine learning [30]. We aggregated all the meteorological datasets into monthly mean values and resampled the data to 1 km spatial resolution.

Pairwise comparison
The pairwise comparison method is based on the concept of 'space as a substitute for time' , which originates from a comparative experiment in ecology research [31]. Currently, this method is increasingly being applied to land cover change research [8,32]. In this study, we defined the sample areas with 3 × 3 km size to satisfy two prerequisites: the climatic condition within each sample area is similar, and the ratio of the pure pixels of green fodder is >10%. Grids were randomly generated on a 3 × 3 km scale as sample areas by utilizing the focal neighborhood function method in ArcGIS Pro software. We performed a neighborhood operation to compute the mean values of parameters within a 3 × 3 range around each pure pixel.
Furthermore, since altitude variation affects the temperature significantly [24], we filtered the generated sample areas with a small elevation range (∼±100 m). In total, 164 sample areas were filtered out for pairwise comparisons of LST and other biophysical factors (i.e. albedo and the energy redistribution factor), as well as for the analysis of the underlying mechanisms for the variations in LST. The significance of these differences was statistically evaluated by a t-test (The P-value less than 0.05 was considered significant).

Temperature response model
As an essential parameter in the temperature response model, the energy redistribution factor characterizes the ability of energy to dissipate from the surface. A larger f value indicates that more energy dissipates through sensible and latent heat fluxes [16]. It is calculated using the following formula: where λ 0 is the monthly mean Planck response to the external radiative forcing at the surface (calculated as 1/4ε s δT s 3 , where ε s is the monthly mean surface emissivity derived from the MYD11A2 dataset, and δ is the Stefan-Boltzmann constant), T s is the monthly mean LST, T a is the monthly mean air temperature, G is the mean value of soil heat flux (calculated as 0.14 (T a,n -T a,n−1 )), and R n is the monthly mean apparent net radiation, which can be computed as where α is the monthly average surface albedo, S is the monthly mean downward shortwave radiation, and L ↓ is the monthly mean downward long-wave radiation.
The revised temperature response model efficiently isolates the contribution to the calculated temperature differences from each biophysical process associated with land cover change [16]. However, the revision is based on many assumptions that may produce a relatively large deviation between calculated results and observations. Since the exact value of ∆f was available, we used (1 + f )(1 + f + ∆f ) instead of (1 + f ) 2 in the revised model [9,17]: In this model, the three right-hand terms are the response to the surface radiative forcing from albedo change (∆T s_alb ), the change in heat conducted by the surface medium (∆T s_G ), and the change in energy redistribution (∆T s_f ), respectively. Predictions with this model were verified to be reasonably robust in a study of local temperature response to land cover change on time series [16]. Here, the land cover types before and after the change were represented by native grassland and green fodder separately, to explore the potential temperature response for increased green fodder area in the study region. The simulated total variation of LSTs (∆T s ) was verified by the observations from MYD11A2.

Seasonal dynamics of LSTs
We generated the average value of daytime and nighttime LSTs during the growing season within each sample area, both for green fodder and native grassland. The seasonal average daytime and nighttime LSTs of green fodder were 24.84 ± 8.36 • C (mean ± SD) and 1.15 ± 2.78 • C, respectively. These LSTs were significantly (P < 0.05) lower than those  On a monthly scale, daytime LST showed an overall downward trend, with the highest value in May (28.58 ± 1.93 • C and 28.57 ± 1.90 • C for green fodder and native grassland, respectively) and the lowest value in September (20.34 ± 1.75 • C and 21.38 ± 1.63 • C for green fodder and native grassland, respectively; figure 2(a)), while the nighttime LST fluctuated with a parabolic shape and peaked in July (4.26 ± 0.84 • C and 4.39 ± 0.79 • C for green fodder and native grassland, respectively; figure 2(b)).
On the sample areas with vegetation cover (green fodder or native grassland), evapotranspiration helps lower the surface temperature during the daytime, especially from June to August when the vegetation grows abundantly [33,34], while the temperature at nighttime is consistent with the overall regional temperature since the vegetation stomata are closed and transpiration gets weaker [35]. The similar trend and magnitude of the nighttime LST in green fodder and native grassland during growing season, as well as the same parabolic shape of nighttime LST and air tempeature in the study area, indicates that regional temperature dominated the nighttime LST changes ( figure 1(d)). In the contrast, significant differences (P < 0.05) were found in monthly mean daytime LSTs (figure 2(a)) due to the intervention of human management activities since the green fodder field was often considered a kind of arable land use type [5,36]. Such interventions further resulted in differences in phenology of green fodder at certain times of the growing season compared with native grassland, reflected as variations of daytime LST, which were most pronounced in the late growth stages (from July to September).

Asymmetric temporal changes in differences of LSTs (∆T s ) during the growing season
We further explored the differences in LSTs (∆T s ) between green fodder and native grassland within each sample area during the growing season ( figure 3(a)) and on a monthly scale ( figure 3(b)) separately. In general, the average temperature of green fodder was lower than that of native grassland during the entire growing season, while the magnitude of ∆T s was greater for the daytime (−0.54 ± 0.98 • C) than for the nighttime (−0.14 ± 0.57 • C; figure 3(a)).
Monthly daytime ∆T s showed an asymmetric variation, with a small magnitude of positive values in May and June (i.e. the early growth stage) but a  Figures (a) and (b) use the same legend, and the colored lines and lightly shaded areas in (a) and (b) represent the mean and SD of all sample areas, respectively. The ∆Ts_ alb and ∆T s_f in (c) demonstrate the LST differences caused by radiative process (albedo) variation and the nonradiative process (energy redistribution factor), respectively. The ∆Ts _cal and ∆Ts_r in (d) denote the LST differences derived from the model calculation and remote sensing data, respectively. larger magnitude of negative differences from July to September (i.e. the late growth stage; figure 3(b)). Specifically, green fodder was detected to be relatively warmer than native grassland in June by 0.42 ± 1.44 • C (P < 0.05), while no significant difference was found in May. The earlier peak of the native grass growth compared to the green fodder led to a higher vegetation coverage on native grassland during the early growth stages, which is likely to increase the magnitude of evapotranspiration on native grassland and decrease the LST by more energy dissipation compared to green fodder [25,37]. However, this condition reversed gradually with increasing vegetation growth in green fodder fields. The maximum cooling effect of green fodder was detected in August at about −1.20 ± 1.68 • C (P < 0.05), and the magnitude of ∆T s in September and July were −1.04 ± 1.43 • C and −0.91 ± 1.89 • C, respectively (both differences were statistically significant at P < 0.05), indicating that conversion from native grassland to green fodder could cool the land surface during late growth stages. The alteration of monthly ∆T s during the nighttime was relatively small and stable during the entire growing season. Since all the differences were detected as negative from May to September, the magnitude of nighttime ∆T s between green fodder and native grassland ranged from −0.13 ± 0.39 • C (in July) to −0.29 ± 0.68 • C (in May).

Mechanism of LST variations between green fodder and native grassland
The mechanism of LST variations between green fodder and native grassland was established by synergizing the effects on radiative and nonradiative processes. Regarding radiative characteristics, the statistical results showed that green fodder had a relatively higher albedo during the entire growing season (figure 4(b)), especially during the late growth period (from July to September) (P < 0.05). The largest deviation of albedo between these two land cover types was found in September, and the albedo of green fodder was higher than that of native grassland by 0.0076 ± 0.0082. The higher albedo of green fodder indicates that more solar radiative energy is reflected away from the underlying surface, which ultimately cools the land surface by reducing the input of the energy budget [38,39]. Thus, a reduced LST during the growing season was found in the radiative process due to the green fodder area increases on native grassland (figure 4(c)).
As for the energy dissipation process (i.e. nonradiative process) on the land surface, a significant divergence of the energy redistribution factor started in July with a difference of 0.14 ± 0.29 and reached a peak in August, when the green fodder showed a higher value than the native grassland by 0.16 ± 0.25 ( figure 4(a)). The larger energy redistribution factor of the green fodder demonstrated a higher capability of energy exchange (sensible and latent heat fluxes) from the land surface to the atmosphere through sensible and latent heat fluxes [15,40], which further induced a lower LST in green fodder fields than in native grassland during the late growth stages ( figure 4(c)).
The quantified mechanism of LST variations between green fodder and native grassland was explored by the temperature response model. Since the magnitude of soil heat flux in equation (3) is negligible in this study, we ultimately compared the monthly ∆T s induced by the albedo change (i.e. radiative process alteration, ∆T s_alb ) with that caused by the energy redistribution factor change (i.e. nonradiative process alteration, ∆T s_f ) during the growing season ( figure 4(c)). The radiative process alteration changed temperature slightly since the ∆T s_alb ranged from −0.23 ± 0.53 • C (in September) to 0.03 ± 0.41 • C (in June), while the magnitude of the ∆T s_f had a fluctuating trend from May to September, which showed positive values during the early growth stages and negative values during late growth stages. In August, the nonradiative process alteration cooled the land surface most dramatically by −1.22 ± 1.72 • C. A much larger magnitude of the ∆T s_f than that of ∆T s_alb, especially during the late growth stage, indicated that the nonradiative process alteration dominated the effects of land cover change on LST.
To verify the performance of the model, we further estimated the similarity between ∆T s derived from remote sensing data (∆T s_r ) and calculated by the temperature response model (∆T s_cal ; figure 4(d)). Based on the statistical test of consistency between the calculated and observed ∆T s , the temperature response model can explain the changes in LST accurately in the northeastern Tibetan Plateau (R 2 ranged from 0.97 to 0.98 during the growing season).
Climatic fragile regions like the Tibetan Plateau are more sensitive to the climate feedback of anthropogenic land cover changes, which means the effects of LST variation may be amplified here [41]. Furthermore, the biophysical effects of the increase in green fodder area that we found in this study are not only relevant for climate fragile areas but are also highly informative for other regions with developed livestock industries, such as the United States, where alfalfa, a major green forage crop, is the third most valuable field crop [42].
Under a global warming background, more scientific attention has been paid to the warming effects caused by biogeochemical processes within the Earth ′ s systems, such as carbon cycling [43]. However, anthropogenic land use changes also significantly alter land-atmosphere interaction via biophysical effects [44], which are vital but not yet fully understood. The cooling effects of increasing green fodder area that we found in this study demonstrated that the biophysical effects of anthropogenic land use change on the climate could be found to be opposite to the warming climate trend, which enhances the necessity of including biophysical effects of anthropogenic land use change in the whole Earth system in future research. Additionally, the larger magnitude of LST variation via energy redistribution illustrates the importance of the nonradiative process (which is related to energy dissipation processes such as evapotranspiration) in this cooling effect. Since the green fodder fields in this study are rainfed, the magnitude of climate feedback via nonradiative processes is potentially larger considering agricultural management like irrigation [45], suggesting the necessity of considering agriculture-induced land cover changes in future anthropogenic climate change research.

Conclusion
This study examined the direction, magnitude, and mechanism of changes in LST due to increased green fodder area on native grassland in the northeastern Tibetan Plateau. Specifically, we provided remote sensing observed evidence of the potential cooling effects of increased green fodder area on native grassland during the growing season, especially in the late growth stage due to differences in land surface phenology between green fodder and native grassland. Pronounced LST differences mainly occurred during the daytime. Further, the mechanism of LST variations between green fodder and native grassland was represented by a synthesis of the effects on radiative and nonradiative processes. Nonradiative process alteration, represented by changes in energy redistribution through sensible and latent heat fluxes, boosted the cooling effects of increased green fodder area in most of the growing season, while radiative process alteration, caused by a relatively higher albedo of green fodder, slightly reduced the LST of green fodder compared to native grassland. Our study suggests the necessity of considering agriculture-induced land cover changes in future anthropogenic climate change research, especially in climatic fragile regions. With the further availability of higher spatial and temporal resolution data, as well as the establishment of more refined modeling of temperature response analysis, climate feedback of anthropogenic land use change could be demonstrated and understood in a more accurate and detailed way.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.