Atmospheric Environment Vulnerability Cause Analysis for the Beijing-Tianjin-Hebei Metropolitan Region

Assessing and quantifying atmospheric vulnerability is a key issue in urban environmental protection and management. This paper integrated the Analytical hierarchy process (AHP), fuzzy synthesis evaluation and Geographic Information System (GIS) spatial analysis into an Exposure-Sensitivity-Adaptive capacity (ESA) framework to quantitatively assess atmospheric environment vulnerability in the Beijing-Tianjin-Hebei (BTH) region with spatial and temporal comparisons. The elaboration of the relationships between atmospheric environment vulnerability and indices of exposure, sensitivity, and adaptive capacity supports enable analysis of the atmospheric environment vulnerability. Our findings indicate that the atmospheric environment vulnerability of 13 cities in the BTH region exhibits obvious spatial heterogeneity, which is caused by regional diversity in exposure, sensitivity, and adaptive capacity indices. The results of atmospheric environment vulnerability assessment and the cause analysis can provide guidance to pick out key control regions and recognize vulnerable indicators for study sites. The framework developed in this paper can also be replicated at different spatial and temporal scales using context-specific datasets to support environmental management.


Introduction
With the booming industrialization and urbanization development, the relationship between environmental quality and human society has come into obvious disharmony [1,2]. The continuous deterioration of atmospheric quality in urban cities, especially in developing countries, is one of the top global crises, engendering severe human health, environmental, social, and economic impacts [3,4]. One of the major future challenges will be to restrain short-and long-term atmospheric environment degradation to fulfil the needs of the fast-growing human population and changes in climate, where effective management and conservation have a key role. In this context, research on efficient recognition of dominant influencing factors and vulnerable areas according to specific circumstance is urgently needed. As an effective means of characterization and identification of the extent to which a system is susceptible to damage from natural and anthropogenic disturbance, atmospheric environment vulnerability assessment is an effective solution to this issue.
To date, most studies about vulnerability assessment have focused on urban groundwater contamination [5,6], climate change [7,8], human ecology [9,10], agriculture [11,12], coastal flooding [13], soil [14,15], socioeconomic issues [16,17], wetland biota [18], transportation networks [19], etc. Berrouet et al. proposed a conceptual framework for the vulnerability assessment of socio-ecological systems to communicate with practitioners and policy makers and identify and and so does the sensitivity index. On the contrary, the adaptive capacity has a negative relationship with atmospheric environment vulnerability, indicating that a higher adaptive capacity leads to lower vulnerability. The spatial change of atmospheric environment vulnerability in this paper can be defined as the change in vulnerability of different cities. The spatial analysis aims to pick out the most stressed area in the study site. The technical route of atmospheric environment vulnerability assessment is presented in Figure 1. correlated with atmospheric environment vulnerability, meaning higher exposure results in increased atmospheric environment vulnerability, and so does the sensitivity index. On the contrary, the adaptive capacity has a negative relationship with atmospheric environment vulnerability, indicating that a higher adaptive capacity leads to lower vulnerability. The spatial change of atmospheric environment vulnerability in this paper can be defined as the change in vulnerability of different cities. The spatial analysis aims to pick out the most stressed area in the study site. The technical route of atmospheric environment vulnerability assessment is presented in Figure 1.

