Investigating Eco-Environmental Vulnerability for China–Pakistan Economic Corridor Key Sector Punjab Using Multi-Sources Geo-Information

Abstract

China-Pakistan economic corridor (CPEC), a critical part of the Belt and Road initiative (BRI), is subjected to rapid infrastructure development, which may lead to potential eco-environmental vulnerability. This study uses multi-source geo-information, and the multi-criteria decision-making (MCDM)-based best–worst method (BWM) to quantify the baseline eco-environmental vulnerability of one key CPEC sector—the Punjab province. The Punjab province is an important connection between northern and southern CPEC routes in Pakistan. In this study, we have established an indicator system consisting of twenty-two influential factors in a geospatial database to conduct eco-environmental vulnerability analysis. The overall setup is supported by a geographic information system (GIS) to perform spatial analysis. The resulting map was categorized into five vulnerability levels: very low, low, medium, high, and very high. The results revealed that the overall eco-environmental health of the Punjab province is reasonably good as 4.64% and 59.45% area of the key sector lies in ‘very low’ and ‘low’ vulnerability categories; however, there also exist highly vulnerable areas, particularly in the proximity of CPEC projects. Although high vulnerability areas constitute a very small percentage, only 0.08% of the Punjab province, still, decision-makers need to be aware of those regions and make corresponding protection strategies. Our study demonstrated that the MCDM-BWM-based EVA model could be effectively used to quantify vulnerability in other areas of CPEC. The findings of the study emphasize that management policies should be aligned with research-based recommendations for ecological protection, natural resource utilization, and sustainable development in regions participating in BRI.

1. Introduction

The China–Pakistan economic corridor (CPEC) is one of the six overland corridors proposed under the Belt and Road Initiative (BRI) [1]. It is a flagship project of BRI, and its successful implementation is vital for social-economic development in the region of Pakistan [2]. It is a massive investment for connecting the southwest of Pakistan (Gwadar) to the northwest of China (Xinjiang province) via an extensive network of rails and roads [3]. This implies that the area of CPEC is huge, including the whole of Pakistan and a part of China. The stakeholders from Pakistan and China are increasingly pressured to build CPEC under the “green” concept to avoid detrimental impacts on the regional eco-systems [4,5,6]. No doubt CPEC will result in the region’s prosperity and socio-economic development [7,8]; however, there is evidence in the growing literature that the operationalization of CPEC might deteriorate the natural eco-system conditions due to increased anthropogenic activities. Wolf [9] has reasoned that CPEC will deteriorate the natural environment of Pakistan on a large scale. Abalakov et al. [10] argued about an increase in the transformation of landscapes within economic corridors and urged for defining the permissible level of the negative impact of anthropogenic activities. Some authors have shown concern over frequent natural disasters, extreme weather events, climate change impacts, and varying geology issues of the CPEC region, which challenges the successful implementation of corridor [11,12,13,14]. Ullah et al. [15] have reported a change in seasonal precipitation in the region of CPEC, which may lead to drought, a threat to food security, loss of biodiversity, and shortage of hydropower [16,17,18]. In addition, researchers have agreed on an increased opportunity of tourism for people of Pakistan and China, but also emphasized sustainable and eco-friendly tourism planning [19,20,21]. One of the main issues in our knowledge is that the interactions among various natural and human systems are not well understood in the CPEC region. This makes CPEC an important study target for eco-environmental vulnerability assessment (EVA), which allows the study of various natural, social, and economic systems in an integrally coupled manner [22].

The concept of “eco-environmental vulnerability” (EV) is generally described as a characteristic of eco-systems that indicates an ecosystem’s tolerance and resilience to external disturbances within a certain region [22]. Research on EVA has remained a growing field of interest among environmental scientists in the past two decades. Numerous studies attempting to investigate the eco-environment conditions and their influencing factors have been conducted in several areas around the globe, such as in China [23,24,25,26,27], Vietnam [28,29], Australia [30], Argentina [31], India [32], Benin [33], and Nepal [34]. Since the eco-environment problems are distinctly distributed in different regions, a uniform model for EVA cannot be constructed. So, the studies have undertaken EVA for their target system by designing an indicator system representing their target system. For instance, Huang et al. [23] conducted EVA for Hainan Island, China, as their target system. Liou et al. [28] established a model based on twelve variables to evaluate spatial–temporal EV in Thua Thien—Hue Province, Vietnam. Chaudhary et al. [34] selected 18 asset-based explanatory variables to target Dordi River Basin as their area of research. Similarly, WANG et al. [35] established a model by combining twelve representative factors to evaluate eco-environment vulnerability in the Yellow River Basin, China as a target system. Furthermore, Huang et al. [36], Xu et al. [37] and Liu et al. [38] are more representative studies that evaluate EV for a watershed, river delta, and mountainous region, respectively, as a target system.

A couple of studies have also addressed the EVA for CPEC [39,40]; however, the results reported in those studies suffered from generalization while attempting to investigate the whole CPEC. More research is needed, which should incorporate local-scale datasets to understand localized eco-environment problems in the CPEC region. Considering this gap, our study chooses Punjab province as a study area for EV investigations. The Punjab province is an essential link of the CPEC route between northern and southern Pakistan and a socio-economically important sector of Pakistan. No previous study, to our knowledge, has undertaken EVA for Punjab province. Conducting EVA of any region provides essential data to policy-making bodies for sustainable development and ecological, environmental management in the corresponding region [41]; therefore, investigations of the CPEC region at a manageable scale are crucial to supplying essential data to decision-makers so that they can make strategies to meet greening and sustainable development targets of the corridor.

