Landsat-observed changes in forest cover and attribution analysis over Northern China from 1996‒2020

ABSTRACT

Forest dynamics provide important information on the ecological environment. The Three-North Shelter Forest Program (TNSFP) is one of the world’s largest reforestation/afforestation programs, however the actual changes in forest cover in the Three-North Regions (TNR) of China resulting from this program are highly uncertain. This study quantified changes in fractional forest cover (FFC) at 30 m using Landsat data from 1996 to 2020. Using the Google Earth Engine platform, more than 40,000 images from Landsat-5, Landsat 7 and Landsat-8 were integrated, and the annual surface reflectance was normalized based on the multi-band least squares regression and maximum normalized difference vegetation index composite method. An ensemble learning model trained using high-resolution Gao-Fen 2 satellite imagery was used to generate the FFC long time-series product. FFC showed an increasing trend with average rates of 0.022/10a in the last 25 years, and 0.03/10a after 2010 largely corresponding to the fourth and fifth phases of the TNSFP. There are significant regional differences in the relationship between FFC and air temperature ($R^2$ = 0.37) and precipitation ($R^2$ = 0.49). The increased air temperature in arid and less rainy areas inhibit the FFC increase, whereas the increase in precipitation had a promoting effect. FFC appeared more sensitive to changes in solar radiation and heat conditions in humid and rainy areas. The attribution analysis revealed that 34% of FFC changes were caused by climatic variables and 66% were caused by non-climatic factors. Among them, afforestation associated with the TNSFP significantly increased FFC, and forest fire is a key factor of forest change in the Greater Khingan Ranges and Lesser Khingan Ranges regions. Planting single tree species caused biological disasters in forests of Xinjiang and Inner Mongolia. Further analysis of the increased FFC using high-level satellite products demonstrated an improvement in environmental conditions with cooler land surface temperature and higher vegetation gross primary production over the TNR.

KEYWORDS:

• Forests
• Northern China
• cover change analysis
• climate variables
• driver contribution rates
• feedback effects

1. Introduction

Forests are critical components of the ecosystem as they play a vital role in regulating global carbon and water cycles, which have significant impacts on the global ecological environment and climate change (Bonan 2008; Cook-Patton et al. 2020; Doughty et al. 2015; Pan et al. 2011; Peng et al. 2014; Stephenson et al. 2014). Changes in forest the growth state and spatiotemporal patterns reflect the feedback and impact of climate factors (Bastin et al. 2019; Ehbrecht et al. 2021; Zellweger et al. 2020), natural factors (Han et al. 2019) and human activities (Liebmann et al. 2016; Sebald, Senf, and Seidl 2021) on ecosystems, including processes such as forest degradation and loss (Curtis et al. 2018; Matricardi et al. 2020; Zhao et al. 2021), forest restoration and recovery (Crouzeilles et al. 2016; Forzieri et al. 2022), and woodland regeneration (Wu et al. 2021). Fractional forest cover (FFC) is a key parameter that represents the percentage of forest canopy coverage in an entire pixel viewed from top to bottom (Jia et al. 2015). It provides a detailed view of the forest canopy closure and detailed stand conditions within the pixel scale and has the potential to reveal forest changes at the sub-pixel level (Hansen et al. 2013; Senf et al. 2020; Yang and Crews 2019), which cannot be obtained through traditional pixel-by-pixel classification methods (Yang and Crews 2019). Time series FFC plays an essential role in forest fire assessment (Yin et al. 2020), forest dynamics and stability monitoring (Amarnath, Babar, and Murthy 2017; Czerwinski, King, and Mitchell 2014), and forest change trend analysis (Sarif, Jeganathan, and Mondal 2017). Despite its importance, there are relatively few studies on forest change and driving forces based on FFC. Therefore, it is crucial to quantitatively analyze forest vegetation change and its driving forces from the perspective of FFC to better understand the changes in forest ecosystems and promote their sustainable management (De Keersmaecker et al. 2015; Liu et al. 2021).

The Three-North Regions of China (TNR) is not only among the most vulnerable regions in China, but also crucial for China’s ecological environment engineering (Suo and Cao 2021; Zhang et al. 2016). The TThree-North Shelter Forest Program (TNSFP), as the largest artificial forest project in the world, has made a positive impact on the protection and restoration of forest resources and ecological environment in the Three North Shelterbelt area. It has played an important role in improving the living conditions of local people and achieving sustainable economic and ecological development (Jia et al. 2015; Suo and Cao 2021; Zhang et al. 2016). Forest change analysis relies heavily on remote sensing data products, which commonly have a relatively coarse spatial resolution and short temporal range. Although these data products can provide insights into general forest change characteristics, they cannot capture the spatial details of changes. Many studies on the TNR use single datasets or small areas and mainly focus on parameters such as vegetation productivity (Zheng et al. 2019), forest biomass (Zhang and Liang 2014), vegetation index (Xu, Yang, and Chen 2015), leaf area index (Zhu et al. 2016) to characterize forest vegetation cover. Forests have unique physiological structure and ecological functions compared to other vegetation types, While the TNSFP covers a vast area over a long period. Given the fragmentation of the landscape and the complexity of vegetation types (Zhang et al. 2022), it is crucial to detect and characterize the TNR forest changes and conduct quantitative analysis using medium-high resolution FFC data products (Gao et al. 2021; Li et al. 2020). Liu et al. (2021) developed an estimation FFC model for Landsat-8 OLI remote-sensing data. Since the OLI sensor has good inheritance from earlier Landsat-5/7 sensors, this provides a theoretical basis for applying the model to other Landsat data. However, different sensors often have inconsistent radiometric information due to varying spectral response and data post-processing methods. Therefore, radiometric normalization of Landsat data is necessary for enhancing the completeness and accuracy of remote sensing information on forest dynamics and enable consistent long-term earth observation.

Forest vegetation is highly sensitive to climate change and external interference (Bastin et al. 2019; Collalti et al. 2020; Mildrexler et al. 2016; Shestakova et al. 2016). Increasing temperatures have the potential to exacerbate the impact of severe drought events, and threats on forest survival (D’Orangeville et al. 2018; Mildrexler et al. 2016), and the abundance of precipitation affects the survival rate of afforestation (Doughty et al. 2015; Wu et al. 2021). Among the non-climatic factors affecting forests, land use-induced forest cover change directly leads to forest ecosystem degradation, This degradation greatly affects the carbon sink capacity of forests. Ecological protection projects such as afforestation support the restoration of forest vegetation and the increase of forest cover area, which improves the productivity of forest ecosystems, increases forest biomass and enhances the carbon storage capacity of forests (Curtis et al. 2018; Danneyrolles et al. 2019; Domke et al. 2020; Tuffour-Mills, Antwi-Agyei, and Addo-Fordjour 2020). Non-climatic natural factors such as forest fires not only destroy the structure of forest ecosystems, but also affect forest biomass, productivity, and biogeochemical cycles (Curtis et al. 2018; Doughty et al. 2015). Forest disturbance and destruction due to pests and diseases caused by factors like climate change or tree species monoculture are also major factors in forest ecosystem degradation. Minor disturbances may not cause significant changes in forest cover, however changes in stand density affect biomass and productivity (Bearup et al. 2014; Coggins et al. 2013).