AHP Method for Index Weights
Multi-criteria decision analysis (MCDA) is an advanced tool to support decision-making that takes several different aspects into account at the same time like AHP, ordered weighted averaging (OWA), weighted arithmetic average (WAA), goal programming (GP), weighted geometric average (WGA), technique for order preference by similarity to an ideal solution (TOPSIS), multi-attribute value theory (MAVT), etc. The purpose of MCDA is to assist a decision maker through the decision-making process via explicit formalized modelling. The AHP-based weight method is the MCDA tool selected in this paper.
AHP was developed by Thomas L. Saaty in the 1970s and has been widely used since then [33][34][35][36]. It is particularly suitable for spatial decision analysis with multiple criteria. AHP offers several advantages over more traditional techniques of weight determination. It can turn multi-objective and multi-criteria problems into a single target with multi-level issues. It is not necessary to ensure statistical generalization and the use of pairwise comparisons can reduce the cognitive burden of index importance. AHP builds a hierarchy (ranking) of decision items using comparisons between each pair of indicators, expressed as a matrix. Paired comparisons produce weighting scores that measure how much importance indicators and criteria have. AHP uses the judgments of decision makers to form a decomposition of atmospheric vulnerability into hierarchies with a combination of qualitative and quantitative methods.
We select 15 indices that are divided into three project layers: exposure, sensitivity, and adaptive capacity, and the three project layers are grouped into a layer for the goal index: atmospheric environment vulnerability. We applied a questionnaire survey to collect the data that are necessary for AHP. The survey was designed to determine the weights of the atmospheric environment vulnerability indices and taken via interview from 10 atmospheric specialists, environment management specialists, and city management experts. The experts were chosen if they (a) have at least 10 years of working/research experience in a related domain, (b) possess sufficient knowledge in atmospheric environment research, and (c) have a relevant educational background. Since it is difficult to gather all experts together for a brainstorming session or free-flowing discussions in a group, we implemented semi-structured interviews to allow more freedom of conversation. The procedure includes: 1. The objective of the interviews: determining the weights for the 15 indices in the atmospheric environment vulnerability assessment.
2. The major questions include: (2.1) What is the contribution rank of exposure, sensitivity, and adaptive capacity for atmospheric environment vulnerability? (2.2) What is the contribution rank of "average annual concentration of PM 2.5 ", "average annual concentration of PM 10 ", "average annual concentration of SO 2 ", "average annual concentration of NO 2 ", "days of AQI equal to or better than grade II", "days of AQI equal to and worse than grade V" for exposure? (2.3) What is the contribution rank of "population density", "percentage of vulnerable groups", "average annual rainfall", "average annual wind speed", "average annual relative humidity" for sensitivity? (2.4) What is the contribution rank of "proportion of secondary industry", "motor vehicle population", "coal consumption", "percentage of urban greenery coverage" for adaptive capacity? (2.5) What is the total contribution rank of the 15 indices for atmospheric environment vulnerability?
After summarizing answers from interviewed experts, the returned data were applied to make pair-wise comparisons and judgement matrix in the AHP process (shown in Table S1). Thus, the contribution of each index to a higher layer is evaluated and the weight of each index in each hierarchy can be achieved (shown in Appendix A). The obtained weights of the 15 indices were then sent to the 10 experts by e-mail to validate the weighting results. The data obtained from the experts were further analysed by manual content analysis. The final weights of the 15 indices are shown in Table 2. All the consistency ratio (CR) values of each matrix are lower than 0.10, which demonstrates the suitability of the calculated weights. The Geometric Consistency Index (GCI) was also used to measure the consistency of the matrices and the results were all less than the threshold determined by Aguarón and Moreno-Jiménez [37,38]. The results verified the consistency achieved from the CR.

Fuzzy Theory for Index Evaluation
The framework relationship between exposure, sensitivity, and adaptive capacity is inexact. Meanwhile, the concept "vulnerability" is fuzzy and not exact; there is no common idea about the exact threshold to distinguish between vulnerable and not vulnerable. For these reasons, it is suitable to use a fuzzy synthesis evaluation method that can handle the uncertainties and aggregate indicators as well as construct a general atmospheric vulnerability index in atmospheric environment vulnerability assessment.
A fuzzy theory is a mathematical theory that values fuzzy set with membership function in the real unit interval [0,1] [39]. A fuzzy set is a pair x, µ A (x): x ∈ X, in which µ A : x→[0,1]. Fuzzy logic operations include fuzzification, combination, and defuzzification. Fuzzification is a process to transfer all input values into fuzzy membership functions. Combination is to execute all applicable rules in the rule base to compute the fuzzy output functions. Defuzzification is to defuzzify the fuzzy output functions to get "crisp" output values (shown in Appendix B).
It is a key step to determine membership function in fuzzy evaluation. We used a truncated trapezoidal membership function and the vulnerability level is divided into a 1-5 scale on which 1 means potential and 5 means highly vulnerable. As the selected vulnerability indices have positive and negative influence, the positive and negative variables are converted to fuzzy sets by means of positive and negative membership functions, respectively. For a positive impact index, its membership functions are calculated by Equations (1)-(3), and the membership functions of a negative impact index are calculated by Equations (4)- (6).
where µ i (x) is the membership of the atmospheric environment vulnerability index at i th (i = 1, 2, . . . ,5) vulnerability class; A i is the threshold value of the i th atmospheric environment vulnerability class. Establishing an appropriate criterion for different atmospheric environment vulnerability levels is essential. As there is no standard to determine the threshold value at each vulnerability class, we built an atmospheric environment vulnerability assessment index system and classification standard ( Table 3) according to the relevant evaluation criteria [40][41][42][43]. The fuzzy membership calculation results for Beijing's atmospheric environment vulnerability in 2015 are shown in Table 4.