Numerous techniques have been adopted in previous studies for EVA, such as spatial principal component analysis (SPCA) [42], the grey theory and fuzzy mathematics [43], analytical hierarchy process (AHP) [36,44,45,46], fuzzy AHP [47], and artificial neural network [24]. The choice of method is guided by the nature of available data and the EVA target system. For example, when a plentiful amount of historical data is available, machine learning techniques (e.g., artificial neural network) can be used by meeting the training and validation requirements of the algorithms [24]. Among all methods, the multiple-criteria decision-making (MCDM) technique AHP has remained by far a commonly used method for performing EVA because of its ability to evaluate a large number of conflicting criteria and alternatives through its organized hierarchy breakdown process; however, AHP is often criticized for its rank reversal problem and complexity. Alternatively, another recently developed MCDM method called the best–worst method (MCDM-BWM) by Rezaei, 2015 [48] has received recognition for its simplicity of use and efficiency [49,50,51]. It has already found its applications in flood risk analysis [52], assessment of solid waste management practices [53], landslide susceptibility mapping [54], and suitable site selection problems, etc. [55]. The advantage of using MCDM-BWM over AHP is that it is a vector-based method and much simpler to use, whereas AHP is a matrix-based method that requires a large number of pairwise comparisons. In addition, MCDM-BWM makes more consistent comparisons; the final weights derived from MCDM-BWM are highly reliable compared to AHP [56]. For this inspirational significance, our study has adopted the use of MCDM-BWM for EVA.

With the rapid development of space and computer technology, the integration of remote sensing (RS) and geographic information systems (GIS) with MCDM has been well-founded and implemented by many researchers and practitioners in a wide range of fields, including nature, society, and economy [57]. The geospatial approach combining MCDM with RS and GIS provides the capacity to conduct investigations at varying spatial and temporal scales. In this paper, MCDM-BWM is integrated with GIS to set up the EVA framework. The main objectives of this work are: (1) propose a flexible EVA framework that would be applicable to other CPEC regions; (2) investigate baseline spatial variation of the EV of key CPEC sector Punjab province in the presence of ongoing and announced CPEC projects [58].

2. Materials and Methods

2.1. The Study Area

Punjab province (located approximately at 30°00 N, 70°00 E) is a key sector of CPEC because it is geographically and the socio-economically significant province of Pakistan (Figure 1). It is a connection between the northern route and southern routes of CPEC. It is Pakistan’s second-largest province in land area, covering 205,344 km². It is also the most populous province of Pakistan. The estimated population is 110,012,442 as of the 2017 census [59]. It can be rightly said as a central province because it borders Islamabad and Khyber Pakhtunkhwa in the north; to the northeast is Azad Kashmir. Baluchistan and Erstwhile Federally Administered Tribal Areas (FATA) are in the west and Sindh in the south. It shares a border with India in the East. Lahore is the capital city of Punjab.

Agriculture is the backbone of the economy of Pakistan and Punjab accounts for 57% of the country’s cultivated area and 69% of cropped area for agricultural production of the country (Government of Punjab, 2013). The Indus River and its four main tributaries Jhelum, Chenab, Ravi, and Sutlej Rivers, traverse Punjab from north to south. This is also one of the reasons that this region has a rich agriculture system irrigated from one of the largest canal-based irrigation systems on Earth [60,61]. During the Monsoon rainfall season, flood is a dominant hazard in this region [62].

Even before of the ideation and launch of CPEC, Punjab already possesses an extensive network of the road (motorways, highways) and railway that serve as the means of trade and transportation. Further, the region is heavily industrialized and consequently undergone rapid urbanization. Given the province’s strategic and socio-economic significance, this region is regarded as one of the key sectors of CPEC.

2.2. Data Sources

The EVA framework adopted in this study is made up of the 22 indicators that best represent the eco-environmental conditions of Punjab; they were chosen after extensive consultation with the relevant literature [36,37,42,44,63,64,65,66]. Once the literature had been reviewed to identify how such frameworks are constructed in similar research, the indicators were chosen per the availability of their data for Punjab province. They are grouped into six broad categories: hydro-meteorology (G1), socio-economics (G2), land resources (G3), topography (G4), hazards (G5), and CPEC projects (G6). These, in turn, are grouped into three overlapping systems: natural, social, and economical. Figure 2 shows that G1, G4, and G5 belong to the natural system, while G3 belongs to both natural and social systems, G2 overlaps both the social and economic systems, and G6 is also multidisciplinary, being important in both the natural and economic systems.

There GIS database consisting of twenty-two indicators established in this study is constructed from six types of data sources:

• Landsat 8 OLI/TIRS imagery;
• Meteorological data;
• Internationally accepted GIS data exchange portals;
• Official reports of the Government of Pakistan (GOP);
• Digital elevation model (DEM) data;
• High-resolution imagery of Google Earth.

The Landsat 8 satellite was launched in February 2013 to ensure a continued supply of data from the Landsat mission. It carries the Operational Land Imager (OLI), which has a spatial resolution of 30 m and a temporal resolution of 16 days [67]. The size of Punjab province is such that 23 Landsat images are needed to cover the entire area. These images correspond to paths 148–152 and rows 37–41 on the standard Landsat Worldwide Reference System (WRS) path/row grid [68].

Table 1. List of Landsat 8 scenes used in the study.

No.	Path/Row	Date (dd/mm/yyyy)	Cloud Cover (%)
1	148/037	10/02/2015	0.94%
2	148/038	10/02/2015	2.66%
3	148/039	10/02/2015	0.59%
4	149/037	16/01/2015	2.31%
5	149/038	16/01/2015	0.03%
6	149/039	21/03/2015	0.2%
7	149/040	21/03/2015	0.01%
8	149/041	05/03/2015	1.89%
9	150/036	07/01/2015	4.81%
10	150/037	08/02/2015	14.98%
11	150/038	08/02/2015	4.71%
12	150/039	08/02/2015	12.08%
13	150/040	08/02/2015	10.03%
14	150/041	08/02/2015	0.21%
15	151/036	19/03/2015	1.7%
16	151/037	30/01/2015	1.73%
17	151/038	30/01/2015	0.06%
18	151/039	14/01/2015	0%
19	151/040	30/01/2015	0.38%
20	151/041	30/01/2015	14.04%
21	152/038	05/01/2015	2.84%
22	152/039	10/03/2015	1.1%
23	152/040	10/03/2015	5.75%

