In terms of the total affected population, the average population affected by rainstorms and droughts in CPSAs from 1984–2019 reached 34 million people/year, accounting for 22.9% of the national population affected by rainstorms and droughts, which is higher than the national share of the population in this region (16.7%) (Figure 5a). From 1984 to 2019, both CPSAs and China's disaster-affected population increased before 2000 and then decreased. The maximum occurred in 2010, reaching 90 million and 326 million, respectively. Moreover, to identify the affected population, the ratio of the average population affected by rainstorms and droughts to total population in CPSAs and China were 16% and 11.7%, respectively. By contrast, the ratio in CPSAs was 1.6% and 6.8% higher than that in China during 1984–1999 and 2000–2019, respectively. Notably, the affected population in CPSAs increased after 2000.

**Figure 5.** (**a**) Rainstorm and drought-affected population (million people) and its proportion in total population (%), and (**b**) direct economic losses caused by rainstorms and droughts (billion yuan) and its proportion in GDP (%) in CPSAs and China during 1984–2019.

As shown in Figure 6, the impact of rainstorms was more serious than droughts. Among the fourteen regions in CPSAs, the Wuling Mountain area (R7) exhibited maximum occurrence of 4.804 million people/year and 2.773 million people/year affected by rainstorms and droughts, respectively. By contrast, the Tibet region (R13) exhibited a minimum occurrence of 0.087 million people/year and 0.048 million people/year, respectively. From 1984 to 2019, the multi-year average of rainstorm and drought-affected population in CPSAs showed different characteristics between the eastern and western sides of the PDL (Figure 6a,b). Divided by PDL, the number of people affected by rainstorms and droughts in the regions to the east of the PDL reached 16.4 million people/year and 12.4 million people/year, which are 8.8 and 3.3 times higher than the population in the regions along and to the west of the PDL, respectively. In particular, the impact of rainstorms and droughts on the population in the regions to the east of the PDL was particularly severe.

#### *4.2. Direct Economic Losses*

In terms of total disaster losses, the average direct economic losses resulting from rainstorms and droughts in CPSAs reached 29 billion yuan/year from 1984–2019, accounting for 12.6% of the total national disaster losses, which is higher than the national share of the GDP of the region (7.7%) (Figure 5b). The direct economic losses due to rainstorms and droughts in CPSAs and China reached a peak of 93.1 billion yuan and 633.9 billion yuan in 1996 and 1998, respectively. The disaster-related losses in 1991, 1994, 1996, 1998, 2010, 2013 and 2016 were relatively high. To identify the extent of economic losses, the ratio of the mean direct economic losses due to rainstorms and droughts to GDP in CPSAs and China were 1.6% and 1%, respectively. Except for 1991, the ratio of disaster losses to GDP in CPSAs was higher than that of in China.

**Figure 6.** Spatial distribution of (**a**,**b**) disasters-affected population (million people) and (**c**,**d**) disasterscaused direct economic losses (billion yuan) due to (**a**,**c**) rainstorms and (**b**,**d**) droughts in CPSAs during 1984–2019.

From 1984 to 2019, the multi-year average direct economic losses resulting from rainstorms in CPSAs also showed significantly distinct characteristics between the eastern and western sides of the PDL. The economic losses caused by rainstorms were higher than those caused by droughts (Figure 6c,d). The region hit hardest by economic losses resulting from rainstorms and droughts was the Wuling Mountain area (R7), with direct economic losses of 7.103 billion yuan/year and 1.445 billion yuan/year resulting from rainstorms and droughts, respectively. The Lvliang Mountain area (R3) and Tibet region (R13) attained the lowest losses resulting from rainstorms and droughts, with losses of 0.098 billion yuan/year and 0.023 billion yuan/year, respectively. Divided by PDL, the direct economic losses due to rainstorms and droughts in the regions to the east of the PDL reached 20.3 billion yuan/year and 5.7 billion yuan/year, which are 9.2 and 7.4 times higher than those in the regions along and to the west of the PDL. Notably, the disaster-related losses in the regions to the east of the PDL were very severe.

As shown in Figure 7a, the trend of the total economic losses attributed to rainstorms and droughts from 1984 to 2019 increased, except in the Dabie Mountain area (R6), which showed a decreasing trend. All 13 other regions showed increasing trends, especially the regions along and to the west of the PDL, where the disaster-related losses increased faster. Among these areas, the Dabie Mountain area (R6), the Wuling Mountain area (R7) and the Luoxiao Mountain area (R11) showed statistically insignificant trends, while the increasing trends in the remaining 11 regions were significant at a confidence level of 95%.

**Figure 7.** Spatial distribution of linear trends of (**a**) direct economic losses caused by rainstorms and droughts and (**b**) disaster resistance capability over CPSAs during 1984–2019. Black dots indicate the results are significant at a confidence level of 95%.

#### *4.3. Disaster Resistance Capability*

To identify the disaster resistance capability in each CPSA, average values of the LGR from 1984 to 2019 were calculated in this study (Figure 8). As shown by the solid line in Figure 8, the annual mean LGR in CPSAs is 1.1%. The LGR in the areas to the east of the PDL was generally higher than the average and the disaster resistance capabilities were relatively low. Among these regions, the disaster resistance capabilities in the Wuling Mountain area (R7) and the Luoxiao Mountain area (R11) were the lowest, with the LGR reaching 3% and 2.9%, respectively. The disaster resistance capabilities in the regions along and to the west of the PDL were relatively high.

