Except global environmental change, a variety of human activities including agricultural production and natural resource exploitation, wildfire, landslide and debris, and invasive species, are placing immense pressure on forests (Bladon, 2018; Erbaugh and Oldekop, 2018; Fischer, 2018; Favero et al., 2020). Large-scale forest restoration and management efforts are required to ensure the forest-related ecosystem services, such as carbon sequestration, biodiversity conservation, and provision contribution (Wang et al., 2007; Chazdon, 2008; Lewis et al., 2015; Ouyang et al., 2016; Lu et al., 2018; Favero et al., 2020). Investments in forest landscape restoration and management play a great role in improving livelihoods and well-being (Kareiva et al., 2011; Ouyang et al., 2016; DeFries and Nagendra, 2017; Erbaugh and Oldekop, 2018; Fischer, 2018). Forest management on large spatial and temporal scales is urgently needed to counteract the effects of environmental change and human activities, which are the result of large-scale ecological patterns and processes that are not limited by administrative borders (Fischer, 2018). Emerging paradigms of forest management emphasize the need to consider decision-making, and how the present and potential decisions may interact with social and ecological conditions and processes that affect landscapes across space and time (Filotas et al., 2014; Messier et al., 2015; Fischer, 2018). Moreover, land management practices influence ecological patterns and processes well into the future with impacts that often go unobserved for long periods of time (Fischer, 2018).

Forest landscapes are typical socio-ecological systems composed of interdependent biophysical components and associated human dimensional factors (Ostrom, 2009; Angelstam et al., 2013; Spies et al., 2014; Guerrero et al., 2018; Wu et al., 2020). A landscape pattern is the key element of a basic landscape structure; therefore, it has a significant role in natural conservation and landscape management (Reed et al., 2017; Zeng et al., 2017; Erbaugh and Oldekop, 2018; Gong and Xie, 2018). The analysis of landscape patterns is a major issue of landscape ecology (Wu and Hobbs, 2002; Zhang et al., 2011; Turner and Gardner, 2015). Changes in landscape patterns can affect the ecosystem’s structure and function, the provision of ecosystem services, through altering surface biophysical parameters (MA, 2005; Ouyang et al., 2010; Kindu et al., 2016). Landscape metrics, the simple quantitative metrics reflecting the composition of landscape structure and attributes of spatial land use allocation (Wu et al., 2002), are useful for studying landscape patterns (Ren et al., 2013), also are useful to reveal relationships between landscape structures and ecological processes (Fry et al., 2009; Ren et al., 2013). On a landscape/regional scale, landscape component directly affects the provision of ecosystem services, while landscape pattern indirectly impacts ecosystem services through changing ecological processes (Fagerholm et al., 2012; Fu et al., 2013; Jia et al., 2014; Rieb and Bennett, 2020). Therefore, understanding the relationship between landscape patterns and ecosystem services is vital to maintaining ecosystem integrity and sustaining Earth (Kozak et al., 2011; Gray and Lee, 2017), and has become more important to the study of landscape ecology and ecosystem management than ever before.

The soil conservation service is a fundamental ecosystem regulation service that is essential to sustainable development. In quantitative terms, it is obtained as the difference between the potential and actual soil erosion (Ausseil et al., 2013; Rao et al., 2014; Liu et al., 2019). The soil conservation service can be obtained with the Revised Universal Soil Loss Equation (RUSLE) (Liu et al., 2019), Water Erosion Prediction Project (WEPP) models (Revuelta-Acosta et al., 2021), Soil and Water Assessment Tool (SWAT) (Cong et al., 2020; Shi and Huang, 2021), and InVEST-SDR (Sediment Delivery Ratio (SDR) in Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST)) (Posner et al., 2016; Sharp et al., 2020), which are all useful tools for quantifying and managing soil conservation (Daily et al., 2009; Kareiva et al., 2011; Fu et al., 2013; Liu et al., 2019; Sharp et al., 2020). SWAT, RUSLE, WEPP, and InVEST-SDR have been extensively validated and used at the watershed and hillslope scales in different environments and countries (Liu et al., 2019; Shi and Huang, 2021; Revuelta-Acosta et al., 2021). InVEST-SDR has been widely used for its reliable feature to predict water erosion in different environments and its integration with ArcGIS and satellite data (Aneseyee et al., 2020; Cong et al., 2020; Ougougdal et al., 2020; Sharp et al., 2020). InVEST-SDR can also be used to obtain customized soil conservation measures and management policy for specific watershed characteristics (Zhang et al., 2019; Aneseyee et al., 2020; Ougougdal et al., 2020; Zhou et al., 2020).