Atmospheric Environment Vulnerability Aggregation
To achieve general atmospheric environment vulnerability results, the first step is to calculate the membership for the three project layers. Taking the exposure layer as an example, for the study site x, the memberships for five scales are calculated by Equation (7), as follows: With the same method as above, the memberships of five scales for the project layer sensitivity (V x,S ) and adaptive capacity (V x,AC ) can be achieved. The memberships of five scales for the goal layer atmospheric environment vulnerability are obtained to aggregate the memberships of three project layers by Equation (8): According to the maximum membership principle, the final atmospheric environment vulnerability result is the maximum value among the memberships for five scales.

Study Area
The BTH region (shown in Figure 2) is located in the northwest part of the North China Plain (36 • 05 -42 • 37 N, 113 • 11 -119 • 45 E), which contains two municipalities (Beijing, Tianjin) and one province (Hebei); the total residential population was 111 million at the end of 2015. The total land area of the BTH region is 218,000 km 2 and occupies 2.3% of the Chinese territory. The BTH region is one of the three hotspot areas of air pollution control according to the Twelfth Five-Year Plan for Air Pollution Control in Key Regions [44]. This area has a semi-humid continental monsoon climate and is in the North Temperate Zone.

Data Resource
The following provides the data and sources used in this study. We collected 15 indices of the 13 cities from 2013 to 2015, yielding a total of 585 values. All the input data used in this paper can be divided into three categories: environmental data, meteorological data, and socioeconomic data. The environmental data are derived from the Environment Quality Bulletin [52][53][54][55][56][57][58][59][60], published by Environmental Protection Bureau (EPB) for 13 cities in the BTH region. The meteorological history data of 13 cities are extracted from the China meteorological data sharing service website.

Data Resource
The following provides the data and sources used in this study. We collected 15 indices of the 13 cities from 2013 to 2015, yielding a total of 585 values. All the input data used in this paper can be divided into three categories: environmental data, meteorological data, and socioeconomic data. The environmental data are derived from the Environment Quality Bulletin [52][53][54][55][56][57][58][59][60], published by Environmental Protection Bureau (EPB) for 13 cities in the BTH region. The meteorological history data of 13 cities are extracted from the China meteorological data sharing service website. Socioeconomic data used in this paper are from the China City Statistical Yearbooks [45][46][47]

Results Description
To better depict the spatial and temporal distribution of atmospheric environment vulnerability in BTH region, each assessment index is divided into five levels from 1 to 5, standing for potential vulnerable, weakly vulnerable, moderately vulnerable, strongly vulnerable, and extremely strong vulnerable, respectively. A higher level presents higher impacts on atmospheric environment vulnerability. The assessment results are shown in

