2.2.1. Identification of Significant Ecological Resistance Indicators for Construction and Development
Under this component, the significant indicators include ecological elements, importance, and resilience that describe the ecological resistance to construction and development in mountainous areas, as suggested by Peng, Ma, Du, Zhang, and Hu [
14], were investigated using the Delphi method. In general, the Delphi method is a method in which the intuitive idea of the participant was used based on a structural survey [
46]. This method is an iterative process that is designed to achieve a consensus among a group of experts on a specific issue, and it is one of the most effective means for participants (experts) to identify criteria or indicators. Barzekar et al. [
47] claimed that the Delphi method is an excellent way to generate a consensus of experts’ opinions when solid scientific data are unavailable. It carries quantitative and qualitative results and has underneath its explorative, predictive, even normative elements [
48].
To implement the Delphi method, a variety of indicators that represent ecological elements, importance, and resilience under three criteria (ecological structural, functional, and dynamics resistance) for construction and development in mountainous areas was first reviewed from a wide range of sources, particularly research papers (
Table 1) and a questionnaire was then prepared for expert interviewing. In this study, 10 experts from the Institute of Geography, Vietnam Academy of Science and Technology, University of Science, Vietnam National University, and Thai Nguyen University of Education were invited to participate in identifying indicators. The selected experts are experienced in different fields of their areas of interest. Dalkey and Helmer [
49] mentioned that a Delphi group possesses the largest confidence when the number of experts is at least 10.
In practice, the prepared first-round questionnaire was sent to all experts for ecological resistance indicators identification and weighting. Each expert was asked to indicate a degree (scale of 1 to 5) to which they agree with indicators as follows:
Indicator is highly irrelevant to the selected criteria,
Indicator is likely irrelevant to the selected criteria,
Indicator is more or less relevant to the selected criteria,
Indicator is likely relevant to the selected criteria, and
Indicator is highly relevant to the selected criteria.
After receiving the first-round response, the information was summarized, collated, categorized, and tabulated into the second questionnaire for the second round. Then, the second questionnaire that incorporated a feedback report was redistributed to experts. They were asked again to give their ratings. The goal of the second round is to attain stability or consensus of experts’ responses. The Delphi procedure is completed when the consensus or stability is reached. The Delphi method finishes when all questionnaire items are either accepted or rejected or over 75% of the questionnaire items have their rating variant values being less than 15% [
61]. The rules for conducting the Delphi method are given in
Table 2.
Where, rating mean (qi) represents the mean of the ratings for questionnaire item qi, rating variant (qi) represents the ratio of experts who change their ratings for qi, and Q is the quartile range.
2.2.2. Ecological Suitability Evaluation for Construction and Development
An integrated ecological resistance (IER) model was used to produce the ecological suitability map for construction and development in this study. The IER model, which is based on the spatial ecology in landscape ecology, was proposed by Peng, Ma, Du, Zhang, and Hu [
14]. Under the IER model, ecological elements, importance, and resilience are used to characterize the resistance of the ecological structure, function, and dynamics for construction and development, respectively, then an IER index is calculated using these three resistances. Actually, the IER index has a negative relationship with the level of ecological suitability, and the highest ecological suitability implies the lowest IER index.
The IER index was calculated using the following equation:
where S is the ecological structural resistance measured by ecological elements, I refers to the ecological functional resistance characterized by ecological importance, and R is the ecological dynamics resistance represented by ecological resilience. In addition, ωs, ωi, and ωr are the weights of three components. Assuming that structure, function, and dynamics are interdependent for the entire ecological environment, equal weight values were applied here to these resistances. The derived IER index was then applied to classify ecological suitability for construction and development into 5 levels using the natural breaks method.
In this study, ecological structural, functional, and dynamics resistances were separately assessed based on the identified indicators from experts as mentioned in
Section 2.2.1. Details of three ecological resistances are summarized below.
(1) Ecological structural resistance assessment
Ecological elements create the ecological structural resistance included (1) elevation, (2) slope, (3) geological hazard frequency, and (4) distance to fracture zones. The weight of each element was determined using the analytical hierarchy process (AHP) from 10 experts. This method is useful for obtaining a single assessment value based on different indicators or criteria [
63]. Herein, the standard numeric scale used for the AHP is on a 1–9 scale that lies between “equal importance” to “extreme importance”, the value 1 indicates equal importance, while value 9 indicates that one indicator is extremely important than others [
64]. To implement the AHP, the Super Decisions software (version, Manufacturer, City, US State abbrev. if applicable, Country) was used to calculate the weight and the inconsistency index. If the inconsistency index is less than 0.10, no correction of judgments is needed [
65].
To prepare ecological elements for ecological structural assessment, elevation and slope were first extracted from DEM and then reclassified into 5 classes using the natural breaks method. By using the natural breaks method, similar values are grouped and therefore the differences between classes are maximized [
66]. The resistance level has a positive relationship with elevation and slope classification. Meanwhile, the spatial characteristics of the geological hazard frequency, which was represented by historical landslides and flash flood records of Hung [
67], was first extracted using the kernel method for surface interpolation under ESRI ArcMap software and reclassified into 5 classes using the natural breaks method. The resistance level of this indicator has a positive relationship with its classification. In the meantime, the distance to fracture zones, which are in crustal strata along where adjacent rocks have been displaced, was extracted based on a Bac Kan geological map using the Euclidean distance method under ESRI ArcMap software (version, Manufacturer, City, US State abbrev. if applicable, Country) and reclassified into 5 zones (0–200 m, 200–500 m, 500–1000 m, 1000–2000 m, and greater than 2000 m). The resistance level of this indicator has a negative relationship with its classifications.
Finally, the ecological structural resistance was assessed using the simple additive weighting (SAW) method (Equation (2)), and it was classified into 5 levels using the natural breaks method. The SAW method is one of the most widely used models for creating composite maps in GIS [
68]:
where
is the total score,
is a normalized weight of ith indicator, and
is the score of the ith indicator (
i = 1, 2, …, n).
(2) Ecological functional resistance assessment
The ecological functional resistance is described by different ecological processes, as well as the relationship between human demands and the natural ecosystem [
69]. The ecological functional resistance for construction and development that was assessed through ecological importance included (1) biodiversity protection, (2) water retention, and (3) soil conservation.
To assess ecological importance for ecological functional resistance assessment, the relative importance of biodiversity protection was classified into 5 levels according to land use types and associated ecosystem services (
Table 3). The functional resistance level caused by this indicator has a positive relationship with its important classification. Meanwhile, the relative importance of water retention was classified based on the experiment result of Hai [
70], which focused on the ability of water retention from different vegetation types and structures. The ability of water retention was assessed through a variety of indicators including soil humidity at different soil depths (
Table 4) and rainfall interception capacity (
Table 5). The final relative importance of water retention that links with vegetation types was classified into 5 levels (
Table 6). Like the previous indicator, the functional resistance level caused by this indicator has a positive relationship with its important classification.
In the meantime, soil conservation was assessed based on the difference between potential and actual soil loss as suggested by [
14,
57]. As the amount of soil conservation increases, the relative importance increases. In this study, actual soil loss was assessed using the universal soil loss equation (USLE) [
71] as:
where:
A is actual soil loss (tons ha−1 year−1),
R is rainfall erosivity (MJ mm ha−1 h−1 year−1),
K is soil erodibility factor (tons h (MJ mm)−1),
LS is topographic factor (dimensionless),
C is cropping management factors (dimensionless), and
P is supporting practice factor (dimensionless).
Brief information on relevant USLE factors was separately summarized below.
(a) Rainfall erosivity factor (R). Rainfall erosivity is defined as the product of the total kinetic energy of the rain and the maximum 30 min rainfall intensity [
71]. Many methods can be used to calculate the annual rainfall erosivity factor [
72,
73,
74]. The process to calculate the R factor is rather complex and requires long-term data collection. In this study, the R factor was estimated based on a local equation that was introduced by Nguyen [
75] using analysis of the rainfall data over 54 years from 253 meteorological stations throughout Vietnam as follows:
where P is the annual precipitation (mm).
In practice, spatial interpolation of annual precipitation based on the existing datasets from 1958 to 2018 was carried out using the inverse distance weighted (IDW) technique to generate the rainfall distribution layer and the derived result was used as the input data to produce rainfall erosivity using Equation (4) as shown in
Figure 3a.
(b) Soil erodibility factor (K). Soil erodibility represents the effect of soil properties and soil profile characteristics on soil loss. Erodibility is a function of soil texture, organic matter content, and permeability [
76]. In this study, K values from different soil types in northern highlands of Vietnam [
77] and Ba Be Lake basin [
78] as a summary in
Table 7 were applied to generate the soil erodibility (
Figure 3b).
(c) Topographic factor (LS). The topographic factor includes slope length (L) and slope steepness (S). The slope has a major effect on soil erosion rate. The higher the slope, the higher the velocity of overland flow, thus increasing the shear stresses on the soil particles. As slope length increases, the overland flow and flow velocity also steadily increase, leading to greater erosion forces applied to the soil surface [
76]. In this study, the equation of Moore and Burch [
79], which was adopted and developed by [
80], was used to calculate the LS factor (
Figure 3c) as:
where
FA is flow accumulation, cell size is the size of the DEM data, slope angle is in degrees (°), and 0.01745 is a parameter to convert degrees to radians.
m and
n are assigned 0.5 and 1.3 respectively as recommended by [
80].
(d) Cropping management factor (C). This factor is the second most important factor that controls soil erosion risk [
81], and it reflects the effect of cropping and management practices on the soil erosion rate [
82]. Generally, the C factor ranges from 0 to 1. The C value equals to 1, it indicates no cover present and the surface is treated as barren land, whereas the C value is close to zero (0), which indicates high plant cover. The C factors were adapted to the Bac Kan land use map from literature reviews for northern Vietnam [
76,
77,
78] (
Table 8 and
Figure 3d).
(e) Supporting practice factor (P). This factor represents the effects of practices to the amount and rate of the water runoff and thus influence the amount of erosion when the supporting practice is high, the value of the P factor is low. It also expresses the agricultural practice in areas such as contour cultivation and arable land terrace [
76]. In this study, the supporting practice factor for A Sap basin [
83] was applied for Bac Kan province, and it was calculated based on group of land use types and the slope as suggested by [
84] (
Table 9 and
Figure 3e).
Meanwhile, the potential soil loss of the landscape was computed as the product of R, K, and LS [
77], and soil conservation was then extracted by subtraction operation between the derived potential and actual soil loss estimation (
Figure 4).
Finally, the ecological functional resistance was assessed using the SAW method with equal weights and then was classified into 5 levels using the natural breaks method.
(3) Ecological dynamics resistance assessment
The ecological dynamics resistance applies the concept of ecological resilience based on ecological circulation theory and taking great human disturbance into account for representing the temporal dynamic characteristics of a natural ecosystem’s self-organization and self-update using the three aspects of resistance, exposure, and interference [
14]. These three aspects correspond to specific indicators include (1) vegetation stability (S
1), (2) ecological sensitivity (S
2), and (3) social disturbance intensity (S
3), respectively. Thus, the ecological dynamics resistance (R) can be expressed as:
Brief information about relevant indicators of ecological resilience for ecological dynamics resistance assessment was separately summarized below.
(a) Vegetation stability (S
1). Vegetation stability shows the vegetation capacity to resist interfering and sustain its condition with environmental variations [
14]. In this study, vegetation stability was quantified through the variation amplitude of vegetation using the normalized difference vegetation index (NDVI) from MODIS product (MOD13Q1) with support of the TIMESAT 3.3 software (Manufacturer, City, US State abbrev. if applicable, Country). Actually, vegetation stability is the reciprocal of variation amplitude of vegetation. The lower the variation amplitude, the greater the vegetation stability, and therefore the higher the resistance to development activities [
14]. In practice, the NDVI dataset of 207 scenes (23 scenes/year) between 2009 and 2018 was used to quantify the variation amplitude of vegetation. The quality of vegetation data was here used to assign the weight for adjusting the NDVI dataset based on the pixel reliability index (
Table 10).
After weighting vegetation data quality, the Savitsky–Golay filtering algorithm [
86] was applied to enhance the quality of the NDVI dataset. At the last step, the variation amplitude of vegetation was extracted from seasonality data and the amplitude value was then classified into 5 classes to represent the vegetation stability using a natural breaks method. The dynamic resistance level caused by this indicator has a positive relationship with its classification.
(b) Ecological sensitivity (S
2). Ecological sensitivity reflects the degree to which human activities and natural changes impact on the ecosystem, as well as the degree to which regional ecological and environmental problems can possibly occur [
87]. In general, ecological sensitivity is assessed using distance from important ecological functional areas or key ecological patches [
88]. A shorter distance from such areas will increase the higher risk for ecosystem degradation when disturbances happen.
In this study, the important ecological function area in Bac Kan province that includes Ba Be National Park, Kim Hy Natural Reserve, and Nam Xuan Lac Species and Habitat Conservation Area from Bac Kan Department of Natural Resources and Environment in 2017 were first chosen to create distance surface data using the Euclidean distance method under ESRI ArcMap software (Manufacturer, City, US State abbrev. if applicable, Country), and it was then classified into 5 levels using the natural breaks method. The dynamic resistance level caused by this indicator has a positive relationship with its classification.
(c) Social disturbance intensity (S
3). Social disturbance intensity refers to the spatial differentiation in the intensity of human activity in the study area, as expressed by spatial population density [
14]. It was generated at the district level using the population data obtained from the People’s Committee of Bac Kan Province in 2018, and it was further classified into 5 classes using the natural breaks method. The dynamic resistance level caused by this indicator has a negative relationship with its classification.
Finally, the ecological dynamics resistance was assessed using Equation (6) and was classified into 5 levels using the natural breaks method.