Land use and landscape patterns affect changes in soil conservation (Fu et al., 2013; da Silva et al., 2015; Turner and Gardner, 2015; Zhang et al., 2019). Landscape metrics have been found useful to show the relationship between landscape patterns and soil conservation services and the main methods of analysis are correlation, regression, residual analysis and spatial econometric models (Shi et al., 2013; Zhang et al., 2017; Jiang et al., 2020). For example, Ouyang et al. (2010) selected landscape metrics to explore the soil erosion dynamics response to landscape patterns and found that contiguous grassland patches reduced soil erosion. da Silva et al. (2015) indicated that landscape metrics, including PLAND (percentage of landscape), NP (number of patches), PD (density of patches), LPI (largest patch index), ED (edge density), LSI (landscape shape index) and FRAC (fractal dimension index) were sensitive to the changes in soil surface micro-morphology. Shi et al. (2013) found that the main factors that influenced watershed soil erosion and sediment yield at the landscape level were SHDI (Shannon diversity index), AI (aggregation index), LPI, CONTAG (contagion) and COHESION (cohesion index). Liu (2017) found that landscape metrics like patch/edge density, shape indices, and diversity indices at the landscape level were effective in linking landscape patterns with soil erosion. Jiang et al. (2020) found that the soil erosion module correlated positively with NP, PD, ED, LSI, SHAPE (mean shape index) and FRAC, while it correlated negatively with LPI, ENN (Euclidean nearest- neighbor distance), IJI (interspersion and juxtaposition index), AI, and PLADJ (percentage of like adjacencies).

Afforestation activities and soil and water conservation measures can reduce erosion to promote soil retention (Wang et al., 2007; Nunes et al., 2011; Alatorre et al., 2012; Feng et al., 2016; McEachran et al., 2018; Jiang et al., 2020; Zhou et al., 2020). For example, Jiang et al. (2020) found that vegetation restoration can improve vegetation cover and soil retention, and that dense forest cover is a major factor in preventing erosion in the Chinese Loess Plateau. Furthermore, soil erosion is affected by various factors like climate, topography, vegetation cover, human interferences, improper land use, agricultural practices, and deforestation, etc. (Zhang et al., 2017; Diwediga et al., 2018). Therefore, it is urgent to understand the relationship between forest landscape patterns and soil conservation and identify landscape metrics to monitor it and benefit forest landscape pattern optimization and management.

Mountainous areas provide diverse goods and services to human societies but are sensitive to rapid land use/cover change, fragmentation of habitats, natural disasters, human activities and climate change (Schröter et al., 2005; Grêt-Regamey et al., 2012; Gong et al., 2014; Xie, 2015; Wang et al., 2017). Protecting fragile mountainous ecosystems through afforestation and ecological restoration to enhance ecosystem services is a huge challenge (Grêt-Regamey et al., 2012). A series of afforestation programs, the most ambitious revegetation and conservation projects, including the Natural Forest Protection Project (NFPP), the Grain to Green Program (GTGP), and soil and water conservation (SWC) programs have been implemented in mountainous areas of western China (Ouyang et al., 2016). The design and management of the forest landscape to improve the effectiveness of afforestation and ecological restoration are of scientific and practical concern (Brunckhorst, 2011; Hainz-Renetzeder et al., 2015). Thus, mountainous ecotones in western China, which undergo the pressures placed by geo-disasters, a large population and intensive human activities, should be paid more attention to forest landscape design and management to enhance their vital ecosystem services (Wang et al., 2017; Fu et al., 2018). In this study, Bailongjiang Watershed (BLJW) in western China, a typical ecotone with fragile environments that experiences frequent landslides, soil erosion, and intensive human activities, is selected to understand how the forest landscape pattern will affect soil conservation services. The objectives of this study were to: 1) assess the soil conservation service on a sub-basin scale via the SDR model of InVEST; 2) identify important forestland landscape indices for soil conservation; 3) characterize the spatial distributions between soil conservation and the selected forestland landscape indices, and 4) establish a robust basis for forestland landscape management to enhance soil conservation services for mountainous areas.

Bailongjiang Watershed (BLJW) (103°00′E–105°30′E and 32°36′N–34°24′N), located in the transitional zone of the Tibetan Plateau, Qinba Mountains region and Loess Plateau, with an area of 18 437.7 km², is vital for soil and water conservation in the western China due to its steep slopes and deep valleys (Fig. 1). BLJW is characterized by a subtropical climate, annual mean temperatures that range from 2℃ to 15℃ and a mean annual precipitation of 450 mm, of which over 80% falls from June to September. The dominant land use types of BLJW are forestland, farmland, and grassland. With the implementation of a series of afforestation/reforestation policies in China, forestland has become the main and most important land-use type of watershed (Zhang et al., 2019). The soils in BLJW are mainly leached-cinnamon soils, brown soils, dark-brown soils and alpine meadow soils. Vegetation types are temperate deciduous broad-leaved forest, evergreen broad-leaved forest, temperate mountain coniferous forest and alpine mountain coniferous forest. The ongoing land use and human activities in the watershed have resulted in severe land degradation and soil erosion.

Figure 1.  Location of the study area: Bailongjiang Watershed and its 37 hydrological sub-basins in Gansu Province of the western China

BLJW was divided into 37 hydrological sub-basins response units via the distributed watershed hydrological model—SWAT—to obtain spatial patches for analysis of landscape pattern and soil conservation (Fig. 1). Based on the sub-watershed unit, the landscape pattern and soil conservation service were quantitatively characterized via the landscape pattern index and SDR model of InVEST. Mathematical and spatial statistics methods were used to reveal the relationship between landscape metrics and soil conservation of the 37 sub-basins in BLJW. Soil conservation of each sub-basin was assessed via the SDR model of InVEST. The landscape metrics of each sub-basin were calculated with Fragstats version 4.2 (McGarigal et al., 2012). Our research framework included several steps: 1) land use classification and mapping via Landsat Thematic Mapper (TM) and Enhanced Thematic Mapper (ETM+) images; 2) calculation of forestland landscape indices (area and edge metrics, shape metrics, and aggregation metrics) via Fragstats version 4.2; 3) assessment of soil conservation service via the SDR model of InVEST; 4) identification of the important forestland landscape indices via mathematical statistics; 5) spatial correlation between the identified forestland landscape indices and soil conservation via spatial statistics methods; 6) relationship between soil conservation and forestland landscape pattern and implications for landscape pattern optimization and watershed management.