Correctly distinguishing and evaluating the impact of different factors on forests, and understanding the response between climate-forest, human-forest, and natural ecology-forest are important prerequisites for formulating sustainable ecosystem development policies (Los et al. 2021; Trumbore, Brando, and Hartmann 2015; Wu et al. 2021); however, the research on TNR forest change and impact factors is still relatively scarce. In the study and analysis of temporal remote sensing data, the non-parametric approach of Sen Slope estimator combined with Mann-Kendall (Sen-MK), which has no requirements on the distribution of serial data, insensitivity to outliers and the absence of assumptions, enables, among others, slope estimation and change significance testing of variables on time series. The method was first widely applied in hydro-meteorology (Fathian et al. 2016; Rauf et al. 2016), and was subsequently gradually applied to spatial and temporal variation and trend analyses of evapotranspiration (Peng et al. 2017), temporal variation of vegetation index (Liang, Liu, and Hashimoto 2020) and differences in vegetation cover change (Hill and Guerschman 2022). Some scholars have explored the contribution rate of driving factors (Gadedjisso-Tossou, Adjegan, and Kablan 2021; Leroux et al. 2017; Pei et al. 2021). The residual trend approach (Evans and Geerken 2004) was used to quantitatively assess the impact of human and climatic factors on vegetation changes based on the assumption that vegetation changes are caused by climatic factors and anthropogenic drivers. Shi et al. (2021) used NDVI to represent vegetation coverage, and precipitation and reference evapotranspiration were identified as key climatic factors. A residual analysis model was established based on the multivariate linear regression method to study the changes in vegetation coverage in the Loess Plateau of China. Gao et al. (2022) A residual analysis model was constructed based on pixel-by-pixel regression analysis to study the contribution rate of drivers of vegetation changes (NDVI) in the Mu Us Sandland of China. Huang and Kong (2016) Used NDVI to represent the vegetation coverage of TNR, and the residual analysis model based on machine learning was combined with correlation analysis to study the contribution rate of human and climatic factors to land degradation in TNR. However, in the context of complex and variable natural environments and anthropogenic activities, analyses and evaluation of FFC changes in the TNR and its driving factors have been relatively rare. To better understand the dynamic process and pattern characteristics of temporal and spatial changes in forest vegetation, and strengthen the understanding of the interaction between climate and human activities on FFC changes, it is necessary to quantitatively evaluate the driving mechanism and contribution rate of driving factors to forest changes in the region.

In response to the current relatively few studies on forest vegetation using continuous variable FFC, and most of the studies on forest change analysis of TNR focus on local areas. In this study, the long time series FFC products based on remote sensing are used to study the forest vegetation of TNR. Specifically, we aimed to answer the following questions: (1) How does forest cover change in the TNR? (2) What are the main factors that drive forest change? (3) What are the consequences of forest cover change on the ecological environment of the TNR?

2. Materials

2.1. Landsat data

Landsat-5 TM, Landsat-7 ETM+ and Landsat-8 OLI have good inheritance, similar imaging ranges and spatial resolution (30 m) in the visible light band, although their spectral response functions differ, however they are generally considered equivalent in research (Flood 2014; Lück and van Niekerk 2016). Since the Google earth engine (GEE) platform provides data for each processing level of the Landsat series, which enhanced the study’s accuracy and reliability. This study utilized the surface reflectance data (T1_SR) of approximately 40,000 images of surface reflectance data (T1_SR) from TM (1996–2011), ETM+(2012) and OLI (2013–2020) covering TNR from 1996 to 2020, Landsat Path/Row number 113/26–151/33 were included, and cloud (including cloud shadow or snow) was removed using Landsat cloud mask data and the quality assessment band (Gorelick et al. 2017; Kumar et al. 2022; Xie et al. 2019).

2.2. The global LAnd surface satellite products

This study utilized several products from the Global LAnd Surface Satellite (GLASS) products suite (Liang et al. 2021), including land-surface temperature (LST) (Li et al. 2021), short-wave daytime land-surface albedo (LSA) (Feng et al. 2016; Liu et al. 2013; Qu et al. 2016), and gross pirmary productivity (GPP) (Zheng et al. 2020). These data offer a a spatial resolution of 1 KM and demonstrate strong spatiotemporal continuity. Since GLASS data started in 2000, this study selected the data covering TNR in the time range of 2000–2020. These products are widely used in various research fields, and their high-quality documentation ensures reliable data quality for analysis.