**Figure 8.** The average GDP per capita (yuan) and loss-to-GDP rate (LGR) (%) in CPSAs during 1984–2019. The serial numbers of CPSAs are the same as Figure 1. (1. Southern area of Daxing'anling Mountains; 2. Yanshan-Taihang Mountain Area; 3. Lvliang Mountain Area; 4. Liupanshan Area; 5. Qinba Mountain Area; 6. Dabie Mountain Area; 7. Wuling Mountain Area; 8. Wumeng Mountain Area; 9. West Yunnan Border Area; 10. Yunnan, Guangxi and Guizhou rocky desertification area; 11. Luoxiao Mountain Area; 12. Three prefectures of Southern Xinjiang; 13. Tibet Region; 14. Tibetan areas in four provinces).

The disaster resistance capabilities in CPSAs increased during 1984–2019, except for the Tibet region (R13) (Figure 7b). Among fourteen regions in CPSAs, the trends in the Yanshan-Taihang Mountain area (R2), the Lvliang Mountain area (R3), the Qinba Mountain area (R5), the Dabie Mountain area (R6), the Wuling Mountain area (R7) and the Luoxiao Mountain area (R11) are significant at a confidence level of 95%. In the regions to the east of the PDL, for example the Dabie Mountain area (R6) and the Wuling Mountain area (R7), the resistance capabilities increased faster than other regions.

#### **5. Summary and Discussion**

#### *5.1. Summary*

The population, GDP and GDP per capita in CPSAs reached 210 million people, 1989.3 billion yuan, and 8943 yuan per person, accounting for 16.7%, 7.7%, and 46.9% of China's total amounts, respectively. To consolidate the achievements in the elimination of absolute poverty in China and prevent rural populations from returning to poverty due to meteorological disasters, this study analyzed the spatiotemporal characteristics of rainstorms and droughts and their socioeconomic impacts in China's CPSAs from 1984 to 2019.

The annual average rainstorm and drought days in CPSAs were approximately 1.9 days/year and 44.6 days/year, respectively. It gradually decreased from southeast to northwest. The number of rainstorm days showed a significant increasing trend of 0.075 days/decade, while the decreasing trend of the number of drought days was not significant. Due to rainstorms and droughts, the average annual affected population and direct economic losses in CPSAs reached 34 million people/year and 29 billion yuan/year, accounting for 22.9% and 12.6% of China's total amounts, respectively. The average annual loss rate due to disasters in this region reached 1.6%, which is 0.6% higher than the national level.

Furthermore, distinct features and socioeconomic impacts of rainstorms and droughts were identified on the southeastern and northwestern sides of the PDL along Tengchong-Aihui in China. Droughts often impacted the regions located along the PDL, while rainstorms and droughts occurred more frequently in the regions to the southeast of the PDL than in the regions to the northwest of the PDL. As a result, the affected number of people and direct economic losses due to rainstorms and droughts in the regions to the southeast of the PDL reached 16.4 million people/year and 20.3 billion yuan/year, and 12.4 million people/year and 5.7 billion yuan/year, which were 8.8, 9.2, 3.3 and 7.4 times higher than those in the regions to the northwest of the PDL. Although there were more affected people and higher direct economic losses due to rainstorms and droughts in the regions to the east of the PDL, the disaster resistance capabilities were significantly improved. In contrast, the regions to the northwest side of the PDL showed a significant increasing trend of losses with relatively low disaster resistance capabilities.

#### *5.2. Discussion*

From the perspective of rainstorms and droughts, drought conditions were more severe in regions along the PDL, and both rainstorms and droughts were relatively severe in the regions to the east of the PDL. In addition, the number of rainstorm days significantly increased in the entire study area. Rainstorm-induced disasters became increasingly serious, and only the number of rainstorm days in the Southern area of Daxing'anling Mountains (R1), the Dabie Mountain area (R6), the West Yunnan Border Area (R9) and the southern part of the Tibetan areas in four provinces (R14) showed a decreasing trend. The number of drought days decreased in the entire study area, especially the drought conditions in the three prefectures of southern Xinjiang (R12), the Tibet region (R13) and the Tibetan areas in four provinces (R14), whereas the drought conditions in the West Yunnan Border Area (R9) and the southern part of the Tibetan areas in four provinces (R14) and the number of drought days in the West Yunnan Border Area (R9), the southern Tibetan areas in four provinces (R14) and the eastern Tibet region (R13) notably increased.

According to a comprehensive analysis of the direct economic losses, economic level and disaster resistance capabilities in CPSAs, the disaster-related losses significantly increased, and at a higher rate than that of the regions to the east of the PDL, while the

disaster resistance capabilities did not increase significantly. The situation in the regions to the east of the PDL was the opposite, with high disaster losses but at a lower rate, and the disaster resistance capabilities were improving. The disaster resistance capabilities in the Lvliang Mountain area (R3), the Liupan Mountain region and the three prefectures of southern Xinjiang (R12) were satisfactory, although the GDP per capita was slightly lower than that in the other regions. It is worth noting that the disaster resistance capabilities in the Yanshan-Taihang Mountain area (R2), the Wuling Mountain area (R7), Yunnan, the Guangxi and Guizhou rocky desertification area (R10) and the Tibet region (R13) were relatively low. The GDP per capita in these 4 regions is higher than that of other regions, their cost in disaster prevention and reduction might be relatively low (Figure 8). These regions should enhance the disaster prevention consciousness and increase investment in disaster reduction and prevention facilities. Meanwhile, the number of affected population due to rainstorms and droughts in CPSAs showed a remarkable increase after 2000 (Figure 5a). This might be ascribed to the increase in rainstorm days and population exposed to the disasters since the population in CPSAs increased after 2000, with a growth of 39.7% relative to 1984–1999 (Figure 2a). Therefore, the government should take the occurrence of rainstorms and droughts into consideration during urban planning and land use planning in CPSAs to reduce the exposure and risk of population and economy to disasters.

For a long time, the international community has paid great attention to the relationship between poverty and natural disasters. For example, the reduction of rainfall has slowed the economic growth of sub-Saharan Africa [42]. Climate change has plunged 6.5% of the population of Latin America into poverty and vulnerability [43]. Previous studies have reported that disasters will aggravate poverty and recovery is not straightforward for poor people [44]. Under global warming, the frequency and intensity of rainstorms and droughts might increase significantly over CPSAs [6,7]. Correspondingly, the average annual direct economic losses due to rainstorms is expected to be 4 times and 17 times higher than it is currently under global warming of 1.5 ◦C and 4.0 ◦C [45]. As the relatively poor area in China, it is necessary to pay close attention to the socio-economic impact and risks of rainstorms and droughts over CPSAs in the future. Hence, the above-mentioned regions with high GDP per capita but low disaster resistance capabilities should strengthen their awareness of disaster prevention and mitigation, enhance early warnings and disaster relief inputs in response to meteorological disasters, and improve the disaster resistance capabilities to meet climate change challenges. This is an important part of consolidating China's gains in poverty eradication and further contributing to the realization of SDGs.

**Author Contributions:** Conceptualization, M.G., T.J. and J.Y.; Data curation, A.L. and H.Z.; Methodology, A.L. and S.G.; Supervision, M.G., T.J. and J.Y.; Validation, A.L. and X.W.; Visualization, A.L. and S.G.; Writing—original draft, A.L. and M.G.; Writing—review & editing, A.L., X.W. and M.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was funded by the National Natural Science Foundation of China (grant no.42071024, 42105158 and 42005126), the project supported by the State Key Laboratory of Earth Surface Processes and Resource Ecology, and the Startup Foundation for Introducing Talent of NUIST.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors would like to thank Ke Tan of Nanjing University of Information Science and Technology for the support in data.

**Conflicts of Interest:** The authors declare that they do not have any conflict of interest.

#### **References**

#### *Article* **Fresh Air–Natural Microclimate Comfort Index: A New Tourism Climate Index Applied in Chinese Scenic Spots**

**Xiaoyan Yang 1, Changshun Li 1,2,\*, Muhammad Bilal 3 and Shaofei Jin 1,4,\***

**Abstract:** Severe air pollution in China has caused significant tourism transformation for pursuing fresh air in microclimate tourism markets. Contemporary practices simply measure the air freshness of destinations and scenic spots using a single index, i.e., primarily negative oxygen ions (O2 −). This index cannot comprehensively reveal scenic spots' air freshness degree and determine the dynamic interactions between air freshness and scenic spots' tourism development, thus inducing an illusion of air freshness for the target scenic spots. Meanwhile, the current fresh air index primarily ignores connections with the microclimate index of scenic spots and cannot provide a multidimensional index for scenic spots to take advantage of both air and microclimate resources for diverse tourism products and service production. Therefore, this study proposes a multidimensional index, the fresh air–natural microclimate comfort index (FAI-NMCI), connecting the fresh air index with the natural microclimate comfort index of scenic spots together from transdisciplinary and multidisciplinary perspectives. This study utilizes FAI-NMCI to measure four scenic spots of Fujian Province, and reveals in-depth results of scenic spots' air freshness and natural microclimate comfort degree together. The results demonstrate that the four scenic spots in Fujian province of China had different levels of air freshness degree and natural microclimate comfort degree in 2018. The natural scenic spots were mostly distributed in Healing Fresh, Very Fresh, and Super Fresh levels of FAI with the most comfortable and comfortable levels of NMCI. The cultural scenic spots were mostly distributed in Relatively Fresh and Healing Fresh levels of FAI with the most comfortable and comfortable levels of NMCI. Meanwhile, the FAI-NMCI of natural and cultural scenic spots also had significant differences within 24 Jieqi, which will promote dynamic and creative utilization of those resources in microclimate tourism development.

**Keywords:** fresh air index; natural microclimate comfort index; fresh air–natural microclimate comfort index; scenic spots; Fujian province

#### **1. Introduction**

Air pollution influences tourists' and residents' (potential tourists') health, which will directly or indirectly affect demand in the tourism market [1]. Exposure to air-polluted environments leads to various health problems [2], such as respiratory infections [3,4], asthma [5,6], stroke, and even lung cancers [7]. Additionally, air pollution, especially fog and haze pollution, also impacts the psychological health [8,9] and mental disorders [10] of tourists and potential tourists, which directly increases healthcare expenditures [11], health costs [12], and other socioeconomic burdens [13,14] and socioeconomic costs [15] for tourists and local residents. Meanwhile, the factors of natural microclimate and the comfort of the natural microclimate as a whole also influence tourists' decision-making for specific destinations or scenic spots, impacting tourists' behavior and the constant and dynamic

**Citation:** Yang, X.; Li, C.; Bilal, M.; Jin, S. Fresh Air–Natural Microclimate Comfort Index: A New Tourism Climate Index Applied in Chinese Scenic Spots. *Sustainability* **2022**, *14*, 1911. https://doi.org/ 10.3390/su14031911

Academic Editor: Andrei P. Kirilenko

Received: 4 December 2021 Accepted: 30 January 2022 Published: 8 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

tourist flow [16] from one region to another because of dynamic natural microclimate changes [17,18]. Under an air pollution background, the increasing fresh air demand of tourists in China has attracted tourists to destinations with fresh air to experience comfort in natural microclimates. Fresh air and natural microclimate together provide friendly, healthy, and comfort natural tourist attractions in the contemporary Chinese tourism market. Air freshness and microclimate comfort degree together as vital environmental information have significantly influenced tourists, destination and scenic spot management, as well as tourism market development, which has seldom been explored by contemporary research and practice.

Actually, the most popular index for measuring air freshness and air quality is the air quality index (AQI), which utilizes primary air pollutant components such as PM2.5, PM10, SO2, NO2, CO, or O3 to measure air quality. AQI studies have developed according to specific research purposes or considering practical utilization; for example, many studies have explored relationships between the AQI and years of life lost [19] and the aggregated effects of internal indicators at a short-term scale [20]. However, AQI reveals the degree of air pollution situations of cities or some air-polluted areas, such as factories or streets, rather than scenic spots. It also neglects the degree of air freshness through adding indicators of O2 − to measure the health and beneficial degree of fresh air [21]. Moreover, current studies of the fresh air index primarily utilize singular indicators, only O2 −, to measure air freshness, without considering air pollutants' influence on the comprehensive air freshness degree [22].

For the natural microclimate of scenic spots, many studies primarily utilize one index or two indexes to analyze and reveal the natural microclimate comfort degree, which can mainly be summarized as the temperature-humidity index (THI) [23], wind-cold index (WCI) [24], index of cloth loading (ICI) [25], and human comfort index (HCI) [26]. However, it is vital to note that current indexes do not consider rain factors to measure the influence of rain volume on the comfort degree, which provides an opportunity to establish a rain comfort index and construct a new comprehensive index including all primary factors of the natural microclimate of scenic spots. Additionally, the natural microclimate comfort degree of scenic spots has seldom been integrated with the fresh air index in contemporary studies due to a lack of comprehensive air freshness data and natural microclimate factor data at the same time, along with the existence of some research barriers.

Integration embraces the power to reveal the actual distribution of real scenic spot air freshness and natural microclimate resources. Based on individual economic development trends of different places, it is essential to point out that the air freshness degree of scenic spots needs to be classified according to WHO standards [27] and actual situations together, which calls for stricter standards of fresh air, with gradual improvement in standards. For more in-depth exploration and development of scenic spots with fresh air and natural microclimate tourism resources, it is indispensable to classify different levels or add more actual levels of air freshness degree and natural microclimate degree in future research because the current classification needs to be developed based on more accurate research.

For sustainable development of destinations and scenic spots with fresh air and natural microclimate comfort resources, this paper establishes a fresh air index (FAI) based on comprehensive identification of fresh air components and the index construction methods and experiences of AQI, constructs natural and cultural FAI according to the different threshold standards, and creatively constructs a natural microclimate comfort index of scenic spots through taking advantage of the current index and adding a rain comfort index, connecting the fresh air index (FAI) and the natural microclimate comfort index (NMCI) together to assess the air freshness degree of destinations and the natural microclimate comfort degree. For in-depth air and microclimate resource distributions in the tourism market, this paper proposes an improved fresh air–natural microclimate comfort index to systematically analyze fresh air resources and natural microclimate comfort resources together in Chinese scenic spots by classifying seven levels of the fresh air index (FAI) and seven levels of the natural microclimate comfort index (NMCI). Based on the comparative advantages of embracing relatively fresher

air, this paper takes four scenic spots of Fujian Province, China, as an example and analyzes the natural and cultural scenic spots' fresh air–natural microclimate comfort degree to promote the sustainable utilization and development of fresh air and natural microclimate comfort resources in the microclimate tourism market.

#### **2. Materials and Methods**

This study utilizes the comprehensive index construction method to construct one fresh air–natural microclimate comfort index (FAI-NMCI), then creatively applies this index to four scenic spots of Fujian Province to reveal the actual air freshness and microclimate comfortable degree at the same time.

#### *2.1. Methodology of Fresh Air–Natural Microclimate Comfort Index Construction and Evaluation* 2.1.1. Definition and Essence of Fresh Air

Different from polluted air, fresh air not only contains basic air components and air structures, but also contains more beneficial air components, such as negative oxygen ions (O2 −), which are also treated as "air vitamins" for promoting human health [28,29] and include no harmful, or less harmful air components, such as PM2.5, PM10, SO2, NO2, CO, or O3, that are harmless standards. Therefore, fresh air has beneficial and no harmful characteristics and can be simply defined and summarized as follows:

$$\text{Fresh air (FA)} = \text{bereficial air components (BAC)} + \text{no harmful air components (NAD)}$$

According to the aforementioned analysis, fresh air could be summarized as having health essences from protecting and promoting human health perspectives, which could be represented by a higher concentration of healthy air components and a lower concentration of polluted air components.

In this study, negative oxygen ions (O2 −) were adopted to represent beneficial air components and were selected to reveal general nonharmful air components that severely caused air pollution in China. SO2, NO2, CO, or O3 was chosen to represent local special air contaminants in different regions of China. Therefore, to illustrate fresh air in specific scenic spots of China, this study constructed a fresh air structure with selected air components, which can be summarized as follows:

#### Fresh Air = Beneficial Air Components (O2 −) + No Harmful Air Components {PM2.5+PM10+(SO2, NO2, CO or O3)}

#### 2.1.2. Improved Fresh Air Index of Natural and Cultural Scenic Spots

We compared the contemporary single-factor fresh air index of Fujian, Hubei, and Zhejiang provinces of China, utilizing only negative oxygen ion (O2 −) degree and PM2.5, to assess scenic spots' air freshness (Table 1). The Meteorological Service Center of Fujian Province, China, and the authors of this paper have creatively proposed a comprehensive fresh air index by considering both the healthy air component factor negative oxygen ion (O2 −) and the low degrees of PM2.5 and PM10 as general polluted air components that exist in different regions of China, and NO2, SO2, CO, and O3 as particular air pollutant factors in specific regions of China, determined by local geological and economic structure features.

The subindices of the beneficial air component IFAIB are expressed as the ratio of the negative oxygen ion (O2 −) concentration to the 24 h average standard O2 − level (1000 ions/cm3) (Table 2) recommended by the World Health Organization Air Quality Guidelines for fresh air [30–32]. That is, IFAIB can be written as:

$$\text{IFAI}\_{\text{B}} = \frac{\text{C}\_{\text{B}}}{\text{S}\_{\text{B}}} \tag{1}$$

where IFAIB represents the subindices of the beneficial air component in scenic spots, CB represents the negative oxygen ion (O2 −) concentration of scenic spots, and SB represents the 24 h average standard O2 − level (1000 ions/cm3) of scenic spots.

**Table 1.** Comparison of methods and classification of air freshness in different provinces of China.

**Table 2.** Standard concentration degree of indicators of fresh air components.

Notes: Data source from Ambient Air Quality Standards of China (GB3095-2012). The category I region used I concentration degree, and the category II region used II concentration degree. The category I region is the place of natural conservation, famous scenic spots, and other special conservation regions that need special protection. The category II region is the residential regions, the mixed region with commercial, transportation, and residents, the cultural regions, industrial regions, and rural regions.

The subindex IFAIpi of nonharmful air components is expressed as the negative ratio of air pollutant concentrations, Cpi, to the recommended short-term Ambient Air Quality Standards of China (GB3095-2012) [33], Spi. That is, IFAIpi can be written as:

$$\text{IFAI}\_{\text{Pi}} = -\frac{\text{C}\_{\text{Pi}}}{\text{S}\_{\text{Pi}}} \tag{2}$$

where IFAIpi represents the subindices of no harmful air components in scenic spots, Cpi represents air pollutant concentrations of scenic spots, and Spi represents the recommended short-term air pollutant concentration standards of scenic spots.

The short-term Ambient Air Quality Standards of China [34] used in natural scenic spots are defined as concentrations with average times of 24 h means of 50 μg/m3 for PM10, 35 μg/m3 for PM2.5, 50 μg/m3 for SO2, 80 μg/m3 for NO2, and 4 mg/m3 for CO, and running an 8 h mean of 100 μg/m3 for O3. The short-term Ambient Air Quality Standards of China used in cultural scenic spots are defined as concentrations with average times of 24 h means of 150 μg/m3 for PM10, 75 μg/m3 for PM2.5, 150 μg/m3 for SO2, 80 μg/m3 for NO2, and 4 mg/m3 for CO, and running an 8 h mean of 160 μg/m3 for O3. The 24 h average standard O2 − level recommended by the World Health Organization Air Quality Guidelines for fresh air is 1000 ions/cm3 in category I and II regions (Table 2).

The natural and cultural fresh air index is the sum of IFAIB and the sum of subindices of no harmful air components, IFAIPi, which can be written as:

$$\text{FAI} = \text{IFAI}\_{\text{B}} + \sum \text{IFAI}\_{\text{Pi}} = \sum \text{IFAI}\_{\text{j}} \tag{3}$$

where FAI refers to the fresh air index, IFAIB refers to subindices of beneficial air components, ∑IFAIpi indicates the sum of subindices of no harmful air components, and ∑IFAIj indicates the sum of subindices of fresh air components j.

The classifications of natural and cultural FAIs of scenic spots and their levels were explored in this paper. According to Table 1, the FAI of natural scenic spots utilizes category I region standards to further explore different threshold levels, while the FAI of cultural scenic spots utilizes category II region standards (Table 2).

According to the China Industrial Technical Regulation on Grade of Air Negative (Oxygen) Ion Concentration (QX/T 380-2017) [32] and Local Technical Regulation on Grade of Air Negative (Oxygen) Ion Concentration of Hubei Province (DB42/T 1198-2016) [30], the concentration degree of negative oxygen ion (O2 −) increase per 500 ions/cm3 will enhance air freshness and have a beneficial influence on human health.

According to the World Health Organization (WHO) air freshness grading standard of the negative oxygen ion index, when 1800 ≤ O2 − ≤ 2100, O2 − has natural healing power, and when O2 − > 2100, it has therapeutic and rehabilitation effects [35]. Therefore, on the basis of a higher standard than the WHO standard for negative oxygen ion therapy, this paper proposes that when 2500 ≤ O2 − ≤ 3000, PM2.5, PM10, O3, NO2, SO2, and CO decreased in proportion, and the air freshness degree is at the Healing Freshness Level, the resulting air freshness level of the scenic spot reaches the healing standard. Then, when 3000 ≤ O2 − ≤ 3500 and when O2 − > 3500, the air could have a more active influence on human health, which could further promote the air freshness degree to very fresh and super fresh levels.

According to the World Health Organization Air Quality Guideline (WHO AQG), a PM2.5 and PM10 decrease of 5 μg/m3 will reduce mortality risk by approximately 3% [31], an O3 decrease of 5 μg/m3 will reduce daily mortality by approximately 0.12–0.25% [31], an SO2 decrease of 5 μg/m3 will reduce daily mortality by approximately 0.09–0.29% [31], a CO decrease of 0.5 mg/m3 will also reduce mortality risk or health risk [31], and an NO2 decrease of 5 μg/m3 will reduce daily mortality by approximately 0.29% [36].

This paper further classifies the FAI threshold range of natural and cultural scenic spots according to the different standards of natural and cultural scenic spot air quality standards in China and the WHO standards, which are summarized in Tables 2 and 3, respectively. The air freshness degree of natural and cultural scenic spots could be classified into seven levels: Not Fresh, General, Fresh, Relatively Fresh, Healing Fresh, Very Fresh, and Super Fresh (Tables 3 and 4).

#### 2.1.3. Improved Natural Microclimate Comfort Index

It is obvious that those indices usually assess two or three factors of natural microclimate, such as THI analyzes temperature and humidity factors, WCI analyzes wind and temperature factors, ICI analyzes wind, temperature, and solar radiation factors, and HCI analyzes temperature, humidity, and wind factors. Therefore, it is a trend in Chinese tourism research to utilize the indices mentioned above to comprehensively analyze the degree of natural microclimate comfort in tourists' destinations [37–40]. However, there is also a lack of consideration of the influence of precipitation factors on the degree of

natural microclimate comfort in contemporary Chinese tourism research. According to climate comfort analysis of scenic spots, when precipitation < 1.0 mm, there is nearly no influence on tourists. When precipitation < 3.0 mm, it is within tourists' precipitation comfort range, which will provide a rainfall microclimate for tourists to experience and appreciate rainfall scenery in scenic spots, while when precipitation > 10.0 mm, tourism will decrease dramatically because large volumes of rainfall make tourists feel uncomfortable and cause them to stop tourism activities and leave scenic spots [41].

**Table 3.** Threshold range of natural scenic spot FAI levels.

**Table 4.** Threshold range of cultural scenic spot FAI levels.

Therefore, this paper proposes a rain comfort index (RCI) to assess the comfort degree to which different volumes of rainfall influence tourists in scenic spots. The RCI refers to the influence of rainfall on the extent of human comfort in a natural microclimate environment. According to the Chinese grade of precipitation standard (GB/T 28592-2012) and considering tourism activities usually within a 12 h range, this paper chose a 12 h precipitation standard to measure rainfall levels; that is, drizzle < 0.1 mm, light rain 0.1~4.9 mm, moderate rain 5.0~14.9 mm, and heavy rain 15.0~29.9 mm. Based on the above analysis, this paper defined the comfort range of precipitation as between 0.1 (drizzle) and 10.00 mm (moderate rain) for tourists in scenic spots and then constructed the RCI as follows:

$$\mathbb{RCI} = \mathbb{R}/\mathbb{R}\_0$$

where RCI refers to the rain comfort index, R refers to actual rainfall, and R0 refers to standard rainfall, that is equal to 0.1 mm (drizzle) in this paper. Then, seven levels of RCI of different comfort degrees were classified, which are summarized in Table 5.

By considering all factors of the natural microclimate, this paper proposes the natural microclimate index (NMCI):

$$\text{NMCI} = \frac{\text{THI} + \text{WCI} + \text{ICI} + \text{HCI} + \text{RCI}}{5}$$

**Table 5.** Threshold range of RCI levels.

By assessing each subindex of the NMCI (THI, WCI, ICI, HCI, RCI) and defining them into their own comfort levels, the results are summarized in Table 6. For subindices of NMCI, such as THI, WCI, ICI, HCI, and RCI, which have their own levels of standards and classification patterns, it is essential to standardize each subindex of NMCI for further calculation of NMCI. This paper summed up and calculated their average comfort levels through reclassifying each subindex of NMCI into seven levels (d, c, b, A, B, C, D) and assigning values (that is, d = 2, c = 3, b = 4, A = 5, B = 6, C = 7, D = 8). Then, the average THI + WCI + ICI + HCI + RCI/5 value was calculated to determine the final levels of NMCI in Table 7.

When 5 ≤ NMCI < 6, the NMCI is at the most comfortable level; when 4 ≤ NMCI < 5, the NMCI is at the comfortable (cool comfort) level; when 6 ≤ NMCI < 7, the NMCI is at the comfortable (warm comfort) level; when 3 ≤ NMCI < 4, the NMCI is at the less comfortable (cold less comfortable) level; when 7 ≤ NMCI < 8, the NMCI is at the less comfortable (hot less comfortable) level; when 2 ≤ NMCI < 3, the NMCI is at the uncomfortable (very cold uncomfortable) level; when NMCI > 8, the NMCI is at the uncomfortable (very hot uncomfortable) level (Table 7).

#### 2.1.4. Combined Fresh Air–Natural Microclimate Comfort Index

After the construction of the fresh air index (FAI) and the natural microclimate comfort index (NMCI), this paper comprehensively constructed a parallel fresh air–natural microclimate index for assessing the air freshness degree and natural microclimate comfort degree of specific scenic spots in China. This approach utilizes the fresh air index and natural microclimate comfort index to analyze scenic spots and then connects them together to explore fresh air and natural microclimate tourism resources in depth (Table 8).

#### *2.2. Data Source*

This study chose four fresh air index monitoring stations within four scenic spots of Fujian Province, China. That is, the Pingnan Baishuiyang-Yuanyang Brook Scenic Spot in Ningde City, the Gutian Conference Site Scenic Spot in Longyan City, the Yongding earth building (Tulou) Scenic Spot in Longyan City, and the Wuyi Mountain Scenic Spot in Nanping City (Figure 1 and Table 9).

Pingnan Baishuiyang-Yuanyang Brook Scenic Spot is a natural scenic spot. It has a total area of 66 square kilometers, integrating streams, peaks, rocks, waterfalls, caves, lakes, and other landscapes. There are more than 10 kinds of national protected animals, such as Silver Pheasant, Python, and Pangolin, which are found in the "Natural Animal and Botanical Garden" (Figure 1 and Table 9).

Gutian Conference Site Scenic Spot is a cultural scenic spot. It is located in the south of Meihua Mountain, with a planned area of 7.6 square kilometers. In October 2015, it was successfully promoted to a national 5A tourist attraction and won the honorary titles of "National Top Ten Excellent Patriotism Education Bases" and "National Red Tourism Classic Scenic Spots" (Figure 1 and Table 9).

**Table 7.** Threshold range of NMCI Levels.

**Table 8.** Threshold range of FAI-NMCI Levels.

**Figure 1.** The four scenic spots and their monitoring stations in Fujian Province, China.

**Table 9.** Four scenic spots and their location information.

Yongding earth building (Tulou) Scenic Spot is a cultural scenic spot. It is a unique magical mountain residential building. It has a long history, unique style, large scale, and exquisite structure. In July 2008, it was successfully listed in the world heritage list (Figure 1 and Table 9).

Wuyi Mount Scenic Spot is a natural scenic spot. It is located at the junction of Jiangxi and northwest Fujian Provinces and at the southeast foot of the north section of Wuyi Mountain vein, with a total area of 999.75 square kilometers. It is a famous scenic spot and summer resort in China (Figure 1 and Table 9).

Since the PM2.5, PM10, and O3 standards in natural and cultural scenic spots have different regulations in practice, this paper attempts to establish natural and cultural scenic spot FAIs to distinguish air freshness evaluations. In addition, natural and cultural scenic spots have different environmental backgrounds and social interaction intensities, which will significantly influence the air freshness degree of scenic spots. Meanwhile, for tourists, the fresher the air in scenic spots, the better the traveling experience, so they could utilize the FAI of natural scenic spots to measure all scenic spots in practice. For destination and scenic spot management, it is vital to distinguish the difference between natural and cultural scenic spots and to evaluate the related and compared freshness within natural scenic spots and cultural scenic spot categories. This paper chose the Wuyi Mountain scenic spot and Baishuiyang scenic spot as natural scenic spots of Fujian Province and the Yongding Tulou scenic spot and Gutian Conference Site scenic spot as cultural scenic spots of Fujian Province.

According to local air characteristics, O3 is a special air component that influences the local air freshness of Fujian Province. Therefore, data on four indicators, O2 −, PM2.5, PM10, and O3, from 1 January to 31 December 2018 were collected. Temperature, relative humidity, wind speed, solar radiation, and precipitation data of the four scenic spots were collected from 1 January to 31 December 2018. Therefore, this paper utilized negative oxygen ion (O2 −), PM2.5, PM10, and O3 data for each day of 2018 of the four scenic spots to calculate their fresh air index (FAI). Then, temperature, relative humidity, wind speed, solar radiation, and precipitation data were utilized on each day of 2018 to calculate the natural microclimate comfort index (NMCI) of the four scenic spots.

The data for natural microclimate factors were collected from the National Meteorological Information Center of China (http://data.cma.cn/, accessed on 1 February 2022) and provided by the Fujian Meteorological Service Center, China (http://fj.cma.gov.cn/, accessed on 1 January 2022), and the data for fresh air index subindices were provided by the Fujian Meteorological Service Center, China (http://fj.cma.gov.cn/, accessed on 1 January 2022). The data were generated from the four scenic spots' microclimate and air freshness index monitoring equipment. More precisely, one set of negative oxygen ion monitoring stations (model HuaTron FR500) and one set of air index stations (model HuaTron AEI365) were built in the above four scenic spots (Figure 1 and Table 8).

#### **3. Results**

#### *3.1. FAI-NMCI of Natural Scenic Spots*

In the Wuyi Mountain scenic spot, the FAI levels were mainly distributed in the upper levels; for example, the Super Fresh level was approximately 164 days (accounting for 47%), the Very Fresh level was approximately 67 days (accounting for 19%), and the Healing Fresh level was approximately 56 days (accounting for 16%), which was approximately

287 days in 2018, accounting for almost 82% (Figure 2A). Meanwhile, the NMCI primarily included three levels. The comfortable level was approximately 246 days (accounting for 67%), which also embraced the cool comfort level (168 days, 46%) and the warm comfort level (78 days, 21%). The most comfortable level was approximately 116 days (accounting for 32%), while the less comfortable level was approximately 3 days (accounting for only 1%) (Figure 2B). Therefore, in the Wuyi Mountain scenic spot, the FAI-NMCI was mainly distributed within the upper fresh-comfort range, which was approximately 246 days of upper freshness levels and 362 days of comfort in 2018 (Figure 2A,B).

**Figure 2.** FAI and NMCI levels of Wuyi Mountain (**A**,**B**) and Baishuiyang scenic spot (**C**,**D**) of Fujian Province, China.

In the Baishuiyang scenic spot, the FAI levels were distributed from the Not Fresh level to the Super Fresh level. More precisely, the Fresh level was approximately 86 days (accounting for 25%), the Relatively Fresh level was approximately 74 days (accounting for 21%), the Super Fresh level was approximately 51 days (accounting for 15%), the Healing Fresh level was approximately 40 days (accounting for 11%), and the Very Fresh level was approximately 36 days (accounting for 10%), which indicated the dynamic fluctuation of FAI in the Baishuiyang scenic spot, and the FAI reached or exceeded the Fresh level for approximately 287 days, accounting for approximately 78.63% (Figure 2C). At the same time, the NMCI was mainly distributed at the comfortable level (185 days, 51%), among which the cool comfort level was 117 days (accounting for 32%) and the warm comfort level was 68 days (accounting for 19%). Meanwhile, the most comfortable level was also in a large proportion (45%) and was approximately 165 days. The less comfortable level was approximately 15 days and accounted for 4% (Figure 2D). Therefore, in Baishuiyang scenic spot, the FAI-NMCI was mainly distributed within the middle fresh-comfort range, which was approximately 160 days of middle freshness and 185 days of comfort in 2018 (Figure 2C,D).

In Wuyi Mountain scenic spot, the FAI mainly fluctuated between −2.69 and 32.7, that is, it fluctuated between Not Fresh and Super Fresh levels. The NMCI mainly fluctuated between 3.6 and 7.2, that is, it fluctuated between cold less comfortable and hot less comfortable levels. Within the whole year of 2018, from January to April, the NMCI gradually increased from 4 to 5 (in the cool comfort level), while the FAI was mostly more than 5 (in the Super Fresh level) in the second half of January, decreased dramatically to −1.87 (in General level) in early February, then slightly increased to Relatively Fresh, Healing Fresh, and Very Fresh levels, and some Super Fresh levels, which also fluctuated in this period. From May to the first half of July, the NMCI increased from 5 to 6 and 6 to 7, which generally entered into the most comfortable level and increased into the warm comfort level in this period, which fluctuated mainly within these two levels; meanwhile, the FAI mainly fluctuated between Fresh and Super Fresh levels and reached the highest FAI of approximately 3.7 in June, which was more than the threshold value of Super Fresh (FAI > 2) by approximately 16.35 times. From the second half of July to the first half of September, the NMCI remained at the warm comfort level, while the FAI decreased gradually from more than 5 to lower than 3, although it was also mainly within the Super Fresh level. Then, the FAI increased to more than 5 in the second half of August and increased to more than 20 in the first half of September. From September to October, the NMCI decreased gradually from the warm comfort level to the most comfortable level, while the FAI also decreased to lower levels, such as the Fresh level and Relatively Fresh level, although some days reached the Healing Fresh level and Very Fresh level. When entering November and December, the NMCI decreased from 5 to lower than 4, which is a decrease from the cool comfort level to the cold less comfortable level, while the FAI increased gradually to the Super Fresh level in the middle of November and then decreased dramatically in the second half of November and the first half of December. When the NMCI decreased to the cold less comfortable level in December, the FAI fluctuated, while most fluctuated between the Healing Fresh and Super Fresh levels (Figure 3A).

**Figure 3.** FAI and NMCI of Wuyi Mountain (**A**) and Baishuiyang scenic spots (**B**) of Fujian Province on each day of 2018.

In Baishuiyang scenic spot, the FAI mainly fluctuated between −2.83 and 34.27, which fluctuated between the Not Fresh level and Super Fresh levels, while the NMCI mainly fluctuated between 3.4 and 7.2, which fluctuated between cold less comfortable and hot less comfortable levels. Within the whole year of 2018, from January to early June, the NMCI gradually increased from 4 to nearly 7, namely, it increased from cool comfort to the most comfortable level and then to the warm comfort level, while the FAI was mostly below 0, namely, it fluctuated between the Not Fresh level and the Relatively Fresh level. From the middle of June to August, the NMCI fluctuated and gradually increased from 6 to more than 7, that is, from the warm comfort level to the hot less comfortable level; meanwhile, the FAI increased from 0 to between nearly 5 or more than 5 on some days, which mainly fluctuated between the General level and the Super Fresh level. From September to November, the NMCI fluctuated mainly between 6.5 and 5, which indicated that the natural microclimate fluctuated from the warm comfort level to the most comfortable level from early September to late November. The FAI slightly and gradually increased in September and then decreased drastically below 0 from late September to early October, namely, decreased from the Super Fresh level to the Not Fresh or General level, mainly

during the National Holiday of China (approximately seven-day long holiday). After the seven days of the National Holiday of China, the FAI increased dramatically and reached the highest level of the year (FAI = 34.27), then fluctuated within the Super Fresh level (mainly between 2 and 14). From late November to December, the NMCI decreased below 5 but higher than 4, namely, at the cool comfort level, and the FAI also decreased to a lower level of the year, which primarily fluctuated between the General level and the Healing Fresh level (Figure 3B).

It is essential to point out that the FAI of the Wuyi Mountain and Baishuiyang scenic spots have some higher FAI values, such as more than 30 or even higher, which is primarily because of the higher O2 − levels and the lower PM2.5, PM10, and O3 levels in scenic spots. When scenic spots have a high O2 − concentration, it will significantly influence the entire FAI value.

More precisely, within 24 Jieqi of Wuyi Mountain scenic spot and Baishuiyang scenic spot, the FAI-NMCI can be summarized as in Figure 4. For example, it is obvious that in Lesser Cold and Greater Cold Jieqi, both the Wuyi Mountain and Baishuiyang scenic spots were at the cool comfort level, although Wuyi Mountain was at the Super Fresh level and was fresher than Baishuiyang, which was only at the Fresh level. At the start of Spring, the FAI of both Wuyi Mountain and Baishuiyang decreased, which was not because of the decrease in negative oxygen ions (O2 −) at the two scenic spots but because of the increase in PM2.5 and PM10 due to the large number of tourists and tourism transportation entering those two scenic spots (Figure 4).

From Rain Water to Spring Equinox Jieqi, although the two scenic spots were also within the cool comfort level, the FAI both experienced a drastic increase from the General to the Super Fresh level in Wuyi Mountain and from the Not Fresh to the Relatively Fresh level in Baishuiyang. In Pure Brightness and Grain Rain Jieqi, Baishuiyang entered into the most comfortable level while Wuyi Mountain remained at the cool comfort level. Meanwhile, for more tourist travel to Baishuiyang, the FAI was at a relatively lower level, that is the General level, than Wuyi Mountain (Super Fresh and Healing Fresh level). At the Start of Summer Jieqi, Wuyi Mountain was at the Super Fresh and warm comfort levels, while Baishuiyang was only at the Fresh level and was still within the cool comfort level. In Grain Full Jieqi, Baishuiyang had increased freshness and became warm. In the Grain in Ear Jieqi, Wuyi Mountain entered the most comfortable level and was in a Super Fresh situation (Figure 4).

#### *3.2. FAI-NMCI of Cultural Scenic Spots*

In the Yongding Tulou scenic spot, the FAI levels were mainly distributed in Healing Fresh (159 days, 38%) and Relatively Fresh (74 days, 18%), while they were also distributed in Not Fresh (115 days, 27%). The Super Fresh (29 days, 7%), Very Fresh (31 days, 7%), and Fresh (14 days, 3%) levels were in small proportions (Figure 5A).

The NMCI was distributed in the comfort domain, that is, the most comfortable domain (149 days, 41%), and the comfort domain (216 days, 59%), which included the cool comfort domain (133 days, 36%) and the warm comfort domain (83 days, 23%) (Figure 5B). In the Gutian Conference Site scenic spot, the highest levels of FAI were at the Healing Fresh level (166 days, 48%) with a large proportion, and the Relatively Fresh level (146 days, 42%), while the Very Fresh (24 days, 7%) and Fresh (12 days, 3%) levels had a small proportion, and the Super Fresh level was not reached (Figure 5C). Meanwhile, the NMCI included three levels: the most comfortable (110 days, 30%), and the comfort (250 days, 69%), which embraced the cool comfort level (160 days, 44%) and the warm comfort level (90 days, 25%) (Figure 5D).

**Figure 4.** 24 Jieqi FAI-NMCI of Wuyi Mountain and Baishuiyang scenic spots.

In the Yongding Tulou scenic spot, the FAI mainly fluctuated between −0.97 and 14.81, which fluctuated between Fresh and Super Fresh levels. The NMCI mainly fluctuated between 3.80 and 6.80, which fluctuated between cold less comfortable and warm comfort levels (Figure 6A). Throughout 2018, the FAI and NMCI of the Yongding Tulou scenic spot experienced dynamic fluctuations. From January to March, the NMCI increased gradually from 3 to above 4, that is from cold less comfortable to the cool comfort level, and the FAI mainly fluctuated between −0.35 and 0.48, that is mainly at the Relatively Fresh level. From late March to the middle of April, the NMCI continued to increase from 4 to above 5, that is from cool comfort to the most comfortable level. Meanwhile, the FAI increased dramatically above 2 and reached the Very Fresh level and fluctuated into the high level of Relatively Fresh. From May to September, the NMCI fluctuated from 5 to above 6, that is from the most comfortable to the warm comfort level, and the FAI mainly fluctuated within the Relatively Fresh level. When entering October, the NMCI decreased from above 6 to below 5, that is from warm comfort to the most comfortable level. The FAI increased drastically and reached the highest level (14.81) of the year, and most of FAI values were higher than 6 in the second half of October, that is it mainly fluctuated at the Very Fresh level. From November to December, the NMCI decreased from above 5 to below 4, that is from the most comfortable to the cold less comfortable level, while the FAI primarily fluctuated within the Relatively Fresh level until in the middle of December, when it increased dramatically

to more than 12, namely, in the Very Fresh level, then decreased quickly and fluctuated between the Fresh and Relatively Fresh levels (Figure 6A).

**Figure 5.** FAI and NMCI levels of Yongding Tulou (**A**,**B**) and Gutian Conference Site scenic spots (**C**,**D**) of Fujian Province, China.

FAI and NMCI of Gutian Conference Site scenic spot

**Figure 6.** (**A**) FAI and (**B**) NMCI of Yongding Tulou and Gutian Conference Site scenic spots of Fujian Province on each day of 2018.

In the Gutian Conference Site scenic spot, the FAI mainly fluctuated between −0.76 and 2.03, that is it fluctuated between Fresh and Super Fresh levels, while the Super Fresh level was lower than the Yongding Tulou scenic spot by approximately 7 times. The NMCI mainly fluctuated between 3.20 and 7.2, that is it fluctuated between cold less comfortable and hot less comfortable levels (Figure 5B). For the first five months of 2018, from January to May, the NMCI increased from more than 4 to more than 6, that is from cool comfort to the most comfortable level, and then to the warm comfort level. The FAI fluctuated between below −0.5 and 1, that is it fluctuated between the Fresh and Healing Fresh levels. From June to August, the NMCI increased from above 6 to above 7, that is from warm comfort to the hot less comfortable level. The FAI was mainly above 0.5 and some days even above 1, that is primarily at the Healing Fresh level. From September to November, the NMCI decreased from above 6 to 5, that is from warm comfort to the most comfortable level, while the FAI mainly fluctuated between 1 and 1.5, that is it fluctuated between Healing Fresh and Very Fresh levels. In December, the NMCI decreased and mainly varied between 4 and 4.2, that is mainly at the cool comfort level, while the FAI also decreased and mainly fluctuated above 0.5 or below 0.5, that is mainly at the Relatively Fresh and Healing Fresh levels (Figure 6B).

More precisely, within 24 Jieqi of the Yongding Tulou scenic spot and Gutian Conference Site scenic spot, the FAI-NMCI could be summarized as in Figure 7 for satisfying tourists' diverse fresh air–natural microclimate comfort preferences and demands through the provision of proper microclimate tourism products and services in each Jieqi, or different Jieqi with approximately similar air freshness and natural microclimate comfort degrees.

**Figure 7.** 24 Jieqi FAI-NMCI Index of Yongding Tulou and Gutian Conference Site scenic spots.

#### **4. Discussion and Conclusions**

Comparing FAI-NMCI in natural scenic spots of Fujian province, the Wuyi Mountain scenic spot had more days of upper-level air freshness and natural microclimate comfort degree than Baishuiyang scenic spot, which was mainly distributed within the middle fresh-comfort range. In cultural scenic spots of Fujian province, the Gutian Conference Site scenic spot had more days of middle-level air freshness and natural microclimate comfort degree than Yongding Tulou scenic spot, which had even more days of air freshness at the Not Fresh level.

Comparing fluctuating ranges of FAI of the four scenic spots in Fujian Province, the natural scenic spots had wider FAI fluctuation ranges, which was primarily determined by the higher concentration of negative oxygen ions and lower concentrations of PM2.5, PM10, and O3, while, the cultural scenic spots had relatively shorter FAI fluctuation ranges, which was profoundly influenced by the relatively lower concentration degree of negative oxygen ions and higher concentration degree of PM2.5, PM10, and O3.

Comparing fluctuating ranges of NMCI of four scenic spots in Fujian province, the NMCI mainly fluctuated between 3.20 and 7.2, which indicated that their natural microclimate comfort degree primarily fluctuated between cold less comfortable and hot less comfortable levels, except the Yongding Tulou scenic spot, which fluctuated between cold less comfortable and warm comfort levels.

In short, the four scenic spots in Fujian Province of China had different levels of air freshness degree and natural microclimate comfort degree in 2018. The natural scenic spots were mostly distributed in Healing Fresh, Very Fresh, and Super Fresh levels of FAI, with the most comfortable and comfortable levels of NMCI. The cultural scenic spots were mostly distributed in Relatively Fresh and Healing Fresh levels of FAI, with the most comfortable and comfortable levels of NMCI. Meanwhile, the FAI-NMCI of natural and cultural scenic spots also had significant differences within 24 Jieqi, which will promote dynamic and creative utilization of these resources in microclimate tourism development.

Therefore, the fresh air–natural microclimate of scenic spots has become a vital tourist attraction in the contemporary Chinese tourism market under the background of severe air pollution. Tourists with different physiological and psychological conditions will be more suited to different levels of air freshness and natural microclimate comfort levels, such as elderly tourists, pregnant women tourists, and children tourists will be vulnerable to polluted air environments and require more fresh and comfortable natural microclimates, such as Super Fresh and the most comfortable level of scenic spots, which has more negative oxygen ions (O2 −) and lower PM2.5, PM10, and O3, with a more comfortable temperature, humidity, wind speed, solar radiation, and precipitation.

Moreover, the dynamic fresh air–natural microclimate fluctuation revealed by FAI-NMCI of scenic spots will properly indicate the air freshness and natural microclimate comfort degree information for both the demand side and the supply side of the tourism market, which will profoundly promote fresh air–natural microclimate tourism resources' distribution and redistribution, accelerate fresh air–natural microclimate in-depth tourism resources combination and integration from a creative microclimate tourism development perspective, and further drive fresh-microclimate tourism development in destinations and scenic spots with comparative advantages of fresh air and a comfortable natural microclimate for triggering their in-depth transformation and improvement.

From a theoretical contribution perspective, this paper constructed the FAI-NMCI to harmoniously combine air resources and microclimate resources together, which could empirically evaluate air and microclimate tourism resource values for tourists, the tourism market, and the society as a whole. Moreover, it could figure out overlapping study areas of air tourism, climate tourism, health tourism, ecological tourism, and so forth, which have significant theoretical innovation through transdisciplinary, multidisciplinary, and interdisciplinary integration research in China's tourism market. This is a new research trend that has seldom been explored by researchers and could become an important research domain for future multidisciplinary, transdisciplinary, and interdisciplinary exploration.

From a socioeconomic practical application perspective, the creation and utilization of FAI-NMCI is significant for fresh air and comfortable microclimate tourism resources' distribution and redistribution in the tourism market. Meanwhile, it could also become a vital health tourism market indicator, which simply reveals the air freshness degree and microclimate comfort degree for tourists and some potential tourists pursuing health tourism experiences and long-term health-improving vacations. FAI-NMCI could indicate a scenic spot's air freshness and microclimate comfort levels specifically on different days of a year, which also could be calculated and summarized into 24 Jieqi, 12 months, 4 seasons, and so forth. Under contemporary China's high-quality tourism development background, this could provide accurate information not only for demand and supply sides of the tourism market but could also become an important and comprehensive index for related governmental management and NGOs.

The limitation of this paper is we only applied the FAI-NMCI to four scenic spots of China because of data collection limitations and the short time construction of the air freshness index infrastructure, with only 1–2 years of data sources. This new index has only been applied in Fujian Province rather than other provinces of China, which need to be explored in the future. This index could reveal not only the natural conditions of tourism resources but could also indicate tourism practices, tourists' and destinations' social behavior patterns, as well as socioeconomic development patterns of some destinations, which could be explored through introducing multidisciplinary methods in future research.

 All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was financially supported by the National Science Foundation of China (42101238), the start-up funding of Minjiang University (32304307), and the Project of Central Leading Local "Fujian Mental Health Human–Computer Interaction Technology Research Center" (2020L3024).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data are not publicly available due to the privacy of data policy of Fujian Meteorological Service Center.

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

#### **References**

### *Article* **Spatial-Temporal Characteristics of Arctic Summer Climate Comfort Level in the Context of Regional Tourism Resources from 1979 to 2019**

**Yutao Huang 1, Xuezhen Zhang 2, Dan Zhang 2, Lijuan Zhang 1,\*, Wenshuai Zhang 1, Chong Ren 1, Tao Pan 1, Zheng Chu 3 and Yuying Chen 1**

**Abstract:** In the context of global warming, a key scientific question for the sustainable development of the Arctic tourism industry is whether the region's climate is becoming more suitable for tourism. Based on the ERA5-HEAT (Human thErmAl comforT) dataset from the European Center for Mediumrange Weather Forecasts (ECMWF), this study used statistical methods such as climatic tendency rate and RAPS to analyze the spatial-temporal changes in Arctic summer climate comfort zones from 1979 to 2019 and to explore the influence of changes in climate comfort on Arctic tourism. The results showed the following: (1) With the increase in the Arctic summer temperature, the universal thermal climate index (UTCI) rose significantly from 1979 to 2019 at a rate of 0.457 ◦C/10a. There was an abrupt change in 2001, when the climate comfort changed from "colder" to "cool", and the climate comfort has remained cool over the past decade (2010–2019). (2) With the increase in Arctic summer temperatures, the area assessed as "comfortable" increased significantly from 1979 to 2019 at a rate of 2.114 × 105 km2/10a. Compared with the comfortable area in the 1980s, the comfortable area increased by 6.353 × 105 km2 over the past 10 years and expanded to high-latitude and high-altitude areas, mainly in Kola Peninsula, Putorana Plateau, and Verkhoyansk Mountains in Russia, as well as the Brooks Mountains in Alaska. (3) With the increase in Arctic summer temperatures, the number of days rated comfortable on 30% of the grid increased significantly from 1979 to 2019 (maximum increase: 31 days). The spatial range of the area with a low level of comfortable days narrowed and the spatial range of the area with a high level of such days expanded. The area with 60–70 comfortable days increased the most (4.57 × 105 km2). The results of this study suggest that global warming exerts a significant influence on the Arctic summer climate comfort level and provides favorable conditions for further development of regional tourism resources.

**Keywords:** Arctic; universal thermal climate index (UTCI); spatial-temporal changes; 1979–2019

#### **1. Introduction**

According to the sixth assessment report of the Intergovernmental Panel on Climate Change (IPCC), the average global surface temperature has increased by 1.09 ◦C over the past century [1]. The report in the Arctic Climate Impact Assessment also pointed out that, in the context of global warming, climate change in the Arctic is mainly manifested in the melting of ice sheets in the Arctic, the continuous decrease in the volume of snow and frozen soil in the Northern Hemisphere, the frequent occurrence of extreme weather events, and the significant changes in the ecosystem in the Arctic [2]. Accordingly, climate change turned into a severe environmental problem, and soon it will have significant impacts

**Citation:** Huang, Y.; Zhang, X.; Zhang, D.; Zhang, L.; Zhang, W.; Ren, C.; Pan, T.; Chu, Z.; Chen, Y. Spatial-Temporal Characteristics of Arctic Summer Climate Comfort Level in the Context of Regional Tourism Resources from 1979 to 2019. *Sustainability* **2021**, *13*, 13056. https://doi.org/10.3390/su132313056

Academic Editor: Oran Young

Received: 17 October 2021 Accepted: 18 November 2021 Published: 25 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

on our lives, our world, and also on tourism [3]. Climate change has affected the length and quality of the tourist season, the choice and consumption of tourist destinations, the decision making of tourist activities, and the operating costs of the tourism industry [4]. Climate change directly causes the redistribution of tourist flows, geographically and seasonally, and climatic conditions that are attractive to tourists may shift to high latitudes and high altitudes [5]. Climate change also indirectly affects tourism by causing water shortages, floods, deterioration of water quality, changes in natural landscapes, reduction in biodiversity, forest fires, and spread of diseases [6,7]. Global warming could lead to the spatial-temporal redistribution of climate comfort [8]. The Arctic is the highest latitude region of the Northern Hemisphere, as well as a region sensitive to global climate change [9]. The heating rate of the Arctic is twice that of the global average [10]; thus, analyzing the characteristics of climate change in the Arctic is of great significance for the development of Arctic tourism resources.

In recent years, tourism in high-latitude regions has become increasingly popular. For instance, the number of tourists in Greenland increased from 26,410 in 1999 to 57,223 in 2008, which was equivalent to the total population of Greenland [11]. The number of summer tourists to Alaska also increased significantly [12]. The average annual growth rate in the number of tourists in Iceland has doubled over the past few years [13]. The increase in the climate comfort level in Norway has led to a huge increase in summer tourism; the number of tourists in the Svalbard Islands, for example, has quadrupled [14]. This significant increase in the number of tourists indicates that the climate comfort level in the Arctic has improved due to global warming, and the need to carry out quantitative studies on these changes is urgent. Roshan et al. [15] used the physiologically equivalent temperature (PET) to identify and evaluate bioclimatic conditions of 40 meteorological stations in Iran. It showed that Iran's comfortable days transferred to early spring and late autumn and an increase in the frequency of thermal comfort condition was observed in almost half of the stations. Roshana et al. [16] also used stations data to show that both SET, PET, and UTCI increased in Iran from 1995 to 2014. Eludoyin et al. [17] have examined the variations in the thermal condition of Nigeria in terms of the temperature, relative humidity, effective temperature (ET), temperature–humidity index (THI), and relative strain index (RSI) and concluded that thermal stress has increased in Nigeria from 2000 at most stations, especially in the south and north-western regions. Feng et al. [18] used the improved classic THI model to analyze the spatiotemporal evolution of annual and monthly average climate comfort in China from 2005 to 2018. The results showed that the annual average comfort level in China changed from cold to comfortable. Wu et al. [19] used 591 stations data to calculate the change of China's climate comfort index values between 1966 and 2016 and came to the same conclusion. In China, the climate comfort levels in high-altitude regions, plateaus, and mountains have improved significantly [20–24]. Zeng et al. [25] investigated the spatial and temporal variations of the universal thermal climate index (UTCI) of the China–Pakistan Economic Corridor (CPEC) from 1979 to 2018 and concluded that the UTCI generally exhibited a positive trend of 0.33 ◦C/10a over the past 40 years, and the seasonal variation characteristics of the UTCI showed an upward trend in all four seasons. The tourist comfort zone in North America moved from low altitude to high altitude and from low latitude to high latitude [26]. For example, southern Florida, Arizona, and particularly Mexico have decreased the warm tourist season and Canada, Georgia, and South Carolina have extended the warm tourist season. Alaska's tourist season also extended, which now starts 10 days earlier than in the 1940s [20]. Some scholars have analyzed the thermal comfort in different regions of Europe. Napoli et al. have assessed the thermal bioclimate of Europe for the summer season and concluded that an increase in heat stress up to 1 ◦C was observed from 1979 to 2016 [27]. Some scholars pointed out that southern Europe (Crete Island in Greece, Eastern Mediterranean, Athens, Madrid (Spain), Novi Sad (Vojvodina, Northern Serbia)) were under high temperature and pressure, and climate warming has led to a significant increase in the frequency and duration of high-temperature heat waves in southern Europe and a decline in climate comfort [28–35]. In central Europe, climate

change has led to an increase in the number of strong and very strong heat stress days in Poland [36]. Tourists from Germany, Ireland, and the UK would spend more holidays in their home country, and the tourism of Ireland and Britain has moved north [37]. In northern Europe, warming has increased the length of nature and wilderness tourism and the "midnight sun" season in Finland and Sweden [5], and increased comfort in Norway has led to an increase in summer visitors [14]. Climate change could considerably enhance northwestern Europe's comfort level [38,39]. It can be seen that climate changes have led to the redistribution of tourism climate resources in time and space. The existing studies are mostly based on station data analysis, with lower spatial resolution. The area of study is mostly concentrated on the low and middle latitudes of the Northern Hemisphere, but little attention has been paid to the changes in comfort level in the Arctic, where there have been the highest rates of temperature increase.

The climate comfort index, which is widely used in public health, urban planning, tourism, and leisure, is an important index for quantifying the climatic suitability of tourism destinations [40] and as a basis for studying the effect of climate change on tourism [41]. With the measurement of the environment variables, direct indices were initially drawn by the researchers. Including wet-bulb temperature [42], Kata temperature [43], wet-bulb globe temperature (WBGT) [44], temperature–humidity index (THI) [45], effective temperature (ET) [46], and wind chill index (WCI) [47], direct indices were widely used because they are simple and easy to use. As the study of climate comfort went further, the researchers realized that besides environmental factors such as temperature and humidity, behavioral variables such as metabolic rate and the thermal insulation and moisture permeability of clothing that affected the heat comfort should also be considered [48]. Thus, some researchers raised experiential indices, such as the predicted four-hour sweat rate (P4SR), that included the behavioral variables [49]. Because of the lack of the scientific simulation of human response in experiential indices, thermodynamic theorem was used to simulate the heat exchange between the human body and environment, by which climate comfort indices such as physiological equivalent temperature (PET) [50] and the universal thermal climate index (UTCI) were raised. Comparing with the PET index, the environmental and behavioral factors that affected the climate comfort indices were more fully considered in the UTCI. The human multi-joint reaction to the external environment was in the model of climate comfort index to simulate the dynamic physiological response of the pedestrian and construct the thermal stress index the human body accepted during exposure outdoors. According to the thermal stress index, the UTCI is divided into 10 comfort levels from hot to cold [51]. The UTCI is considered the most advanced, comprehensive, and universal human comfort index [52,53], and it is applicable to various climates, seasons, temporal, and spatial scales, as well as being widely used in urban zoning and tourism planning [54,55]. The ERA5-HEAT (Human thErmAl comforT) dataset from the European Center for Mediumrange Weather Forecasts (ECMWF) was used in this study. The objectives of the present study were to (1) determine the spatial pattern of Arctic summer climate comfort, (2) reveal the spatial-temporal characteristics of Arctic summer climate comfort, and (3) discuss the influence of climate change on Arctic summer climate comfort. It provides a scientific basis for sustainable development of tourism in the region.

#### **2. Materials and Methods**

#### *2.1. Study Area*

In this study, 800 km2 of the northern land in the Arctic Circle (66◦34 N) was taken as the study area, including parts of Greenland (territory of Denmark), Canada, Alaska (the United States), Russia, Norway, Sweden, Finland, and Iceland. Located in the North Frigid Zone, the study area is covered by sea ice year-round and experiences polar weather (Figure 1). The average temperature in the coldest month reaches −40 ◦C–−20 ◦C, so the area is an important cold source for the world. The weather is characterized by a long, extreme cold winter, and a short, cool summer; the annual average precipitation is 100–250 mm, which reaches 500 mm in the Greenland sea area, and summer precipitation is usually in the form of rain. Due to its climate characteristics, the Arctic region has unique tourism resources, and current tourism in the region mainly focuses on glaciers, auroras, gorges, endemic species, skiing, underwater diving, and boating on the ice sea.

**Figure 1.** Location of the study area (shaded area) and the average temperature in the Arctic from 1979 to 2019.

#### *2.2. Datasets*

1. Thermal comfort data: The climate comfort data were drawn from the ERA5-HEAT dataset of ECMWF from 1979 to the present, with a spatial resolution of 0.25◦ × 0.25◦ and a time interval of one hour. The UTCI from 1979 to 2019 was analyzed in this study, and the data were downloaded from the ERA5-HEAT website Available online: https://cds.climate.copernicus.eu/cdsapp#!/dataset/ (accessed on 28 May 2020). According to the thermal physiological response of the human body, the UTCI thermal stress was divided into 10 levels (Table 1).

**Table 1.** UTCI equivalent temperatures categorized in terms of thermal stress and thermal perception.

2. Arctic climate data: To analyze the climate background of the changes in Arctic comfort and analyze the changing factors of the UTCI, we chose temperature, wind speed, relative humidity, and sunshine duration for the analysis. These climate elements were obtained from the ERA5 dataset. Again, the time range drawn from ERA5 was from 1979 to the present, the spatial resolution was 0.25◦ × 0.25◦, and the time interval was one hour. The temperature and wind speed data were used directly, whereas the sunshine duration and relative humidity required further processing: sunshine duration was replaced by total cloud cover (TCC) [56] and relative humidity was calculated by the average temperature and dew-point temperature. The specific equation is as follows:

$$f = \frac{E}{EW} \times 100\% \tag{1}$$

$$E = 6.11 \exp 5417.753 \langle (1/273.16) - [1/(\; 273.16 + Td)] \rangle \tag{2}$$

$$Ew = 6.11 \exp \text{5417.753} \langle (1/273.16) - [1/(\ 273.16 + T)] \rangle \tag{3}$$

where *f* is relative humidity, *E* is actual vapor pressure, *Ew* is saturated vapor pressure, *Td* is dew point, and *T* is temperature.

#### *2.3. Methods*

2.3.1. Trend Analysis Method

A univariate linear regression equation of the cold region area variable (y) and the corresponding time (x) was established:

$$y = a\mathbf{x} + b \quad \text{( $i = 1, 2, \dots, n$ )}\tag{4}$$

where *a* is the linear regression coefficient indicating the rate of change in the spatial extent of the cold region. The positive or negative value of *a* indicates that the spatial extent of the cold region is increasing or decreasing over time. The positive value of *a* indicates the spatial extent of cold region is increasing while the negative value of a indicates the spatial extent of cold region is decreasing.

#### 2.3.2. Rescaled Adjusted Partial SUMS (RAPS) Method

The rescaled adjusted partial sums (RAPS) visualization approach can be used for preliminary visual inspection of a time series, to gain a feel for the data, and/or to guide and focus subsequent statistical tests and analyses [57]. The RAPS of the *Xk* are defined as follows:

$$X\_k = \sum\_{t=1}^k \frac{Y\_t - \overline{Y}}{S\_Y}; k = 1, \dots, n \tag{5}$$

where *Y* = [*Yt; t =* 1*,...,n*] represent a time series of a climatic variable, *Y* is sample mean, *SY2* is variance, *n* is number of values in the time series, and *k* is counter limit of the current summation.

#### 2.3.3. Correlation Analysis

Correlation analysis is a statistical method used to discover if there is a relationship between two variables/datasets, and the relatedness and negative/positive correlation of this relationship. The Pearson correlation method is adopted in this paper to analyze the relationship between comfortable area, comfortable days, and climate indicators, quantitatively. In each pair, the Pearson's correlation coefficient (*r*) is calculated as:

$$r = \frac{\sum\_{i=1}^{n} (x\_i - \overline{x})(y\_i - \overline{y})}{\sqrt{\sum\_{i=1}^{n} (x\_i - \overline{x})^2} \sqrt{\sum\_{i=1}^{n} (y\_i - \overline{y})^2}} \tag{6}$$

where *x* represents the average comfortable areas or days and *y* is the average of a given climate variable.

The study based on the UTCI data in the ERA5-HEAT dataset provided by ECMWF used MATLAB software for spatial data calculation and overlay analysis (Figure 2). We used the Pcolor function of MATLAB to achieve graphical visualization.

**Figure 2.** The framework of the study.

#### **3. Results**

#### *3.1. Change of UTCI in the Arctic Summer (1979–2019)*

The average UTCI value in the Arctic summer from 1979 to 2019 was 0.025 ◦C and the comfort level was "cool". The UTCI value increased significantly at a rate of 0.457 ◦C/10a (*p* < 0.01), for a total increase of 1.874 ◦C, and the thermal stress level also increased (Figure 3a). The results of RAPS analysis showed (Figure 3b) that there was an obvious change in the UTCI value in the Arctic summer in 2000, which indicates that a transformation occurred in the comfort level in the early 21st century. The UTCI anomaly map (Figure 3c) also shows that the negative anomaly changed to a positive anomaly in the early 21st century (2001) and the UTCI changed from a relatively low value stage to a relatively high value stage. Further analysis reveals that the UTCI value was −0.501 ◦C and the comfort level was "colder" before the abrupt change (1979–2000). After the abrupt change (2001–2019), the UTCI value was 0.634 ◦C and the comfort level was "cool", which indicates that there was an obvious change in the comfort level of the Arctic summer in the 21st century.

**Figure 3.** Change in Arctic summer UTCI from 1979 to 2019: (**a**) time series, (**b**) RAPS, (**c**) anomaly series.

#### *3.2. Change in the Area Rated "Comfortable" in the Arctic Summer*

It can be seen from Section 3.1 that the spatial distribution of comfortable areas expanded with the increase in the UTCI value for the Arctic summer from 1979 to 2019 (Figure 4a). The area of summer comfortable zones in the Arctic from 1979 to 2019 was 1.541 × 106 km2, covering 19.152% (i.e., nearly 1/5) of the total Arctic land area. The comfortable area in the decade before and after was 1.218 × 106 and 1.864 × 106 km2, respectively. The summer comfortable area was largest in 2019, at 2.44 × 106 km2, and smallest in 1996, at only 0.921 × 106 km2, with an annual variation coefficient of 23.304%. The increase rate for the summer comfortable area in the Arctic from 1979 to 2019 was 2.114 × 105 km2/10a, which was significant (*p* < 0.01). The total increase area was 8.667 × 105 km2, which indicates that the summer comfortable area in the Arctic has been expanding. The RAPS analysis results indicate (Figure 4b) that there was an obvious change in the size of the Arctic summer comfortable area in 2000, which was consistent with the change in the UTCI value. The comfortable area before and after the abrupt change was 1.307 × 106 and 1.811 × 106 km2, respectively. This finding was further verified by the area anomaly (Figure 4c), which indicates that there was a significant expansion in the Arctic summer comfortable area in the early 21st century.

**Figure 4.** Change in Arctic summer comfortable area from 1979 to 2019: (**a**) time series, (**b**) RAPS, (**c**) anomaly series.

The spatial distribution of the areas with different comfort levels in the Arctic over the past ten years was further analyzed (Figure 5). The Arctic summer comfortable area from 2010 to 2019 was 1.864 × 106 km2, accounting for 23.2% of the Arctic land area. These comfortable areas were mainly distributed in Eurasia and North America (66.5◦–70 ◦N), mid-eastern Russia, Norway, Sweden, and Finland in northern Europe, Alaska and the Yukon territory in Canada, as well as scattered areas along the coast of Greenland and islands of Canada, north of 70 ◦N (Figure 5a). Compared with the ten years from 1979 to 1988, the Arctic summer comfortable area increased by 6.353 × 105 km2, accounting for 7.9% of the Arctic land area (Figure 5b). In terms of spatial distribution, the comfortable area mainly expanded to high-latitude and high-altitude regions (66.5◦–70 ◦N) and, specifically, in the Kola Peninsula, Putorana Plateau, Verkhoyansk Mountains in Russia, and Brooks Mountains in Alaska (Figure 5c). It is, therefore, possible for tourism to be developed in the Arctic regions rated comfortable, especially those newly designated as comfortable. Sightseeing in gorges, fishing, and watching rare species can be carried out in the peninsula regions and adventures, such as polar mountain climbing and outdoor expeditions, can be developed in the mountainous and plateau regions.

**Figure 5.** Spatial distribution and variation of Arctic summer comfortable area: (**a**) 2010–2019, (**b**) 1979–1988, (**c**) spatial variation of the decade before and after.

#### *3.3. Change in the Number of Arctic Summer Comfortable Days from 1979 to 2019*

There was an average of 23 Arctic summer days rated comfortable from 1979 to 2019, accounting for a quarter (25%) of the entire summer, while in the decade before and after there were 21 and 26 such days, respectively. The largest number of comfortable days was seen in 2019 (30 d), and the smallest was 18 d in 1992. The annual variation coefficient was 10.955%. The increase rate in comfortable days from 1979 to 2019 was 1.458 d/10a, which was significant (*p* < 0.01), and the total increase was 6 d, which indicates that the season for Arctic summer tourism activities has been significantly extended (Figure 6a). The RAPS analysis results showed that there was an obvious change in the number of comfortable days around 2000, which is consistent with the changes in UTCI values (Figure 6b). The number of comfortable days before and after the abrupt change was 22 d and 25 d, respectively. This finding was further verified by the anomaly of days (Figure 6c), indicating that a significant change occurred in the number of comfortable days in the early 21st century.

**Figure 6.** Change in the number of Arctic summer comfortable days from 1979 to 2019: (**a**) time series, (**b**) RAPS, (**c**) anomaly series.

Figure 7 shows the spatial distribution of Arctic summer comfortable days from 1979 to 2019. The number of days in all grid cells rated comfortable has clearly been increasing, the area with a statistically significant increase accounted for 30% of the total Arctic area, and the increase rate ranged −1.972–7.498 d/10a. The largest increase in the number of comfortable days was in the Qeqertarsuaq region of western Greenland (31 d) and the second largest increase was in northern and eastern Siberia of Russia (8–24 d). The summer tourism seasons were extended in the areas with a large increase in the number of

comfortable days, and the tourism structure could be adjusted appropriately to invest in tourism infrastructures to meet the needs created by this extension.

**Figure 7.** Rate of spatial change for the number of Arctic summer comfortable days from 1979 to 2019.

These Arctic summer comfortable days not only appear in the increase in grid cells but also in the change of space areas with the same level of comfortable days. Table 2 displays the area distribution of comfortable days at various frequency levels between the period of 1979–1988 and 2010–2019. The areas with 0–10, 10–20, 20–30, and 40–50 comfortable days decreased from 1979 to 1988 and from 2010 to 2019, with regions experiencing only 0–10 comfortable days declining most (by 2.96 × 105 km2). Areas with 30–40, 50–60, 60–70, 70–80, and 80–90 comfortable days increased, with areas having 60–70 comfortable days increasing most (4.57 × 105 km2).

**Table 2.** Areas of comfort days at different levels in 1979–1988 and 2010–2019.

In the past ten years (2010–2019), there was an average of 26 comfortable Arctic summer days, which were distributed in North America and Eurasia, excluding Greenland (Figure 8a). These basically followed a zonal distribution and gradually decreased with an increase in latitude. Mid-eastern Russia, northern Europe, Canada, and Alaska (66.5◦–70 ◦N) had a large number of comfortable days (maximum: 90 d) and the regions north of 70 ◦N had a smaller number of comfortable days (average: 20 d), accounting for 21.74% of the Arctic summer. Compared with the early 1980s (Figure 8b), the number of comfortable days increased by 5 d.

In summary, there was an obvious increase in the range and number of comfortable days (between 10 d and 30 d) in the Arctic, which were mainly distributed in mid-eastern Russia (66.5◦–70 ◦N), Greenland, and Ellesmere Island, north of 80 ◦N (Figure 8c). The increase in the number of comfortable days extends the summer tourism season in the Arctic, providing tourists more time to enjoy the natural scenery and creating the potential for higher economic benefits.

**Figure 8.** Spatial distribution and change in the number of Arctic summer comfortable days from 1979 to 1988 and from 2010 to 2019: (**a**) 2010–2019, (**b**) 1979–1988, (**c**) Spatial transformation.

#### **4. Discussion**

As an area sensitive to global warming, the climate comfort level in the Arctic has changed with the warming climate, which has had positive effects on the regional tourism industry. Therefore, it is of great significance to study the changes in the areas of the Arctic that can be rated comfortable during the summer, with a particular focus on the influence of climate comfort level on tourism in the context of global warming. In this study, a quantitative analysis was carried out on the spatial-temporal evolution of the Arctic summer climate comfort zones from 1979 to 2019, as well as the influence of the change in the climate comfort levels on Arctic tourism. The results of this analysis may be helpful in optimizing the Arctic tourism industry, promoting the development of tourism destinations, and providing scientific and quantitative support for the development of Arctic tourism. The following aspects were discussed in this study.

Factors that influence the climate comfort level include temperature, relative humidity, wind speed, and sunshine duration. The relationship between meteorological elements and climate comfort level was analyzed based on the ERA5-HEAT dataset. Figure 9 shows the interannual changes in Arctic summer average temperature, relative humidity, wind speed, and TCC from 1979 to 2019, which had averages of 3.809 ◦C, 79.461%, 4.128 m/s, and 70.249%, respectively. The temperature rose significantly from 1979 to 2019 (*p* < 0.01), with a rate of 0.418 ◦C/10a and a total increase of 1.714 ◦C. The relative humidity, wind speed, and TCC fluctuated greatly, but the change was not significant (*p* > 0.05), which indicates that the temperature increase was the main climate change factor in the Arctic.

**Figure 9.** Changes in the average temperature, relative humidity, wind speed, and total cloud cover in the Arctic summer from 1979 to 2019.

The correlation between the comfortable area and the number of comfortable days and various meteorological factors was also analyzed (Table 3). The comfortable area and number of days were positively correlated with the average temperature (*p* < 0.01) and negatively correlated with TCC; the positive correlation coefficient was much larger than the negative correlation coefficient. Thus, it appears that the comfortable area and days are mainly affected by the temperature and TCC, with the temperature being the main factor.

**Table 3.** Correlation coefficients between comfortable area, comfortable days, and meteorological factors in arctic summer from 1979 to 2019.

Note: "\*", "\*\*": Significance at 0.05 and 0.01 levels.

The spatial correlation coefficients between the number of comfortable days and temperature, relative humidity, wind speed, and TCC in the Arctic summer are analyzed (Figure 10). It showed that there are significant spatial correlations between comfortable days and temperature, relative humidity, wind speed, and TCC, which are basically distributed between 66.5–70 ◦N and can reach 80 ◦N in areas such as Canada's Queen Elizabeth Islands and Ellesmere Island, indicating that the Arctic summer comfortable days of 66.5–70 ◦N is related to temperature, relative humidity, wind speed, and TCC. However, there are differences between different elements in the relevant regions. Statistics showed that the areas where summer comfortable days are significantly correlated with temperature, relative humidity, wind speed, and TCC are 5.719 × 106, 3.24 × 106, 3.974 × 106, and 5.412 × 106 km2, respectively. The average significant correlation coefficients are 0.776, −0.447, −0.506, and −0.609, respectively, indicating that the Arctic summer comfortable days have a more significant correlation with temperature, and temperature is the most important factor affecting summer climate comfort.

**Figure 10.** Spatial correlation coefficients between the number of comfortable days and temperature, relative humidity, wind speed, and TCC in Arctic summer from 1979 to 2019 (R: correlation coefficient).

The spatial pattern of and change in climate comfort levels in the peri-Arctic nations and regions has been previously studied. For instance, the spatial pattern of the summer UTCI index in Russia from 2001 to 2015 was analyzed by Vinogradova et al. [58] based on data from 500 meteorological stations in Russia. The results showed that the comfortable zones were distributed in most areas of Russia in summer, and a moderate thermal stress level was observed in daytime in southern Europe. In a study by Yu et al., the hourly weather data from meteorological stations were used to calculate the modified climate index for tourism [20]. The authors suggested that the weather conditions in parts of regions in Alaska had improved greatly from 1942 to 2005, which extended the tourism season. Based on the temperature, precipitation, TCC, and visibility data, Forland et al. found that "warm" weather in Norway, north of the Arctic Circle, increased from 1981 to 2010, the summer tourism season was extended, and the improvement in the summer climate comfort level led to an increase in the number of tourists, quadrupling in the Svalbard Islands [14]. According to statistics, from 1990 to 2005, the annual growth rate of Iceland's Arctic cruise tourists was 19.19%, from 7952 in 1990 to 266,070 in 2015, an increase of more than 30 times; in 2019, the number of overnight tourists in Greenland increased to 264,711 people, equivalent to 3.47 times the total population of the island; from 1996 to 2018, cruise tourism in Svalbard also increased greatly, and cruise tourists increased

by 1.94 times [59]. These studies, to some extent, verify the findings of the present study, which showed that the Arctic summer comfortable zones have expanded and the tourism areas increased, while the number of comfortable days has significantly increased and the summer tourism season has been extended.

The results of this study showed that the Arctic summer comfortable zones have expanded and there are more comfortable days, which has created favorable conditions for tourism. In recent years, Norway, Iceland, Russia, and Canada have set up Arctic tourism routes and promoted the development of ecological tourism in the region [60]. The Arctic is the region least affected by human activities in the world, and unplanned or inadequate development of tourism resources could eventually lead to the degradation of the natural environment. Tourists are also likely to bring new species and germs into the region [61]. To achieve the sustainable development of Arctic tourism, the various countries involved should adhere to the concept of global integration and cooperate to complete legal systems for Arctic tourism development to ensure ecological balanced and sustainable growth in the region.

#### **5. Conclusions**

The average summer UTCI value in the Arctic was 0.025 ◦C from 1979 to 2019 and the comfort level was cool. The UTCI increased significantly at a rate of 0.457 ◦C/10a and there was an abrupt change in 2001. The Arctic summer comfort level changed from colder to cool after this abrupt change.

The average comfortable area in the Arctic summer from 1979 to 2019 was 1.541 × 106 km2, accounting for a fifth of the Arctic land area. These comfortable zones were mainly distributed in mid-eastern Russia, Norway, Sweden, and Finland in northern Europe, as well as Alaska and the Yukon territory in Canada. The comfortable area increased significantly at a rate of 2.114 × 105 km2/10a from 1979 to 2019. Compared with the 1980s, the comfortable area has increased by 6.353 × 105 km2 in recent decades, mainly expanding to high-latitude and high-altitude regions, specifically distributed in the Kola Peninsula, Putorana Plateau, Verkhoyansk Mountains in Russia, and Brooks Mountains in Alaska.

The average number of Arctic summer comfortable days was 23 from 1979 to 2019, accounting for a quarter of the entire summer. These were distributed zonally and gradually decreased from low to high latitudes. The number of comfortable days increased significantly from 1979 to 2019, at a rate of 1.458 d/10a. Compared with the early 1980s, the number of comfortable days has increased by 5 d over the past decades. The range and number of summer comfortable days has increased significantly, and the largest number of comfortable days was seen in the Qeqertarsuaq region of western Greenland (31 d), while the area seeing 60–70 summer comfortable days expanded the most (by 4.57 × 105 km2).

This study is only limited to the Arctic region and does not conduct a detailed analysis of the Northern Hemisphere or typical Arctic tourism cities. It has not conducted an in-depth study on how the Arctic climate comfort changes under the background of future climate change. It will be subject to special research in the future.

The improvement of Arctic summer climate comfort and the extension of suitable tourist seasons will enhance the attractiveness of Arctic tourist destinations while promoting the development and utilization of Arctic natural resources. Sightseeing in gorges, fishing, and watching rare species can be carried out in the peninsula regions and adventures, such as polar mountain climbing and outdoor expeditions, can be developed in the mountainous and plateau regions. Despite vigorously developing Arctic tourism activities, it is necessary to protect and manage Arctic tourism resources to ensure the sustainable development of Arctic tourism.

**Author Contributions:** Y.H. analyzed the data and drafted the manuscript; L.Z., D.Z. and X.Z. completed the manuscript and made major revisions; W.Z. completed data visualization; T.P. and C.R. downloaded data and searched references; Y.C. and Z.C. checked and proofread the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by National Natural Science Foundation of China (Grant No. 41771067. Supported by the National Key R&D Program of China (Grant No. 2017YFA0603301).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Acknowledgments:** Thanks to ECMWF for providing the ERA5-HEAT dataset.

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

#### **References**

### *Article* **How the Updated Earth System Models Project Terrestrial Gross Primary Productivity in China under 1.5 and 2** *◦***C Global Warming**

**Chi Zhang 1, Shaohong Wu 1, Yu Deng 2,\* and Jieming Chou 3**

**Abstract:** Three Earth system models (ESMs) from the Coupled Model Intercomparison Project phase 6 (CMIP6) were chosen to project ecosystem changes under 1.5 and 2 ◦C global warming targets in the Shared Socioeconomic Pathway 4.5 W m−2 (SSP245) scenario. Annual terrestrial gross primary productivity (GPP) was taken as the representative ecological indicator of the ecosystem. Under 1.5 ◦C global warming, GPP in four climate zones—i.e., temperate continental; temperate monsoonal; subtropical–tropical monsoonal; high-cold Tibetan Plateau—showed a marked increase, the smallest magnitude of which was around 12.3%. The increase was greater under 2 ◦C of global warming, which suggests that from the perspective of ecosystem productivity, global warming poses no ecological risk in China. Specifically, in comparison with historical GPP (1986–2005), under 1.5 ◦C global warming GPP was projected to increase by 16.1–23.8% in the temperate continental zone, 12.3–16.1% in the temperate monsoonal zone, 12.5–14.7% in the subtropical–tropical monsoonal zone, and 20.0–37.0% on the Tibetan Plateau. Under 2 ◦C global warming, the projected GPP increase was 23.0–34.3% in the temperate continental zone, 21.2–24.4% in the temperate monsoonal zone, 16.1–28.4% in the subtropical–tropical monsoonal zone, and 28.4–63.0% on the Tibetan Plateau. The GPP increase contributed by climate change was further quantified and attributed. The ESM prediction from the Max Planck Institute suggested that the climate contribution could range from −12.8% in the temperate continental zone up to 61.1% on the Tibetan Plateau; however, the ESMs differed markedly regarding their climate contribution to GPP change. Although precipitation has a higher sensitivity coefficient, temperature generally plays a more important role in GPP change, primarily because of the larger relative change in temperature in comparison with that of precipitation.

**Keywords:** GPP; climate change; CMIP6; ESM

#### **1. Introduction**

Gross and net primary productivity (GPP and NPP, respectively) are representative indicators that reflect ecosystem production capacity [1–3]. Many previous studies have considered future GPP/NPP change. For example, Huang et al. [4] evaluated NPP variations in the 21st century under various climate scenarios using the Lund–Potsdam–Jena dynamic global vegetation model. They found that total NPP in China is projected to increase continuously under different scenarios, with CO2 concentration playing the dominant role. Using a machine learning model to constrain the spatial uncertainty in GPP projections, Schlund et al. [5] predicted a higher increase in GPP in northern high latitudes

**Citation:** Zhang, C.; Wu, S.; Deng, Y.; Chou, J. How the Updated Earth System Models Project Terrestrial Gross Primary Productivity in China under 1.5 and 2 ◦C Global Warming. *Sustainability* **2021**, *13*, 11744. https:// doi.org/10.3390/su132111744

Academic Editors: Xiaodong Yan, Jia Yang, Shaofei Jin and Adriana Del Borghi

Received: 18 July 2021 Accepted: 9 October 2021 Published: 24 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

over the 21st century under the Representative Concentration Pathway [6] 8.5 W m−2 (RCP8.5) in comparison with regions closer to the equator. Under 1.5 ◦C of global warming, the GPP in China is expected to increase by 15.5% ± 5.4% on a stabilized pathway and by 11.9% ± 4.4% on a transient pathway [7]. Zhang et al. [3] explored the trend features of GPP/NPP in the 21st century under the Shared Socioeconomic Pathway [8] 24.5 W m−2 (SSP245) with the Beijing Climate Model. Their results predicted the overall trends of increase in both the near-term and long-term terrestrial GPP/NPP. However, in certain districts, the trend of GPP/NPP showed an initial increase followed by a decrease. Wang et al. [9] investigated the variation in NPP over the 21st century using the Earth system models (ESMs) of the Coupled Model Intercomparison Project phase 5. The results obtained under all four RCP scenarios suggested an increasing trend of NPP over China, especially in western areas.

In summary, GPP/NPP in China under different scenarios is expected to show a trend of increase in the 21st century. However, large uncertainties exist in the various ESMs [5,9]. Under the global warming targets of 1.5 and 2 ◦C above preindustrial levels set by the Paris Agreement, many regional impacts wait to be assessed. In particular, as the Coupled Model Intercomparison Project enters into the 6th phase (CMIP6), more and more ESMs have distributed their latest climate simulation under the Shared Socioeconomic Pathways (SSPs). How the latest ESMs will project the future ecosystem change in China and the corresponding climate attribution remains to be determined and revealed. In a recent study on performance in presenting historical terrestrial GPP in China, three out of seven ESMs evaluated were found to perform well in terms of climatological GPP, spatial pattern, and the ecosystem–climate relationship [10]. Consequently, these three ESMs were chosen in this study to predict ecosystem change under the warming targets. The ecological indicator of annual GPP was applied to measure the general state of the ecosystem. Changes in annual GPP predicted using the different ESMs were quantified with respect to the different climate zones in China. Furthermore, the relationship between the ecosystem and climate variables was tested and built through linear correlation and multiple regression. Relying on the model-specific parameters of the ecosystem response to the climate, the climate-related GPP changes were revealed and quantified.

#### **2. Materials and Methods**

#### *2.1. Data*

In this study, three ESMs that performed well in historical GPP reproduction [10] were chosen to project future GPP in China: (1) the Beijing Climate Center Climate System Model (BCC-CSM2-MR) [11], (2) the Euro-Mediterranean Centre on Climate Change coupled climate model (CMCC-CM2-SR5) [12], and (3) the Max Planck Institute for Meteorology Earth System Model version 1.2 (MPI-ESM1.2-HR) [13]. Specifically, BCC-CSM2-MR and MPI-ESM1-2-HR, out of seven ESMs, gave the best estimation of climatological GPP in China from 1980 to 2013. MPI-ESM1-2-HR performed best in characterizing the spatial structure. BCC-CSM2-MR and CMCC-CM2-SR5 best captured the response of the ecosystem to the climate [10]. The land surface models used for the three ESMs were BCC-AVIM2.0, CLM4.5, and JSBACH3.2. Major improvements or parameterizations have been made to these models in comparison with their predecessors; they make use of new scientific understanding to better simulate vegetation phenology [12–14]. These ESMs could provide not only the monthly GPP, but also the monthly surface air temperature and precipitation. In CMIP6, new SSPs were employed for climate modelling. The SSPs included five narratives describing alternative socio-economic developments, such as sustainable development, fossil-fueled development, etc. [8]. The middle of the road development—i.e., the SSP2 scenario—featured a continuation of the current fossil fuel-dominated energy mixed with intermediate challenges for both mitigation and adaptation, which resembled the historical pattern most [8]. SSP245, as the sole scenario of SSP2 implemented in CMIP6, was thus chosen to represent the most possible future world. Historical data (1980–2013) were applied to determine the ecosystem–climate relationship—i.e., correlation and multiple linear

regression. Data from the BCC-CSM2-MR, CMCC-CM2-SR5, and MPI-ESM1.2-HR ESMs were output as 1.125◦ × 1.125◦, 0.9375◦ × 1.25◦, and 0.9375◦ × 0.9375◦ grids, respectively. Because the grids were not uniform, they were first transformed to a 1◦ × 1◦ grid through bilinear interpolation for comparative purposes.

A climate division map of China was applied for regional analysis. It divided the country into four climate zones—i.e., temperate continental, temperate monsoonal, subtropical– tropical monsoonal, and high-cold Tibetan Plateau, as in He et al. [15] and Zhang et al. [16].

#### *2.2. Methods*

#### 2.2.1. Bilinear Interpolation

Bilinear interpolation can produce a smoother interpolation than that achieved using the nearest neighbor method [17]. Thus, it was applied to transform fields from various grids of the ESMs into the formal 1◦ × 1◦ grid. In this approach, g (*n*1, *n*2) is defined as a linear combination of the values of its four nearest neighbors.  Area Weighting

Regional and global mean variables—e.g., temperature, precipitation, and GPP—on the 1◦ × 1◦ grid are calculated through area weighting:

$$V\_{\text{reg}} = \frac{\sum\_{i} V\_{i} \cdot \frac{\pi R}{180} \cdot \frac{\pi R}{180} \cos(\theta\_{i})}{\sum\_{i} \frac{\pi R}{180} \cdot \frac{\pi R}{180} \cos(\theta\_{i})} = \frac{\sum\_{i} V\_{i} \cos(\theta\_{i})}{\sum\_{i} \cos(\theta\_{i})},\tag{3}$$

where *θ* represents the latitude of the grid, *R* is the Earth's radius, and *V* is the variable.

#### 2.2.3. Linear Correlation and Multiple Regression

The correlation coefficient *r* is used to test the relationship between ecosystem productivity and climate factors. The formula can be expressed as follows:

$$r = \frac{n\sum\_{i=1}^{n} E\_i \mathbf{C}\_i - \sum\_{i=1}^{n} E\_i \sum\_{i=1}^{n} \mathbf{C}\_i}{\sqrt{n\sum\_{i=1}^{n} E\_i^2 - \left(\sum\_{i=1}^{n} E\_i\right)^2} \cdot \sqrt{n\sum\_{i=1}^{n} \mathbf{C}\_i^2 - \left(\sum\_{i=1}^{n} \mathbf{C}\_i\right)^2}}\tag{4}$$

where *E* and *C* represent ecosystem productivity and climate factors, respectively.

The interannual variation in GPP reflects year-to-year differences attributable mainly to climate variations [18,19]; therefore, the relationship between GPP and climate—i.e., precipitation and surface air temperature—was explored using a multiple regression approach [20]:

$$y = a\mathbf{x}\_T + b\mathbf{x}\_P + \varepsilon,\tag{5}$$

where *y* is the detrended anomaly of the carbon flux GPP, variable *xT* is the detrended annual temperature anomaly, and *xP* is the detrended annual precipitation anomaly. The fitted regression coefficients *a* and *b* define the apparent carbon flux sensitivity to interannual

variations in temperature and precipitation, and *ε* is the residual error term. The use of the detrended time series instead of the original nonstationary time series in the above linear correlation and regression analysis provides a robust estimate of their relationship [21–23].

Some definitions set 1986–2005 as a reference period when the global surface air temperature was 0.61 ◦C warmer than preindustrial levels [24,25]. We adopted this definition and defined the 1.5 and 2 ◦C warming periods as the first time when the 20-year-moving-average global temperature was 0.89 and 1.39 ◦C warmer, respectively, than that from 1986–2005 in the models. The corresponding changes in ecosystem and regional climate were based upon the reference period of 1986–2005. It is also necessary to point out that 20 years is a duration that is commonly applied to represent a climate state in the scientific world [26–28].

#### **3. Results and Discussion**

#### *3.1. GPP Distribution and Projected Changes*

The climatological GPP distribution produced by each of the three ESMs from 1986 to 2005 is shown in Figure 1. The ESMs all produce a similar spatial pattern of GPP, showing high (low) values in the southeast (northwest) of China. Regionally, the GPP in the subtropical– tropical monsoonal zone is largest, followed in descending order by the temperate monsoonal zone, Tibetan Plateau, and temperate continental zone. The three ESMs produced comparable estimates in the climate zones except in the monsoonal regions, where CMCC-CM2-SR5 produced larger estimates, especially in the subtropical–tropical monsoonal zone.

**Figure 1.** Climatological GPP distribution (unit: gC m−2 yr−1) from 1986 to 2005 from (**a**) BCC-CSM2-MR, (**b**) CMCC-CM2- SR5, and (**c**) MPI-ESM1-2-HR. The four delineated regions represent the climate zones temperate continental (TC), temperate monsoonal (TM), subtropical–tropical monsoonal (STM), and high-cold Tibetan Plateau (TP). (**d**) Regional GPP from the ESMs (PgC yr−1: petagram of carbon per year) (BCC: BCC-CSM2-MR; CMCC: CMCC-CM2-SR5; MPI: MPI-ESM1-2-HR). The bars and lines represent the standard deviation during the 20 years.

The GPP change under 1.5 ◦C of global warming is shown in Figure 2. Throughout China, the GPP of all three ESMs showed a positive anomaly except over certain individual grid points. In the subtropical–tropical monsoonal with MPI-ESM1-2-HR, the negatively changed grids tend to concentrate (Figure 2e). The areal GPP reduction may be related to reduced local precipitation. However, many factors could contribute to the GPP change in addition to precipitation and temperature, such as land-use change, soil moisture, wind speed, humidity, solar radiation, nitrogen deposition, etc. Thus, it is really hard to be conclusive. Moreover, the aggregated negative grids tend to dissipate under 2 ◦C global warming (Figure 3e). The GPP change patterns differ among the models (Figure 2). For example, the largest anomaly in the output of BCC-CSM2-MR appears over the southeastern Tibetan Plateau, whereas the largest anomalies in the output of CMCC-CM2-SR5 and MPI-ESM1-2-HR appear in the central and southern parts of the subtropical–tropical monsoonal zone, respectively. The absolute GPP change is largest in the subtropical–tropical monsoonal zone; however, the relative change is rather small and stable among the models (Figure 2 right). The absolute change is smallest in the temperate continental zone owing to its low base value in GPP. The relative change is large over the Tibetan Plateau, and there are strong differences in the magnitude of the GPP increments among the different models.

The GPP anomaly under 2 ◦C of global warming shows a spatial pattern similar to that found under 1.5 ◦C global warming but with a stronger intensity (Figure 3). Regional statistics indicate that the regional GPP changes will be larger under 2 ◦C of global warming. The projected increment of GPP in China under the different warming targets is consistent with previous findings [5,7,29]. This suggests that from the perspective of GPP, there is no ecological crisis in the projected future climate within the studied domain [3]. As with 1.5 ◦C warming, the subtropical–tropical monsoonal zone with the highest GPP value contributed the most to the increment in China's GPP under 2 ◦C of warming. However, the increase rate does not show much difference in magnitude when compared with that of other regions. It is worth noting that the rate of increase in GPP is substantial on the Tibetan Plateau—i.e., the increase is nearly 63% with regard to BCC-CSM2-MR. Thus, the Tibetan Plateau would appear to be the region most susceptible to the effects of climate warming, although the influence could be considered positive and beneficial.

The seasonal GPP anomalies under the 1.5 and 2 ◦C warming scenarios are shown in Figures 4 and 5, respectively. The spatial modes between the two warming scenarios are similar, noting that the magnitude in the 2 ◦C warming is much larger than in the 1.5 ◦C warming. In spring and summer, the GPP anomalies are the most prosperous, as they correspond to the growing season in China, while they drop to become the weakest in winter. In spring and winter with all ESMs, the GPP all over China generally shows a positive anomaly, with only sporadic negative points. Some negative changes occur in summer and autumn, especially with BCC-CSM2-MR and MPI-ESM1-2-HR. For BCC-CSM2-MR in summer, the negative GPP anomalies concentrate in the Huaihe River, which divides the subtropical–tropical monsoonal and temperate monsoonal regions. However, there were no negative changes in GPP at the zone scale. For MPI-ESM1-2-HR in summer and autumn, we observed some negative changes over the grassland in the temperate continental, which is similar to the results of Ma et al. [30]. They found that large areas in Northern China showed a decreasing trend in NPP under global warming, although the overall NPP increased significantly. The fact that only one ESM obtained similar results also indicates the large inter-model spread in representing the future GPP change. The negative changes in the temperate continental were weakened under the 2 ◦C warming scenario (Figure 5).

**Figure 2.** (Left) Spatial distribution of GPP anomaly (gC m−2 yr−1) and (right) the respective regional GPP changes under 1.5 ◦C of global warming. (**a**,**b**) BCC-CSM2-MR, (**c**,**d**) CMCC-CM2-SR5, and (**e**,**f**) MPI-ESM1-2-HR. Numbers in the right-hand panels represent the relative change in comparison with regional GPP from 1986 to 2005.

**Figure 3.** (Left) Spatial distribution of GPP anomaly (gC m−2 yr−1) and (right) the respective regional GPP changes under 2 ◦C of global warming. (**a**,**b**) BCC-CSM2-MR, (**c**,**d**) CMCC-CM2-SR5, and (**e**,**f**) MPI-ESM1-2-HR. Numbers in the right-hand panels represent the relative change in comparison with regional GPP from 1986 to 2005.

#### *3.2. Climate Attribution*

The variation in GPP is closely related to climate [20,31], and the correlation parameters within the studied ESMs are provided in Table 1. It can be seen that GPP is correlated significantly with at least one climate variable. There are cases in which GPP correlates negatively with temperature, such as in the temperate continental zone with MPI-ESM1- 2-HR and in the temperate monsoonal zone with BCC-CSM2-MR, which imply inherent differences in ecological modeling between the different ESMs [9,10]. On the Tibetan Plateau, it is unanimous within the ESMs that GPP is most closely related to temperature rather than to precipitation. Because the ESMs substantially overestimate precipitation over the Tibetan Plateau [32–34], it is possible that precipitation is not the primary climate factor constraining the regional ecosystem.

**Figure 4.** Seasonal GPP anomaly distribution under 1.5 ◦C of global warming (gC m−2 yr−1). The four rows from up to down represent spring, summer, autumn, and winter, respectively. The three columns from left to right represent BCC-CSM2-MR, CMCC-CM2-SR5, and MPI-ESM1-2-HR, respectively.

The corresponding climate changes—i.e., precipitation and surface air temperature— under the warming targets are shown in Tables 2 and 3, respectively. The mode of temperature change in the four regions is consistent among the models. The hottest region—i.e., the subtropical–tropical monsoonal zone—increases least under the effects of warming. Both CMCC-CM2-SR5 and MPI-ESM1-2-HR produced similar estimates of temperature change, whereas the estimates from BCC-CSM2-MR were larger, especially over the monsoonal regions.

**Figure 5.** Seasonal GPP anomaly distribution under 2 ◦C of global warming (gC m−2 yr−1). The four rows from up to down represent spring, summer, autumn, and winter, respectively. The three columns from left to right represent BCC-CSM2-MR, CMCC-CM2-SR5, and MPI-ESM1-2-HR, respectively.

There is a greater model variety regarding the change in precipitation. Under 1.5 ◦C of warming, there are negative changes—e.g., in the temperate monsoonal zone with CMCC-CM2-SR5 and in the subtropical–tropical monsoonal zone with both BCC-CSM2-MR and CMCC-CM2-SR5. Conversely, under 2 ◦C of warming, there are no negative changes, but the incremental differences for one certain region are huge. These results indicate the large uncertainty in the precipitation projections made by the ESMs. It is also worth noting

that in comparison with their variabilities, the change in temperature under the warming scenarios is reasonably large, whereas the precipitation change is rather limited [10].

The apparent sensitivity of climate to the ecosystem of each of the ESMs is shown in Table 4. The response of the ecosystem to climate varies strongly among the models. For some ESM regions, climate plays a crucial role, such that the variation in climate explains more than half of the variation in GPP. However, for certain other ESM regions, the degree of explanation attributable to climate is rather small—e.g., MPI-ESM1-2-HR in the subtropical–tropical monsoonal zone and BCC-CSM2-MR on the Tibetan Plateau. For one particular region, the same climate factor might affect the ecosystem differently in the various models. Taking the subtropical–tropical monsoonal zone as an example, precipitation is the major influencing factor and affects the ecosystem positively with BCC-CSM2-MR and CMCC-CM2-SR5. However, with MPI-ESM1-2-HR, the correlation between precipitation and GPP is insignificant and negative (Tables 1 and 4). Moreover, the overall climate contribution to ecosystem variation with MPI-ESM1-2-HR is very small (Table 4). These features further reflect the inherent differences of ecological modeling within the ESMs.

**Table 1.** Correlation parameter *r* between GPP and climate variables during the historical period of 1980–2013. \* denotes correlation that is significant at the 0.1 level; \*\* denotes correlation that is significant at the 0.05 level; \*\*\* denotes correlation that is significant at the 0.01 level.

**Table 2.** Regional temperature change (◦C) in comparison with that of the reference period of 1986–2005 under 1.5 and 2 ◦C (in parentheses) of global warming.

**Table 3.** Regional precipitation change (mm yr−1) in comparison with that of the reference period of 1986–2005 under 1.5 and 2 ◦C (in parentheses) of global warming.

**Table 4.** Sensitivity parameters of GPP to climate variables during the historical period of 1980–2013 through the multiple linear regression approach. The final row shows the explained variance of climate to GPP through multiple regression.

Based on the ecosystem–climate relationship (Table 4) and the known climate changes (Tables 2 and 3), the GPP change over the climate zones with the different ESMs is attributed quantitatively in Figure 6. As mentioned before, some contributions from the climate factors are negative, for which there are two major reasons. First, the climate change is negative—e.g., the negative precipitation change leads to a negative contribution to GPP. Second, the correlation between the climate factor and GPP is negative—e.g., temperature and GPP in the temperate continental zone in MPI-ESM1-2-HR. A positive anomaly in temperature could also lead to GPP reduction. The climate contribution to the variation in GPP changes among the studied ESMs. Even with the same model—e.g., MPI-ESM1- 2-HR—it can be −12.8% in the temperate continental zone and 61.1% on the Tibetan Plateau. Generally, under the effects of global warming, the influence of temperature on the ecosystem is larger than that of precipitation. This is mainly because the relative change in temperature is much larger than that in precipitation. On the Tibetan Plateau, where temperature is the most constraining factor (Tables 1 and 4), temperature plays a more dominant role than precipitation in the increase in GPP (Figure 6). In addition, it is observed that the relative GPP increase in the Tibetan Plateau is much larger than that in other regions. On the one hand, it is related to the low baseline value of GPP in the Tibetan Plateau (Figure 1d). A light increase in the GPP of the Tibetan Plateau is salient in relative values compared to the respective change in the subtropical–tropical monsoonal region. On the other hand, this may be related to the vegetation structure on the Plateau. Demonym plants can be divided into three types based on their photosynthesis patterns (i.e., C3, C4, and crassulacean acid metabolism). C3 photosynthesis produces a threecarbon compound during the Calvin cycle, while C4 photosynthesis makes an intermediate four-carbon compound that splits into a three-carbon compound for the Calvin cycle. They favor different conditions of nature. The conditions on the frigid Tibetan Plateau are unsuitable for the growth of C4 plants [35]. Consequently, the plateau is dominant by C3 plants [36]. C3 plants are more efficient in vegetative growth than C4 plants in response to the increasing air CO2 [37]. As a result, GPP increases more rapidly with increased air CO2 in the Tibetan Plateau than in other regions containing both C3 and C4 plants.

**Figure 6.** Climate contributions to GPP change with the different ESMs under (**a**) 1.5 and (**b**) 2 ◦C of global warming. Bars represent the absolute climate contributions (PgC yr−1). Percentages represent the relative climate contributions to the overall GPP changes.

#### **4. Conclusions**

To evaluate the GPP change under 1.5 and 2 ◦C of global warming, this study selected three CMIP6 ESMs (i.e., BCC-CSM2-MR, CMCC-CM2-SR5, and MPI-ESM1-2-HR) that performed well in historical GPP modeling; the principal conclusions derived are as follows:

Similarly, GPP is projected to increase by 12.3–16.1% in the temperate monsoonal zone, 12.5–14.7% in the subtropical–tropical monsoonal zone, and 20.0–37.0% on the Tibetan Plateau. Under 2 ◦C warming, the increase in GPP is projected to be even greater—i.e., 23.0–34.3% in the temperate continental zone, 21.2–24.4% in the temperate monsoonal zone, 16.1–28.4% in the subtropical–tropical monsoonal zone, and 28.4–63.0% on the Tibetan Plateau.

3. Climate change is projected to contribute positively to GPP change, except in the temperate continental zone with MPI-ESM1-2-HR. Although precipitation has larger sensitivity parameters, temperature generally plays a more important role in GPP change because of the larger change relative to its own variability in comparison with that of precipitation.