Effects of Land Use/Cover Changes and Urban Forest Configuration on Urban Heat Islands in a Loess Hilly Region: Case Study Based on Yan’an City, China

In this study Yan’an City, a typical hilly valley city, was considered as the study area in order to explain the relationships between the surface urban heat island (SUHI) and land use/land cover (LULC) types, the landscape pattern metrics of LULC types and land surface temperature (LST) and remote sensing indexes were retrieved from Landsat data during 1990–2015, and to find factors contributed to the green space cool island intensity (GSCI) through field measurements of 34 green spaces. The results showed that during 1990–2015, because of local anthropogenic activities, SUHI was mainly located in lower vegetation cover areas. There was a significant suburban-urban gradient in the average LST, as well as its heterogeneity and fluctuations. Six landscape metrics comprising the fractal dimension index, percentage of landscape, aggregation index, division index, Shannon’s diversity index, and expansion intensity of the classified LST spatiotemporal changes were paralleled to LULC changes, especially for construction land, during the past 25 years. In the urban area, an index-based built-up index was the key positive factor for explaining LST increases, whereas the normalized difference vegetation index and modified normalized difference water index were crucial factors for explaining LST decreases during the study periods. In terms of the heat mitigation performance of green spaces, mixed forest was better than pure forest, and the urban forest configuration had positive effects on GSCI. The results of this study provide insights into the importance of species choice and the spatial design of green spaces for cooling the environment.


Introduction
It is well known that urbanization is one of the most powerful and visible anthropogenic forces on Earth [1][2][3]. The most obvious aspect of urbanization is that the natural landscape is increasingly being replaced by anthropogenic land use/land cover (LULC) types, which can lead to many ecological and environmental problems, such as urban heat islands [4,5]. In China, cities are expanding rapidly as the economy grows, but land suitable for development is in short supply, particularly in mountainous areas, where about one-fifth of the population lives. Thus, in the last decade, local governments have begun removing the tops of mountains to fill valleys and create land suitable for building. In cities such as Chongqing, Shiyan, Yichang, Lanzhou and Yan'an, tens of square kilometers of land have been created. However, in the Yan'an city area, the relatively warmer urban land surface/atmosphere compared with the rural surroundings is very obvious in the summer, which is usually referred to as an urban heat island (UHI).