The DEM dataset with a spatial resolution of 30 m was downloaded from the Geospatial Data Cloud (http://www.gscloud.cn). Land-use maps of BLJW in 1990, 2002, and 2014 were obtained via Landsat TM/ETM satellite data in the summer (July and August) with a spatial resolution of 30 m downloaded from the United States Geological Survey (http://glovis.usgs.gov). In this study, land use types were divided into six categories, including farmland, forestland, grassland, water, construction land, and unused land, according to the Land-Use and Land-Cover Change (LUCC) classification standard and the current situation of BLJW. The overall accuracies in 1990, 2002, and 2014 were 89.02%, 90.71%, and 91.08%, respectively, which are acceptable for land use change analysis (Gong and Xie, 2018). The soil map in BLJW was obtained from the 1∶1 000 000 soil type map of Gansu Province; the soil particle and soil organic carbon data were obtained via field sampling and laboratory analysis in July 2012. Meteorological data, including the monthly precipitation of the 17 meteorological stations in BLJW and its adjacent area from 1990 to 2014, were obtained from the Chinese Meteorological Data Service Center (http://cdc.cma.gov.cn).

The soil conservation of ecosystems can best be described as soil erosion prevention (Ausseil et al., 2013; Sharp et al., 2020). Based on the InVEST-SDR model, the soil conservation service was obtained by the sediment retention index (SEDRET) with the revised universal soil loss equation (RUSLE). An index of sediment retention on pixel i was calculated by the model as follows:

where SEDRET_i is an index of sediment retention weighted by the SDR factor. SDR_i is the SDR ratio for pixel i, R_i is the rainfall erosivity ((MJ·mm) / (ha·h)), K_i is the soil erodibility ((t·ha·h) / (MJ·ha·mm)), LS_i is the slope length-gradient factor, C_i is the crop management factor, P_i is the practice factor (Renard et al., 1997), SDR_max is the maximum theoretical SDR and is set to an average value of 0.8 (Vigiak et al., 2012), and IC₀ and k are calibration parameters that define the shape of the SDR-IC relationship (see more details in the InVEST User’s Guide). The procedure for the calculation of R_i, K_i, C_i and P_i refers to the research work of Zhang et al. (2019).

Fragstats (version 4.2) was used to calculate the landscape pattern indices from the grid files of land use in 1990, 2002, and 2014. More than 110 landscape pattern metrics can be obtained for each land use grid. However, if all the possible landscape pattern metrics had been used as effective variables for regression analysis, it would have caused workload redundancy, also introduced the curse of dimensionality due to the linear relations between the various metrics (Zhang et al., 2017). Therefore, it is vital to select the most representative metrics that are closely correlated with soil conservation. There are three kinds of landscape pattern metrics at the class level: area and edge metrics, shape metrics and aggregation metrics (McGarigal et al., 2012). In this study, 11 landscape indices were selected based on their suitability to characterize the forest landscape pattern: PLAND (percentage of landscape), LPI (largest patch index), PD (patch density), AREA_MN (mean patch size), ED (edge density), LSI (landscape shape index), SHAPE_MN (mean patch shape), COHESION (patch cohesion index), DIVISION (landscape division index), SPLIT (splitting index), AI (aggregation index). More details on the exact information and calculation can be found in McGarigal et al. (2012).

Pearson correlation analysis between soil conservation and forestland landscape metrics was carried out to select the landscape pattern metrics that significantly correlated with soil conservation (Table 1). The preliminary analysis showed that most of the forestland landscape metrics were significantly related to soil conservation at the 0.05 level including PLAND, LPI, ED, LSI, SHAPE_MN, COHESION, DIVISION, SPLIT, and AI. Furthermore, it was determined that soil conservation was affected by the combination of area and edge metrics, shape metrics and aggregation metrics. The correlation coefficient between PLAND for forestland and soil conservation was the largest one, which indicated that the ratio of forestland played an important role in watershed soil conservation in western China. Finally, nine forestland landscape metrics were used as the independent variables for the multiple linear regressions.

Regression analysis was performed to determine the correlations between two or more variables with cause-effect relations. Multiple linear regression models were used to describe how a single response variable depends linearly on some predictor variables. The multiple linear regression equation has the form:

where y is the observed value of the soil conservation, $\hat y$ is the predicted value of the soil conservation. β₀ is a constant term, and β₁, ..., β_p (p = 1, 2, ..., 9) are partial regression coefficients. x_PLAND, x_LPI, ..., x_AI are the values of nine landscape metrics for forestland including PLAND, LPI, ED, LSI, SHAPE_MN, COHESION, DIVISION, SPLIT, and AI. ε is the residual error.

For multiple linear regression analysis, the significance test of the equations and variables were used to estimate the accuracy of the regression model. The regression equation and the contribution of each independent variable should be statistically significant to ensure that the multiple linear regression equation consistent with the data characteristics (Zhang et al., 2017). Therefore, hypothesis testing was needed for the established regression equation.

(1) Hypothesis test for the regression equation

The F test was selected as the hypothesis test for the regression equation. The null hypothesis is H₀: β₁ = β₂ = … = β_p = 0, that is, no regression model can be established. The F statistic is formulated as the following (Zhang et al., 2017):

where S_reg is the sum of squares of the regression, S_res is the sum of squares of the residuals, $\hat y$ is the predicted value of the dependent variable, $\bar y$ is the mean value of the dependent variable, y is the observed value of the dependent variable, m is the degrees of freedom of the regression, n is the number of sample data and (n−m−1) is the degrees of freedom of the residuals. If F > F_{α = 0.05 (m, n−m−1)}, the null hypothesis is rejected, which shows that there is a linear relationship between the dependent and the independent variable.

(2) Hypothesis test of the partial regression coefficient

When a regression model is tested for significance, the t-test is used to determine which models are important explanatory variables. The null hypothesis of the t-test is H₀: β_i = 0 (i = 1, 2, …, m; m = 9). The t-statistic is determined as:

where b_i is the partial regression coefficient of the independent variable, S_bi is the standard error of b_i.

Stepwise regression is a useful method for selecting independent variables with no significant collinearity for multiple linear regression models (Neter et al., 1996). Specific steps of stepwise regression were:

(1) The simple linear regression models (only include one variable) were fitted with nine candidate independent variables (PLAND, LPI, ..., AI). In addition to this, k (k ≤ 9) linear regression models satisfied with the F test, and PLAND was first introduced into the model because of the smallest p value.

(2) After the PLAND variable was introduced into the model, the independent variables were combined as PLAND and LPI, PLAND and ED, and so on. Since the P value of PLAND and AI was the smallest and satisfied with the F test, the AI variable was also added to the model.

(3) Then, the PLAND and AI were incorporated in the model. The t-test was used to determine whether the PLAND variable still had statistical significance. If the PLAND did not pass the t-test, it was removed from the model. If the regression models of the other seven variables including AI were fitted well, then the combination of independent variables with the lowest P value was introduced into the model. The computing process was terminated until all statistically significant independent variables were included in the model.

(4) If PLAND passed the t-test, it was incorporated in the model. The regression analysis continued to fit the model with the remaining seven variables, beginning with the PLAND and AI variables. The rest of the variables were tested for statistical significance, and the variable with the smallest P value was gradually incorporated into the model. This kind of computing process was repeated until the independent variables without statistical significance were excluded from the model.

The bivariate Moran’s I, including the global bivariate Moran’s I (I_sl) and the local bivariate Moran’s I (I’_sl), is an effective way to explore spatial clustering (positive spatial correlation) and spatial dispersion (negative spatial correlation). In this study, bivariate LISA was used to reveal the spatial correlations between soil conservation and forestland landscape metrics at the sub-basin unit based on GeoDA software (Anselin and Rey, 2014). The I_sl and I’_sl are the global and local bivariate Moran’s I for soil conservation and each forestland landscape metric (EFLM), respectively. The calculation equations can be found in Anselin and Rey (2014). The values of I_sl /I’_sl range from −1 to 1. A positive I_sl /I’_sl value indicated a positive spatial correlation between a soil conservation value and EFLM, which signifies that the sub-basin with a high soil conservation value is likely to be surrounded by a sub-basin with high EFLM values; a negative value indicates that there is a negative spatial correlation. The greater the absolute value of I_sl /I’_sl, the closer of the spatial correlation between soil conservation and EFLM.

The bivariate LISA method can be used to visualize the spatial correlations by generating cluster maps. The cluster map generated by the bivariate LISA represents four types of local spatial autocorrelations, 1) the high-high type (shorted as HH) shows high values of soil conservation surrounded by high values of the forestland landscape index; 2) the high-low type (shorted as HL) shows high values of soil conservation surrounded by low values of the forestland landscape index; 3) the low-high type (shorted as LH) shows low values of soil conservation surrounded by high values of the forestland landscape index; and 4) the low-low type (shorted as LL) shows low values of soil conservation surrounded by low values of the forestland landscape index.