Sensitivity
Figure 4c shows that the third-level (53.84%) and fourth-level (30.77%) sensitivity indexes account for a major proportion of the BTH region. This suggests that the atmospheric environment of most cities is influenced by perturbations or stressors. The longitudinal comparisons of the three years showed that the sensitivity index also had an obvious spatial discrepancy. The sensitivity degree for Tianjin and Hengshui improved during the research period. On the contrary, the sensitivity levels of Langfang, Handan, Chengde and Zhangjiakou were aggravated. showed that the sensitivity index also had an obvious spatial discrepancy. The sensitivity degree for Tianjin and Hengshui improved during the research period. On the contrary, the sensitivity levels of Langfang, Handan, Chengde and Zhangjiakou were aggravated.  Figure 4c shows that the third-level (53.84%) and fourth-level (30.77%) sensitivity indexes account for a major proportion of the BTH region. This suggests that the atmospheric environment of most cities is influenced by perturbations or stressors. The longitudinal comparisons of the three years showed that the sensitivity index also had an obvious spatial discrepancy. The sensitivity degree for Tianjin and Hengshui improved during the research period. On the contrary, the sensitivity levels of Langfang, Handan, Chengde and Zhangjiakou were aggravated.  Figure 5c shows that the adaptive capacity status in the BTH region is mainly in the first level (30.77%) and third level (30.77%). This indicates that the general adaptive capacity index in the BTH region is relatively positive, with conspicuous spatial heterogeneity generally declining from central cities to southern cities and then to northern cities. The longitudinal comparison from 2013 to 2015 showed that the adaptive capacity situation in the three years presented almost similar patterns in eight cities, while five cites (Shijiazhuang, Tangshan, Handan, Xingtai, and Zhangjiakou) improved. This implies that changes in related adaptive capacity attributes were minor in the three years and half of the cities presented a relatively lower adaptive capacity.

Aggregated Atmospheric Environment Vulnerability
The assessment results of aggregated atmospheric environment vulnerability are shown in Figure 6. Figure 6c indicates that, in terms of the general atmospheric environment vulnerability, the fourth level (53.84% of 13 cities) took up the largest proportion of the BTH region in 2015. In 2015, the atmospheric environment vulnerability exhibited obvious spatial discrepancy. The general trend shows the vulnerability level increasing from northern to southern regions. In the northern regions, the atmospheric environment in Zhangjiakou, Chengde and Qinhuangdao, at the second level, can be described as stable, with a high capacity to recover from disturbance. In the central regions, Beijing, Tianjin and Cangzhou represented moderately vulnerable atmospheric vulnerability, which indicated the medium anti-interference capacity of the atmospheric system.

Aggregated Atmospheric Environment Vulnerability
The assessment results of aggregated atmospheric environment vulnerability are shown in Figure 6. Figure 6c indicates that, in terms of the general atmospheric environment vulnerability, the fourth level (53.84% of 13 cities) took up the largest proportion of the BTH region in 2015. In 2015, the atmospheric environment vulnerability exhibited obvious spatial discrepancy. The general trend shows the vulnerability level increasing from northern to southern regions. In the northern regions, the atmospheric environment in Zhangjiakou, Chengde and Qinhuangdao, at the second level, can be described as stable, with a high capacity to recover from disturbance. In the central regions, Beijing, Tianjin and Cangzhou represented moderately vulnerable atmospheric vulnerability, which indicated the medium anti-interference capacity of the atmospheric system. The widely distributed area with fourth level in the southwestern regions (Langfang, Baoding, Shijiazhuang, Hengshui, Xingtai and Handan) and Tangshan can be treated as an unstable atmospheric environment with a low anti-interference ability for disturbance. This suggests that the atmospheric environment of these cities is more potentially affected by minatory disturbance. For a temporal comparison, the general atmospheric vulnerability of 12 cities improved with declining level and Zhangjiakou remained at a weakly vulnerable level with little change.