Related Research
UHIs were identified for the first time in London during 1833 [6]. An UHI is more obvious under calm, cloudless conditions, where it depends greatly on the urban morphology and land use type [7][8][9]. UHIs are now considered major environmental problems in the 21st century [7,10,11], because the higher temperatures in UHIs can lead to increases in urban energy consumption [12][13][14][15], raised pollution levels [16][17][18], and they may even affect the habitability of cities and increase mortality [12,15,19]. Therefore, methods for mitigating UHIs are now a major research focus. UHIs exhibit different forms associated with at least four spatial regions: boundary-layer UHI (BLUHI), canopy layer UHI (CLUHI), surface UHI (SUHI), and subsurface UHI (SubUHI). The urban canopy layer extends upward from the surface to approximately the mean building height, whereas the urban boundary layer is located above the canopy layer [20]. CLUHI and BLUHI are atmospheric heat islands caused by warming of the urban atmosphere, whereas SUHI refers to the relatively higher warmth of urban surfaces compared with surrounding rural areas. SubUHI is part of the overall UHI, which denotes the relative warmth of urban ground temperatures compared with the rural background [21]. It is known that atmospheric UHIs are larger at night whereas surface UHIs are larger during the day [22]. SUHIs have stronger effects on the daily outdoor activities of people. Moreover, SUHIs can be measured conveniently at different spatial-temporal scales and simulated based on remote sensing data. The effects of SUHIs can be evaluated using air temperature measurements and satellite land surface temperature (LST) measurements. The air temperature measurements are obtained either by traversing a city or by comparing point temperature measurements [23,24]. In general, SUHI measurements based on air temperature have high temporal resolution with extensive time coverage, and they can effectively describe the temporal variation in the effects of UHIs [25]. However, spatially continuous analysis is often difficult because of the sparse distribution of observation stations. Fortunately, the LST is a universal and important parameter for analyzing SUHIs. To address these problems, many studies of the effects of SUHIs have considered LST data, which were mainly measured using two approaches. Traditionally, LST data are measured via ground-based observations obtained from automobile transects and weather station networks [20]. However, at present, due to the development of thermal remote sensors, satellite-based imaging technology is now employed widely to detect UHIs remotely and regionally because it facilitates the straightforward and consistent determination of the spatial-temporal LST distribution [26]. Thus, AVHRR and ATSR [27,28], Moderate Resolution Imaging Spectoradiometer (MODIS) [29], Landsat Thematic Mapper (TM) [30,31], Landsat Enhanced Thematic Mapper Plus (ETM+) [32], Landsat Operational Land Imager and Thermal Infrared Sensor (OLI & TIRS) [33,34], Chinese HJ-1B Infrared Multispectral Scanner (IRMSS) [35][36][37] alone or combined with high-resolution satellite images, such as SPOT [38], Gaofen-1 [29], or Quickbird [25,[39][40][41], can be used for the rapid retrieval of LST data. In addition, hyperspectral thermal infrared data, such as those acquired by the Infrared Atmospheric Sounding Interferometer (IASI) and Cross-track Infrared Sounder (CrIS), have been used for retrieving LST data [42]. Therefore, LST data derived from thermal infrared remote sensors are among the most commonly used indicators for heat island analysis. In the last two decades, many studies have focused on identifying the factors that contribute to SUHIs [7,10,24,[43][44][45][46][47], which have shown that urban built-up land and the impervious surface area have positive effects on SUHIs [48], whereas water bodies and green landscape have negative effects on SUHIs [24,43,49,50]. In addition, previous studies have demonstrated that there are complex relationship between the landscape composition [25,39,44], land use and cover changes (LUCC) [3,36,51,52], and SUHIs.
Urban green spaces (UGS), including urban parks, green belts, attached green space, and productive plantation areas, are considered to be important components of urban vegetation, where they are cooler than their surrounding built-up areas and they can form a "green space cool island" (GSCI).
A previous study suggested that the GSCI intensity is affected by the size, shape, seasonal changes (summer/autumn), and the forest structure (stem density, diameter, tree height, basal area, leaf area index (LAI), and canopy density) for urban parks in Shenyang, China [24]. The urban landscape configuration (LULC types and landscape pattern metrics) also influences SUHIs [43,53]. However, the relationships between the GSCI intensity and other aspects of urban forest structures (e.g., species composition, three-dimensional green biomass (TGB), and health) in UGSs have rarely been studied and they are not yet fully understood. Furthermore, previous studies of the spatiotemporal distribution and factors that affect SUHIs have mainly considered municipalities, provincial capitals, or coastal cities, whereas few studies have investigated cities at the prefecture level in the loess hilly region of China.

Study Objectives
In this study, based on field surveys and Landsat series (5/7/8) satellite remote sensing data (1990-2015, 5-year intervals) acquired for Yan'an City, China, we investigated the effects of the LUCC and UGS tree-layer structure on the SUHI intensity. The aims of this study were: (1) to explore the quantitative relationships between the GSCI intensity and UGS structure, and to determine whether the UGS structure significantly affects the GSCI intensity; (2) to examine the spatiotemporal trends in the LST and LUCC during the last 25 years (1990-2015); and (3) to identify the main LULC types that significantly affected the LST during six research periods. The results of this study should be useful for urban planners and designers by facilitating the design of UGSs that maximize the GSCI intensity and mitigate UHIs.

Study Area
The case study area, in the Baota District of Yan'an City (36 • 22 44"-36 • 45 53" N, 109 • 14 7"-109 • 46 9" E) in the middle reaches of the Yellow River, where it is located in an arid and semiarid region of the Loess Plateau in Northwestern China, has been affected by severe soil erosion which makes this district a key region for ecological restoration. The case study area covers 1174.54 km 2 (comprising 1600 × 1432 grid cells at 30-m resolution) and the elevation ranges from circa 900 m to circa 1466 m above sea level ( Figure 1). In 2013, this district had a population of 478,500 (Shaanxi Statistical Yearbook in 2014), 63% of whom were rural. This region has a typical semiarid continental climate, with average annual rainfall of approximately 470 mm (with high variability in recent years), where over 65% occurs between June and September, mainly in the form of heavy rain (Baota Meteorological Observatory). The study area contains 260 villages according to the Second National Land Survey of China. This area has undergone long-term soil and water conservation, and it was selected as one of the pioneer demonstration areas for the large-scale ecological restoration project known in China as "Grain to Green" [4,54,55].
Petroleum exploitation has been conducted throughout the study area by the Shaanxi Yanchang Petroleum (Group) Co. Ltd. (Yan'an, China). During the past three decades, and road construction on forestland or farmland in order to transport the crude oil from exploitation sites to oil refineries has destroyed large amounts of vegetation.
One of the largest construction projects in the study area, which started in April 2012 (Figure 1), doubled the city's current area by creating 78.5 km 2 of flat ground [56], on which the new Yan'an City was built. Most of the vegetation was removed from more than 30 hills.

Data Sources
In this study, three types of data were utilized: (1) Baota district meteorological data (Baota district bureau of meteorology), (2) Baota district land use data (provided by Baota District Branch of Yan'an Municipal Bureau of Land and Resources), and (3) Landsat satellite images (path/row: 127/35), downloaded from the United States Geological Survey [58]. Further details of these data are given in Table 1.

Technical Details
The present study of the changes in the fraction of suburban vegetation and urban green spaces (UGSs), and their effects on the thermal environment considered four main issues ( Figure 3): (1) the fraction of suburban vegetation and the characteristic changes in UGSs were studied using spatial analysis and statistical analysis based on land use and cover change (LUCC) and associated remote sensing indicators; (2) the characteristic changes in the thermal environmental were examined using the land surface temperature (LST) single-channel algorithm and the landscape assessment method; (3) the effects of LUCC on the thermal environment were analyzed using a spatial linear regression method for the LST data versus land surface remote sensing indicators; and (4) the mitigating effects of UGSs on surface urban heat islands (SUHIs) were assessed based on urban forest configuration-related indicators.

Technical Details
The present study of the changes in the fraction of suburban vegetation and urban green spaces (UGSs), and their effects on the thermal environment considered four main issues ( Figure 3): (1) the fraction of suburban vegetation and the characteristic changes in UGSs were studied using spatial analysis and statistical analysis based on land use and cover change (LUCC) and associated remote sensing indicators; (2) the characteristic changes in the thermal environmental were examined using the land surface temperature (LST) single-channel algorithm and the landscape assessment method; (3) the effects of LUCC on the thermal environment were analyzed using a spatial linear regression method for the LST data versus land surface remote sensing indicators; and (4) the mitigating effects of UGSs on surface urban heat islands (SUHIs) were assessed based on urban forest configurationrelated indicators.

Derivation of the Normalized Difference Vegetation Index (NDVI), Index-Based Built-Up Index (IBI), and Modified Normalized Difference Water Index (MNDWI), and LULC Classification
Six remote sensing images acquired in the summer (Table 1, Serial Number = 4, 5, 6) were used to produce multi-temporal sets of LULC maps of the study area for 1990, 1995, 2000, 2005, 2010 and 2015, with the assistance of the auxiliary data (Table 1, Serial Number = 1, 2, 3) and using maximum likelihood classification in ENVI TM 5.1 (Exelis Visual Information Solutions, Inc., Boulder, CO, USA) [59]. The detailed image processing techniques, including radiometric calibration, geometric correction, image classification, and accuracy assessment, were described previously [60]. The "Present Situation of Land Use Classification (GB/T221010-2007)" and the actual details of the study area were employed to extract the characteristics of the remote sensing data. Images of the study area were used to classify the LULC into six types: construction land, farmland, forest, grassland, water, and unused land ( Table 2). Manual correction was then applied to ensure the accuracy of the classification. The accuracies of the classified products were assessed by manual interpretation using Google Earth Pro ® . In total, 5% of the patches measuring more than 15 ha were selected randomly as samples. Using the spatial join function in ESRI ArcGIS TM version 10.0 (ESRI, Redlands, CA, USA) [61], the manual interpretation results and the original results for these selected samples were compared to produce confusion matrixes [62]. The results of the Jeffries-Matusita distance separability inspections were all above 1.8000. The accuracies of the maps were above 90% in all 6 years. Finally, the kappa indices were calculated and the results were above 0.85 in all 6 years. These values satisfied the accuracy requirements for land-use change analysis [62]. All of the spatial data were transformed to a uniform coordinate system (datum: Beijing_54, ellipse: Krasovsky, projection: Transverse Mercator, zone: 19N). The NDVI is a proxy of vegetation cover, which is frequently used in ecological and environmental studies: The MNDWI can enhance open water features while efficiently suppressing and even removing built-up land noise as well as vegetation and soil noise, and thus MNDWI is suitable for enhancing and extracting water information for a water region with a background dominated by built-up land areas [63]. We used MNDWI to represent water areas.
IBI was proposed for the rapid extraction of construction land features in satellite imagery [64]: where ρ Green , ρ Red , ρ N IR , and ρ MIR represent the surface reflectivity of the green band, red band, near-infrared (NIR) band and, mid-infrared band, respectively, of the Landsat satellites images after Fast Line-of-sight Atmospheric Analysis of Hypercubes (FLAASH) atmospheric correction with ENVI 5.1.

LST Inversion
In this study, the LST was inverted from Landsat thermal infrared data (Table 1, Serial Number = 4, 5, 6, resampled to a spatial resolution of 30 m) in terms of TM-6 or ETM+-6 using the single-channel algorithm [65], and for TIRS-10/11 using the generalized single-channel method [66] combined with the split-window covariance-variance ratio technique (SWCVR) [67]. The main steps in this process comprised: (1) correcting the radiometric and geometrical distortions; (2) converting calibrated digital numbers (DNs) into absolute units of at-sensor spectral radiance; (3) converting the at-sensor spectral radiance into the at-sensor brightness temperature; and (4) correcting for the spectral emissivity of different land cover types (NDVI values) to generate the LST data [34,[65][66][67]. More details of the computational process as follows: Step 1: At-sensor radiance (L sen ).
Gain and offset are the slope and intercept of the radiance/digital number (DN) conversion function, respectively, which were obtained from the metadata files.
where the two correction constants K 1 and K 2 , whose unite is in W/(m 2 ·sr·µm) and K, respectively. More details value were shown in Table S1 (Supplementary Materials).
(2) Atmospheric parameters (ψ 1 , ψ 2 , ψ 3 and w) The values of constants (a 1 , a 2 , a 3 , b 1 , b 2 , b 3 , c 1 , c 2 and c 3 ) were shown in Table S1. The w is water vapor content in g/cm 2 , which can be obtained from atmospheric profile (TM-6 and ETM+-6) or SWCVR method (TIRS-10): where FV is the vegetation fraction, NDVI S and NDVI V correspond to the bare soil and fully-vegetated NDVI, respectively, which can be extracted from the NDVI histogram according to the cumulative percentage (5% sand 95% respectively) in the corresponding NDVI data. Refer to TM-6 and ETM+-6: τ is transmissivity; L up and L down are up-welling and down-welling atmospheric radiance, respectively, which were calculated using a web-based tool (http://atmcorr.gsfc.nasa.gov).
Refer to the band 10 of Landsat 8: where T 10,k represents the brightness temperature of TIRS band 10 for pixel k; T 11,k represents the brightness temperature of TIRS band 11 for pixel k; T 10 and T 11 represents the mean brightness temperature of TIRS band 10 and 11 over N pixels, respectively.
Measurement of the Relative SUHIs

• Green Space Cooling Island Intensity
The GSCI intensity is calculated as: where T gs is the daily average LST m for a certain UGS interior and T u is the daily average land surface temperature measured by temperature and humidity probes (LST m ) in the external 10 m buffer (excluding other green spaces and water) of the corresponding UGS, which is sufficiently wide to include the neighboring urban thermal information for roads, residential or business buildings, and parking spaces. In this context, the unit for temperature is degrees Celsius ( • C). The thermal field variance index (LST grade/heat island intensity) was graded into five levels ( Table 3). The spatial distribution of heat islands in the study area is shown in Figure 3. Table 3. Thermal landscape classification obtained using the mean-standard deviation method.

Thermal Landscape Category (LST Grade/Heat Island Intensity) LST Division
Hot/extremely strong y) is the surface temperature at location (x, y), and m and std are the mean and standard deviation for the LSTs, respectively.
• Thermal Landscape The thermal landscape in the urban and suburban areas was divided into five categories using the LST mean-standard deviation method [68,69] (Table 3).

Landscape Pattern Analysis
To assess the changes in the structural characteristics of the land and thermal landscapes at a scale of 100 m from 1990 to 2015 (5-year interval), the FRAGSTATS 4.2 program (University of Massachusetts, Boston, MA, USA) [70] was used to calculate the following six landscape metrics: fractal dimension index (FRAC), percentage of landscape (PLAND), aggregation index (AI), division index (DI), Shannon's diversity index (SHDI), and expansion intensity (EI) ( Table 4). These metrics have been used frequently to assess the structural characteristics of landscapes and to monitor changes in land use [44,53,[71][72][73]. Table 4. Evaluation indices for LST and LUCC types [72,73].

Evaluation Index Description Formula
Fractal dimension index (FRAC) Ranging between 1 and 2, where a greater value indicates more complex characteristic of the plaque and landscape.
A is the total area and P is the perimeter of a patch.
Characteristic of a certain class area relative to the proportion of the total. A is the total area, a is the plaque area, and n is the number of patches.
Characterization of the degree of plaque accumulation, ranging between 0 and 100, where a lower value indicates a greater degree of dispersion for the representative. g ii is adjacent to a number of patches relative to a class plaque.
Measure of the plaque distribution, ranging between 0 and 1, where a value closer to 1 represents a more severe split. A is the total area, a i is the area of the ith plaque, and n is the number of patches.
Shannon's diversity index (SHDI) Diversity measure that increases with the number of patch types and as the proportional distribution of the area among patch types becomes more equal.

Expansion intensity (EI)
Measure of the intensity of spatial expansion. A i + j and A i are the areas in years i + j and i, respectively.
FRAC, PLAND, AI, DI and EI were calculated for LULC and SHUI in the whole area, suburban area and urban area at the landscape level [74].

Surveying and Measurement of SUHI-Related Indicators at the Plot Level
Size and Shape of UGSs We used the Google Earth Professional online platform (Version 7.1) and images captured on 25 July 2015 under the same conditions (perspective elevation of ca. 1500 m and a parallel projection state). We delineated boundary polygons of 34 representative green spaces (area larger than 1 ha; Figure 2) with relatively even distributions in the urban core area, which ranged in size among large, medium, and small [75,76]. Each boundary polygon was saved separately in the form of a KML file, which was then transformed into a shape file (shp) in Global Mapper V10.01 and employed for extracting area, perimeter, external buffer, and subset/clipping-related spatial data. The UGS area, perimeter, and landscape shape index (LSI) were used to describe the UGS size and shape, and the perimeter/area ratio described the complexity and the edge effect for an urban park. A larger LSI indicated a more complex urban park shape [24]. LSI was calculated using ArcGIS 10.0 as follows: where P t is the total perimeter around a green space and A is the area of the green space.

Surveys of UGS Forest Structure and Temperature Measurements
We determined the urban forest structures based on the tree-layer species number (TSN), TGB and LAI for forests in the UGSs according to our preliminary investigations. The number of plots sampled for each green space is shown in Figure 2. Each of the 34 sampling quadrats was defined as a 20 m × 20 m square (0.04 ha).
TSN: Number of species (diameter at breast height ≥5 cm, with a certain number) comprising the tree layer of a forest in an UGS. LAI: One half the total leaf surface area per unit ground area [77]. Every tree in each sampling quadrat was measured directly with a canopy analyzer (fish-eye camera) and adjustable tripod. Each single tree's LAI was calculated using analysis software (WinSCANOPY Reg v2005a) and the mean value was taken as the LAI for the whole sampling quadrat.
TGB: The space volume occupied by living plant stems and leaves [78]. TGB was evaluated for each tree in a sampling quadrat using an empirical formula according to corresponding shape of the tree canopy (Table 5).  The near-ground (approximately 1.5 m above the ground) air temperature (AT 1.5 ) and land surface temperature (LST m ) were measured using temperature and humidity probes (Yangling Qiantai Electronic Science and Technology Co. Ltd., Yangling, China) outside (30 m from the center, without trees) and inside 34 different types of green space during 26-31 July 2015, from 08:00 to 18:00, under roughly the same weather conditions (sunny, cloudless and wind-free). The data recording interval was 20 min and the daily average values were used in the statistical analyses.

Statistical Analysis
Ordinary least squares multiple linear regression models and correlation coefficients were used to determine the effects of the UGSs configurations on LST m at a scale of 100 m. In the analysis of the relationships between the linear regression models, GSCI/LST was used as the dependent variable, and the amount of urban vegetation (LAI and TGB) and shape (LSI) of the UGSs were used as independent variables. Spatial linear regression analysis of LST versus NDVI, FV, IBI and MNDWI in urban and suburban areas were performed by using SAS 9.2 (SAS Institute Inc., Cary, NC, USA) in order to analyze their quantitative relationships during the six study periods based on over 500 random sampling points with distance ≥200 m, which were evenly distributed in the study area.

Characteristics of the Mean Annual and Monthly Air Temperature, and Summer Heat Islands
June, July and August were the warmest months of the year (Figure 4) and we used Landsat satellite images, including the thermal infrared bands, for LST retrieval. According to the meteorological station at Yan'an (36 •  was attributable to the strong global EI Nino climatic effect as well as severe local deforestation and construction activities due to construction of the new Yan'an city area, which caused dry and warm weather. These effects are readily discernible in the annual mean air temperature charts shown in Figure S1.  Figure S2 showed the 20-min variations in the SUHI intensity at 34 open impervious surface sites located close to the centre of Yan'an City in July. The SUHI intensity values ranged from −6.9 • C to 10.1 • C, and the standard deviation ranged from −2.7 • C to 9.4 • C during the day (08:00-18:00), where the SUHI occurred predominantly during the nearly noon-time hours and it reached a maximum by 11:40 (Beijing Time).

Variations in the LST among Different Land Use Types
The average LST distributions during the six research periods for urban area of Yan'an city (UA), and suburban area of Yan'an City (SUA) are showed in Figure 4. In the two regions (UA and SUA), the lowest average LST always occurred in forestland and water area among the six LULC types during the six research periods. By contrast, the highest average LST occurred in construction land (1990, 2000, 2005, 2010 and 2015) or unused land (1995) because of the bare soil in those two LULC types caused by soil erosion, which mainly affected by local anthropogenic activities, such as land reclamation, oil exploitation, and urbanization. A comparison of the average LST among the two regions (UA and SUA) showed that the average LST for each LULC type in the descending order by UA and SUA. Thus, these were significant spatial gradients in the LST from the city center to the surrounding suburban area. In Figure 4, the length of the error bar represents the spatial difference in the LST, where a longer error bar indicates the higher heterogeneity of the corresponding thermal landscape between the two regions (UA and SUA) for each LULC type. In general, the length of the error bar for each LULC type in the ascending order by UA and SUA, from 1995 to 2015, which also suggested that the heterogeneity and fluctuations in the average LST decreased gradually along a suburban-urban gradient in Yan'an City.

Relationships between the Spatial Distributions of SUHI and LULC
In this study, we selected six Landsat images (Table 1), where the distribution of urbanized areas and SUHIs (medium-hot and hot areas) changed from a sparsely dotted pattern in 1990 to a chain/areal pattern in 2015 due to the gradual increase in the urban area of Yan'an city during the study period. A comparison of the LULC maps derived from images acquired on 29 August 1990 and 1 July 2015 shows the dramatic expansion of the urbanized area by 2015 ( Figure 5). In 1990, the hot areas only occurred on construction land at the center of Yan'an city (Figure 5a According to the spatial-temporal analysis of the variations in the SUHI intensity at 5-year intervals, construction land had the highest LST. There were significant spatial gradients in the LST from the city center to the surrounding suburban and rural areas. Based on land use/cover changes trajectory tracking, from 1990 to 2015, the amount of construction land in the urban area increased by 56.72 km 2 , from 11.53 km 2 to 68.25 km 2 , where 33.32 km 2 of UGS was transformed into construction land. Therefore, 42.42% (17.61 km 2 ) of the new land use types (41.52 km 2 in all) were transformed from UGS. Table S2 shows the time series of the six types of land use/cover area in percentage.

Relationship between SUHI and LULC
At the landscape scale, the spatiotemporal changes in the six landscape metrics (FRAC, DI, AI, SHDI, PLAND and EI) for the thermal landscape (classified LST) were similar to those in the corresponding urbanization landscape metrics for LUCC (especially construction land) during the last 25 years. The landscape metrics showed that long-term anthropic intervention (petroleum exploitation and rapid urbanization) produced a highly fragmented and diversified landscape, as shown by the decrease in FRAC for LULC and thermal landscape because of the decline in the natural landscape ( Figure 6). The decreases in DI of LULC and thermal landscape of urban area, except for thermal landscape of suburb area (Figure 6b), and the increases in AI, PLAND and EI for urban area (construction land) because of urbanization (Figure 6c,e,f). The decrease in SHIDI (Figure 6d) showed that there was a less dispersed and even distribution of land-use types, especially in the new Yan'an city, which was confirmed by the landscape spatial distribution analysis ( Figure 5).

Effects of UGS Size, Shape, and Tree-Layer Structures on GSCI
Based on GSCI order, after considering the diversity of tree species and the cooling effect, the top three tree species compositions for UGSs in the urban area of Yan'an city were: (1) 6PT + 4RP + 3PO + 3SJ + 3PCa + 2SC + 2PH + 1PA (2) 20PB + 36SJ + 7SM + 4As + 2GB + 1PH + 1PT, and (3) 15GB + 5UP + 3JF.At the sampling plot scale, the linear regression analysis based on ordinary least squares showed that there were significant strong positive relationships between GSCI and LAI, and the logarithm of TGB, and there was a negative correlation (p < 0.05) between GSCI and LSI when LSI between 1.062 and 1.717, while GSCI LSTm had a significant positive correlation with GSCI AT1.5 according to their corresponding independent variables (Table 7).  Figure 7 shows the following results: (1) There was a significant strong positive correlation between LST and IBI, and inverse correlations between LST and FV, NDVI, and MNDWI in both the urban and suburban areas during the six periods. (2) The correlation coefficients (Pearson's r) differed for the ground remote sensing indexes (FV, NDVI, IBI, MNDWI) versus LST, where in urban areas, IBI was the key positive factor related to the increases in LST, whereas NDVI and MNDWI were the crucial factors related to decreases in LST during the six periods.

Main Reasons for the Changes in Vegetation and LST
In response to the national "Grain-to-Green" program, Yan'an Municipal People's Government issued the "Closing the mountain and planting trees, feeding animals in shed" act on 16 October 1999. Subsequently, deforestation was replaced by afforestation. Moreover, the establishment of Yan'an city as a "National Forest City" was affirmed by the State Forestry Bureau of China on 29 January 2012. Yan'an Municipal People's Government initiated the "National Forest City" establishment project and issued related operational documents, such as "The overall planning of construction for national forest city in Yan'an" and "Overall planning of Yan'an City (2012-2030)." However, specific weather conditions and large-scale development projects affected the urban forest and its construction. For example, torrential rain occurred during 2013, where the total precipitation in 17 days (from 28 May to 7 October) was 725 mm, while the annual precipitation was 929.1 mm, which was nearly twice of the normal year. Petroleum exploitation and rapid urbanization, including the leveling of mountains to build the new city, meant that all the vegetation was removed from over 30 hills, i.e., a total area of 78.5 km 2 . The changes in vegetation were significantly influenced by these land use policies and human activities. In the suburban area, due to ecological restoration, the average vegetation cover (forest land) increased by 56.54% (from 34,172 ha to 53,493 ha) between 1990 and 2015, which was similar to the trend in Baota District [66]. By contrast, in the urban area, the average vegetation cover (forest land) decreased by 11.15% (from 17,661 ha to 15,692 ha) between 1990 and 2015. The leveling of mountains to build the new city area was the major cause for the decrease in the urban vegetation coverage and the corresponding increase in the LST [56].

Correlations between Different Land Surface Indicators and LST
Odindi et al. [33] noted that impervious surfaces are heat sources, whereas green spaces are the major heat sinks in Ethekwini Municipal Area, South Africa. Impervious surfaces and bare land are high temperature zones [31,74], and construction land, farmland, and unused land contained high proportions of impervious surface and bare land in our study area. Our results showed that construction land, farmland and unused land had higher average LST values than the other types of land use during the six study periods (Figure 4). Our results also showed that the four ground indexes were all significantly correlated with the mean LST, i.e., a positive correlation with IBI and negative correlations with FV (vegetation fraction), NDVI (normalized difference vegetation index) and MNDWI (modified normalized difference water index), shown in Figure 7. These findings are consistent with those of other studies. For example, Xu [64] found a significant positive relationship between the average LST (land surface temperature) and IBI (an index-based built-up index), while other studies reported significant negative relationships between the mean LST and the amount of impervious surface [79,80], NDVI [43,51,81] and MNDWI [64]. Related studies had shown that the key factors affecting urban LST are not only land cover patterns, but also other anthropogenic forces, especially land use. Therefore, the explanation of urban LST or SHUI (surface urban heat island) by land cover alone is inadequate. Especially at fine spatial scales, information on land use is more meaningful than that of land cover to indicate the impacts of urbanization on ecosystems [3,5,36,41,47,51,52,81].

Spatial Characteristics of LULC and Variations in Vegetation vs. LST Along the Urban-Rural Gradient
We found that the zones closest to the city center did not have the highest average LST values in all six periods ( Figure 5). Considering that the percentage of construction land and forest (Table  S2) had a strong significant correlation with the mean LST, then PLAND (percentage of landscape)for construction land and forest helped to explain why the zones close to the city center did not always have the highest mean LST (Figure 5c). The changes in the spatio-temporal distribution of the high LST zone were clearly divided into two stages: in the first stage from 1990 to 2005, the high LST zones had a plaque-like distribution in farmland (Figure 5a-h); and in the second stage from 2010 to 2015, the high LST zones were distributed mainly in construction land (Figure 5i-l). Previous studies have shown that some screened landscape pattern metrics can be used to delineate the effects of uraban heat islands (UHIs) [44,70,72,82]. In the present study, six urban landscape evaluation indexes. i.e., FRAC (fractal dimension index), PLAND, AI (aggregation index), DI (division index), SHDI (Shannon's diversity index) and EI (expansion intensity) and previous research conclusions were employed to calculate the thermal and land-type characteristics of the two regions (urban area and suburban area) in the study area. The PLAND indicator showed that the area of impervious surface (construction land) and land surface temperature increased in the study area. The EI indicator demonstrated that the intensity of the increase in the construction land and land surface temperature were greater during 2005-2015.
The FRAC indicator showed that the regional complexity of LST and LULC were decreasing, and revealed that the reduction in patch boundary complexity and degeneration of natural landscape. The AI and DI indicators demonstrated that the spatial division and aggregation of the different LST grades changed consistent with that of the land use/cover through the study periods ( Figure 6). These findings highlight the relationship between the LULC pattern and LST at the 100 m landscape scale.

Major Factors That Influenced the GSCI Intensity
Our results demonstrate that the shape of UGSs as well as features of the urban forest structure (TSN, LAI and TGB) in UGSs significantly affected the magnitude of the GSCI (green space cooling island) intensity (Tables 6 and 7). We found significant positive relationships between GSCI and TSN, LAI and GSCI AT1.5 . Ren et al. [24] also found that the park cooling island intensity had significant positive correlations with features of the urban forest structure in parks (i.e., canopy density, LAI, basal area, height, diameter and stem density), while Vidrih and Medved [83] showed that the optimal length for a park with LAIsp = 3 (specific dimensionless coefficient of the leaf area) to achieve the best cooling intensity is 130 m. Tree volume had the highest impact on the nocturnal UHI intensity in Amsterdam within 40 m and a one degree reduction in temperature was predicted for an increase in the tree canopy volume of 60,000 m 3 in its 40 m buffer [84].

Conclusions
Urban heat islands (UHIs), especially surface urban heat islands (SUHIs), are mainly influenced by rapid local urbanization. The spatial variations in the LST levels in Yan'an City were related to SUHIs and LUCC (land use and cover change) from 1990 to 2015 with five-year intervals were revealed in this studies. The higher LST values were usually distributed in low vegetation cover land types (construction land, farmland and unused land). Land surface remote sensing indexes, an index-based built-up index (IBI), normalized difference vegetation index (NDVI), vegetation fraction (FV) and modified normalized difference water index (MNDWI) and landscape pattern metrics, i.e., Fractal dimension index (FRAC), percentage of landscape (PLAND), aggregation index (AI), division index (DI), Shannon's diversity index (SHDI) and expansion intensity (EI) described the spatiotemporal relationships well between SUHI and LUCC at the landscape scale (100 m) in our study area. Urban green spaces (UGSs) can mitigate SUHIs at the sample scale (20 × 20 m) and reduce the LST by 0.67 • C to 8.57 • C. Our results demonstrated that several factors affected the GSCI (green space cooling island) intensity, i.e., positive effects of the tree canopy (LAI) and TSN (number of species comprising the tree layer of a forest in an UGS, diameter at breast height above 5 cm, with a certain number), a negative effect of UGS shape (LSI), logarithmic positive effects of the spatial green mass of trees and shrubs, i.e., the space volume occupied by living plant stems and leaves (TGB), and complex relationships with the configuration of tree layer species. We found that mixed forest was better than that pure forest for mitigating SUHIs. According to the important value (VI) of tree species in the 34 sample plots (Table S3) Table S1: The constants for LST retrieval with different types of thermal infrared data, Table S2: The time series of the area percentage of the six land use/cover (%), Table S3: The importance value of tree species in the 34 sample plots of the UGSs (%), Figure S1: Fluctuations in the monthly (a) and yearly (b) mean air temperature from 1990 to 2015, where the mean value was calculated using data obtained from Yan'an city weather stations, Figure S2: The 20-min variations in the land surface heat island intensity during July 2015 under sunny weather at 34 open impervious surface sites located close to various green spaces in the core urban area of Yan'an City. The mean (blue line) and standard deviation (gray line) for each 20 min period are also shown.