provides a list of Landsat 8 scenes, path/row detail, their acquisition date, and cloud cover percentage. The time for data acquisition is between January to April 2015. In this study, data from 2015 were taken as a priority for baseline evaluation because CPEC was officially launched in 2015. Data for a few variables were not available for that year, so data from later years were included, although special attention was paid while selecting them to ensure that the variables had not undergone a rapid change; examples of these variables are dominant hazards in the region, population, land use and land cover (LULC).

Meteorological data were obtained from the Pakistan Meteorological Department (PMD) by following official procedure [69]. Upon request, PMD supplied figures for monthly mean temperature (°C) and monthly total rainfall (mm) for the period from January 2015 to April 2015.

The third data source includes internationally accepted GIS data exchange portals. They include: (i) a humanitarian data exchange website where population data for various developing countries are hosted and distributed [70], (ii) OpenStreetMap (OSM) data available from the Geofabrik portal [71], and (iii) the regional level railway network and inland water (hydrological network) dataset from DIVA-GIS platform [72].

The fourth source is a credible GOP report titled Integrated Context Analysis (ICA): On Vulnerability to Food Insecurity and Natural Hazards—Pakistan [73]. To review hazards in the region, several other reports from the National Disaster Management Authority (NDMA) of Pakistan are also considered. The reports are accessible at http://cms.ndma.gov.pk/publications (accessed on 4 January 2021).

The fifth source, the DEM, is the digital representation of the land surface elevation with respect to any reference datum [74] and is used to extract topographical variables for the EVA model. The National Aeronautics and Space Administration (NASA) has provided Shuttle Radar Topography Mission (SRTM) digital elevation data for over 80% of the globe. These data are currently distributed free of charge by the U.S. Geological Survey (USGS) earth explorer.

The sixth and last source takes advantage of the computer-assisted GIS digitization technique, in which basemap imagery, such as Google Earth, is used to trace features of interest.

2.3. Data Processing