Aggregated Atmospheric Environment Vulnerability
The assessment results of aggregated atmospheric environment vulnerability are shown in Figure 6. Figure 6c indicates that, in terms of the general atmospheric environment vulnerability, the fourth level (53.84% of 13 cities) took up the largest proportion of the BTH region in 2015. In 2015, the atmospheric environment vulnerability exhibited obvious spatial discrepancy. The general trend shows the vulnerability level increasing from northern to southern regions. In the northern regions, the atmospheric environment in Zhangjiakou, Chengde and Qinhuangdao, at the second level, can be described as stable, with a high capacity to recover from disturbance. In the central regions, Beijing, Tianjin and Cangzhou represented moderately vulnerable atmospheric vulnerability, which indicated the medium anti-interference capacity of the atmospheric system. The widely distributed area with fourth level in the southwestern regions (Langfang, Baoding, Shijiazhuang, Hengshui, Xingtai and Handan) and Tangshan can be treated as an unstable atmospheric environment with a low anti-interference ability for disturbance. This suggests that the atmospheric environment of these cities is more potentially affected by minatory disturbance. For a temporal comparison, the general atmospheric vulnerability of 12 cities improved with declining level and Zhangjiakou remained at a weakly vulnerable level with little change.

Discussion
The atmospheric environment vulnerability of the BTH region showed diverse patterns in space. There are two patterns in the northern cities (Zhangjiakou, Chengde and Qinhuangdao). The weakly vulnerability of Zhangjiakou was caused by weakly exposure, sensitivity, and adaptive