The spatial regression method is used to reveal the spatial dependence of soil conservation on forestland landscape metrics for 1990, 2002, and 2014. Spatial Regression Models (SRMs), including Spatial Lag Model (SLM) and Spatial Error Model (SEM), were developed via Ordinary Least Squares (OLS) by introducing a spatial weight matrix into the regression (McMillen, 2004; Chi and Zhu, 2008). SLM depicts spatial dependency by dependent variables, whereas SEM accounts for spatial error dependency (Anselin et al., 2006). The equations of SLM and SEM can be found in Anselin and Rey (2014). The model performance comparison of SLM and SEM is carried out by the Lagrange Multiplier (LMLAG and LMERR) and the Robust Lagrange Multiplier (R-LMLAG and R-LMERR). If the LMLAG is more statistically significant than the LMERR, and the R-LMLAG is significant while the R-LMERR is not, then SLM is more appropriate; conversely, SEM is more appropriate. If the model performance comparison cannot be established, statistical parameters like R², Log likelihood (LogL), Akaike information criterion (AIC) and Schwarz criterion (SC), can be used. Normally, the improved model performance is indicated by the increased R² and LogL, as well as by decreased AIC and SC (Anselin et al., 2006).

The soil conservation of BLJW initially decreased and then increased from 1990 to 2014 (Fig. 2). The minimum value of soil conservation decreased continuously, which indicated that soil conservation declined in some sub-basins. Spatially, there were few substantial changes in soil conservation over the three years. High values of soil conservation were mainly distributed in the southeastern parts of BLJW (such as in sub-watershed No. 1, 2, 3, 4, 5, 9, 10 and 11) (Fig. 2), which mostly consist of nature reserves with less human activities. Low values of soil conservation were distributed in the northern and middle parts of BLJW (such as in sub-basins No. 30, 31, 14, 20, 21 and 22), where agricultural activities are more intense (Fig. 2).