2.3. Topographical, climatic, and statistical data

Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model (ASTER GDEM) (http://www.gscloud.cn/sources/accessdata)) data were selected as the data source for the terrain factor. The spatial resolution of the data was comparable to the spatial scale of Landsat pixels (1 rad, approximately, about 30 m) (Hermance, Augustine, and Derner 2015; Satgé et al. 2015). Terrain analysis was performed on the data, and the corresponding slope, aspect, and elevation data were calculated (Connette et al. 2016; Peng et al. 2019).

The temperature and precipitation data were obtained from the National Qinghai-Tibet Plateau Scientific Data Center (https://www.tpdc.ac.cn/zh-hans/). The original data set has a spatial resolution of 1 km, covering the land area of China from 1901 to 2020 (Peng et al. 2017, 2018, 2019). This study extracted the data covering TNR from 1996 to 2020, which was then aggregated into annual precipitation and mean temperature datasets. In addition, this study utilize data from the China Statistical Yearbook of the National Bureau of Statistics (http://www.stats.gov.cn/tjsj/ndsj/) and the China Forestry Statistical Yearbook of the China Forestry Information Network (http://www.lknet.ac.cn/), which provide statistics on the annual afforestation area of each province in northern China, as well as the number of fires that occur and the forest area affected by fire and pests.

3. Methods

This study developed and applied a main process of evaluating and analyzing TNR forest changes which is shown in Figure 1. Firstly, time-series annual synthetic remote sensing data were generated using Landsat remote sensing data on GEE platform by employing the M-OLS and MVC method. Next, long time-series FFC datasets were constructed based-on machine learning models. Secondly, meteorological, statistical, and remote sensing data products were jointly used to assess and analyze the spatio-temporal variability characteristics, the contribution of drivers of change and feedback effects of TNR forests using Sen-MK, CV, VD, residual trend analysis, and other methods. Finally, the ecological effect assessment of TNR forests and ecological engineering management are discussed.

Figure 1. Flowchart of the analysis of forest cover changes and driving factors in the three-north of China.

3.1. Long time series FFC dataset

This study used GEE to generate long-term FFC datasets according to the following process: (1) To improve the spatial and temporal representativeness of the data and reduce data discrepancy caused by non-sensor factors, image with a sensor imaging time difference or less than 2 days and the containing different surface types were selected for pre-processing. Surface reflectance image pairs from Landsat-7 ETM + and Landsat-8 OLI were uesd, and only good quality pixels were retained. (2) A multi-band joint least squares (M-OLS) approach was used to construct a conversion model (Equation (1)(1) $OL{I}_i=\sum _{i=1}^n{\rho }_i\times ETM{+}_i+{\theta }_i$(1) ) from ETM+ to OLI. This approach achieves consistent conversion of the reflectance between the two sensors while considering the interaction between bands. Although there is no time overlapping imaging between TM and OLI, previous research suggests that the visible light bands of TM and ETM+are equivalent (Flood 2014), so the conversion model built by ETM+and OLI can be used to convert between TM and OLI. (3) To reduce the amount of data processing and improve the processing efficiency, the consistent converted reflectance data were stacked by year to form a year-by-year time series remote sensing image collection; Among them, due to the lower NDVI values of abnormal pixels such as clouds or fog, NDVI is used as an indicator to remove noise elements. The maximum normalized difference vegetation index (NDVI) composite (MVC) was adopted to remove clouds, shadows, or edge pixels in time series Landsat data, and synthesize the annual stacked data into the annual composite data (Myroniuk et al. 2020; Piao et al. 2021). (4) In our previous study (Liu et al. 2021), we developed the FFC estimation model for Landsat 8 imagery based on an ensemble learning algorithm using Gaofen-2 data as a reference, but we did not produce a long-time series FFC data over the entire TNR because of the massive data volume. This study applies this model to produce the TNR long time series FFC dataset using more than 40,000 Landsat images. Composite data and other feature parameters are combined into a data cube consistent with the model input to generate a time series FFC dataset.

(1) $OL{I}_i=\sum _{i=1}^n{\rho }_i\times ETM{+}_i+{\theta }_i$(1)

where:${\rho }_i$ is the conversion coefficient, $ETM{+}_i$ denotes the ith visible band of ETM+ (i = 1,2,3,4,5,6), $OL{I}_i$ denotes the band of the OLI sensor corresponding to the $ETM{+}_i$ band, and θ is a constant term.

3.2. FFC change trend detection and analysis

3.2.1. Sen-Mann-Kendall trend test

The Theil-Sen slope estimator is an nonparametric trend estimate method widely utilized in remote sensing and ecological studies. Unlike paremetric methods, it does not assume any particular data distribution and can provide accurate regression slopes even in the presence of outliers or non-normality. The Mann-Kendall (MK) is a non-parametric test, also known as a distribution-free test. It is insensitive to outliers, and does not require time series data to follow a certain distribution (Pei et al. 2021; Yang et al. 2019; Zhang et al. 2016). Moreover, the Sen-MK trend test allows us to statistically test whether there is a monotonically increasing or decreasing trend within a time-series FFC dataset, providing valuable insight into long-term environmental changes.

3.2.2. Quantification of forest cover change

The coefficient of variation (CV) is a statistical measure that quantifies the ratio of the standard deviation to the mean of a time series data, reflecting the degree of variability or stability of the data over time. In this study, the CV is utilized to characterize the fluctuation level of FFC. Vegetation Dynamics (VD) was used to evaluate the extent of changes in vegetation area between the end and base period of the study. To quantify the degree of change in forests at different time intervals, the FFC is divided into 11 gradients, starting from 0% and with a gradient interval of 10%. The characteristics of forest changes for each FFC gradient are analyzed. Additionally, the land use transition matrix was employed to describe the transition scale of land use types during different periods. It indicates the numbers of different land-use types that remained constant and changed during the study period. In this study, forests with different FFC levels were regarded as different types. According to the 20% interval, FFC was divided into low (L; FFC: 0–20%), medium low (M-L; FFC: 20–40%), medium (M; FFC: 40–60%), medium high (M-H; FFC: 60–80%) and high (H; FFC: 80–100%) (Liu et al. 2018; Zhang and Liang 2014).

3.3. Attribution analysis of FFC changes

3.3.1. Correlation analysis

We assessed the degree of influence of temperature and precipitation on forest cover change by calculating the Pearson correlation coefficient (CC) and partial correlation coefficient (PCC) of FFC with temperature and precipitation. The CC reflects the degree of association between two variables; while the PCC enables the correlation analysis of two variables (FFC and temperature or FFC and precipitation), while controlling for the effect of the other variable (precipitation or temperature), which can provide a better understanding of the relationship among the three (Shi et al. 2021; Zhao et al. 2020).

(2) $r=\frac{{\sum }_{i=1}^n\left(x_i-\bar{x}\right)\left(y_i-\bar{y}\right)}{\sqrt{{\sum }_{i=1}^n\left(x_i-\bar{x}\right)}\sqrt{{\sum }_{i=1}^n\left(y_i-\bar{y}\right)}}$(2)

(3) $r_{Ft,p}=\frac{R_{Ft}-R_{Fp}\ast R_{tp}}{\sqrt{\left(1-R_{Fp}^2\right)\ast \left(1-R_{tp}^2\right)}}$(3)

(4) $r_{Fp,t}=\frac{R_{Fp}-R_{Ft}\ast R_{tp}}{\sqrt{\left(1-R_{Ft}^2\right)\ast \left(1-R_{tp}^2\right)}}$(4)

Where:$\hspace{0.17em}{x}_i$ and $y_i$ represent the values of the two variables $x$ and $y$ involved in the correlation analysis in the ith year; $\bar{x}$ and $\bar{y}$ are the mean values of the variables x and y, respectively; R is the CC of the two elements; $R_{Ft,p}$ represents the PCC between FFC and temperature; $R_{Fp,t}$ represents the PCC between FFC and precipitation; $R_{Ft}$ represents the CC between FFC and temperature; $R_{Fp}$ represents the CC between FFC and P; $R_{tp}$ represents the CC between temperature and precipitation.

3.3.2. The analysis of residual trend

The residual trend method, proposed by Evans and Geerken in 2004, is a widely used approach for analyzing vegetation changes. It is based on the assumption that changes in vegetation cover are closely related to climatic factors (primarily precipitation and temperature) and that the effects of non-climatic factors, such as human activities, can be represented by model residuals – by unexplained changes between vegetation cover indicators and climatic factors (Evans and Geerken 2004; Shi et al. 2021; Wang et al. 2015). In this study, we utilized residual trend method to investigate forest vegetation changes. In the first step, to establish a model that describes the relationship between forest vegetation and climate factors, based on this model, the estimated value of FFC under ideal conditions can be calculated from climate data; however, due to the influence of other factors, there may be differences between the estimated and the actual observed value of FFC. Therefore, in the second step, we calculated the FFC residual, which represents the difference between the estimated and actual values of FFC. In the third step, we conducted a trend analysis on the residuals. The trend of residuals is considered to be unrelated to climate, mainly reflects changes in forest vegetation caused by non-climate factors. By jointly analyzing the trends of observed, estimated and residual values, it is possible to identify the contribution of both climatic and non-climatic factors to forest vegetation change. A positive trend of the residuals indicates that non-climatic factors contribute to the improvement of forest vegetation, while the opposite indicates that non-climatic factors lead to the degradation of forest vegetation (Burrell, Evans, and Liu 2017; Wessels, van den Bergh, and Scholes 2012). However, there are some limitations in the conventional use of this method: Due to the complex response mechanisms between climate factors and vegetation, linear regression cannot accurately model the relationship between climate factors and forest vegetation (Huang and Kong 2016; Shi et al. 2020).

This study is grounded on the assumption that the changes in FFC are impacted by both climatic and non-climatic factors. Firstly, we analyzed the temporal variation trend of temperature, precipitation, and FFC data from 1996 to 2020 and select the pixels with high stability. During this process, we constructed a dataset containing 10,094 samples with temperature and precipitation as predictors and FFC as response variables. Then, considering that random forest (RF) is a machine learning method for multi-decision integration and has obvious advantages in fitting the nonlinear relationships among the covariates (Zhang and Yang 2020), we selected RF to construct an estimation model of FFC and temperature and precipitation. Thus, 70% of the sample data was used as the training set, and 30% of the data is used as the validation set. RMSE served as the evaluation index of the model, which attained a final validation accuracy of RMSE = 0.09. Finally, the value of FFC (C_FFC) under the influence of climate factors is calculated by using the constructed model combined with the temperature and precipitation of TNR. The reference FFC value $\left(\mathrm{O}\mathrm{\_}\mathrm{F}\mathrm{F}\mathrm{C}\right)$ minus the FFC value under the influence of climate change$\mathrm{C}\mathrm{\_}\mathrm{F}\mathrm{F}\mathrm{C}$, the obtained residual ($\mathrm{R}\mathrm{\_}\mathrm{F}\mathrm{F}\mathrm{C}$) can be considered as the influence of non-climate factors including human factors on the FFC ($\mathrm{H}\mathrm{\_}\mathrm{F}\mathrm{F}\mathrm{C}$):

(5) $\mathrm{H}\mathrm{\_}\mathrm{F}\mathrm{F}\mathrm{C}=\mathrm{O}\mathrm{\_}\mathrm{F}\mathrm{F}\mathrm{C}-\mathrm{C}\mathrm{\_}\mathrm{F}\mathrm{F}\mathrm{C}$(5)

To analyze whether the influence of non-climatic factors on the FFC change was positive or negative, this study conducted a pixel-by-pixel trend analysis of the obtained residual sequence, and obtained the spatial distribution pattern of the residual change trend. The trend can reflect the influence of non-climatic factors on the FFC variation. Combined with the FFC change trend $\mathrm{O}\mathrm{\_}\mathrm{F}\mathrm{F}\mathrm{C}$, under the condition that the FFC change trend was significant, the contribution rates of climatic and non-climatic factors to FFC change could be calculated by referring to Equations (6)(6) $\mathrm{C}\mathrm{c}=\frac{\mathrm{S}\mathrm{l}\mathrm{o}\mathrm{p}{\mathrm{e}}_{\mathrm{C}\mathrm{\_}\mathrm{F}\mathrm{F}\mathrm{C}}}{\mathrm{S}\mathrm{l}\mathrm{o}\mathrm{p}{\mathrm{e}}_{\mathrm{O}\mathrm{\_}\mathrm{F}\mathrm{F}\mathrm{C}}}\times 100\mathrm{\%}$(6) and (7(7) $\mathrm{C}\mathrm{h}=\frac{\mathrm{S}\mathrm{l}\mathrm{o}\mathrm{p}{\mathrm{e}}_{\mathrm{H}\mathrm{\_}\mathrm{F}\mathrm{F}\mathrm{C}}}{\mathrm{S}\mathrm{l}\mathrm{o}\mathrm{p}{\mathrm{e}}_{\mathrm{O}\mathrm{\_}\mathrm{F}\mathrm{F}\mathrm{C}}}\times 100\mathrm{\%}$(7) ) (denoted by Cc and Ch, respectively) (Shi et al. 2021; Yang et al. 2017).

(6) $\mathrm{C}\mathrm{c}=\frac{\mathrm{S}\mathrm{l}\mathrm{o}\mathrm{p}{\mathrm{e}}_{\mathrm{C}\mathrm{\_}\mathrm{F}\mathrm{F}\mathrm{C}}}{\mathrm{S}\mathrm{l}\mathrm{o}\mathrm{p}{\mathrm{e}}_{\mathrm{O}\mathrm{\_}\mathrm{F}\mathrm{F}\mathrm{C}}}\times 100\mathrm{\%}$(6)

(7) $\mathrm{C}\mathrm{h}=\frac{\mathrm{S}\mathrm{l}\mathrm{o}\mathrm{p}{\mathrm{e}}_{\mathrm{H}\mathrm{\_}\mathrm{F}\mathrm{F}\mathrm{C}}}{\mathrm{S}\mathrm{l}\mathrm{o}\mathrm{p}{\mathrm{e}}_{\mathrm{O}\mathrm{\_}\mathrm{F}\mathrm{F}\mathrm{C}}}\times 100\mathrm{\%}$(7)

where, $\mathrm{S}\mathrm{l}\mathrm{o}\mathrm{p}\mathrm{e}$ represents the slope of the corresponding time series data with the change of the year, the Sen method is used to calculate the slope, which can improve the robustness and reduce the influence of outliers. Based on the significance of the $\mathrm{O}\mathrm{\_}\mathrm{F}\mathrm{F}\mathrm{C}$ value of FFC and equations (6)(6) $\mathrm{C}\mathrm{c}=\frac{\mathrm{S}\mathrm{l}\mathrm{o}\mathrm{p}{\mathrm{e}}_{\mathrm{C}\mathrm{\_}\mathrm{F}\mathrm{F}\mathrm{C}}}{\mathrm{S}\mathrm{l}\mathrm{o}\mathrm{p}{\mathrm{e}}_{\mathrm{O}\mathrm{\_}\mathrm{F}\mathrm{F}\mathrm{C}}}\times 100\mathrm{\%}$(6) and (7(7) $\mathrm{C}\mathrm{h}=\frac{\mathrm{S}\mathrm{l}\mathrm{o}\mathrm{p}{\mathrm{e}}_{\mathrm{H}\mathrm{\_}\mathrm{F}\mathrm{F}\mathrm{C}}}{\mathrm{S}\mathrm{l}\mathrm{o}\mathrm{p}{\mathrm{e}}_{\mathrm{O}\mathrm{\_}\mathrm{F}\mathrm{F}\mathrm{C}}}\times 100\mathrm{\%}$(7) ), the calculation method for the contribution rate of driving factors is shown in Table 1 (Leroux et al. 2017):

Table 1. Residual trend method to calculate the contribution rate to FFC changes.

$\mathrm{S}\mathrm{l}\mathrm{o}\mathrm{p}{\mathrm{e}}_{\mathrm{O}\mathrm{\_}\mathrm{F}\mathrm{F}\mathrm{C}}$	$\mathrm{S}\mathrm{l}\mathrm{o}\mathrm{p}{\mathrm{e}}_{\mathrm{C}\mathrm{\_}\mathrm{F}\mathrm{F}\mathrm{C}}$	$\mathrm{S}\mathrm{l}\mathrm{o}\mathrm{p}{\mathrm{e}}_{\mathrm{h}\mathrm{\_}\mathrm{F}\mathrm{F}\mathrm{C}}$	Climate Contribution%	Non-climate contribution%
significant increase (p<0.05)	>0	>0	Cc	Ch
	>0	<0	100	0
	<0	>0	0	100
significantly Decrease (p<0.05)	<0	<0	Cc	Ch
	<0	>0	100	0
	>0	<0	0	100

3.4. Evaluating impacts of FFC change

We implemented a spatial-temporal approach (Li et al. 2015; Peng et al. 2014; Yang, Huang, and Tang 2020) to compare the differences between LST, GPP, and land-surface albedo (LSA) of high-FFC forests and adjacent low-FFC forests to quantify the impact of forests on physical processes at the surface. First, to eliminate the effects of different climatic conditions and ensure that the selected pixels had similar climatic backgrounds, we established a square window with an initial edge length of 40 km, centered on a forest pixel with high-FFC. Second, we searched for forest pixels with low-FFC within this window. To mitigate the effect of elevation changes, pixels with an elevation difference >50 m from the central image element were excluded. We calculated the annual average LSA, LST, and GPP values. To account for the possible impact of human activities on the results, the data extracted in this study were mainly pixels with different FFC. The LST, GPP, and LSA differences ($\mathrm{\Delta }LST\mathrm{\Delta }LSA\mathrm{\Delta }GPP$) of forest pixels with different FFC were compared.

(8) $\mathrm{\Delta }LST=LS{T}_{FF{C}_{High}}-LS{T}_{FF{C}_{low}}$(8)

(9) $\mathrm{\Delta }LSA=LS{A}_{FF{C}_{High}}-LS{A}_{FF{C}_{low}}$(9)

(10) $\mathrm{\Delta }GPP=GP{P}_{FF{C}_{High}}-GP{P}_{FF{C}_{low}}$(10)

Taking LST as an example, a positive (negative) $\mathrm{\Delta }LST$ indicates a warming (cooling) effect of increasing forest cover, and the differences in $\mathrm{\Delta }LSA\mathrm{\Delta }GPP$ were similarly defined.

4. Result

4.1. Fractional forest cover spatiotemporal variations

4.1.1. Fractional forest cover mapping and time series trends

The OLI and ETM+ band conversion parameters (Table 2) of the M-OLS regression were obtained. This was applied to ETM+ to obtain the adjusted transformed reflectance data. The performance of the reflectance conversion model was evaluated using independent validation data.

Table 2. ETM+ and OLI multi-band least squares regression (M-OLS) conversion parameters.

Band name	Conversion factor/constant term
OLI_Blue	[0.1495, 0.4948, 0.0579, −0.0271, −0.0112, 0.0184]/0.0001
OLI_Green	[−0.5206, 1.1124, 0.1615, 0.0055, −0.0513, 0.0316]/0.0069
OLI_Red	[−0.5780, 0.4115, 0.9665, −0.0388, −0.0225, 0.0157]/0.008
OLI_Nir	[−0.4978, 0.2105, 0.1033, 0.9571, 0.1265, −0.1269]/0.0101
OLI_SWIR1	[−0.3482, 0.1683, 0.1363, −0.0382, 0.9591, −0.0606]/0.0169
OLI_SWIR2	[−0.3342, 0.1959, 0.1087, −0.045, 0.0127, 0.886]/0.0162

Figure 2 shows the reflectance scatter plot of the ETM+ and OLI raw bands used as the validation data. The horizontal axis of the plot represents the ETM+ band, and the vertical axis represents the OLI band. There was a difference between the bands corresponding to the two sensors, and the blue band had the weakest correlation with a slope of only 0.82. The red band had the strongest correlation. The RMSE values between the two sensor bands were between 0.016 and 0.026. The slopes of the blue band and SWIR1 deviated the most from the 1:1 line, with both < 0.9. Figure 3 shows the conversion result of the M-OLS model. The consistency between several bands showed a high improvement, the slope of the blue band also increased (to > 0.91), and the correlation (R²)between the bands was significantly improved (i.e. all were > 0.8) and their deviations reduced.

Figure 2. Comparison of raw reflectance between ETM+ (LE07) and OLI (LC08).

Figure 3. Comparison of ETM+ (C_LE07) and OLI (LC08) reflectance after coefficient adjustment.

Figure 4 shows the FFC dataset of the TNR from 1996 to 2020. The dataset was generated based on the GEE platform and had the characteristics of spatiotemporal continuity and high spatial resolution (30 m). Meanwhile, considering that in the definition of forest, there are different criteria for canopy cover, in orser to clarify the temporal characteristics of forest cover changes, we set 0, 10%, 20%, 30%, and 40% as gradient thresholds, respectively, obtained forest cover at different thresholds from 1996 to 2020, and performed linear fitting on the time-series forest cover (Figure 5). The results showed that the overall forest area of the TNR fluctuated and increased, with an average annual growth rate of 0.02/10a. Taking the statistical results of the 20% threshold as an example, the forest area in 2003 was the lowest, and the forest cover area in 2018 was the highest over the past 25 years. In Figure 5, 2003 is the first year with the least forest cover since 1996 when the threshold of FFC > 20% is used. Meanwhile, 2010 is the end of the fourth phase of the TNSFP, so the time series are divided into 1996–2003, 2003–2010, and 2010–2020 by using two years, 2003 and 2010, as the time nodes. There was a decreasing trend from 1996 to 2003, and the forest cover increased after 2003. Before 2010, there was an increasing trend in forest cover of 0.022/10a. After 2010, forest cover increased more rapidly, i.e. at 0.03/10a. According to the afforestation area statistics in each province (Tab S12), after 2006, the annual increase in the afforestation area in each province was the factor for the increase in the forest area in TNR after 2010.

Figure 4. The distribution of FFC in the TNR from 1996 to 2020.

Figure 5. Time series variation characteristics of FFC (10% is the gradient threshold to represent different fraction cover).

4.1.2. Fractional forest cover spatial distributions and changes

To analyze the relationship between terrain factors and forest spatial distribution, the overall characteristics of the terrain were obtained by calculating the elevation and slope, and the FFC map was overlaid with the elevation map to obtain the spatial characteristics of the forest distribution under different terrain conditions. Figure 6(a) shows that the distribution of the mean FFC of the TNR from 1996 to 2020 had obvious spatial heterogeneity. Figure 6(b) shows the distribution box plot of FFC at different altitudes. In the range of altitudes <200 m, the forest was dominated by low-FFC. At altitudes <800 m, FFC increased with increasing altitude. At altitudes >2000 m, FFC decreased with increasing altitude. According to statistics on the distribution of TNR forests (Figure 6(c)), approximately 2.5% of the forests were distributed at altitudes <200 m and approximately 6.5% were distributed at altitudes >2000 m.

Figure 6. Average spatial distribution of FFC from 1996 to 2020 (a), Boxplot of FFC distribution with elevation, representing the coverage distribution range (b), Characterizing the proportion of forest pixels at different elevations (c).

FFC of the TNR showed an overall increasing trend. As shown in Figure 7(a), the regions with significant increases were mainly distributed in the western part of the Tianshan Mountian (TSM), Loess Plateau (LP), and Greater Khingan Mountains (GKM), and there was a significant decreasing trend in the southern region of the Changbai Mountains (CBM) and Lesser Khingan Mountains (LKM). The FFC slope was divided into five categories: very significant increase (slope ≥0.005), significant increase (0.002 ≤ slope < 0.005), no significant change (−0.002 ≤ slope < 0.002), significant decrease (−0.002 ≤ slope < −0.005), or very significant decrease (slope < −0.005) (Figures 7b and 8(a,b)). The area with a very significant increase accounted for 55% (p < 0.05) of the total area. Whereas, approximately 29% showed no significant change. Approximately 15% of the forest cover showed a very significant downward trend (Table S3). Among them, GKM mainly showed a significant increasing trend (84%). However, there was also degradation of approximately 14% in the northern part of GKM, and significantly degraded areas of LKM and CBM accounted for > 30%. The LP and CBM are stable forests that occupy a large proportion. In general, FFC of these areas has mainly been increasing, and forests showing a downward trend exhibited discrete distribution characteristics. From Figure 8(c), it can be seen that approximately 64% of the forest cover in LP exhibited a higher fluctuation, and approximately 55% in TSM; GKM, LKM, and CBM areas were dominated by medium and low fluctuations (i.e. 50%, 57%, and 55%, respectively).

Figure 7. MK change trend of FFC.

Figure 8. Statistics of FFC change characteristics in TNR sub-region. Trend in sub-regions (a), Slope in sub-regions (b), CV in sub-regions (c).

4.1.3. Dynamic change characteristics of fractional forest cover

Based on FFC, Table 3 and Tab S1 show the degradation and restoration characteristics during different periods. From 1996 to 2000, the forests with FFC of 0, 10–20%, 20–30%, and 70–80% were positive, and the other FFC level of forests were negative, only approximately 2% of low-FFC level forests changed to a higher FFC, and the degraded proportions of M-L and M forests were lower than their restoration proportions. The degradation ratio of M-H and H forests was higher than that of restoration, indicating that forest degradation was more serious during this period. From 2000 to 2005, the VD of the low-FFC group was negative, and that of the high-FFC group was positive. Forest degradation of the M-L FFC, M FFC, and M-H FFC showed slowing trends compared with the previous period, and the proportion of restoration increased significantly. This showed that forest growth had improved and been restored to a certain extent, especially the forests with FFC in the ranges of 80‒90% and 90‒100%. From 2005 to 2010 and 2010 to 2015, the degradation trend slowed, the restoration trend accelerated, and the overall forest coverage increased. The VD from 2005 to 2010 showed a positive change in M and H FFC. From 2010 to 2015, the VD value of 0 showed a negative value, indicating that the area of non-forest land decreased and the area of newly added forestland increased, the decrease in M FFC and the increase in H FFC also indicated improved forest cover during this period. From 2015 to 2020, the overall increasing trend slowed, and the VD of forests with H FFC levels was negative, indicating that the degradation ratio of forests with M FFC levels and above increased during this period. The negative VD value for non-forest land indicated that continuous forestland addition had been taking place. Comprehensive analysis showed that there has been continuous new forest land in this area over the past 25 years. Forests with M FFC levels fluctuated greatly, whereas, forests with H FFC levels showed a trend of initially increasing and then degrading with time. Among them, 2000‒2015 was a period of increasing trends in TNR forest restoration, however, the forest restoration rate decreased after 2015.

Table 3. The dynamic degree of FFC change in the TNR.

	1996–2000	2000–2005	2005–2010	2010–2015	2015–2020
0	0.0024	−0.0009	0.0001	−0.0022	−0.0019
(0,10]	−0.0259	0.0187	−0.0395	0.0155	0.0572
(10,20]	0.0072	−0.0014	−0.0345	−0.0547	0.0147
(20,30]	0.0049	−0.0031	−0.0207	−0.0442	−0.0012
(30,40]	−0.0035	−0.0018	−0.0146	−0.0336	0.0068
(40,50]	−0.0149	−0.0071	0.0041	−0.0237	0.0047
(50,60]	−0.0221	−0.0046	0.0202	−0.0269	−0.0001
(60,70]	−0.0031	−0.0109	0.0359	−0.0316	0.0094
(70,80]	0.011	−0.0008	0.0294	0.0039	0.012
(80,90]	−0.0406	0.0371	0.0062	0.119	−0.0016
(90,100]	−0.0875	0.0463	0.0054	0.6083	−0.0283

The gradient trend of the time series can represent the temporal characteristics of forest changes, and when combined with statistical data, the reasons for forest changes can be explored and analyzed. FFC of LP showed an overall upward trend (Figure 9a), and the forests with L FFC (0‒30%) showed a trend of initially increasing and then decreasing, forests with the FFC of 30‒50% fluctuated greatly, increasing from 2000 to 2010 and decreasing after 2010. Forests with M and H FFC (50‒100%) showed overall upward trends, with a higher rate of increase after 2010 than before 2010, indicating that the forest vegetation in this area greatly improved after 2010. The inter-annual slope of FFC in the CBM region over the past 25 years was from 0.02% to 0.2%, and the overall upward trend was not obvious (Figure 9b). Forests with an FFC levels of 40‒50% showed a significant decline after 2010, and according to the forest fire statistics in this area (Tab S13 Tab S14), the area of forests affected by fires after 2010 was significantly higher than that before 2010. Combined with the forest changes, it can be considered that fires in this area have had a greater impact on forests with M and L FFC. According to Figure 9(c), forest cover in the LKM exhibited a relatively gentle upward trend. Combined with the statistics (Tab S13 Tab S14) of forest change data in the region, fire was the main factor for fluctuations in the region in 2003. Forests with M and L FFC showed greater instability, forest change with an FFC levels of 60‒80% had obvious time nodes, whereby it showed an upward trend before 2010 and subequently decreased. Forests with H FFC (80‒100%) showed an upward trend. The above analyses showed that forests with H FFC in this area had good stability, whereas, forests with M and L FFC were easily disturbed and fluctuated greatly. The GKM analysis showed an overall upward trend. Figure 9(d) shows that large fluctuations occurred from 1998 to 2003. Combined with the forest change statistics (Tab S13 Tab S14), the higher number of fires and burned areas in the region in 2003 were the main disturbance factors. Since 2015, Inner Mongolia, Jilin, and other regions have gradually stopped commercial logging of natural forests, and the increasing trend of forests with different FFC level grades in the region after 2015 shows that this measure has had a positive effect on improving forest coverage. Overall, the TSM forests showed a fluctuating and increasing trend (Figure 9e); however, there was a significant decrease in 2005. Combining the data (Tab S12, Tab S13, Tab S14, Tab S15), it can be seen that the number of fires in this area increased since 2004, and that 2008‒2016 was a period of high incidence of insect pests in this area. The area maintained a relatively high afforestation level after 2010, the number of fires gradually decreased after 2012, and the area affected by biological disasters gradually decreased after 2016. All of these factors have been beneficial to forest restoration in this area and have contributed to the increase in the forest area.

Figure 9. Temporal variation characteristics of FFC at different levels in sub-regions.

4.2. Relative contributions of climate and anthropogenic activities

4.2.1. Correlation analysis of fractional forest cover with temperature and precipitation

The spatial distribution of the annual average precipitation and air temperature of the TNR over the past 25 years is shown in Figure 10(a,c). We overlayed FFC with annual precipitation and air temperature maps to analyze the relationship between FFC, precipitation, and air temperature. Figure 10(b) shows a boxplot of the relationship between precipitation and FFC distribution. The forests with L FFC levels were dominant in areas with precipitation <370 mm. Forests with precipitation between 370 and 410 mm showed an increasing trend toward high FFC, and were dominated by M and H FFC in areas with precipitation >410 mm. According to statistics, 30% of the forest pixels were distributed below the 370 mm precipitation line and 70% of the forest pixels were above the 370 mm precipitation line. Figure 10(f) shows a boxplot of the relationship between air temperature and the FFC distribution. It can be concluded that the forest was sparsely distributed in the area where the air temperature was < −5°C. When the temperature was in the range of −5 to 0°C, combined with Figure 6, the forest was mainly dominated by M and H FFC. With the increase in temperature, the distribution of forests tended to be stable, whereas, the FFC levels of forests decreased in areas with temperatures > 10°C. Based on the spatial distribution patterns of temperature and precipitation and their distribution law with FFC, it can be concluded that the growth and distribution of forest vegetation were inseparable from water and heat conditions.

Figure 10. Spatial pattern of annual precipitation and temperature in TNR. (a) Mean annual precipitation, (b) The relationship between FFC and precipitation distribution, (c) Mean annual temperature, (d) The relationship between FFC and air temperature distribution.

Partial correlation (Pc) analysis was used to explore the independent correlations between precipitation, temperature, and FFC (Figure 11). Among them, the Pc was positive between FFC and precipitation in the forest cover area of the LP, the northern and eastern areas of the GKM, the central and northern areas of the CBM, and the TSM, indicating that precipitation had a promoting effect on FFC in these areas when temperature was not considered. FFCs in the central and southern LKM, southern CBM, and southeastern GKM regions showed negative partial correlations with precipitation, indicating that precipitation negatively affected forest growth in these regions. Figure 11(b) shows the distribution of the Pc coefficient between FFC and temperature without considering the precipitation factor. 70.7% of the forest cover areas showed a positive correlation. These areas were mainly distributed in the east of LP, north of GKM, west and south, east of LKM, south of CBM Region, and in the TSM region. Approximately 28.5% of the forest cover showed a negative correlation (Tab S9 and Tab S10), distributed in the western part of the LP, in the middle of the GKM, and in the border area between the LKM and the CBM. The Pc between precipitation and FFC was higher than that between temperature and FFC. However, forests had different responses to precipitation and temperature due to different hydrothermal conditions in different regions. In LP, there was aridity and minimal rainfall and the hydrothermal conditions were sufficient. Precipitation could promote forest vegetation growth in this area, and the temperature increase had an inhibitory effect on FFC. In the southern part of CBM, it was humid and rainy, and the change in hydrothermal conditions affected vegetation growth by affecting plant photosynthesis, thus affecting the FFC change. However, in the TSM, the hydrothermal and water vapor conditions were lower than in other regions. Therefore, precipitation and temperature could simultaneously promote forest vegetation growth in the TSM region.

Figure 11. Partial correlation between FFC and precipitation (a) and temperature (b).

4.2.2. Impact of climate and anthropogenic activities on fractional forest cover

Figure 12(a) shows the potential FFC (C_FFC) variation trend under the influence of climate simulated by the model established between climate factors and FFC. Figure 12(b) shows the residual (R_FFC) between O_FFC and C_FFC, which can be considered as the impact of non-climatic factors on FFC, including human activities (logging, planting, etc.), biological disasters (diseases, pests, species invasion, etc.), and forest fires. Calculation of the change trends of C_FFC and R_FFC showed a significant difference between C_FFC and O_FFC: > 56% of the areas showed a very significant decline, including the northern and southern parts of GKM, most of the CBM, and the southern part of LP,only 41% of the areas showed a significant increase (LKM and TSM), R_FFC showed an upward trend overall, indicating that non-climatic factors played a positive role in the increase of TNR forests, and most areas of LP and GKM showed the most obvious trend, and most areas of the LKM showed a downward trend. Owing to the high population density and intensive human activities in the areas where the mountains and plains meet in the southern part of the GKM and the northern part of CBM, the forests in these areas are more vulnerable to human activities, especially urban expansion, and the increase in agricultural area has had negative effects on forest change.

Figure 12. Variation trend of FFC under the influence of climate (a) and non-climate (b) from 1996 to 2020.

Figure 13 shows the spatial distribution of the contribution rates of the influencing factors obtained using the contribution rate estimation model. The main area where climatic factors were influential was in the northeast (Figure 13(a,c,e)), including the southern part of the GKM, LKM, and part of the CBM, which mainly reflected a positive effect, whereas, the northern part of the GKM and the southern part of the CBM showed a negative effect with discrete distribution characteristics. The non-climatic factors mainly affected a wider area (Figure 13(b,d,f)), primarily in the central part of the GKM, LKM, and Yanshan to the LP and other areas. According to statistics (Figure 14), the relative contribution rate of climatic factors to FFC changes was 34%, of which the proportion of positive effects was 65% and that of negative effects was 35%. Overall, climate factors positively affected changes in FFC, and the relative contribution rate of non-climatic factors was 66%, of which the positive and negative effects were 74% and 26%, respectively, which also had a positive effect on the change in FFC. Based on the contribution rate analysis, it can be concluded that both climatic and non-climatic factors have had a positive effect on FFC, and that non-climatic factors have played a leading role in the change in FFC in the past 25 years.

Figure 13. Distribution of contribution rate. climate factor (a), Non-climate factor (b), Positive effect of climate factor (c), Negative effect of climate factor (d), Positive effect of non-climate factor (e), Negative effect of non-climate factor (f).

Figure 14. Statistics of the contribution rate of influencing factors in the change of FFC.

The periods from 2001 to 2010 and 2011 to 2020 were the fourth and fifth phases of the TNSFP project, respectively. Taking 2010 as the time node, the changes in FFC in the two periods of 1996‒2010 and 2010‒2020 were analyzed. As shown in Figure 15, the growth trends of the northern GKM and LP were significantly higher than that of other regions from 1996 to 2010, and that after 2010, the growth trend of FFC of the TNR began to decline. Except for the southern part of the GKM, the growth areas showed a discrete distribution. As shown in Figure 15(b,e), the influence of non-climatic factors on FFC of the two periods was mainly positive feedback, especially before 2010, most of the forests showed an increasing trend, however the northern part of the GKM showed a downward trend after 2010. Combined with Figure 15(c), it can be concluded that the downward trend was mainly caused by non-climatic factors. Overall FFC of the TNR showed an increasing trend in both periods, and the overall growth trend after 2010 was higher than that before 2010. However, it was also found that the growth trend of LP before 2010 was higher than that after 2010.

Figure 15. The changing trend of FFC in the TNR. (a) 1996–2010 O_FFC, (b) 1996–2010 C_FFC, (c) 1996–2010 R_FFC, (d) 2010–2020 O_FFC, (e) 2010–2020 C_FFC, (f) 2010–2020 R_FFC.

4.3. Consequences of FFC changes

We evaluated the environmental consequences of FFC changes. An increase in FFC can significantly reduce the surface albedo, which is a warming effect. Figure 16 shows a strong negative correlation between FFC and LSA in each region of the TNR, i.e. the higher the FFC, the lower the LSA. This was because, under normal circumstances, the brightness of surfaces such as grassland and bare soil is higher than that of forest, and the forest with low FFC level was affected by grass or bare soil in the background and had a higher albedo, theus, albedo decreased when FFC increased. Figure 16(f) shows the quartile map of the effect of FFC on LSA (ΔLSA) in different regions of the TNR. The black horizontal line in Figure 16(f) represents the mean value, the red dotted line represents the zero value line, and a value < 0 indicates that ΔLSA was negative. It can be concluded that the increase in FFC had a reducing effect on the albedo of the region overall, among which FFC of the TSM had the most obvious reduction effect on LSA (−0.104[0.077]), and the forest in the LP area had the lowest reduction effect on LSA (−0.015[0.013]).

Figure 16. Relationship between FFC and LSA.

Increases in FFC have a net cooling effect, largely due to the increased evapotranspiration despite the albedo warming effect. Figure 17 shows the relationship between FFC and annual average LST, which showed different degrees of negative correlation in different regions of the TNR, indicating a cooling effect with an increase in FFC. The effects of forests in the GKM, LKM, and CBM regions on LST were lower than those in the LP and TSM regions, i.e. the increase in FFC in the LP and TSM regions had a more obvious cooling effect. Figure 17(f) shows a quartile map of the effect of forests in different regions on LST (ΔLST). The forests in the LP area had the most obvious cooling effect. Although there are some warming samples in other areas, the average value was negative, therefore, the overall impact was a cooling effect. We further analyzed the seasonal temperature changes in these areas (Tab S11), and FFC had the most obvious cooling effect in summer, during which the ΔLST of the first three quarters of GKM, LKM, and CBM were all negative; however, these were positive in winter, indicating that the increase in FFC manifested as a warming effect in winter. The LP and TSM regions exhibited cooling effects.

Figure 17. Relationship between FFC and LST.

FFC is a useful indicator of vegetation productivity. Figure 18 shows the correlation between forest coverage and GPP, and that the increase in FFC had a positive effect on the increase in GPP. There were obvious positive correlations between FFC and GPP in CBM, GKM, LKM, and LP, and FFC in the TSM region was positively (but not obviously) correlated with GPP. Figure 18(f) shows a quartile map of the forest effect on GPP (ΔGPP). Overall, the increase in FFC had a positive effect on GPP; however, there were significant differences in different regions. The change in FFC in the LP region caused the largest change in GPP, whereas, in the GKM region, the increase in GPP caused by FFC was the smallest (ΔGPP = 0.168 gC/m³).

Figure 18. Relationship between FFC and GPP.

5. Discussion

This study utilized FFC data to investigate the forest changes in the TNR, which demonstrated the strengths of the FFC dataset in representing forest change information and providing clear insight into the characteristics, status, and trends of forest distribution in the TNR.

5.1. Discussion of impact factors

5.1.1. Climatic factors

In the TNR, temperature and precipitation have had significant impacts on FFC change over the past 25 years, with their effects exhibiting spatial heterogeneity (Peng et al. 2021). Specifically, in the midwest of the TNR, precipitation has had a stronger positive effect than temperature’s negative impact on FFC change. which can be attributed to the climatic conditions and geographical location of the TNR. The eastern part of the TNR is characterized by humid and rainy conditions, where forest growth has a higher demand for light. In contrast, the central and western parts of the TNR generally experience less precipitation but sufficient sunlight, leading to excess light and heat that result in soil evaporation exceeding precipitation levels. High temperatures accelerate plant transpiration and increase the respiration rate of vegetation, thereby accelerating dry matter consumption. Plants are subjected to water stress and experience physiological drought, which affects plant growth, resulting in reduced FFC level. Precipitation plays a crucial role in alleviating water stress and increasing soil moisture, making vegetation in these areas more sensitive to the precipitation level (Jiang et al. 2020; Zhang et al. 2021). In the eastern region, annual average precipitation is abundant; however, light and heat are lacking, resulting in higher sensitivity of forest vegetation to temperature in this region.

5.1.2. Anthropogenic drivers

In addition to climatic factors, anthropogenic activities are considered to be an important factor in VD (Wei et al. 2022). In the TNR, on one hand, through afforestation activities, the forest coverage rate increased from 5.05% in 1977 to 13.57% in 2020, and has continued to increase (Tab S12). Several ecological projects in northern China are of great significance for increasing the total amount of forest resources, controlling land desertification, reducing soil erosion, and increasing forest carbon sequestration. On the other hand, the implementation of ecological protection policies and the need for production development have affected the restrictive effect of natural geographical conditions on the change in the FFC pattern, i.e. flat, low altitude (<200 m), and adequate water and heat conditions suitable for forest growth areas. As a result of increased human activities, forest vegetation in these areas was sparse or even non-existent, and the rate of forest growth was higher in mountainous areas with steeper slopes than in areas with flat terrain. For example, in the LP area, the slopes of the main reforestation areas are > 15°, and the relative proportion of forest cover increased with the increase in slope. This is because as the human population has increased, urbanization and agricultural cultivation activities have expanded, and flat terrain areas have been used more as areas for human production-related activities. In areas with large slopes, the terrain fluctuates violently, soil erosion becomes more intense, and soil barrenness gradually deepens, making it increasingly unsuitable for agricultural production and cultivation, therefore, people occupy flat land for development. At the same time, high-slope areas are also fully utilized to carry out suitable afforestation and grass-planting projects, which can improve the vegetation coverage in these areas and further improve soil and water conditions, resulting in sparse forest coverage in low-altitude and flat areas, and a gradual increase in forest cover in high-altitude and large sloped areas.

5.1.3. Natural factors other than climatic factors

Forest degradation is driven by several factors, including aging trees whose protective capacity declines, particularly in the early stage of the shelterbelt project that is over 30 years old. In the TNR, drought, water shortage, and poor soil and water conditions prevail, causing instability and degradation of stand structure. These unfavorable conditions also pose challenges for afforestation efforts. The frequency and extent of forest fires, prevention and control of forest diseases, insects, and rodents, and the forest resources of the TNR were counted (Tab S12, Tab S13, Tab S14, Tab S15). TNR is an area with a high incidence of forest fires. Comparing the relationship between the temporal changes of FFC and forest fire, it was found that they had obvious temporal similarity, indicating that fire has been a major natural factor causing forest degradation or loss. Among them, the areas of forest fire damage in Heilongjiang and Inner Mongolia accounted for > 97%, and the four periods of 2003‒2004, 2007‒2008, 2011, and 2018 were periods of high incidence of forest fires in this region. Pests have also had a significant impact on FFC of the TNR, i.e. more than the impact of fire. The areas where pest disasters occurred were mainly concentrated in Inner Mongolia and Xinjiang regions, and have been due to the relatively single species of trees. For example, in Xinjiang region, spruce woodland accounted for 45% of the entire forest area. Dominance of a single tree species could easily cause the spread of pests and diseases, and the forest loss caused by factors such as forest fires could further increase the risk of pests and diseases. In the eastern part of the TNR, the area of disasters was also higher than in other areas because of the larger forest area. The area of biological disasters in Inner Mongolia’s forests peaked in 2014‒2015, and the Yanshan area had a relatively high area of biological forest disasters from 2011‒2017. In general, disturbance of the TNR due to biohazards has been one of the main natural factors for the reduction in FFC.

5.2. Feedback effects on forests

Changes in forest cover affect the climate by altering biogeochemical (e.g. atmospheric CO2 concentration) and biogeophysical (surface albedo, surface temperature, roughness length, and evapotranspiration) processes (Godinho et al. 2016; Peng et al. 2014; Zhang et al. 2021). LST has a strong coupling relationship with the energy radiation and thermodynamic properties of the Earth’s surface, and a higher sensitivity than air temperature to changes in vegetation density (Li et al. 2021; Yang, Huang, and Tang 2020), which were used to characterize the changes in climate caused by changes in forest cover (Peng et al. 2014). LSA determines the amount of solar radiation absorbed by the surface. The consumption of these energy as latent heat and sensible heat fluxes is controlled by vegetation activities and soil water conditions. Therefore, LSA is the main biogeophysical driving force that causes climate change. The change of forest vegetation cover plays a key role in climate change by affecting the energy and water balance (Alibakhshi et al. 2020; Godinho et al. 2016; He et al. 2018). GPP determines the initial material and energy of terrestrial ecosystems, and changes in forest cover directly or indirectly affect ecosystem productivity.

Analysis of the feedback effects of FFC changes in this study showed the same characteristics caused by changes in FFC in different regions of the TNR, such as the same trends in annual LSA, LST, and GPP. Due to differences in climatic conditions, latitude, and other factors, FFC changes have had different degrees of impact on LSA, LST, and GPP in different regions. FFC of GKM, LKM, and CBM was mainly characterized by warming in winter. This was because the radiative and non-radiative effects of afforestation on surface temperature differ at different latitudes. In winter, the proportion of non-radiative effects at high latitudes is small and the radiation effect is dominant, therefore, the effect of afforestation on the surface temperature is a warming effect. Enhanced evapotranspiration during the growing season results in stronger non-radiative effects and, therefore, stronger cooling effects. Forest vegetation productivity is also related to the forest type (Feng et al. 2022); therefore, the changes in GPP caused by changes in FFC showed differences in various regions of the TNR.

5.3. Limitations and shortcomings

Complex interactions exist between human activity and climate change. This study conducted an independent analysis of their impacts on forest change, ignoring their mutual influence, which may have impacted the accuracy of the estimation of the contribution rate of forest change. The response of vegetation, climate, and ecosystems is a complex and interactive process. In this study, only temperature and precipitation were considered, thereby simplifying the limiting factors that affect forest vegetation growth. By considering the effects of topography (slope and slope direction), economy, population, policy, and other climatic factors (CO₂, solar radiation, sunshine hours, etc.) on forest vegetation, a deeper and more precise explanation of the causes of forest change can be obtained. The resolution of the FFC data used in this study was 30 m, whereas the resolution of other data was 500 m or 1 km. To maintain consistency between the data, FFC was resampled to 1 km, which would have caused a lack of data, and at the same time would have reduced the ability to extract the change information on a fine spatial scale. Thus, there is uncertainty in the estimation results. As the resolution of various data types increases, data with consistent resolution will reduce the uncertainty and improve the accuracy of. In addition, this study used the inter-annual synthesis of Landsat to construct annual FFC. Forest vegetation in the TNR is of multiple types which have different growth characteristics in different seasons, and the monsoon-type climate of the TNR also causes obvious hydrothermal differences in different seasons. Thus, they will have different interactions and responses in different seasons. Conducting research on seasonal forest cover changes is conducive to revealing the response patterns of forest and climate change as well as mutual driving, revealing the cycle and operation mechanism of the ecosystem in the region, and further guiding the construction of scientific ecological engineering implementation plans.

6. Conclusions

Based on long time series, 30-m resolution fractional forest cover data combined with multi-source data and integrated with multiple analysis methods, this study revealed the spatiotemporal variation in forests of the three-north regions based on pixel and regional scales.

The trend analysis by time period demonstrated the growth status and change characteristics of the forests in three-north regions. Fractional forest cover showed a decreasing trend before 2003 and an increasing trend from 2003 to 2017. The increasing trend after 2010 was stronger than that before 2010. The stability analysis showed that 56% of the forests were in a state of low-to-medium volatility, mainly in the northeast. The Loess Plateau and Tianshan Mountian areas mainly showed high or higher fluctuations, among which the forest vegetation of Loess Plateau was more unstable owing to the influence of climatic factors, whereas, the forest cover fluctuation in the Tianshan Mountian area was caused by pests and diseases. Vegetation dynamics were combined with a transition matrix to jointly analyze the forest changes of each fractional cover level. After 2017, the transition of high-density forests to medium-density forests indicated that high-density forests were degrading. Medium-density forests did not shift to low-density forests, indicating that high-density forests were in the early stages of degradation. It may be that forest aging led to a decline in the forest, resulting in a decrease in fractional cover.

The climate data trend showed that the temperature in the three-north regions has been increasing. The correlation analysis showed that the correlation between precipitation and fractional forest cover in this region was higher than that between temperature and fractional forest cover. Based on the residual estimation concept, a method for estimating the contributions of climatic and non-climatic factors was developed using machine-learning methods. The driving factors affecting fractional forest cover were analyzed, and the results showed that non-climatic factors (including human factors) played a dominant role in the impact of fractional forest cover (contribution rate = 66%), and that the contribution rate of climatic factors was 34%. Among non-climatic factors, forest fires were one of the main factors affecting forest fluctuations in the Greater Khingan Mountains and Lesser Khingan Mountains regions. The use of a single tree species for artificial afforestation is likely to cause harmful biological disasters in the forests of Xinjiang and Inner Mongolia, thus causing fractional forest cover fluctuations in these areas. The impact of human factors on the distribution and change of forest cover is complex, among which afforestation activities have been the most direct factor in the increase in forest cover in the three-north regions. However, human intervention has restricted the effects of natural geographical conditions on changes in the fractional forest cover pattern, the most obvious of which is that the increase rate of forests in mountains with higher slopes is higher than that of areas with flat terrain. This is because people occupying flat land for production and living activities also make full use of high-slope areas for suitable afforestation and grass-planting projects, hence, the forest coverage in high-altitude and high-slope areas gradually increases. By analyzing the feedback effect of fractional forest cover, it can be concluded that the increase in fractional forest cover has had a positive effect on the overall ecology of the three-north regions, that is, the improvement of forest vegetation coverage through ecological projects such as afforestation has increased vegetation productivity, and can also have a cooling effect on the surface.

Data and codes availability statement

Data available on request from the authors

Acknowledgments

We thank the Global LAnd Surface Satellite (GLASS) team for providing the products to the public (http://glass.umd.edu/index.html, http://glass-product.bnu.edu.cn/index.html). Thanks for the data support provided by the National Qinghai-Tibet Plateau Scientific Data Center (http://data.tpdc.ac.cn). Thanks to the National Bureau of Statistics and China Forestry Information Network for the statistical yearbook data (http://www.stats.gov.cn/tjsj/ndsj/, http://www.lknet.ac.cn/).

Disclosure statement

No potential conflict of interest was reported by the authors.

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/15481603.2023.2300214.

Additional information

Funding

This work was supported by the The Open Research Program of the International Research Center of Big Data for Sustainable Development Goals [No. CBAS2022ORP01]; The National Natural Science Foundation of China [No. 42090011].