Discussion
The atmospheric environment vulnerability of the BTH region showed diverse patterns in space. There are two patterns in the northern cities (Zhangjiakou, Chengde and Qinhuangdao). The weakly vulnerability of Zhangjiakou was caused by weakly exposure, sensitivity, and adaptive capacity. In the other two cities, the potential adaptive capacity relieved the moderate exposure and sensitivity, resulting in weakly atmospheric environment vulnerability. In the three central cities (Beijing, Tianjin and Cangzhou), the moderate vulnerability was principally influenced by moderate exposure and sensitivity. In the other seven cities with strong vulnerability, there were two possible patterns. One was for Langfang and Hengshui, where the potential and weakly adaptive capacity could not relieve the strong vulnerability to exposure. The other was for Shijiazhuang, Tangshan, Handan, Xingtai and Baoding, showing that the exposure, sensitivity and adaptive capacity were all at moderate to strong vulnerable levels.
The diverse patterns of atmospheric environment vulnerability in space are reasonable. As the sensitivity in most cities was at a high degree with tiny temporal and spatial differences, the spatial change of atmospheric environment vulnerability in the BTH region was mainly influenced by differences in exposure and adaptive capacity indices. In the exposure indices, air quality is gaining increased attention in China, especially the PM 2.5 and PM 10 pollution. According to the air pollution source and transport rationale, the inhomogeneous air pollution distribution is mainly caused by diversity in topographic and meteorological condition, coal consumption, vehicle emissions, industrial emissions, and atmospheric transportation between southern and northern cities. In terms of the wind rose of the BTH region (shown in Figure S1), when the prevailing northerly clean wind from Inner Mongolia blows into the northern cities of the BTH region, the pollutants in the atmosphere are blown away, which lowers the exposure level. Another factor influencing exposure differences was the fact that the southern cities had received many heavy industrial companies from Beijing. Spatial comparison demonstrated that the adaptive capacity was generally lower in the northern areas than in central and southern areas. This circumstance primarily resulted from spatial differences in coal consumption and secondary industry index distribution. As for the coal consumption index diversity, cities in the first and second vulnerable classes of the coal consumption index are mainly distributed in northeastern areas; meanwhile, central and southwestern areas are mainly in the third and fourth classes. For secondary industry index distribution, southwestern cities show a higher proportion of secondary industry than northeastern cities.
The longitudinal comparisons highlight that the exposure index had an obvious downtrend from 2013 and 2015. This is mainly because, since 2012, people and the government have been more sensitive to the severe particulate air pollution that threatens human health. China's Ministry of Environmental Protection also issued Ambient air quality standards [39] on 29 February 2012, which officially induce PM 2.5 into normal atmospheric environment quality assessment. The Ministry of Environmental Protection also set Implementation Rules of Air Pollution Prevention Plan for the Beijing-Tianjin-Hebei region and its adjacent areas [61] in September 2013, aiming to decrease the PM 2.5 concentration to 25% of the 2012 level by 2017. These two important decisions decreased the exposure level for most cities. Although the atmospheric environment vulnerability in the BTH region over the three years has generally improved, the atmospheric environment vulnerability is still unacceptable. This implies that decreasing air pollution is still a major long-term task for China's government. This result is consistent with those of Jin et al. and Guo et al., who found that strenuous effort is required to improve air quality [2,62]. The longitudinal comparison also indicated that the sensitivity and adaptive capacity experienced a tiny change from 2013 to 2015.
The atmospheric environment vulnerability patterns could further impact the atmospheric environment protection management [63,64]. For example, the generally poor situation of the exposure index in the BHT region was caused by severe PM 2.5 and PM 10 pollution and bad AQI day indices. These exposure indices of most cities were mainly at the fourth or fifth level in the three years studied. This suggests that decision makers and managers should pay more attention to air pollution control. More specifically, all the central and southern cities need to implement comparatively strict policy limitations on air pollution emissions, such as implementing Odd-Even Day Vehicles Prohibition, shutting down highly polluting industries, and temporary suspension of kindergartens, primary and middle schools when the urban AQI is equal to or worse than grade V. These cities should concentrate on controlling PM 2.5 and PM 10 , particularly in southwestern cities and Tangshan. Zhangjiakou, Shijiazhuang and Tangshan should improve SO 2 pollution regulations. Shijiazhuang, Tangshan, Handan, Baoding, Xingtai, Hengshui and Langfang should reinforce NO 2 control. Meanwhile, for the sensitivity index, the whole situation of population density and percentage of vulnerable groups is very negative in the BTH region, except for some northern cities (Zhangjiakou and Chengde). Effective and appropriate population management is an urgent task for decision makers and managers. Moreover, improving the adaptive capacity in most cities is a long-term job to decrease the general atmospheric environment vulnerability. The southern cities need to reinforce the developing of the third industry so as to accelerate economic restructuring and shift development mode. Beijing, Tianjin, Shijiazhuang, Tangshan and Baoding should be more concerned about motor vehicle population control. Tianjin, Xingtai, Baoding and Cangzhou should make more effort to strengthen the urban greening rate. The central and southern cities should increase their usage of clean energies and control coal consumption.
To provide support for atmospheric environment vulnerability management, three partitions are proposed according to the assessment results of aggregated atmospheric environment vulnerability (shown in Figure 6). The three district partitions are presented as follows: I-Districts for general control, II-Districts for key control, and III-Districts for strict control.
(1) Districts for general control. These districts consist of the regions with "potentially vulnerable" and "slightly vulnerable" vulnerability classes. These districts are mainly in the northern areas (Zhangjiakou, Chengde and Qinhuangdao), accounting for 38.75% of the area of the 13 cities in the BTH region. (2) Districts for key control. These districts are in the vulnerability class "moderately vulnerable", lying mainly in the central cities (Beijing, Tianjin and Cangzhou), with 16.72% of the area of the 13 cities in the BTH region. (3) Districts for strict control. The cities with "significantly vulnerable" and "extremely vulnerable" vulnerability classes generate these districts, lying mainly in the southwestern (Baoding, Shijiazhuang, Hengshui, Xingtai, Langfang and Handan) and southeastern cities (Tangshan), accounting for 44.53% of the area of the 13 cities in the BTH region.
The developed atmospheric environment vulnerability assessment framework is demonstrated to be valid at identifying districts where natural and anthropogenic indices are accountable for current atmospheric environment vulnerability. Taking the time dimension into account makes this framework able to evaluate the influence of historical policies on atmospheric environment vulnerability and provide guidance for further atmospheric environment vulnerability management. The proposed framework allows for a comparison of atmospheric environment vulnerability if the status quo is maintained or if new decisions are introduced and allows for representation of natural and anthropogenic change under intrinsic or extrinsic stresses.