Figure 2.  Soil conservation service of the sub-basins in Bailongjiang Watershed in 1990, 2002, and 2014

The independent variables were the nine forestland landscape indices and the dependent variable was soil conservation. Upon completion of the regression analysis, four regression models were obtained (Table 2). The Sig. F indicated the significance of the regression equation at the 0.05 level, and the Sig. F values of the four models were 0.000, which indicated that the models provided a good fit to the data. The R² values reflected the proportion of the variance of the original data explained by the regression equation. The R² of the four models were 0.183, 0.231, 0.298 and 0.394, respectively, indicating that the fourth model explained the largest amount of variance in the data. However, the R² coefficient has some limitations when evaluating the model. For example, the value of R² may increase even if the introduced variable is not statistically significant. Thus, the adjusted R² was used as an additional statistical parameter to evaluate the model. In this study, the R² and adjusted R² both increased as the number of independent variables increased.

Five variables were excluded by model 4 and its test results are listed in Table 3. In the study, the t-test was applied to the independent variables at the 0.05 level. As shown in Table 3, the Sig. T values of the five variables added to model 4 were greater than 0.05, which indicated that the five variables were not statistically significant for regression analysis. Moreover, the variance inflation factor (VIF, which is the reciprocal of tolerance) of the collinearity analysis was greater than 15 for the variable COHESION, indicating significant collinearity between variables. Therefore, this variable was excluded from the regression model.

Table 4 indicated the coefficients of the variables in model 4, including the standardized coefficients of the model variables, the probability of the t-test, and the collinearity statistics (tolerance and VIF). The probabilities of the t-test for the constant and the four independent variables were less than 0.05, which indicated that the variables were statistically significant. The VIF values of the four variables were <10, which indicated that there was no obvious collinearity between the four independent variables. The mean value of the standardized variable was 0 and the standard deviation was 1, so the con-stant term of the fitted regression in this study was 0. The final multiple linear regression equation of soil conserva-tion is 1.079PLAND−1.050AI−0.356LSI−0.463SPLIT.

Based on stepwise regression analysis, four variables (PLAND, AI, LSI, and SPLIT) were determined to be more suitable for soil conservation without significant collinearity. The spatial distributions of soil conservation and each of the four forestland landscape metrics within the sub-basin unit for the three years is shown in the bivariate LISA maps (Fig. 3). We found a clear similarity in the spatial clustering patterns for the three years of each of the four forestland landscape metrics and soil conservation. As for the bivariate LISA map on the relationship between soil conservation and PLAND of forestland, LL areas were mainly concentrated in the middle of BLJW (sub-watersheds No. 14, 20, 21, 22 and 23), which indicated that areas with lower soil conservation were surrounded by areas of forestland with lower values of PLAND (Fig. 3). Compared with the clustering pattern above, the LL areas of soil conservation and the AI of forestland occupied less area were mainly distributed in sub-watersheds No. 14 and 21. According to the bivariate LISA map on the relationship between soil conservation and LSI of forestland, HL areas were mainly concentrated in the south of BLJW (such as in sub-watershed No. 1 and 5) (Fig. 3). According to the bivariate LISA map on the relationship between soil conservation and SPLIT of forestland, LH areas were distributed in sub-watersheds No. 14, 21 and 22. These maps revealed that soil conservation in sub-basins No. 1, 5, 14, 20, 21, 22 and 23, was closely related to the forestland landscape pattern of 1990, 2002, and 2014 (Fig. 3).

Figure 3.  Bivariate LISA maps of soil conservation and forestland landscape indices of the sub-basins in Bailongjiang Watershed in 1990, 2002, and 2014. PLAND, percentage of landscape; AI, aggregation index; LSI, landscape shape index; SPLIT, splitting index

Although there was a clear similarity of the clustering pattern in 1990, 2002, and 2014 between spatial distributions of each of the four forestland landscape metrics and soil conservation, there was also non-negligible dissimilarity. The relationship between soil conservation and PLAND of forestland showed larger HH areas distributed in sub-watershed No. 1, 3, 5, 6, 26 and 34 in 2002. In addition, the spatial clustering pattern between soil conservation and LSI of forestland was different in each of the three years, and was mainly concentrated in sub-watersheds No. 33, 34, 35, 36 and 37. This could be caused by the larger variation in soil conservation in these sub-basins.

