Do Protected Areas Improve Ecosystem Services? A Case Study of Hoh Xil Nature Reserve in Qinghai-Tibetan Plateau

Although there is a consensus that protected areas (PAs) can provide various ecosystem services, it is unclear whether protected areas effectively contribute to the preservation and enhancement of ecosystem services. We conducted a case study of the Hoh Xil Nature Reserve (HXNR) in Qinghai-Tibetan Plateau, China, in order to examine the effectiveness of PA in the conservation of ecosystem services. First, the dynamics of land use/land cover (LULC) were analyzed based on remotely sensed data sets. Then, the ecosystem service value (ESV) in the PA and non-PA were evaluated using a modified benefit transfer method that had been adjusted using normalized difference vegetation index (NDVI). Finally, hotspot analysis was implemented to reveal the ESV changes for the different districts of the PA. The results of the comparison experiment indicate that: (1) The ESV of the HXNR has considerably increased after it was designated as protected, which had been in decline in the previous stage. The ESVs in a near-by non-PA showed opposite results where the values initially increased but then dropped due to urban expansion and desertification. (2) The areas in HXNR with increased ESV significantly outnumbered the areas that had declining values from 1980 to 2018. For the non-PA, the areas that had increased ESV in 1980–1995 saw a decline in value in 1995–2008; moreover, new areas with decreasing ESV emerged in 2008–2018. (3) The HXNR was found to be more effective than non-PA in improving ecosystem services. (4) The core zone of the nature reserve demonstrated better effectiveness in ecosystem service preservation.


Introduction
Protected areas (PAs; e.g., national parks, nature reserves, wilderness areas, landscape protected areas) are geographical spaces with unique natural and cultural resources, where specific regulations are implemented to reduce anthropogenic pressure and to achieve long-term environmental conservation alongside their associated ecosystem services and cultural values [1,2]. Since the tenth meeting of the Conference of the Parties to the Convention on Biological Diversity (CBD) in 2010, the goals of PAs have been extended from mainly biodiversity purposes into preserving ecosystem services (including maintaining biodiversity) [3][4][5]. This development incorporates ecosystem services into veritable index [52] (e.g., the net primary productivity (NPP) of natural vegetation [55,56]) while overlooking spatial heterogeneity. However, NPP is not always positively correlated with ecosystem services [57]. In consideration of spatial heterogeneity of ecosystem services, we propose the use of the normalized difference vegetation index (NDVI) to adjust the equivalent factors, which would be used to estimate the ESV in the HXNR. In order to test the effectiveness of PAs spatially, hotspot analysis on ESV changes is employed, which is vital in integrating ESV into the implementation of conservation policies [29,58].
This study assesses the dynamics of LULC in the interior and exterior of the HXNR, before and after its formal establishment (periods of 1980-1995, 1995-2008, and 2008-2018) using land-use change indicators. The changes in the different types of ESV have been evaluated using an NDVI-based modified benefit transfer method. Hotspot analysis of ESV change was employed to examine the spatio-temporal characteristics of ESV. Finally, the ESV dynamics inside and outside of the HXNR were compared to determine the effectiveness of PA towards ecosystem service conservation. This study aims to address the following questions: (1) What were the spatio-temporal features of ESV in the HXNR due to LULC change from 1980 to 2018?; (2) Did the ESV in PA improve significantly compared to non-PA?; and, (3) Which district of PA achieved more effectiveness of ecosystem services preservation?

Study Area
The Hoh Xil Nature Reserve (34 • 19'-36 • 16'N and 89 • 25'-94 • 05'E) is located in the hinterland of the Qinghai-Tibetan Plateau, with an official size of 4.5 × 10 4 km 2 and includes part of the Zhiduo country in Qinghai Province, China ( Figure 1). It extends to the Kunlun Mountains at the south and the Ulan Ula Mountain at the north and runs from the Qinghai-Tibet Railway in the east to Tibet Autonomous Region in the west. The average altitude of the HXNR is above 4600 m and is geographically composed of gentle and low undulating lake planes, plains, and terraces [59]. The climate in the region is a typical alpine climate characterized by low temperatures and low precipitation [60]. The annual average temperature of the Hoh Xil Reserve is −10.0 to 4.1 • C, and the lowest recorded temperature is −46.4 • C; the average annual precipitation is 173.0 to 494.9 mm, which is mainly concentrated in summer. These conditions make living in Hoh Xil region extremely difficult, and given the region's inhospitable climate, the HXNR has preserved its original landscape and rare special species [45,46], providing unique and vital ecosystem services and contributing to ecosystem conservation and human well-being at the regional and global scales. However, while the region has been designated as a protected area since 1995, maintaining and improving the ecological environment remains to be a major challenge due to the changing climate [60] and pressures from human activities (e.g., grazing and infrastructure construction) [61].
For comparison, we selected Zhiduo county as control, which is a region not designated as a protected area and is located in the southeast of the HXNR with an area of approximately 1.08 × 10 4 km 2 . The area between the control area and HXNR is part of the Sanjiangyuan Nature Reserve, and will not be discussed in this study due to the limitation of data.