Conclusions
In conclusion, this paper evaluated the atmospheric environment vulnerability in the BTH region with spatial and temporal comparison by integrating AHP, fuzzy theory and GIS spatial analysis into an ESA framework. In the ESA framework, the logical relationships between atmospheric environment vulnerability and indices of exposure, sensitivity and adaptive capacity were elaborated upon. The vulnerability assessment effectively transforms complicated conceptual insights into operational quantification methodologies, procedures and guidance.
The results indicate that the atmospheric environment vulnerability in the BTH region exhibits obvious spatial discrepancy, with southwestern and southeastern cities showing higher vulnerability and lower vulnerability mainly in the central and northern cities. The diverse spatial patterns indicate that exposure and adaptive capacity indices mainly influence the spatial change of atmospheric environment vulnerability and sensitivity in most cities with the same pattern, with tiny temporal and spatial differences. As the vulnerability patterns could further impact the atmospheric environment vulnerability management, managers should preferentially focus on recognized vulnerable indices.
Future improvements of our framework should rely on better spatial databases. For decision makers in BTH regions, the policy challenge now will be to detect how to balance economic development with atmospheric environment protection in the long term.
Supplementary Materials: The following are available online at www.mdpi.com/1660-4601/15/1/128/s1: Figure S1: The wind rose of the BTH region; Table S1: The judgement matrix and consistency ratio (CR) for atmospheric vulnerability assessment. Author Contributions: The study was designed by Yang Zhang and Jing Shen, and the data were collected by Yang Zhang. Yang Zhang and Yu Li contributed to the overall planning of the project, analysing and interpreting the data, as well as the main draft of the manuscript. All authors read and approved the final manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A . Weights Determination
In our study, the AHP technique is applied to determine the weights of each attribute. Six steps are contained in this method, as follows: (1) Deconstruct the decision-making problem into a hierarchical structure with "atmospheric environment vulnerability" as decision objective, which is deconstructed into the second layer with exposure index, sensitivity index, and adaptive capacity index. Each group is further deconstructed into a third layer according to its own attributes. (2) Make decision tables for each layer of the hierarchical decomposition with pair-wise comparisons.
A preference scaling approach is executed in the pair-wise comparisons with 17 scales numbers: 9, 8, 7 . . . , 1/7, 1/8, 1/9 where 9 means the most important, 8 the second most important, and so on down to 1/9, the least important. (3) Make a judgement matrix with the scales numbers from step 2 to establish the relationships between attributes and group index, as shown in Table S1. In the judgement matrix, 1 are the values of diagonal. If the i th row is more important than the j th column, then the value of (i, j) is more than 1, otherwise the j th column is more important than the i th row. Meanwhile, the value of (j, i) in the matrix is the reciprocal of that in (i, j). (4) Determine the weights for each attribute and group index with the largest eigenvalue (λ max ) of the judgement matrix as shown in Equation (A1): (5) Estimate the consistency of the judgement matrix by consistency ratio (CR). The CR is defined as the ratio between consistency index (CI) and random index (RI), as shown in Equation (A2): where the consistency index (CI) is the consistency of a given evaluation matrix and calculated by Equation (A3): and the RI is a random matrix defined as the average of resulting consistency index depending on the order of the matrix. The value of CR is suggested to be less than 0.1, which means that the judgement matrix has reasonable consistency. All the CR values of each matrix in our paper are lower than 0.10, which demonstrates the suitability of the calculated weights. (6) Aggregate the local weight of each parameter to achieve the weight of corresponding group index, as shown in Table 2.