Spatial regression analysis was used to reveal the relationship between soil conservation and the four variables retained in model 4 by the GeoDA software. The OLS results showed that the PLAND, AI, and SPLIT of forestland were highly significant at the 0.05 level for all soil conservation regressions in 1990, 2002, and 2014 (Table 5), while the LSI of forestland was not significant in 2002. Compared to the adjusted R² values of the three regressions, the adjusted R² value was highest (0.650) in 2002, which indicated that soil conservation was explained by 65% of the four forestland landscape metrics (PLAND, AI, LSI, and SPLIT).

Due to the spatial dependence of soil conservation within the sub-basin unit, it was necessary to diagnose the estimation results of OLS to verify whether there was a more reasonable spatial regression model in the three years. The diagnostics of OLS showed that apparent spatial dependence existed in the regression in 1990 (Table 6). Statistical diagnostics showed that LMLAG was more statistically significant than the LMERR, and that R-LMLAG was significant while R-LMERR was not significant at the 0.05 level. That is, SLM is more appropriate. However, both R-LMLAG and R-LMERR were not significant for the regression in 2014, which indicated that it is impossible to directly determine which regression model (i.e., SLM, SEM, and OLS) was better. Four statistical parameters can be considered as R², LogL, AIC and SC. In general, the larger R² and LogL and the smaller AIC and SC, the more reasonable the OLS. Because the P value of all the tests was far greater than 0.05 and the adjusted R² value of OLS was larger for the regression in 2002 (Table 4 and Table 5), the spatial regression analysis was not carried out.

The estimation results of the spatial regression models showed that SLM was better than OLS for the regression in 1990, as R² and LogL increased and AIC and SC decreased (Table 5 and Table 7). R² and LogL are higher when the regression model better explains the variance of the dependent variable. Compared to the R² and LogL values of OLS, SLM, and SEM of the regression in 2014, the results showed that R² and LogL values in SLM were the highest, followed by the SEM and OLS (Table 5 and Table 7). The lower the AIC and SC values, the simpler the model, therefore, SLM was used as a more reasonable regression model for soil conservation in 2014, even if the SC value in SLM was slightly higher. The coefficients of the ρ variable in SLM were significantly positive at p < 0.05 of the soil conservation regressions in 1990 and 2014, which indicated that the soil conservation of 1990 and 2014 was positively influenced by the surrounding environmental system. The coefficient of the λ variable in SEM was not significant at the 0.05 level of the soil conservation regression in 2014, which showed that the soil conservation of 2014 was mainly influenced by the four forestland landscape metrics from the perspective of landscape patterns.

Pearson correlation analysis shows that the forestland landscape indices, including PLAND, LPI, ED, LSI, SHAPE_MN, COHESION, DIVISION, SPLIT, and AI, played a crucial role in soil conservation. Soil conservation was positively correlated with PLAND and AI, but negatively correlated with LSI and SPLIT (Table 1). The value of PLAND ranged from 0 to 100, indicating the ratio of the landscape occupied by patch type. The correlation coefficient showed that there was a strong positive correlation between soil conservation and PLAND, indicating that the proportion of forestland in the sub-basin played a dominant role in soil conservation. Soil conservation was also positively correlated with AI, which indicated that the decrease of forestland patch aggregation resulted in reduced soil conservation in the sub-basin. The negative correlation between soil conservation and LSI means that soil conservation was reduced as the shape of forestland became more complex. Soil conservation was negatively correlated with SPLIT, which indicated that human disturbance may have a negative effect on soil conservation.

The stepwise regression analysis also indicated that soil conservation was negatively correlated with LSI and SPLIT, but positively correlated with PLAND (Table 4). However, in contrast to the results of the Pearson correlation analysis, the stepwise regression results showed that soil conservation was negatively correlated with AI (Table 4). This situation was acceptable because the Pearson correlation analysis focused on the relationship between two variables without consideration of the mutual effects of other variables, while the stepwise regression placed emphasis on the mutual effects of independent variables, which may change the direction of the Pearson correlation coefficients due to the collinearity between the independent variables. Although stepwise regression is an effective method for eliminating the independent variables that cause obvious collinearity, it is still impossible to completely avoid collinearity. In this study, the VIF values of the four variables PLAND, AI, LSI, and SPLIT were <10 (Table 4), indicating that there was no apparent collinearity between the four independent variables. However, the stepwise regression only changed the direction of the Pearson correlation coefficients for AI, which may be because the VIF of AI was slightly higher than that of the other three variables. It did not affect the reasonableness of the stepwise regression analysis for watershed soil conservation.

The stepwise regression analysis also showed that four regression models were acquired with adjusted R² values of 0.175, 0.217, 0.278 and 0.371, respectively (Table 2). Apparently, the fourth model was the best one of the models. The fourth model showed that the forestland landscape pattern was closely related to watershed soil conservation, not only to the forestland area but also to the forestland shape and aggregation. Because forestland is one of the main land-use types in BLJW, optimizing forest landscape patterns can improve watershed soil conservation (Zhang et al., 2017, 2019; Gong and Xie, 2018).