Data Sources and Processing
The LULC data for the period 1980-2018 were retrieved from the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (RESDC) (http://www.resdc.cn) [62,63]. The LULC data initially had seven categories and 25 sub-categories (see Table 1). The LULC data for 1980, 1995, and 2008 were generated using Landsat Thematic Mapper™/Enhanced Thematic Mapper (ETM) remotely sensed imagery, while the LULC data for 2018 were updated based on the 2015 LULC data from Landsat 8 remotely sensed imagery and generated by manual interpretation. Using field surveys and records for comparison, the overall identified accuracy of the LULC data sets for the seven land use/cover types reached 94.3%, while for the accuracy for the 25 subclasses types reached 91.2% [62,63]. The spatial resolution was 30 × 30 m, and the land use/land cover were reclustered into six primary types: forest, grassland, water bodies, wetland, urban, and desert.
The NDVI data set, which was used to revise the equivalent factors in this study, was downloaded from the Geospatial Data Cloud site of the Computer Network Information Center at the Chinese Academy of Sciences (http://www.gscloud.cn) [64].

Data Sources and Processing
The LULC data for the period 1980-2018 were retrieved from the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (RESDC) (http://www.resdc.cn) [62,63]. The LULC data initially had seven categories and 25 sub-categories (see Table 1). The LULC data for 1980, 1995, and 2008 were generated using Landsat Thematic Mapper™/Enhanced Thematic Mapper (ETM) remotely sensed imagery, while the LULC data for 2018 were updated based on the 2015 LULC data from Landsat 8 remotely sensed imagery and generated by manual interpretation. Using field surveys and records for comparison, the overall identified accuracy of the LULC data sets for the seven land use/cover types reached 94.3%, while for the accuracy for the 25 subclasses types reached 91.2% [62,63]. The spatial resolution was 30 × 30 m, and the land use/land cover were reclustered into six primary types: forest, grassland, water bodies, wetland, urban, and desert.
The NDVI data set, which was used to revise the equivalent factors in this study, was downloaded from the Geospatial Data Cloud site of the Computer Network Information Center at the Chinese Academy of Sciences (http://www.gscloud.cn) [64].

LULC Dynamic Analysis
The spatiotemporal characteristics of the LULC in HXNR and the control area were analyzed using land-use dynamic degree index (LUDD), comprising the single LUDD (SLUDD) and comprehensive LUDD (CLUDD) [65]. SLUDD calculates the change rates of certain single land-use types, while CLUDD determines the overall change rate of LULC. SLUDD and CLUDD can be calculated using Equations (1) and (2), respectively: where A i,t1, and A i,t2 represent the areas of land use type i at the start (t 1 ) and end (t 2 ) of the study period, respectively; ∆A i−j is the area of land use type i converted to land use type j (j = 1,2, . . . n, i j) during the study period; n is the number of land-use types in the study area; and, T is the study period.

Calculation of ESVs
In this study, we adopted the value transfer method developed by Xie et al. [28] to estimate the ESVs inside and outside of the HXNR. The benefit transfer model proposed by Costanza et al. [15] and Xie et al. [28] has been commonly used to evaluate ESVs using monetary values designated to the different LULC types. Considering the inapplicability of the ecosystem services classification and equivalent coefficients proposed by Costanza et al. [15] for China, Xie et al. [28] proposed modifications to the ecosystem services classification and the equivalent factors table based on a survey of more than 700 ecological experts. In the modified version, the ESV for food production of croplands was set to 1, while the equivalent ESV coefficients for the other ecosystems were based on their relative importance to croplands. The value of each equivalent ESV factor was identified as 1/7 of the market value of the annual regional average grain yield, equal to 58.5 USD/ha in China. The equation to calculate ESV is given as follows: where VC i,k is the kth type of ESV coefficient for land use type i; A i represents the ith type of LULC; m and n indicate the number of sub-categories of ecosystem services and the LULC types, respectively. However, the value transfer method proposed by Xie et al. [28] overlooked the spatial heterogeneity [52,57]. Previous studies have enhanced Xie's model by considering biomass and vegetation cover index; however, they still have deviations on ESV evaluation [55,56]. In contrast, the NDVI index is shown to have a strong correlation with meteorological condition and local resource background, especially in areas with a large number of different vegetation covers [66], which is believed to reflect the spatial patterns of ecosystem background and associated with ecosystem services. Thus, we use NDVI to adjust the ESV coefficients to localize the equivalent factors in the Hoh Xil region compensating for the limitations of the previous model from neglecting spatial heterogeneity. The corrected ESVs were evaluated by Equations (4) and (5). where NDVI f, is the average annual normalized difference vegetation index of grid f ; NDVI max and NDVI min represent the maximum and minimum of the average annual NDVI in China; V f is the normalized NDVI f ; and, V is the average of normalized NDVI in China.

Hotspot Analysis of ESV Change
Hotspot analysis is a widely used to reveal the spatial locations of statistically significant high-value clusters (hotspots) and low-value clusters (cold spots). Theoretically, a hotspot indicates that its surrounding areas are statistically in high values, while a cold spot is surrounded by low value spots. Therefore, the effects of abnormally high and low values can be excluded [57,58]. In this study, we apply hotspot analysis to identify the distribution of hotspots and cold spots showing ESV changes in the HXNR. As the hotspots and the cold spots of ESV change illustrate the significant points of ESV increase and decrease, it would reveal the effectiveness of PA spatially. Based on the calculation of adjusted ESV for each LULC patch in different years, the change in ESV was assigned to each patch. To enhance feasibility and flexibly, the patches with modified ESV were selected and vectorized and were then applied to the Hotspot Analysis (using the Getis-Ord Gi* statistic) to reveal the statistically significant ESV change areas. The z-scores and p-values within the analysis indicate standard deviations and probabilities. The typical probabilities are 0.01, 0.05, and 0.1, and the z-scores at < −1.65 or > 1.65, < −1.96 or > 1.96, and < −2.58 or > 2.58 are critical for 90%, 95%, and 99% confidence levels. A positive and significant z-score indicates that the change in ESV in the vicinity is relatively in high value (hotspots of the ESV changed points). In contrast, a z-score that is negative and significant suggests that the change in ESV in the surrounding area is relatively in low value (cold spots of the ESV changed points).

LULC Change from 1980 to 2018
For the HXNR, due to its high altitude and harsh climate, the region is considered as not suitable for human habitation. As a result, there were no croplands or construction lands, as shown in Table 2. The main LULC types from 1980 to 2018 were desert and grassland ( Figure 2). Desert accounted for 55.98%, 55.58%, 56.03%, and 55.08% in 1980, 1995, 2008, and 2018, respectively; grassland constituted 35.31%, 36.23%, 35.27%, and 35.91% for the same periods. The average proportions of waters and wetlands in the HXNR were about 7.1% and 1.4%, while dispersed patches of forest had a total area of approximately 7.20 km 2 . For the control area, there were no cultivated lands similar to the HXNR, as shown in Figure 2. The area of grassland was the largest among the six LULC types, comprising around 90% of the total extent from 1980 to 2018 (Table 2). Desert in the control area was significantly smaller than in the HXNR but still was the second-largest LULC type, covering approximately 10% of the control area. And despite the control area being the capital of Zhiduo county where the population and industries are located, the urban area still accounts for less than 0.1%. During the period 1980-2018, the areas of grassland, wetland, and waterbodies in the HXNR increased significantly by 293.36 km 2 , 110.43 km 2 , and 37.42 km 2 , respectively, with the SLUDD of 0.05%, 0.03% and 0.44% (see Table 3). Meanwhile, the SLUDD of the desert in the HXNR is −0.05%, with a total reduction of 441.22 km 2 . In contrast, the grassland in the control area had been reduced considerably, which diminished by 99.35 km 2 (the SLUDD was −0.03%) from 1980 to 2018. The urban and desert areas showed relatively significant expansion, increasing by 5.36 km 2 and 94.39 km 2 and having the SLUDD of 0.24% and 19.08%, respectively.  For the control area, there were no cultivated lands similar to the HXNR, as shown in Figure 2. The area of grassland was the largest among the six LULC types, comprising around 90% of the total extent from 1980 to 2018 (Table 2). Desert in the control area was significantly smaller than in the HXNR but still was the second-largest LULC type, covering approximately 10% of the control area. And despite the control area being the capital of Zhiduo county where the population and industries are located, the urban area still accounts for less than 0.1%.
During the period 1980-2018, the areas of grassland, wetland, and waterbodies in the HXNR increased significantly by 293.36 km 2 , 110.43 km 2 , and 37.42 km 2 , respectively, with the SLUDD of 0.05%, 0.03% and 0.44% (see Table 3). Meanwhile, the SLUDD of the desert in the HXNR is −0.05%, with a total reduction of 441.22 km 2 . In contrast, the grassland in the control area had been reduced considerably, which diminished by 99.35 km 2 (the SLUDD was −0.03%) from 1980 to 2018. The urban and desert areas showed relatively significant expansion, increasing by 5.36 km 2 and 94.39 km 2 and having the SLUDD of 0.24% and 19.08%, respectively. For the period 1980 to 1995 (Figure 3a), the CLUDD of the HXNR was estimated to be 6.05%, compared to 8.43% in the control area. During this period, the conversion of desert to green spaces (i.e., forests, grasslands, and wetlands) was substantial, with a total area of 1508.27 km 2 (Table 4). However, newly desertified lands had an aggregated area of 1312.6 km 2 , mainly due to grassland degradation. Similarly, the area of desert greening in the control area was also pronounced. The desert area converted into grasslands, waterbodies, and wetlands reached 778.86 km2, while the area of newly desertified lands was relatively small, with a total area of 70.89 km 2 . Another noticeable change in the control area was the conversion of waterbodies and wetlands into grasslands with converted areas of 79.12 km 2 and 44.36 km 2 , respectively. At this stage, urban expansion was slow, expanding only by 0.96 km 2 .
Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 20   From 1995 to 2008, large areas of land in the HXNR were converted from and into deserts (Figure 3b). While 1326.43 km 2 of desert lands were reclaimed through desert greening, this area was exceeded by desertified land, which reached 1548.07 km 2 (Table 4). Previous conversions of waterbodies and wetlands into grasslands had been reversed to a certain extent at this period. The area of HXNR grasslands converted into waterbodies and wetlands reached 84.83 km 2 . As for the control area, desertified lands (781.54 km 2 ) were substantially more extensive than those that have been reclaimed through desert regreening (63.32 km 2 ). Similar to the HXNR, the grasslands converted into waterbodies and wetlands (44.98 km 2 ) were also significantly larger than the waterbodies and wetlands being converted into grasslands (5.11 km 2 ). The CLUDD of the HXNR and the control area were 6.14% and 8.39%, respectively.
From 2008 to 2018 (Figure 3c), protection policies and actions had been greatly improved, and human activities have been effectively monitored and regulated in the HXNR (e.g., poaching activities in the HXNR have been eliminated since 2008). The CLUDD in the HXNR was 3.55%, exceeding that in the control area (2.34%). The reclaimed area from desert greening (i.e., transformed into forests, grasslands, wetlands, Remote Sens. 2020, 12, 471 9 of 19 and waterbodies) totaled 839.86 km 2 , which is much larger than the desertified regions (368.16 km 2 ) ( Table 4). In the same period, the area of wetlands also increased substantially. The area converted between wetlands and waterbodies reached 400.55 km 2 , which includes water-to-wetland (271.81 km 2 ) and wetland-to-water (128.72 km 2 ). In contrast, the control area still reflected a distinct trend of desertification. The desert expanded by 159.84 km 2 , while desert greening was only able to reclaim 73.21 km 2 . Urban expansion was significant at this stage, with the SLUDD of 66.51% (Table 1) and an expansion area of 5.51 km 2 , which mainly encroached on grasslands, wetlands, and deserts.

ESV Change in Different Periods
From 1980 to 2018, the total ESV of the HXNR showed an overall growth trend, rising from 764.11 million USD in 1980 to 777.92 million USD in 2018, which is an increase of 13.81 million USD (Table 5). In contrast, the total ESV of the control area declined, with a total reduction of 5.33 million USD and an average annual decrease of 0.14% during the study period. The values of HXNR's provisioning, regulating, supporting, and cultural services increased by 0.59, 9.54, 2.48, and 1.19 million USD from 1980 to 2018, while the values for the control area decreased by 0.38, 2.77, 1.84, and 0.33 million USD respectively. Even though the ESVs in the HXNR have shown upward trends, their changes varied in the different stages (Table 5). From 1980 to 1995, the total ESV decreased annually by 0.01%, while the values of provisioning and supporting services increased by 0.15% and 0.13%, respectively. During 1995-2008, the total ESV also declined, but the yearly rate of decline was reduced to 0.004. The values of the provisioning and supporting services decreased by 0.18% and 0.16%, while the values of regulating and cultural services increased by 0.07% and 0.05%. In the period of 2008-2018, the ESVs of the HXNR generally expanded, with the total ESV increasing annually by 0.20%. The values of provisioning, regulating, supporting, and cultural services grew by 0.17%, 0.24%, 0.13%, and 0.19%.In contrast to the HXNR, the ESVs in the four ecosystem categories simultaneously increased or decreased for each period. From 1980 to 1995, the total ESV rose by 0.19%, with the provisioning, regulating, supporting, and cultural ecosystem services increasing in value by 0.34%, 0.11%, 0.30%, and 0.10%, respectively (Table 5). For 1995-2008, the total ESV showed an annual decline of 0.22%, with all the ecosystem service types decreasing in value. For 2008-2108, the downtrend trend in the total ESV continued having a yearly decline of 0.07%, and all ecosystem service types registering negative rates of change.

Spatial Change of ESV Per Unit
In both the HXNR and in the control area, the total ESV per unit have been spatially high in the south while low in the north (see Figure 4), which is similar to the spatial configuration of the average NDVI. The highest value per unit in the HXNR reached 2785 USD/ha between 1980 and 2018, while the control area had the highest value at 4797 USD/ha. High-value zones of the HXNR were mainly distributed in the southeast of the reserve, which were mostly situated around Trache Lake, Rigachi Mountain, and Dongduoqu River, while the rest were located along the route of Lake Trache-Dorge Co-Lake Hoh Xil. Due to high NDVI values, the adjusted ESV per unit in the control area was approximately four times higher, on average, than in the HXNR (Table 6). There were still low-value zones in the control area, mainly located along the 313 and 408 Provincial Highway. Table 6 shows the change in average ESV per unit from 1980 to 2018. Although the average value per unit slightly dropped for two periods (1980-1995 and 1995-2008     Spatially, the changes in the total ESV per unit at the HXNR for 1980-1995 ( Figure 5a) and 1995-2008 (Figure 5b) were mainly situated around the Kunlun Mountain Pass, south of the Rigachi Mountain, and the line from Dorge Co Lake to Chuma'er River, primarily due to the conversion of deserts and grasslands. In some of these places, value losses occurred in 1980-1995 but were regained in 1995-2008 (e.g., the areas around the Kunlun Mountain Pass and Dorge Co Lake), while other areas had the reverse trend. From 2008 to 2018 (Figure 5c), the areas that increased of the total ESV per unit in the HXNR were significantly more than those that declined, and are mainly concentrated around Xijir Ulan Lake, north of Lake Hoh Xil and north of Daxue Peak. The ESV increases in these areas can be attributed to the progress of desert greening and the transformation of lakes and wetlands. In contrast, areas with increased ESV could scarcely be found in the control area, while those with decreased ESV are situated around the towns of Jiajiboluoge and Lixin due primarily to urban growth and desert expansion.

Hotspots of ESV Change
Using the Hot Spot Analysis (Getis-Ord Gi*) tool in ArcGIS Pro 2.1, the statistically significant high values (hotspots) and low values (cold spots) of ESV change in the HXNR were identified. Table  7 summarizes the statistical results of hotspot analysis at 95% confidence level. From 1980 to 1995, the hotspots (49% of the total ESV changed patches) outnumbered the cold spots (44%), but there were more cold spots (31%) in the core zone than hotspots (9%), which indicates that at this stage, the decline of ESV in the core zone had played a crucial role in reducing the total ESV in HXNR. From 1995 to 2008, while cold spots (50%) exceeded the number of hotspots (44%) in total, there were more

Hotspots of ESV Change
Using the Hot Spot Analysis (Getis-Ord Gi*) tool in ArcGIS Pro 2.1, the statistically significant high values (hotspots) and low values (cold spots) of ESV change in the HXNR were identified. Table 7 summarizes the statistical results of hotspot analysis at 95% confidence level. From 1980 to 1995, the hotspots (49% of the total ESV changed patches) outnumbered the cold spots (44%), but there were more cold spots (31%) in the core zone than hotspots (9%), which indicates that at this stage, the decline of ESV in the core zone had played a crucial role in reducing the total ESV in HXNR. From 1995 to 2008, while cold spots (50%) exceeded the number of hotspots (44%) in total, there were more hotspots found in the core area (31%) than cold spots (10%). Since the total ESV was still in a downward trend at this period, this means that the increase in ESV in the core zone was not enough to offset losses in the buffer zone. From 2008 to 2018, the proportion of hotspots was significantly higher than that of cold spots, accounting for 52% and 28% of the total changed patches, respectively, prompting a notable increase in HXNR's aggregated ESV. For the entire period of 1980-2018, the hotspots in the HXNR (51%) were far greater than cold spots(30%). Hotspots (28%) in the core zone outnumbered those in the buffer zone (23%), while cold spots in the core zone (11%) were less than those found in the buffer zone (19%). This suggests that the HXNR's core zone had been more effective in improving ecosystem services. In terms of spatial configuration (Figure 6), for 1980-1995, the hotspots of ESV change were mainly concentrated in the south of Fengxue Beach (located in the southwest corner of the core zone of HXNR), the south of Trache Lake-Rigachi Mountain, and the north of Chuma'er River in the buffer zone. Around Dorge Co Lake, there were also some scattered hot spots. Meanwhile, cold spots were mainly located around Cuodarima Lake-Dorge Co Lake, Xijir Ulan Lake, and Lianhu Lake in the core zone, as well as in the Chuma'er River section and Lubei Beach. From 1995 to 2008, the hotspot and cold spot regions were noticeably in reverse of the preceding period, and both distributions expanded significantly, either in the core area or in the buffer zone. For the succeeding period of 2008-2018, the distribution of hotspots and cold spots was more widespread and uniform. The regions with relatively concentrated hotspots were the south of Trache Lake-Rigachi Mountain and the areas surrounding Hoh Sai Lake, Xijir Ulan Lake, Lexie Wudan Lake, Lianhu Lake, and Moon Lake. The cold spots were more concentrated in the areas around Haidingnuo'er Lake and Yanhu lake and the line of Daxue Peak-Buka Daban Mountain Peak. significantly, either in the core area or in the buffer zone. For the succeeding period of 2008-2018, the distribution of hotspots and cold spots was more widespread and uniform. The regions with relatively concentrated hotspots were the south of Trache Lake-Rigachi Mountain and the areas surrounding Hoh Sai Lake, Xijir Ulan Lake, Lexie Wudan Lake, Lianhu Lake, and Moon Lake. The cold spots were more concentrated in the areas around Haidingnuo'er Lake and Yanhu lake and the line of Daxue Peak-Buka Daban Mountain Peak.

Characteristics of ESV Change Due to LULC Dynamics.
Based on the LULC data sets, we used NDVI to enhance the equivalent ESV coefficients proposed by Xie et al. [28] in order to evaluate the ESVs from 1980 to 2018 in the study area. In the HXNR, although declining slightly in the first two periods (1980-1995 and 1995-2008), the total ESV exhibited a general upward trend, increasing from 764.11 million USD in 1980 to 777.92 million USD in 2018. Previous studies have explored the ESV change of the Sanjiangyuan area (including parts of the HXNR) and found that the ESV decreased initially and then increased considerably [67]. Our findings support this conclusion that there is a downward trend from the 1980s to the mid-1990s.
Spatially, the ESV change in HXNR can be found in the Kunlun Mountain Pass, south of the Rigachi Mountain, and the line from Dorge Co Lake to Chuma'er River. The areas with increased ESV can be found in the north of Hoh Xil Lake, north of Daxue Peak, and the surrounding regions of Xijir Ulan Lake, which could be attributed with the conversion of desert to grasslands in north of Hoh Xil Lake and north of Daxue Peak and in the expansion of lakes and wetlands around the Hoh Xil Lake and Xijir Ulan Lake. Yao et al. [50] have reported the continuous expansion of Hoh Xil Lake and Xijir Ulan Lake from 1970 to 2010 due to the increase in precipitation and temperature, which is

Characteristics of ESV Change Due to LULC Dynamics
Based on the LULC data sets, we used NDVI to enhance the equivalent ESV coefficients proposed by Xie et al. [28] in order to evaluate the ESVs from 1980 to 2018 in the study area. In the HXNR, although declining slightly in the first two periods (1980-1995 and 1995-2008), the total ESV exhibited a general upward trend, increasing from 764.11 million USD in 1980 to 777.92 million USD in 2018. Previous studies have explored the ESV change of the Sanjiangyuan area (including parts of the HXNR) and found that the ESV decreased initially and then increased considerably [67]. Our findings support this conclusion that there is a downward trend from the 1980s to the mid-1990s.
Spatially, the ESV change in HXNR can be found in the Kunlun Mountain Pass, south of the Rigachi Mountain, and the line from Dorge Co Lake to Chuma'er River. The areas with increased ESV can be found in the north of Hoh Xil Lake, north of Daxue Peak, and the surrounding regions of Xijir Ulan Lake, which could be attributed with the conversion of desert to grasslands in north of Hoh Xil Lake and north of Daxue Peak and in the expansion of lakes and wetlands around the Hoh Xil Lake and Xijir Ulan Lake. Yao et al. [50] have reported the continuous expansion of Hoh Xil Lake and Xijir Ulan Lake from 1970 to 2010 due to the increase in precipitation and temperature, which is consistent with the findings of this study. Moreover, given the almost opposing change trends in ESVs before and after the establishment of PA, the restoration and maintenance of ecosystem services in the HXNR have clearly benefited with the establishment of the PA after 1995.

Effectiveness of the HXNR on Ecosystem Services Protection
In this study, we chose a non-protected region as a control area, situated in the southeast of the HXNR. By comparing the ESV changes between PA and non-PA, We determined the effectiveness of the HXNR in promoting ecosystem services. Compared with the HXNR, the ESV in non-PA demonstrated a contrary trend change, which increased by 2.85 million USD in 1980-1995, decreased sharply in 1995-2008 (decreased by 2.87 million USD), and continued dropping by 0.66 million USD in 2008-2018, mainly due to the impact of urban expansion and desertification. Without protection, the non-PA suffered more ESV losses, even with an initial upward trend. Additionally, the results of hotspot analysis showed that most of the sites with increased ESV belong to the HXNR's core zone, significantly contributing to the growth of the total ESV. Given the characteristics of ESV change through the different phases in establishing the PA and the contrasting developments between PA and non-PA, the establishment of the HXNR has clearly mitigated the degradation of the ecosystem and has been effective in enhancing ecosystem services. While there has been limited research evaluating the conservation effectiveness in the HXNR, several studies have observed improvements in habitat suitability and biological population [46][47][48], and the effective regulation of human activities in the HXNR after the 1990s [61], which are consistent with the conclusions of this study.
Even though our findings suggest that the HXNR was effective in improving ecosystem services, it is still less clear whether all PAs can successfully contribute to the enhancement of ecosystem services. Palomo et al. [68] demonstrated that high levels of PAs, such as the National Parks and Natural Reserves, did not ensure the provision of multiple ecosystem services due to the disconnect between protected areas and society. Castro et al. [36] suggested that the current PA networks supply only slightly higher levels of regulating services than non-PAs at the landscape level. On the other hand, based on nine individual cases, Eastwood et al. [5] assert that PAs deliver higher levels of ecosystem services than non-PAs. Given the variations in research scales, methods, and settings, there is still considerable uncertainty whether PAs are effective in promoting ecosystem services. The differences in outcomes would be largely dependent on geographical conditions, targets, and policies of PAs [68]. Thus, to ensure the effectiveness of PAs in promoting ecosystem services, ecosystem assessments and environmental evaluations should integrate more the use ecosystem services, particularly at the local scale.

Application of ESV on the Effectiveness Assessment of PAs
Based on LULC data from remotely sensed imagery, this study evaluated the ESV using an NDVI-based modified version of the benefit transfer method. The ESV evaluation was integrated into the effectiveness assessment of PAs by comparing the ESV changes at different phases between PA and non-PA. This approach exhibited feasibility and practicality in various aspects that could be implemented in other PAs. First, with the use of remotely sensed imagery, which has wide observation range, fast acquisition speed, and short update period, a quick evaluation of ESV can be implemented on a large-scale, long-term basis, and at low costs [69], particularly in remote and data-scarce areas.
Second, this approach assesses the various types of ecosystem services in monetary units using benefit transfer method, which demonstrates flexibility in analyzing the levels of different ecosystem services and practicability in comparing economic benefits of PAs. Such an approach that analyzes benefits in pecuniary terms can be integrated into evaluation studies of PAs to provide a more effective platform that would be easily understood by decision-makers and other stakeholders [9,70]. For example, based on the ESV change in monetary units, it is workable for the government and environmental protection organizations to evaluate the income on the investment of ecosystem protection projects, and also provides a basis for the improvement of ecological compensation standards for local residents. Since the benefit transfer method proposed by Costanza et al. [15] and Xie et al. [28] overlooked spatial heterogeneity [52,57], our revised approach uses NDVI to modify equivalent coefficients in order to overcome this limitation and lessen deviations inherent in other techniques (e.g., modification based on NPP) [52,55,56]. In vegetation monitoring, the NDVI varied in space with the meteorological condition and local resource background [66], especially in the regions with numerous different vegetation covers [71]. Thus, NDVI is adequate to represent the spatial heterogeneity of ecosystems and to adjust the benefit transfer method.
Additionally, the ESV evaluation based on grid cells and the hotspot analysis of ESV change provide more spatial details of ESV and allows the spatial evaluation of ecosystem service improvements in PAs. This offers additional value in the decision-making of PAs' protection and governance [53] since strategies and policies promoting ecosystem services need to be adjusted in accordance with the varying spatial conditions.

Limitations and Future Directions
This study reclassified land use/land cover into six primary categories based on remotely sensed LULC data and evaluated the ESV of the HXNR by using an NDVI-revised benefit transfer method. This approach is preferred when quick ESV evaluation is required but the available data are limited. However, one major limitation in this approach is that the differences in areas having the same LULC type were largely ignored, which may influence the accuracy of the evaluation results. In fact, the LULC primary type can be further subdivided into multiple subclasses (e.g., dry lands and paddy fields in croplands, broad-leaved woodlands, and shrublands in the forest), such that the ecosystem services provided by each LULC subclass could significantly vary [16,55]. Considering the diversity of bio-species and landscape in PAs, more detailed research on equivalent coefficients in benefit transfer models is essential, particularly those focusing on variations in LULC subclasses. Moreover, this method is also limited by ground truthing as the study focused on a long-term dynamic process covering a large area. Hence, incorporating reliable field verification into the ESV assessment needs to be fully considered in future studies.
Considering the strong correlation with meteorological condition and local resource background [66], the NDVI index was used to modify the benefit transfer method, which make up for the shortage of ignoring the spatial heterogeneity of ESV. However, the uncertainty of the NDVI index was overlooked. For instance, the changeable relationship between NDVI and temperature in different regions [71] may affect the accuracy of ESV evaluation. Accordingly, the NDVI-based localized benefit transfer method needs to integrate multiple-indicates in future ESV evaluation studies. Previous studies developed various satellite-based vegetation indices applied in monitoring the vegetation [69,71]. Such indices (e.g., the Vegetation Condition Index (VCI), the Vegetation Health Index (VHI), the Perpendicular Vegetation Index (PVI), and the two-band version of the enhanced Vegetation Index (EVI2) et al.) representing a range of spectral responses of vegetation conditions [69], would be appropriate to adjust the value transfer method. However, each vegetation index varies in performance due to environmental changes (e.g., weather) [71]. Therefore, there is still a call for studies on the application scenarios of different vegetation indices in modifying value transfer method.
In addition, by comparing the ESV change in PA and non-PA at different periods, our findings demonstrated the effectiveness of PA in improving the different types of ecosystem services. Although the spatial effectiveness was evaluated using hotspots analysis, it remains incomplete as spatial conflicts and trade-offs within different stakeholders on ecosystem services were not considered. Trade-offs of ecosystem services, arising from human preferences for specific sets of ecosystem services (e.g., human preferences often prioritize provisioning services over regulating services), have shown significant influence on policies and strategies in environmental management [70,72]. Thus, future research focusing the effectiveness of PAs would need to consider ecosystem service trade-offs, in order to provide a comprehensive understanding of ecosystem services and to ensure that the target of ecosystem conversation is not dominated by short-term needs [70], which could be more valuable in policy-making and governance.

Conclusions
We analyzed the dynamics of LULC in the HXNR and the control area using remote sensing techniques and evaluated the ESV change due to LULC using the NDVI-based modified value transfer method. Comparing ESV changes between PA and non-PA at different periods, we confirmed the effectiveness of the HXNR in improving ecosystem services and identified the spatial effectiveness at the different districts. The main findings of this study are as follows: (1) The ESV of the HXNR initially decreased in value and then increased considerably over time, with the total ESV of 764.11 million USD in 1980 and increasing to 777.92 million USD in 2018. The four categories of ecosystem services (provisioning, regulating, supporting, and cultural) all increased in value, largely due to the conversions of desert-to-green and waterbodies-to-wetland. In contrast, the ESVs of the control area improved in 1980-1995 and declined in the succeeding years, with the total ESV decreasing by 5.33 million USD mainly due to urban expansion and desertification.
(2) The change of ESV per unit in the HXNR were concentrated in regions around the Kunlun Mountain Pass, Xijir Ulan Lake, north of Hoh Xil Lake, south of the Rigachi Mountain, north of Daxue Peak, and the line from Dorge Co Lake to Chuma'er River. Generally, the areas that achieved ESV growth significantly outnumbered areas with declining values. As for the control area, the regions that saw an increase in ESV in 1980-1995 experienced a decline in values in 1995-2008; also, the areas around the towns of Jiajiboluoge and Lixin experienced a value decline in 2008-2018.
(3) Comparing the ESV changes between the HXNR and the control area at different periods, our results show the HXNR has been more effective in improving the value of ecosystem services than the control area.
(4) The core zone of the nature reserve achieved a higher level of effectiveness in improving ecosystem services, given that the hotspots of ESV change were more prominent in the core zone, while the cold spots were much concentrated in the buffer zone.
The methodology used in this study can be widely applied for prompt and effective evaluation of PA in enhancing ecosystem services, especially for long-term conservation assessment in areas with limited data. In order to comprehensively understand the effectiveness of PAs on ecosystem services conservation, more research and case studies from around the world have to be undertaken, particularly at the local scale.