Appendix B . Fuzzy Theory for Index Evaluation
Fuzzy logic is one of the most effective techniques to handle uncertainties due to vague human understanding and imprecise information. Human understanding or expert opinion is best represented with degrees of certainty (fuzzy logic) rather than absolute certitude (classical logic). Zadeh (1965) introduced the fuzzy set theory in his pioneering work to deal with imprecise information. A fuzzy set is a pair x, µA (x):x ∈ X, in which µA: x→[0,1]. Equation (B1) represents a mathematical form of the fuzzy set theory.
where F is the fuzzy set; U is the universe of discourse; x is any variable; and µ (x) is the membership operatory on x ranging from 0 to 1. The set "F" is a combination of the series of subsets represented by a variable value and a membership value. Membership grades or values represent the belongingness level of a variable value on a scale of 0-1 on the linguistic scales that are more comprehensible by human intelligence. The "degree of confidence" indicates the membership grade shown in Equation (B1). Essentially, the fuzzy set theory recognizes the partial truth concept and assigns a degree of truth values between 0 and 1 to the linguistic scales.
Fuzzy synthetic evaluation based on the concept of fuzzy logic has been used in this study to develop fuzzy membership functions for each indicator. The fuzzy membership functions quantitatively model uncertainties in indicator benchmarks and inputs. Fuzzy synthetic evaluation involves fuzzification, aggregation and defuzzification to evaluate the system performance. These steps are described below.
(1) Step 1: Fuzzification Fuzzification is a process to transfer all input values into fuzzy membership functions. In this process, defining benchmark values is the most important step. Benchmarking for atmospheric vulnerability contains indefinite, inexact, and vague information. In this study, the benchmarks have been fuzzified to handle such uncertainties. The major steps for identifying and fuzzifying the benchmarks of indicators are: Step 1.1: Define atmospheric vulnerability indicators for each criterion under a sub-project. A list of criteria and indicators for the project of "atmospheric vulnerability" developed in this study is presented in Table 1.
Step 1.2: Establish a common linguistic scale for developing benchmarks for each indicator. Indicator values are arbitrarily assigned to each of the linguistic scales based on experience and the data available in the different rating systems stated above. The benchmark of atmospheric environment vulnerability assessment indices under five scales is presented in Table 3. Table 3 illustrates the five levels of performance (or linguistic scale) with different membership functions. The triangular fuzzy functions for each linguistic scale are developed such that the most likely values (membership = 1) of one linguistic scale coincides with the minimum value (membership = 0) of next linguistic scale and the maximum value (membership = 1) of the preceding linguistic scale (2) Step 2: Aggregation Aggregation process synthesizes information (e.g., fuzzy sets) from the bottom of the hierarchy to the elements present at higher hierarchical levels. In this study, the aggregation process is used to develop the fuzzy sets for atmospheric vulnerability objectives by aggregating the fuzzy sets of underlying atmospheric vulnerability indicators. All indicators under each layer have been given equal importance with equal weights.
A direct assignment procedure is applied to evaluate the relative weights of atmospheric vulnerability objectives to calculate the atmospheric vulnerability at the component level. This procedure uses a linear and discrete 1-5 Likert scale represented by potential vulnerable, weakly vulnerable, moderately vulnerable, significantly vulnerable, and extremely vulnerable, respectively, and allows decision makers to prioritize objectives using simple linguistic descriptors. When one of the linguistic scales is selected, corresponding numbers on a 1-5 scale are generated. The generated numbers are normalized across all the objectives under each key component. The normalization process results in the relative weights of indicators. The fuzzy sets for atmospheric vulnerability objectives are obtained from the aggregation of related indicator fuzzy sets using the following equations: · · · · · · · · · · · · · · · · · · · · · µ P In µ W In µ M In µ S In µ ES where F P is the sub-project corresponding to the atmospheric vulnerability objective, F A is the sub-index corresponding to the sub-project of the atmospheric vulnerability objective, and W In are the weights of the indicators, F I is the set of fuzzified values of the corresponding indicators, and µ P In , µ W In , µ M In , µ S In , µ ES In are the corresponding fuzzified values of the respective indicators. The fuzzy sets of atmospheric vulnerability objectives are further aggregated to determine an overall atmospheric vulnerability fuzzy set, "F AEV ", for a complete atmospheric environment system: where W A are the weights of the project layers (sub-project) and F A is the set of fuzzified values of corresponding indices under each project layer. (3) Step 3: Defuzzification The aggregated fuzzy sets are converted to crisp numbers using defuzzification methods in order to represent atmospheric vulnerability. In this study, to obtain crisp output of fuzzified inputs and outputs, the following centre of gravity (COG) method is used for defuzzification: In Equation (B5), the centre of gravity of the fuzzy set is located and a corresponding x-scale value is returned as a crisp number. Subsequently, the defuzzified values of atmospheric vulnerability objectives are plotted on a radar diagram.