The bivariate LISA maps indicated that the number of sub-basins with significant spatial correlations between soil conservation and forestland landscape indices including PLAND, AI, LSI, and SPLIT was 15, 19, and 20 in 1990, 2002, and 2014, accounting for 40.54%, 51.35%, and 54.05% of the total number of sub-basins, respectively (Fig. 3). This showed that there is a scientific basis for promoting watershed soil conservation by optimizing the forest landscape pattern of the sub-basins. Furthermore, the SLM was superior to the OLS for the regression in 1990 and 2014, which indicated that soil conservation within the sub-basin unit was spatially dependent, while OLS was more appropriate for the regression in 2002, which indicated that the forestland landscape pattern played a greater role in watershed soil conservation.

The soil conservation service in BLJW initially decreased and then increased from 1990 to 2014. There is higher soil conservation in the western parts and the lower reaches of BLJW (Diebu, western Zhouqu and Wenxian). The areas with lower soil conservation service were mostly located along the Bailongjiang valley and eastern parts of BLJW. The results also showed that the sub-watersheds with higher soil conservation services were mostly located in the lower reaches of BLJW with higher forest vegetation coverage, and the sub-watersheds with lower soil conservation service were mostly distributed along Bailongjiang valley with lower forest coverage and intensive human activities (e.g., industry and agricultural production, urbanization, road construction). This is similar to the research results of Xie (2015), Gong and Xie (2018) and Gong et al. (2021). The sub-watershed with high soil conservation service was also characterized by single landscape types composition, non-uniform distribution among different landscape types, dominant patches, and a low degree of landscape separation. The low degree of landscape separation was also found by Zhang et al. (2019), who indicated that the landscape pattern design and planning would play a key role to prevent soil erosion, especially in the afforestation/ecological restoration projects in the future.

The soil conservation service can control soil erosion and retain sediment, and therefore generate positive effects on land productivity and many other ecosystem services that ensure human welfare (Fu et al., 2013; Liu et al., 2019; Zhang et al., 2019). Studies found that afforestation and vegetation restoration play a vital role in controlling soil loss and enhance the soil conservation service by decreasing runoff and soil erosion (Long et al., 2006; Xu et al., 2006; Mehri et al., 2018; Jiang et al., 2020; Wang et al., 2020). Some previous studies have also revealed vegetation as the most important influence factor against soil erosion (Jiang et al., 2020; Zhang and Wei, 2021). For example, Mehri et al. (2018) highlighted that forest cover plays a vital role in reducing rainfall erosivity and subsequent soil erosion. Sun et al. (2014) found that the erosion rate and susceptibility of soil loss were obviously decreased by forest and dense grassland due to the increasing of vegetation cover.

The landscape pattern has an important role in landscape management, nature conservation, ecosystem service provision, and human well-being (Wu and Hobbs, 2002; Chazdon, 2008; Ouyang et al., 2010; Turner and Gardner, 2015; Kindu et al., 2016; Erbaugh and Oldekop, 2018; Rieb and Bennett, 2020). Jiang et al. (2020) found that properly fragmented landscapes can effectively reduce soil erosion in the Loess Plateau. Landscape heterogeneity and landscape structure also generally lead to variations in ecological and hydrological conditions which alter sediment production, transportation, and connectivity (Ai et al., 2015; Li et al., 2019; Sun et al., 2019). Our study showed that the four forestland landscape indices, including PLAND, AI, LSI, and SPLIT, were vital for soil conservation in the typical ecotone watershed, BLJW, in western China. That is, landscape types’ diversity and evenness were the most important indicators affecting soil conservation in the mountainous area. Therefore, the soil conservation service can be enhanced by forestland landscape pattern optimization and management (Fu et al., 2013; Liu, 2017; Gong and Xie, 2018; Zhang et al., 2019; Rieb and Bennett, 2020).

Since the late 1970s, China has carried out four national key afforestation/reforestation programs to protect the environment and restore the degraded ecosystems, especially affecting western China. These projects include 1) the Three-North Shelter Forest Program, which was initiated in 1978 and is known as the ‘Green Great Wall’, 2) the Yangtze River Shelter Forest Projects, launched in 1989 to combat floods and soil loss, 3) the Natural Forest Protection Project, which was initiated in 1998 for biodiversity conservation, reduction of soil loss and flood risk, and combat the natural disasters associated with deforestation, 4) the Grain for Green Program, also known as China’s Sloping Lands Conversion Project, was initiated in 1999 and has enhanced the change of croplands in hilly areas to forests. Recent studies have identified that the implementation of the national ecological restoration projects and soil and water conservation programs has improved ecosystem services like soil conservation, biodiversity protection, water retention (Ouyang et al., 2016; Wang et al., 2016; Lu et al., 2018; Jiang et al., 2020; Wang et al., 2020). However, there are still challenges to implementing the key afforestation/reforestation projects more productively and efficiently for not only local, but also global ecological restoration and sustainability (Ma, 2005; Xu et al., 2006; Chazdon, 2008; Elmqvist et al., 2015). What has been especially challenging is determining how to carry out nature-based afforestation design, planning, and management, with particular consideration of the ecotones with a fragile natural system, intensive human activities and frequent landslide and debris events in western mountainous areas. Such design, planning, management policy, and interventions with specific solutions require a better understanding of the local environmental and social and human dimensions of mountainous contexts. The concepts of nature-based solutions (NBS) (Almenar et al., 2021) and multi-functional landscapes (Peng et al., 2015) have recently emerged with the potential to supply multiple ecosystem services (Peng et al., 2015; Almenar et al., 2021). Silviculturist and landscape ecologists must therefore integrate a landscape science perspective when spatially planning a contemporary strategy, designing new forestry socio-ecological systems adapted to local areas, and addressing the development and sustainability challenges at the regional scale (Liu et al., 2007; Dale et al., 2013; Chopin et al., 2017; Gong and Xie, 2018; Wu et al., 2020). Moreover, agricultural expansion is inevitable in most mountainous areas with fragile environments and heavy populations for food. Consequently, spatial land-use trade-offs may be inevitable among reforestation, agricultural expansion, and other land uses (Gong and Xie, 2018; Liu et al., 2020). Therefore, it is obligatory to maintain crop yields and conserve lands through rational agricultural management. Conservation agriculture has been proven to be a win-win way to meet this goal (Liu et al., 2020). The future planning and management of ecological restoration, including afforestation and revegetation, should focus more on scientific planning to strengthen the nature-related aspects in spatial planning and integrated landscape management (Frank et al., 2012; Zeng et al., 2017; Mann et al., 2018; Jiang et al., 2020; Hersperger et al., 2020; Rieb and Bennett, 2020). The strategic spatial plans and implementation refer to landscape functioning, soil and water conservation agriculture, economic forestry and fruit industry, and contribution of forest landscapes to well-being (Zeng et al., 2017; Gong and Xie, 2018; Jiang et al., 2020; Liu et al., 2020; Gong et al., 2021).