Figure 3 illustrates the technical flow for the processing of data sources. The processing of the data has been carried out as follows. Google Earth Engine (GEE) (accessible at: https://code.earthengine.google.com/ accessed on 4 January 2021) was used to compute four indicators, namely normalized difference moisture index (NDMI), normalized difference water index (NDWI), normalized difference built-up index (NDBI), and normalized difference vegetation index (NDVI) from Landsat 8 OLI dataset. These four indices are well-defined and reliable in the field of RS [75,76,77,78]. GEE code editor was used to fetch the Landsat 8 TOA collection in the GEE environment “LANDSAT/LC08/C01/T1_TOA/”. In the next step, cloud-free (<20%) scenes from the collection were selected for the defined dates, mosaicked, and clipped according to our area of study. In the later step, selected bands according to indices formulae were processed to obtain NDMI, NDWI, NDBI, and NDVI.

PMD provided the temperature and precipitation data on request. Kriging interpolation was applied to generate raster surfaces of temperature and precipitation variables [79]. This interpolation weighs surrounding measured values to derive a prediction for an unmeasured location. It can be presented mathematically as:

$$\widehat{Z}\left(s_0\right)={{\displaystyle \sum }}_{i=1}^N{\lambda }_{i}Z\left(s_i\right).$$

(1)

where Z(s_i) is the measured value of the i-th location, λ_i is the unknown weight for the measured value of the i-th location, s₀ is the prediction location, and N is the number of measured values.

The population indicator layer was constructed after downloading data from the humanitarian data exchange website. The data were clipped to the extent of the study area. Primary and secondary roads networks were extracted from OSM data after downloading data from the Geofabrik portal. Similarly, railway network and inland water (hydrological network) datasets were downloaded from the DIVA-GIS platform. The Euclidean distance tool was used on the above datasets in ArcMap 10.2.2 to prepare distance from the road network, distance from railway, and distance from hydrological network variables.

Four important indicators, including dominant land cover, flood hazard, landslide hazard, and drought hazard maps, were digitized from the ICA report mentioned in the data sources section. This report is prepared by NDMA and World Food Program, Pakistan (WFP). The dominant land cover is sourced from ESA GlobeCover 2009 and analyzed by FAO Pakistan, whereas flood hazard, landslide hazard, and drought hazard data have been contributed by NDMA.

Furthermore, additional socio-economic variables (which depict anthropogenic and economic activities) such as major cities and CPEC projects were created by the GIS digitization. Google Earth imagery was used as a base visualization to trace features and CPEC official website was used to identify the scope of projects [58]. Euclidean distance tool was applied to convert major cities data to distance from major cities, and buffer distance was computed for CPEC projects variables.

Topography-related variables such as elevation, slope, and aspect are derived from SRTM DEM [80].

Finally, we converted all datasets into a 30 m × 30 m resolution and prepared a GIS database for EVA model implementation. The resolution of Landsat imagery is 30 m, which is taken as reference resolution in this work. A significant amount of data was in a vector format such as roads, railways, major cities, and CPEC projects; however, when variables such as distance from roads, distance from railway, and distance from major cities are developed, the output format was a raster dataset. Special attention has been paid to restrict export data to 30 m raster datasets. This is how the 30 m × 30 m resolution GIS database was prepared for EVA model implementation.

Table 2. The allocation of indicators in each group and their significance in the EVA framework.

Groups	Indicators		Significance in EVA Framework	Maps/1234 Figures
Hydrometeorology (G1)	Normalized difference moisture index (NDMI)	j1	NDMI evaluates the different contents of humidity from landscape elements, especially for soils, rocks, and vegetation. It is one of the important indicators adopted in EVA. The assessment of moisture content is essential to incorporate the health of vegetation perspective in the natural eco-environment of the region [81]. NDMI is calculated using the equation below [78]: NDMI = (NIR − SWIR1)/(NIR + SWIR1) (2)	Figure A1a
	Normalized difference water index (NDWI)	j2	NDWI is also called leaf area water-absent index, which implies the water content within vegetation. The index is derived from the reflectance properties of green vegetation, dry areas (lack of vegetation), and soils by using two near-IR bands. NDWI is calculated using Equation below. ([82]): NDWI = (G − NIR)/(G + NIR) (3)	Figure A1b
	Distance from hydrological network	j3	Surface water is an important driver of natural eco-systems. According to experts, the areas that lack surface water are prone to desertification, forest fire, and landfill sites; therefore, areas farther from the hydrological network are more vulnerable and assigned more weightage in our EVA system.	Figure A1c
	Temperature	j4	Temperature is a vital variable for EVA because it represents the climate conditions of the region. As the greenhouse gases in the atmosphere increase, the temperature will also increase. As a result, there is an increase in the melting of glaciers, which will first lead to extreme flooding followed by a scarcity of irrigation water [83,84]. This affects the vegetation growth of that area; therefore, the higher temperature implies higher vulnerability, and the lower temperature implies low vulnerability.	Figure A1d
	Precipitation	j5	Precipitation is a vital variable for EVA because it represents the climate conditions of the region. Climate change is likely to involve changes in the amount and intensity of precipitation [15]. Changes in precipitation impact the hydrological, ecological, and biogeochemical processes, either directly or indirectly [85]; therefore, precipitation is incorporated in the EVA indicator system established in this study. Areas receiving huge rainfall annually are likely to be hit by extreme events such as floods and thus more vulnerable to eco-system health.	Figure A1e
Socioeconomics (G2)	Normalized difference built-up index (NDBI)	J6	This variable is beneficial for the identification of human settlements and some elements of surrounding construction. It is a rapid and accurate method of mapping urban areas [86]. It can be calculated as: NDBI = MIR − NIR/MIR + NIR (4) Where NIR is the reflectance in the near-infrared band (Landsat OLI-TIRS band5 (0.845–0.885 μm)) and MIR represents the reflectance of the middle infrared band (OLI-TIRS band6 (1.560–1.651 μm)).	Figure A2a
	Population	j7	The population is an important variable for the quantification of vulnerable populations, in particular, vulnerable eco-systems. The high density of population reflects the pressure on natural resources and the intensity of economic activities [87].	Figure A2b
	Distance from the road network	j8	Similar to population, the distance to roads is an indicator in our EVA system to reflect the degree of human involvement. To avoid overestimating the influence of roads, only major roads such as highways, motorways are incorporated in the analysis. The regions nearest to the road network are assigned more weightage and areas farther away are less weightage.	Figure A2c
	Distance from cities	j9	CPEC is primarily a corridor or connection between small and big cities for trade and economic growth; therefore, distance from the city center is indirect consideration of anthropogenic activity and urbanization. In general, cities stand for increased anthropogenic activity, industrial pollution, and fossil fuel exhaust from motor vehicles. Hence, the areas in close vicinity of cities are assigned more weightage in our EVA system.	Figure A2d
	Distance from the railway network	j10	The construction of a new railway as well as upgrading existing railway lines are essential aspects of CPEC. The network of rail and road will boost trade and economic activities in the region. Thus, considering the distance from railway lines is vital for conducting EVA. Areas in closer proximity to railway lines will usually have higher vulnerability.	Figure A2e
Land resources (G3)	Normalized difference vegetation index (NDVI)	j11	The normalized difference vegetation index is generated from the red and near-infrared (NIR) bands by the following simple expression: NDVI= (NIR-Red)/(NIR+Red). (5) It is a popular indicator used in EVA frameworks for incorporating the influence of vegetation and land use. High values of this index are obtained for areas covered by green vegetation.	Figure A3a
	Land use land cover (LULC)	j12	Land use/land cover is an important indicator of how human activities alter the natural land. The land surface characteristics have a direct relationship with surface energy balance, hydrological cycle, and eco-system services [88,89,90,91,92]. As revealed in literature, LULC is the essential indicator used extensively in EVA research [29]. Our EVA model also uses the dominant LULC layer as one of the inputs in the model.	Figure A3b
Topography (G4)	Elevation	j13	The elevation indicator in the EVA model represents the relief and terrain properties of the region. The elevation is a most important topographic characteristic that cannot be neglected in modeling-based research [93].	Figure A4a
	Slope	j14	Slope direction and degree is a basic factor for eco-environmental vulnerability assessment because it impacts on climate conditions of the region. The slope influences hazards such as floods.	Figure A4b
	Aspect	j15	The aspect is a popular indicator used in EVA research by various scientists. It plays a significant role in depicting surface characteristics in environmental research.	Figure A4c
Hazards (G5)	Drought hazard	j16	Drought is a cycling recurring natural event that affects environmental, economic, and social conditions [94]. It has been widely incorporated in vulnerability modeling research. In our EVA framework, more weightage is assigned to an area where drought hazard is “very high”.	Figure A5a
	Flood hazard	j17	Flood is the most prevalent and devastating natural disaster among all-natural disasters that adversely impact human health and natural and artificial environments [95,96]. It has been widely incorporated in vulnerability modeling research. In our EVA framework, more weightage is assigned to the area where flood hazard is “very high”.	Figure A5b
	Landslide hazard	j18	Landslides are the downslope movements of debris, rocks, or earth material under the influence of the force of gravity. The areas where landslides occur frequently are highly vulnerable in terms of deaths, infrastructure damages, and environmental losses [97]; therefore, it is a widely used indicator in vulnerability modeling research in literature. In our EVA framework, more weightage is assigned to the area where landslide hazard is “very high”.	Figure A5c
CPEC projects (G6)	Distance from polygon CPEC projects	j19	The group of CPEC projects represents a truly multidisciplinary system that signifies environmental concerns and economic growth. A total of eight projects are included in the study. Five out of eight can be represented as polygons. These projects include coal power plants, special economic zones, etc. Buffer distance is used to assign more weightage to the area very close to the project.	Figure A6a
	Distance from linear CPEC projects	j20–22	Three out of eight projects can be represented as line features in the GIS environment, such as transmission line, railway, and orange line in-city mass-transit project. Their influence in the EVA model is different from those that can be represented as polygons; therefore, the polygon and linear projects are separated in our EVA model.	Figure A6b–d

presents the group-wise systematics organization of indicators, their significance in the EVA model, and reference to the figures/maps of the indicators’ GIS data layers.

2.4. Methodology

As illustrated in Figure 4, the data gathering and processing in eco-environment research are followed by the determination of the evaluation unit to allow all indicators for fair participation in the EVA model. It is imperative to consider that all the indicators had their own dimension and distribution, making it difficult to compare and operate directly. To overcome this, each data layer of indicators was divided into several classes using manual-interval classification [98] and reclassified by assigning a rating value on a numeral scale from 1, 2, 3, …, 8; where one refers to a positive relationship with eco-environment and eight corresponds to a negative effect on eco-environment. In other words, one can be termed as best and eight can be thought of as worst. For example: in the ‘distance from surface water’ variable, according to experts, the lower vulnerability should be located near a river/water body, and higher vulnerability should be found in the area farther away. Similarly, for the population data layer, as the vulnerability directly relates to anthropogenic activity, the lower vulnerability should be located in the least populous areas, and the higher vulnerability should be located in densely populated areas.

Table A1. Classification and index value assigned to twenty-two indicators of six groups.

Hydrometeorology (G1)					Socio-economics (G2)					Land Resources (G3)		Topography (G4)			Hazards (G5)			CPEC Projects (G6)
j1	j2	j3	j4	j5	j6	j7	j8	j9	j10	j11	j12	j13	j14	j15	j16	j17	j18	j19	j20	j21	j22	Index
value	value	meter	°C	mm	value	people.sqkm	meter	meter	meter	value	description	° angle	meter	degree	description	description	description	meter	meter	meter	meter	no units
No Data	No Data	No Data	No Data	No Data	No Data	No Data	No Data	No Data	No Data	No Data	No Data	No Data	No Data	No Data	No Data	No Data	No Data	4000.01–261,643.781	1000.01–312,381.25	3000.01–215,535.672	500.01–215,535	0
0.6–0.942	−0.768–0.024	0–400	14.95–16.906	14.5–26.896	−0.98–(−0.52)	0.006–2.743	30,000.01–97,583.297	130,000.01–140,910	17,000.01–138,078	0.466–0.811	Irrigated croplands and built-up areas	0–1.61	0–1.61	Flat (−1–0)	Very low	Very low	Very low	3500.01–4000	500.01–1000	2500.01–3000	400.01–500	1
0.332–0.6	0.024–0.209	400.01–1200	16.906–18.091	26.896–37.676	−0.52–0.378	2.743–9.131	18,000.01–30,000	100,000.01–1,300,000	10,000.01–17,000	0.329–0.466	Rainfed croplands and irrigated croplands	1.62–5.79	1.62–5.79	North (0–22.5)	-	-	-	3000.01–3500	-	2000.01–2500	300.01–400	2
-	-	1200.01–2500	18.091–19.082	37.676–52.767	-	9.131–20.08	14,000.01–18,000	70,000.01–100,000	7000.01–10,000	-	Irrigated croplands and bare areas	5.8–12.22	5.8–12.22	Northeast (22.5–67.5)	Low	Low	Low	2500.01–3000	-	1500.01–2000	200.01–300	3
-	0.209–0.301	2500.01–4000	-	52.767–70.554	−0.378–0	20.08–36.505	10,000.01–14,000	20,000.01–70,000	5000.01–7000	0.18–0.329	Irrigated croplands and vegetation	12.23–21.51	12.23–21.51	East (67.5–112.5)	-	-	-	2000.01–2500	-	1000.01–1500	100.01–200	4
0.107–0.332	-	4000.01–8000	-	70.554–86.723	-	36.505–60.229	6000.01–10,000	15,000.01–20,000	3000.01–5000	−0.01–0.18	Rainfed croplands and vegetation	21.55–32.15	21.55–32.15	Southeast (112.5–157.5)	Medium	Medium	Medium	1500.01–2000	-	-	-	5
-	0.301–0.393	8000.01–10,000	19.082–19.909	86.723–101.275	0–0.4	60.229–90.34	4000.01–6000	10,000.01–15,000	1500.01–3000	-	Bare areas & cropland vegetation	32.16–43.08	32.16–43.08	South (157.5–202.5)	-	-	-	1000.01–1500	-	-	-	6
−0.017–(0.107)	-	10,000.01–130,000	-	101.275–119.6	-	90.34–136.876	2500.01–4000	5000.01–10,000	500.01–1500	−0.706–0.01	Bare areas and sparse vegetation	43.09–54.97	43.09–54.97	Southwest (202.5–247.5) West (247.5–292.5)	High	High	High	500.01–1000	-	500.01–1000	-	7
−0.568–(−0.017)	0.393–0.7	13,000.01–153,633.078	19.909–21.975	119.6–151.939	0.4–0.824	136.876–232.685	0–2500	0–5000	0–500	-	-	54.98–81.98	54.98–81.98	Northwest (292.5–337.5) North (337.5–360)	Very High	Very High	Very High	0–500	0–500	0–500	0–100	8

summarizes the division of classes and rating assigned.

When all the sub-indicators, their groups, and target problem of EVA are organized in a structural order, it forms a three-level hierarchy (

Table 3. Weights of six groups variables and twenty-two indicators used for EVA of Punjab province.

First Grade (EVA)	Second Grade (G)	No.	Weights-G	Third Grade (ji)	No.	Weights-i
EVA	Hydrometeorology	G1	0.138	Normalized difference moisture index (NDMI)	j1	0.068
				Normalized difference water index (NDWI)	j2	0.052
				Distance from hydrological network	j3	0.454
				Temperature	j4	0.245
				Precipitation	j5	0.181
	Socio-economics	G2	0.197	Normalized difference built-up index (NDBI)	j6	0.063
				Population	j7	0.591
				Distance from Road Network	j8	0.097
				Distance from Major Cities	j9	0.113
				Distance from Railway Network	j10	0.136
	Land Resources	G3	0.09	Normalized difference vegetation index (NDVI)	j11	0.643
				Land use and land cover	j12	0.357
	Topography	G4	0.176	Elevation	j13	0.597
				Slope	j14	0.182
				Aspect	j15	0.091
	Hazards	G5	0.137	Drought hazard	j16	0.236
				Flood hazard	j17	0.451
				Landslide hazard	j18	0.313
	CPEC Projects	G6	0.262	Dist. from Polygon CPEC Proj	j19	0.347
				Dist. from Transmission CPEC Proj	j20	0.193
				Dist. from Railway CPEC Proj	j21	0.284
				Dist. from Orangeline CPEC Proj	j22	0.176

). The top-grade layer in the hierarchy is EVA itself as a final target; the second-grade layers consist of groups (G1 to G6), while the third-grade layers consist of sub-indicators that belong to each group. In the next step, the MCDM-BWM was used to determine the weights of all variables in each group at the third level of the hierarchy. In this step, each group is taken as a separate case for MCDM-BWM. The weights are given as ‘weights-i’ in the table. Similarly, in the next step, comparative weights were determined for all six groups at the second level of the hierarchy using MCDM-BWM. They are given as ‘weights-G’ in the table.

Explanation of Best–Worst Method (BWM)

The computation steps of the BWM are presented below:

Step 1: Identify a relevant list of criteria. In our study, twenty-two relevant indicators and six main groups were identified to constitute EVA.

Step 2: Choose best (B) and worst (W) criteria for groups and indicators. Indicators belonging to each group are treated as separate cases in BWM. For instance, the first five indicators belonging to G1 are subject to identifying best and worst and weighting exercises. Similarly, the set of indicators belonging to G2 are subject to BWM weighting exercise separately.

Step 3: Using a scale of 1 to 9, conduct a pairwise comparison between best criteria B over all the other criteria. This will result in vector $A_B=a_{B1},a_{B2},\dots \dots a_{Bn}$.

Step 4: Similar to the above, the next step is to conduct pairwise comparison ratings of all the other criteria with the worst criterion (W). This will also result in vector $A_w={(a_{1w},a_{2w},\dots \dots a_{nw})}^T$.

Step 5: Next is to obtain the optimized weights $\left(w_1^{*},w_2^{*},\dots \dots ,w_n^{*}\right)$ for all the criteria.

That is, we obtain the weights of criteria so that the highest absolute variations for all j can be minimized for $\left\{\left|w_B-a_{Bj}{w}_j\right|,\left|w_j-a_{jw}{w}_w\right|\right\}.$ The following minimax model will be determined:

$$\min \max \left\{\left|w_B-a_{Bj}{w}_j\right|,\left|w_j-a_{jw}{w}_w\right|\right\}.$$

(2)

$$s.t. {\sum }_{j }{w}_j=1, w_j\ge 0, \mathrm{for} \mathrm{each} \mathrm{criterion}.$$

(3)

Model (3) is transformed to a linear model and shown as ξ^L,

Subject to:

$$\left\{\begin{array}{c}\left|w_B-a_{Bj}{w}_j\right| \le {\mathsf{\xi }}^{\mathrm{L}*}, all j\\ \left|w_j-a_{jw}{w}_w\right| \le {\mathsf{\xi }}^{\mathrm{L}*}, all j\end{array}\right..$$

(4)

$${\sum }_{j }{w}_j=1,w_j\ge 0, \mathrm{for} \mathit{all} \mathit{j}.$$

(5)

Model (5) can be solved to obtain optimal weights $\left(w_1^{*},w_2^{*},\dots \dots ,w_n^{*}\right)$ and optimal value ξ^L* Consistency (ξ^L*) of attribute comparisons close to “0” is desired [48].

The rest of the methodology includes the following steps:

1. Firstly, the weights of third-grade layers (weights-i) are multiplied with their respective data layers to produce one layer for each group. Mathematically, it is presented as below:

$$G_i={{\displaystyle \sum }}_{j=1}^n{w}_j*j.$$

(6)

where $G_i$ represent the groups from 1 to 6, whereas $w_j$ represent the weight (weights-i) assigned to each sub-indicator in the group. $j=1$ and $j=n$ represent first and last indicators in the group, respectively.

2. In the second step, the linear transformation was applied to the result of Equation (6) as a method of standardization before further processing.

$$X_s=\left(\frac{X_i-X_{min}}{X_{max}-X_{min}}\right)\ast 100.$$

(7)

where $X_s$ is a standardized value, $X_i$ is the original value, $X_{min}$ and $X_{max}$ are lowest and highest original values, respectively. The transformed data were magnified 100 times for ease of processing.

3. Third, the standardized layers were again multiplied with their second-grade respective weights (weights-G) to synthesize the EVA layer. The below equation represents the process mathematically.

EVA = αG₁ + ßG₂ + γG₃ + δG₄ + ηG₅ + θG₆.

(8)

where α, ß, γ, δ, η and θ are weights assigned for groups 1 to 6, respectively. The weights for indicators and weights for groups are presented in Table 3.

We used ArcMap 10.2.2. software for the above model implementation and algebraic computations. The final EVA map was again subject to linear transformation for ease of interpretation of results.

3. Results

3.1. Overall Interpretation of EVA Results

Based on the technical flow adopted, six group variables were combined computationally, and a resultant map of eco-environmental vulnerability was produced. The final map (Figure 5) is categorized into five levels representing very low, low, medium, high, and very high eco-environmental vulnerability regions based on the equal interval classification method [99]. The results show that the overall baseline eco-environment health of Punjab province is good, with 4.64% and 59.45% of the total area is in ‘very low’ and ‘low’ EVA zones, respectively. The very low and low vulnerability zones mostly extend along rivers and green croplands. The rest of Punjab comprises 33.32% medium, 2.51% high, and 0.08% very high vulnerability zones (Figure 6). The medium-level vulnerability primarily extends in the region between major cities, while areas close to major cities bore high vulnerability. Finally, the very high vulnerability zone is very small, covering only 166.90 km² of the total area, and that is mostly the areas in very close proximity of CPEC projects. Given that influence of eight big CPEC projects was incorporated in the analysis, the overall results favor the fast progress of CPEC and negate the false criticism on the life-changing project for the people of Pakistan and China.

3.2. District Wise Comparison of EVA Results

The interpretation of EVA results based on its district-wise distribution in Punjab reveals that Rawalpindi, Multan, Sahiwal, and Lahore are prominent districts that contain large values of ‘very high’ vulnerability zone (Figure 7). One of the reasons is that these districts are an economic hub and connected directly via a route of CPEC; however, neither of these regions are spatially connected, nor exist in similar climatic conditions, so the reasons for high vulnerability are complex. For instance, Rashid et al. [100] reported unchecked urbanized expansion in Rawalpindi, where ecological footprints are higher than national standard values. Further, Shabbir and Ahmad, 2016 [101] found water resources in the Rawalpindi region vulnerable to climatic and socio-economic factors. Similarly, studies for Multan, Sahiwal, and Lahore have reported their problems (e.g., crops productivity, floods, CO₂ discharge, land-use change, and pollution) [102,103,104,105,106,107]. The rest of the 32 districts have zero or negligible values in the ‘very high’ vulnerability zone.

Similarly, Rawalpindi, Dera Ghazi Khan, Multan, and Lahore districts are among those who present large areas of ‘high’ vulnerability zone. Moreover, Dera Ghazi Khan and Rahim Yar Khan are two districts with peak values from the ‘medium’ vulnerability zone. In contrast, Mandi Bahauddin is the district with the least value from ‘medium’ vulnerability, constituting only 0.26% of its total area (Figure 8). The ecological problems in Dera Ghazi Khan, such as rangeland degradation, groundwater pollution, and drastic land use patterns change, are already reported by individuals in literature [108,109,110]. Similarly, floods and droughts are reported problems of Rahim Yar Khan [62]; therefore, our EVA result is similar to a superset containing an essence of all eco-environmental influencing factors.

Bahawalpur region shows a peak under the ‘low’ vulnerability zone, with ‘low’ vulnerability covering 91.27% of its total area. The reason for this result is that the Bahawalpur region has the largest desert in Pakistan, most of the area being barren and uncultivated [111]; therefore, all the chosen indicators collectively did not contribute significantly, and the region yields low values.

Finally, Okara, Hafizabad, and Kasur regions show dominance of ‘very low’ vulnerability zones, covering 26.11%, 42.38%, and 26.30% of their total areas respectively. The reason is that these districts are located far from the main CPEC. Further, the economic activity is not on the same scale as Lahore and Rawalpindi.

To summarize, the overall results show that the eastern part of Punjab province mostly has ‘low’ and ‘very low’ vulnerability zones with very few footprints of ‘high’ and ‘very high’ vulnerabilities. The northern part has portrayed a dominantly ‘medium’ vulnerability zone. Moreover, the northwest and western regions predominately have ‘medium’ level eco-environment severity. The ‘high’ and ‘very high’ vulnerability extends from northwest to north–south, and then along the province’s center to the southern region in a uniform pattern, following the distribution of CPEC projects closely. The result now provides evidence to local district governments to formulate specific strategies to protect and restore the eco-environment in their districts.

4. Discussion

4.1. Significance of the Study and EVA Framework

The eco-system provides the service functions in the form of natural resources and living environments for human beings. The scientific and logical evaluation of the EV offers a theoretical basis for the sustainable development of resourced pressure eco-system. The CPEC aims for regional connectivity through massive infrastructural developments; however, climate change and ecological problems represent potential threats to the smooth progress of CPEC. Few studies concerning the EVA of CPEC are conducted. This research uses twenty-two indicators to construct the EVA model with the support of MCDM-BWM and GIS for analyzing the EV situation in the Punjab region of CPEC. Our study is an addition to scientific literature concerning EVA for CPEC, such as Li et al. [39] and Wu et al. [40]; however, their studies do have some downsides. The paper from Li et al. [39] might have been more convincing if ground truth data had been combined with RS data for better accuracy. Currently, the study has solely relied on RS data. The latter study by Wu et al. [40] mainly focused on indexes relevant to natural systems and failed to incorporate data related to socio-economic systems. Only one index named as human disturbance index used in the EV model is insufficient to extract a representation of socio-economics results. Compared with these studies, our study is superior in terms of adopting a balanced approach of combining indicators from natural, social, and economic systems. Furthermore, limiting our study area to Punjab helped us to derive more targeted results instead of generalizing it for whole CPEC area.

The successful implementation of the framework in this study will reduce the dominance of commonly adopted methods such as AHP and PCA. In addition, our work encourages researchers to consider more sophisticated weighting mechanisms for multi-criteria decision problems. More importantly, the devised framework can be applied to other regions of CPEC, and even other regions of BRI so that decisions and policies can be based on sound information. The results from EVA research provide useful information for policymakers to formulate effective eco-environment protection strategies.

4.2. Validation

Validation of the EVA map of Punjab province is necessary to promote the usefulness of the framework adopted in this study. The detailed validation of the results requires tremendous ground survey efforts, which is beyond the scope of this paper. Alternatively, we have attempted a visual validation scheme using very high spatial resolution Google Earth satellite imagery following the recommendation of Tolche et al. [112] (Figure 9). Two randomly selected sites in each EVA zone are used for validation. When analyzed carefully, Figure 9(i,ii) shows that area designated by our EVA model as very low is abundant in natural resources and lacks anthropogenic disturbance. Figure 9(iii,iv) corresponds to low vulnerability areas. It is visible that natural vegetation is less dense as compared to the very low zone. Further, Figure 9(v,vi) shows medium vulnerability areas. It is visible in satellite imagery that the medium vulnerability zone characterizes areas with bare soil combined with greenery and human settlement. Next, the high vulnerability zone generated by our EVA model has pollution from industries and dense population (Figure 9vii,viii). Finally, the sample area from a very high vulnerability zone (Figure 9ix,x) exhibits the scarcity of natural resources and is subject to urban heat island effect or heatwave due to the high density of impervious surface [113]. The confirmation of EVA results through very high-resolution satellite imagery supports the reliability of our results.

4.3. Proposed Actions for Protecting the Eco-Environment of the Region

The results show a large ‘medium’ eco-environmental vulnerability zone that might be transferred to ‘high’ or ‘very high’ zone amid poor management practices and policies; therefore, favorable ecological protection policies should be practiced in the region to prevent the eco-environment from further degradation. China already changed BRI to Green BRI and established ‘Green Belt and Road Initiative center’ to develop research-based policies, advisory, and actions to build ecologically friendly and green BRI corridors [114]. Similarly, recognizing that overlooking the sustainability of eco-system services, the goals of CPEC investment cannot be met in isolation. Pakistan–China governments have decided to turn CPEC into a green model BRI to protect and conserve the region’s natural environment [6]. Accordingly, the GOP has launched the “Clean Green Pakistan” project to ramp up efforts to combat climate change and restore natural habitats. In the first phase, 13 cities from Punjab are also included in the project. There is a need for strengthening the knowledge and practices among the communities about eco-environmental protection. There should be more effective green policies such as promoting the use of electric vehicles and planting forests whenever trees in a region are cut due to economic development such as road construction. Eco-environment preservation rules must be defined for industries, and authorities should ensure the stringent implementation of policies.

4.4. Prospects of the Study

Although our research has strong policy implications, there are still some limitations in the current study, and some suggestions on future related research are proposed. The indicator system developed for EV evaluation can be improved by considering more suitable indicators. For example, more types of socio-economic data, as suggested by Batar et al. [115], can be included in the study. Their suggested parameters are GDP, per capita income, employment ratio, sex ratio, electric power usage, employability, human development index, and funds for science and technology. Similarly, other datasets such as net primary productivity (NPP), detailed LULC, and updated facts and figures about CPEC should be part of future studies to achieve more actionable results in the region; however, these types of datasets are either not produced on scale or are closed source, so their inclusion is subject to accessibility of these datasets. In future studies, the indexes can be further screened and applied via a multicollinearity test [116]. Further, this study will be expanded to other parts of the CPEC project. We also plan to extend our study to different years to investigate the spatial and temporal variation of the EV. The participating governments should encourage data-sharing initiatives such as Digital Belt and Road (DBAR) program and the Big Earth Data Science Engineering Project (CASEarth) for facilitating the research and development in the region [117,118].

5. Conclusions

This research concentrates on the EVA of Punjab province, one of the important regions in CPEC. By fully considering the eco-environmental conditions and availability of datasets in Punjab, the present study chose twenty-two indicators that are broadly categorized into six groups to construct its EVA model. The model is supported by GIS for spatial analysis and MCDM-BWM for calculating weights of relative importance for indicators. The indicators are constructed from RS data, meteorological data, published reports, DEM, and CPEC projects data. The indicators and their weights are computationally combined to synthesize the EV map of Punjab province. The resultant map of EV in Punjab is classified into five levels, namely very low, low, medium, high, and very high vulnerability levels. The results reveal 4.64%, 59.45%, 33.32%, 2.51%, and 0.08% of the total area in each vulnerability level, respectively. Furthermore, there are 10 districts where the ‘medium’ vulnerability zone is dominant and 26 districts where ‘low’ vulnerability zone is dominant. Similarly, there are 12 districts where ‘very low’ vulnerability zone portrays high values in parallel with other vulnerability zones values. The EVA of Punjab province reflects the comprehensive effect of natural factors and anthropogenic activities booming due to CPEC realization.

This study gives an opportunity to researchers by providing ideas and a basic framework to research EVA in BRI countries while looking at a wide range of scales from local, regional, and country levels. The potential future extension of this study is to update eco-environmental vulnerability status after every five years and integrate more detailed information on CPEC projects in the EVA framework.

The countries participating in BRI in general and countries advancing CPEC (Pakistan and China) in particular should devise a mechanism to publicly share datasets such as high-resolution satellite RS data for construction and validation of indicators; socio-economic datasets such as financial self-sufficiency rate, human development index, per capita income, scientific research funds, etc.; climate data such as temperature, precipitation, etc. to allow scientists to expand the EVA framework by incorporating more indicators. Naboureh et al. [119] already reported the challenges due to the non-availability of detailed reference LULC data in the BRI region. The timely availability of the dataset will not only ensure progress in EVA research, but other fronts such as hydrological, geological, agriculture, and the political regime will also benefit. This way, the funds allocated for research can be utilized fully, and more actionable insights can be sought from research activities. The continuous monitoring of the eco-environment status of all the regions of BRI will ensure that the success of the BRI continues to be sustainable and that the world can benefit from this massive investment.