Although the implementation of the national ecological restoration programs can improve ecosystem services (Wang et al., 2016, 2020; ;Ouyang et al., 2016; Lu et al., 2018; ), the approach of optimizing landscape patterns to enhance ecosystem services cannot be ignored because it is easier to implement and less costly (Zhang et al., 2019; Rieb and Bennett, 2020). Increasing the proportion of forestland, simplifying the shape of forestland, and reducing human disturbance on forestland will improve the soil conservation service in BLJW. Integrated landscape management has become a new strategy for landscape governance to address the growing land-use conflicts in response to the multifunctional management of landscapes worldwide (Freeman et al., 2015; Mann et al., 2018).

The soil conservation service is a vital component of regional ecosystem services and human welfare. Our study results provide the scientific basis for afforestation planning and integrated landscape management, especially for the fragile ecotones facing the dilemma of ecological restoration and economic development due to limited land resources. However, there inevitably are some uncertainties due to the complexity of soil erosion, the intricate relationship between soil conservation, and the landscape pattern. First, model outputs, like InVEST-SDR and landscape metrics, always carry uncertainties related to input data due to spatial variability, data availability and many other factors (Diwediga et al., 2018; Sharp et al., 2020). This might be pertinent for the BLJW, especially with the recurrent scarcity of measured data and reference studies. Second, uncertainties in the data sources and methods, especially the InVEST-SDR models and landscape metrics, are still worth discussing, and future applications of soil conservation assessments and landscape pattern design and management should be emphasized. With long-term in situ monitoring, the growing availability of remote sensing data, advanced analytical tools (machine learning, coupled ecological process-hydrological modeling) within a systematic context could improve the quantitative analysis of soil conservation and landscape in BLJW (Gong and Xie, 2018; Zhang and Wei, 2021). Third, it is evident that large-scale afforestation in water-limited areas will inevitably exacerbate water scarcity, while afforestation in energy-limited environments will help to decrease flood risk (Calder, 2007). For example, large-scale reforestation programs in the semi-arid Loess Plateau in China caused substantial streamflow reductions that led to approach water resource limits (Feng et al., 2016). Therefore, afforestation and management must be designed and planned with full consideration of future climate, water consumption and availability, and hydrological regime, as well as topography and watershed size. There are still many challenges that afforestation and multifunctional landscape planning and management face in enhancing ecosystem services and resilience of the mountainous socio-ecological system. More attention should be paid to these challenges in future research.

Forest landscape patterns played a vital role in watershed soil conservation in Chinese western mountainous areas. Nine forestland landscape indices, including PLAND, LPI, ED, LSI, SHAPE_MN, COHESION, DIVISION, SPLIT, and AI, were significantly correlated with soil conservation. Moreover, PLAND, AI, LSI, and SPLIT of landscape, especially for forestland, were more important for soil conservation.

In 1990 and 2014, SLM was superior to OLS for soil conservation regression, which indicated that soil conservation within the sub-basin unit was spatially dependent. However, in 2002, OLS was more appropriate for the regression, which indicated that the forestland landscape pattern had a strong effect on soil conservation. Furthermore, the number of sub-basins with significant spatial correlations between soil conservation and forestland landscape indices, such as PLAND, AI, LSI, and SPLIT, was considerable, that is, there is a scientific basis for promoting soil conservation by optimizing the forestland landscape pattern of the sub-basins. Integrated landscape management can enhance both soil conservation and regional sustainable development, especially in mountainous areas. Future planning and management of ecological restoration, including afforestation and revegetation, should focus more on strategic spatial planning and integrated landscape management to strengthen nature-related aspects with full consideration of future climate, water consumption and availability, and hydrological regime, as well as topography and watershed size.