Global Land Surface Temperature Influenced by Vegetation Cover and PM_2.5 from 2001 to 2016

Abstract

Land surface temperature (LST) is an important parameter to evaluate environmental changes. In this paper, time series analysis was conducted to estimate the interannual variations in global LST from 2001 to 2016 based on moderate resolution imaging spectroradiometer (MODIS) LST, and normalized difference vegetation index (NDVI) products and fine particulate matter (PM_2.5) data from the Atmospheric Composition Analysis Group. The results showed that LST, seasonally integrated normalized difference vegetation index (SINDVI), and PM_2.5 increased by 0.17 K, 0.04, and 1.02 μg/m³ in the period of 2001–2016, respectively. During the past 16 years, LST showed an increasing trend in most areas, with two peaks of 1.58 K and 1.85 K at 72°N and 48°S, respectively. Marked warming also appeared in the Arctic. On the contrary, remarkable decrease in LST occurred in Antarctic. In most parts of the world, LST was affected by the variation in vegetation cover and air pollutant, which can be detected by the satellite. In the Northern Hemisphere, positive relations between SINDVI and LST were found; however, in the Southern Hemisphere, negative correlations were detected. The impact of PM_2.5 on LST was more complex. On the whole, LST increased with a small increase in PM_2.5 concentrations but decreased with a marked increase in PM_2.5. The study provides insights on the complex relationship between vegetation cover, air pollution, and land surface temperature.

1. Introduction

Climate change has widely impacted human society and the natural environment, and an increasing number of studies have researched global warming [1]. According to the Intergovernmental Panel on Climate Change (IPCC) Fifth Assessment Report (AR5), the average global surface temperature increased by approximately 0.85 °C from 1880 to 2010 [2]. Ji et al. [3] found that air temperature had an obvious warming trend over the global land surface, particularly in northern mid-latitudes. Meanwhile, the population of the world is increasing gradually, which may result in more consumption of fossil fuel and aggravates the global warming.

Land surface temperature (LST), the temperature of the Earth’s skin, plays an important role in heat and energy exchange between land surfaces and the atmosphere [4,5]. Knowledge of LST is a cornerstone of others earth system science [6,7]. As such, it has been used in soil moisture, urban heat island, evapotranspiration [8,9,10,11], and other studies. Nevertheless, because of the strong heterogeneity of land surface pattern, the accuracy of LST must be retrieved at 1K or better [7]. Thus, a large number of studies have focused on retrieval algorithms [7,12,13] and new validation schemes of LST [14,15,16]. Some studies indicated that location, land cover, vegetation fraction, and elevation had important influences on LST [17,18]. For example, compared with other land cover types, bare land has the highest mean LST and obvious seasonal characteristics [19]. Carlson et al. [20] analyzed surface biophysical parameters through the relationship between vegetation cover and LST. The results showed that LST could be used to monitor ecosystem development. In addition, forest fires [21], urban heat islands [22] and climatic anomalies such as La Niña and El Niño all affect LST.

The distribution and variation of vegetation is becoming increasingly important in the ecosystem, and they can also be used to estimate the effects of global warming, biodiversity reduction or habitat degradation [23]. Numerous vegetation indices have been designed to reflect vegetation growth by quantifying foliage density and canopy properties [24]. The normalized difference vegetation index (NDVI), with a range from −1 to 1, is one of the most widely used parameters to estimate land surface roughness and monitor vegetation cover and growth [4] due to the positive relation with characteristics of plant coverage [21]. Moreover, many different ecological and environmental issues have been researched through the NDVI [23,25,26]. It is a reliable reference to assess the potential connection between the variation in regional ecosystems and biocenosis due to the spatiotemporal changes in climatic conditions [27,28,29]. Several studies have analyzed the relationships between LST and NDVI [20,29,30,31]. An opposite pattern of the spatial distributions of NDVI and LST was found in a typical karst area, the northeastern part of Guizhou Province in China [18]. Seasonal difference of the relationship between LST and NDVI was identified in the North America [32,33]. The positive relationship occurred in cold seasons but negative relationship appeared in warm seasons [32]. Other studies researched drought by monitoring LST and NDVI [34]. For example, an increasing trend of aridity was found in arid and semi-arid areas of Southern Europe by monitoring LST and NDVI via Harmonic Analysis of Time Series software from 1982 to 1999 [31].

Fine particulate matter (PM_2.5) is a kind of air pollutant, with tiny air droplets or particles with a diameter of 2.5 microns or less [35]. Due to rapid industrialization and fossil fuel burning, PM_2.5 increased rapidly in many countries in the world, and seriously damaged human health and environment [36,37]. In China, due to the coal combustion, the emission amount of air pollutants and days of poor air quality increased markedly in mega cities [38,39]. Direct damage on human health such as lung cancer, heart disease, and acute respiratory infections have been reported [40]. The variations of LST and surface urban heat island (SUHI) have drawn extensive attentions. The change of land surface properties can cause the variation in albedo and energy storage in the urban surface [41]. The increasing anthropogenic activities produced more air pollutants including particular matter/aerosols [42]. High aerosol optical depth occurred mainly in city central areas, much higher than the surrounding rural areas [43]. By scattering and absorbing solar radiation, aerosols reduced urban surface temperature [44]. Since the industrial revolution, the cooling effect of the increase in aerosol concentrations on the atmosphere has been found and it may mask some of warming effect due to the growing emissions of greenhouse gas [45]. Previous studies indicated that air pollution might produce a fog-like haze which could reduce the solar radiation and might affect the precipitation and microphysical properties of clouds [42,46,47]. For example, haze caused a 0.7 K temperature difference between urban and rural area at nights in semi-arid cities of China [48]. However, most of current studies about PM_2.5 researched the influence on meteorological conditions at a regional scale. Only a few studies have researched the spatio-temporal characteristics of PM_2.5 at a global scale.

Ground-based measurements could provide the most reliable values of land surface temperatures. However, due to limitations caused by sparsity and the inhomogeneous distribution of observation stations, these measurements are insufficient for global-scale research. For example, in some regions, such as Africa and the Antarctic, where very few stations are available and therefore this will prevent us from understanding the changing surface temperatures and extreme weather [17]. New observation methods are still needed. With the development of remote sensing technology, remote sensing images have been increasingly used on assessing changes in large spatiotemporal scales. The modern satellite techniques can provide continuous and temporally repetitive images and long-term observations at the global scale, and it provides the only possibility to quickly measure really information of the entire globe [7,24,49]. The moderate resolution imaging spectroradiometer (MODIS) is an important sensor that is installed on the Terra and Aqua satellites of the Earth Observing System (EOS) and operated by the National Aeronautics and Space Administration (NASA). Due to combination of the global coverage, intermediate spatial resolution, higher temporal resolution, and abundant products, MODIS data have been widely used in environmental research [50,51,52].

In short, previous studies on LST mainly focused on retrieval algorithm, accuracy verification and urban heat island at the regional scale. There is a lack of comprehensive understanding among LST, vegetation cover and PM_2.5 at a global scale. Therefore, it is necessary to explore the influence vegetation cover and PM_2.5 on LST. In this study, Terra MODIS LST and NDVI data were obtained from the Land Processes Distributed Active Archive Center (LP DAAC) managed by the NASA Earth Science Data and Information System (ESDIS) project. PM_2.5 data were provided by Dalhousie University Atmospheric Composition Analysis Group. Multitemporal data were used to analyze the spatiotemporal variations of LST. The effects of vegetation cover and PM_2.5 on LST were estimated during the period of 2001–2016 at a global scale. The objectives of this study are (1) to determine the spatiotemporal changes in LST, vegetation cover and PM_2.5 concentrations; (2) to explore the relationship between LST and vegetation cover changes measured by the seasonally integrated normalized difference vegetation index (SINDVI) and atmospheric contamination variations measured by PM_2.5 concentrations; and (3) to analyze the factors influencing LST.

2. Materials and Methods

2.1. Data and Processing

2.1.1. LST

LST data were collected by the MODIS sensor, a key instrument of the Terra and Aqua satellites. The LST products were extracted from the MOD11C3 collection 6 dataset from 1 January 2001 to 31 December 2016, which provided monthly composited and averaged LST at a 0.05° spatial resolution (5.6 km × 5.6 km at the Equator), and obtained from Earthdata Search (https://search.earthdata.nasa.gov/). There are 17 scientific data sets (SDS) in this product. The band named LST_Day_CMG was extracted as day time LST to analyze the variation of LST. The standard deviation of the MODIS V006 LST errors were within 1 K [53]. The mean value composite (MVC) method was used to calculate annual average values of LST in the period of 2001–2016.

2.1.2. NDVI

The Terra MODIS NDVI products were extracted from MOD13C2 collection 6 datasets, which were composed of 16-day global composite images from January 2001 to December 2016, and provided by the Earthdata Search, with a spatial resolution at 0.05°. It is derived from the red and near-infrared reflectance ratio as follows [23]:

NDVI = (NIR − R)/(NIR + R)

(1)

where NIR and R are the amounts of near-infrared (the reflectance of band 2 of MODIS) and red light (the reflectance of band 1 of MODIS).

Although the product has eliminated some clouds and bidirectional reflectance of ground objects [54], MODIS NDVI time series data are still influenced by the clouds and atmospheric water vapor. A simple but efficient method, the mean-value iteration filter, was used to reduce the noise and reconstruct a high quality of NDVI time series by using the following equation [55]:

Δ_i = |NDVI_i − (NDVI_i−1 + NDVI_i+1)/2|

(2)

where NDVI_i represents the ith monthly NDVI, with range from 1 to 192 (16 years). A small percentage of the multiyear NDVI average of each pixel was set at a threshold (Δ). The results of several test showed that 10% was the optimal threshold in this study. When Δ_i > Δ, the NDVI_i was replaced by (NDVI_i−1 + NDVI_i+1)/2. The procedure is iterated until all Δ_i are within the given threshold.

The seasonal characteristics of plants can cause substantial differences in vegetation cover between seasons. Thus, the seasonally integrated normalized difference vegetation index (SINDVI), sum of NDVI values exceeding the threshold (commonly defined as NDVI > 0.1) [56,57,58], was used to monitor the interannual change of vegetation. It can identify vegetation distribution and growth state simultaneously. Due to the large influence of intermittent climatic factors such as precipitation on vegetation growth, compared with NDVI, SINDVI could better describe the interannual variation [59].

2.1.3. PM_2.5

PM_2.5 data were significantly consistent (R² = 0.81) with the monitoring results of PM_2.5 concentrations [60] obtained from the Atmospheric Composition Analysis Group (http://fizz.phys.dal.ca/~atmos/martin/) with a spatial resolution of 0.1°. To match the resolutions of the LST and SINDVI datasets, the data were resampled to 0.05° latitude × 0.05° longitude using the nearest neighbor method. The unit for these data was μg/m³. These data only included the region spanning from 69.8°N to 54.8°S and 180°E to 180°W. The objective of this study was to detect global change trends and impact of PM_2.5 on LST. Because of lower PM_2.5 concentrations in the Arctic and Antarctica, these data could meet the target of this study.

2.2. Methodology

2.2.1. Inter-Annual Variation Analysis

To determine the spatiotemporal differences, the annual change tendency was quantified at a global scale. Monthly time series data were processed to form yearly time series data; in this way, for each pixel the time series has a length of 16, spanning from 2001 to 2016, and the number of pixels in each band are 3600 × 7200. Linear regression and ordinary least squares methods were used to estimate the rate and range of changes at regional and global scales [9,57,61,62]. Random noises within 1K were added to LST data to test the sensitivity of constant biases. It was worth nothing that the slope and range were not affected by possible constant biases present in the data. Therefore, the slope and range were used to assess the interannual variations in LST, SINDVI, and PM_2.5. The slope represents the change rate of each pixel from 2001 to 2016 calculated by the ordinary least-squares estimation via linear regression as follows:

$$slope=\frac{n\times {\sum }_{i=1}^n\left(i\times A_i\right)-{\sum }_{i=1}^{n}i\times {\sum }_{i=1}^n{A}_i}{n\times {\sum }_{i=1}^n{i}^2-{\left({\sum }_{i=1}^{n}i\right)}^2}$$

(3)

where n represents the length of the time series (in this research, n = 16); A_i represents the ith year for SINDVI, PM_2.5 or the average LST; and the value of i ranges from 1 to n.

The range can evaluate the magnitude of change. Each pixel value represents the total variation during the period of 2001–2016.

$$range=slope\times \left(n-1\right)$$

(4)

where n is similar to the definition of slope and represents the length of the time series.

If slope > 0, it means that there is an increasing trend during the past 16 years, whereas the opposite means that it is decreasing. A large absolute value indicates a dramatic change, while the absolute value close to 0 indicates little variation.

2.2.2. Correlation Analysis

For clear descriptions of the relationships between LST with SINDVI and PM_2.5, Pearson’s correlation coefficients (r) was used to indicate the correlation. It was conducted as follows:

$$r_{xy}=\frac{{\sum }_{i=1}^n\left(x_i-\overline{x}\right)\left(y_i-\overline{y}\right)}{\sqrt{{\sum }_{i=1}^n{\left(x_i-\overline{x}\right)}^2}\sqrt{{\sum }_{i=1}^n{(y_i-\overline{y})}^2}}$$

(5)

where x_i and y_i represent LST and either SINDVI or PM_2.5 in the ith year, respectively; and $\overline{x}$ and $\overline{y}$ represent the corresponding averages over 16 years.

At the same time, partial correlations were conducted to assess the dominant factor of each pixel. The partial correlation analysis refers to the process of removing the influence of the third variable when two variables are related to the third variable at the same time, and it only analyzes the correlation between the other two variables [56]. When the degree of correlation between LST and SINDVI was estimated, the influence of PM_2.5 could be removed. The partial correlation coefficient was calculated as follows:

$$r_{xy\left(z\right)}=\frac{r_{xy}-r_{xz}{r}_{yz}}{\sqrt{\left(1-r_{xz}^2\right)\left(1-r_{yz}^2\right)}}$$

(6)

where r_xy(z) is the partial correlation coefficient, representing the degree of correlation between x and y without z; and r_xy, r_xz and r_yz are simple correlation coefficients.

The C program was used to calculate the average LST for each year from 192 scene images, reconstruct the NDVI, summarize the NDVI during growing season in order to obtain the annual SINDVI and analyze the correlation. The Mollweide projection, an equal-area pseudo-cylindrical map projection, was used to count regional areas.

3. Results

3.1. Spatiotemporal Changes in LST, SINDVI, and PM_2.5 Concentrations

The global land surface temperature presented an increasing trend, with a global mean rise of 0.17 K from 2001 to 2016 (

Table 1. Changes of LST, SINDVI, and PM_2.5 at different latitudinal bands from 2001 to 2016.

	Global	60°–90°N	30°–60°N	0°–30°N	0°–30°S	30°–60°S	60°–90°S
LST	0.17	1.10	0.23	0.07	0.47	0.31	−0.56
SINDVI	0.04	0.03	0.12	0.05	0.01	0.01	0.00
PM2.5	1.02	−0.20	−0.19	4.03	1.20	0.47	--

). LST increased at various degrees across 62.4% of the land area, with the exception in some regions such as the Antarctic, central North America, Central Asia, and eastern Australia (Figure 1a). An LST increase of 15 K was found in the Aral Sea, where water level decreased by 23 m, surface area decreased by 74%, and lake volume decreased by 90% [63], experiencing extreme desiccation [62]. The temperature of most regions increased within the range between 0 and 3 K, which accounts for 61.9% of the global land area. The Arctic experienced a marked warming trend, especially in northern Russia; however, the Antarctic, an area of approximately 14 million km², showed a cooling trend. In the mid-latitudes of the Northern Hemisphere, such as southeast Canada, northern United States, eastern Kazakhstan, Mongolia, and Northeast China, dramatic decreases in LST were detected.

Variable solar radiation causes different climate zones and LST along the latitude. As shown in Figure 1a and Table 1, LST had an increasing trend, except for the areas at high latitudes in the Southern Hemisphere. In the high latitudes of the Northern Hemisphere, the change range increased firstly before a sharp decline. LST increased by an average of 0.30 K in the mid-latitudes of northern hemisphere, and it fluctuated slightly in the low-latitude areas. Compared with the Northern Hemisphere, there were two distinct peaks both at low- and mid-latitudes, and the larger one was an increase of 1.85 K at 48°S.

SINDVI presented a greening trend all over the world [57], with an average value of 0.04 during the period of 2000–2016. Approximately 58.4% of the global area showed an increasing trend. Arctic, Antarctic, and high mountains were coved by snow for most time of the year. Approximately 9.8% of the global land was covered with snow all the year round, and SINDVI had no change in these regions. About 31.8% of other areas experienced a downward trend. Approximately 54.5% of the area’s increase range of SINDVI was between 0 and 1, and 40.5% ranged from −1 to 0. The vegetation cover in central and southeastern China, northwestern India, and eastern Australia increased markedly. In addition, an increase appeared in Europe and most parts of North America. However, the SINDVI decreased significantly in eastern Brazil, northern Argentina, northern Kazakhstan, southeastern Africa, and western Australia.

The latitudinal SINDVI varied from −0.11 to 0.33, with an obvious increasing trend in the Northern Hemisphere, which indicated the great increase in the vegetation coverage in the Northern Hemisphere (Figure 1b and Table 1). SINDVI in Southern Hemisphere (between 15°S to 30°S) also increased, with fluctuation within ±0.15. However, a dramatic decrease (from 0.32 to −0.09) was found between 34°S and 40°S, and there were more rapidly latitudinal fluctuations between 40°S and 55°S.

PM_2.5 concentrations also had an increasing trend (1.02 μg/m³). The global distribution of PM_2.5 was shown in Figure 1c. The air quality in America, Europe, and eastern Russia improved greatly. PM_2.5 concentrations in Africa as a whole showed an increasing trend, and most regions had increased by 10 μg/m³. PM_2.5 in Asia, such as Saudi Arabia, India, and China, increased noticeably, indicating serious air pollution. Slight fluctuations were found in Australia.

In the Northern Hemisphere, PM_2.5 concentrations decreased at mid-latitudes, but increased markedly at low latitudes. In the Southern Hemisphere near the tropics, the concentrations of PM_2.5 increased slightly.

3.2. Correlation Analysis

The Pearson correlation coefficients were calculated based on each pixel, as shown in Figure 2 (left). The results indicated a significant positive correlation between the LST and SINDVI in most areas of the Northern Hemisphere. In contrast, a highly significant negative correlation appeared in most areas of the Southern Hemisphere, especially Brazilian plateau, Pampas steppe, southeastern Africa, and Australia.

Latitudinal patterns of Pearson correlation were shown in Figure 2 (right). In the Northern Hemisphere, a positive correlation between LST and SINDVI was found at high latitudes, but a negative correlation was found in other latitudinal zones. In the Southern Hemisphere, the correlation between LST and SINDVI increased then decreased along with a rising latitude. Strong negative correlations between LST and SINDVI occurred at 15–35°S (i.e., mid-to-low latitudes in the Southern Hemisphere), and the negative correlation gradually weakened from latitude of 40°S.

Except for low latitudes in the Northern Hemisphere, a positive correlation between LST and PM_2.5 appeared around the world. Compared with regions in north of the equator, those in south of the equator had stronger positive correlations. Positive correlations between LST and PM_2.5 existed in central Eurasia, South America, and eastern Australia, where PM_2.5 fluctuated slightly. The negative correlations between LST and PM_2.5 were mainly concentrated in Europe, West Asia, and south of Sahara Desert. In these areas, PM_2.5 experienced dramatic changes, increasing in West Asia especially in Saudi Arabia and decreasing in Europe such as Poland, Germany, and Romania. Therefore, the effect of PM_2.5 on LST was complex, and it was related to the changes of PM_2.5 concentrations. A small amplitude increasing of PM_2.5 concentration increased LST due to the heat retention effect of PM_2.5 on global surfaces, with a result of a positive correlation. However, a marked increase in PM_2.5 concentration prevented solar radiation from reaching the surface, reducing surface energy and causing a decrease in LST.

3.3. Dominant Factor Analysis

The partial correlation coefficients (i.e., r _{LST SINDVI (PM2.5)} and r _{LST PM2.5 (SINDVI)}) were calculated based on pixels to analyze the main factors causing surface temperature changes (Figure 3). If r _{LST SINDVI (PM2.5)} was larger than r _{LST PM2.5 (SINDVI)}, the SINDVI was considered to be the main factor. The opposite result indicated PM_2.5 had a larger effect on LST.

In most parts of the world, SINDVI was found to be the dominant factor (Figure 3). SINDVI increased markedly in the mid-latitude of the Northern Hemisphere, especially over the East European Plain, Labrador Plateau, and India. It decreased obviously in east of the Brazilian Plateau and southeast of Africa. The decrease in vegetation cover caused by human activities led to the increase in LST.

3.4. Significant Decreases in LST on Regions of Interest (ROIs)

Four regions of interest (ROIs) where LST decreased sharply were found in Saudi Arabia (ROI-a), India (ROI-b), Northeast China (ROI-c), and Australia (ROI-d). The spatial positions of ROIs were shown in Figure 1a in yellow and the locations were shown in

Table 2. Locations of four regions of interest (ROIs).

Region Code	Bottom Left Corner	Top Right Corner	Pixel Numbers
ROI-a	25.85°N, 47.30°E	27.30°N, 48.50°E	696
ROI-b	21.30°N, 73.25°E	25.05°N, 76.20°E	4425
ROI-c	46.10°N, 125.20°E	47.30°N, 127.20°E	960
ROI-d	34.15°S, 145.05°E	30.10°S, 148.70°E	5913

.

3.4.1. Saudi Arabia

LST in Saudi Arabia (ROI-a) has decreased by 2.09 K over the past 16 years (Figure 4a). SINDVI in this region also decreased with a slope of −0.03, but PM_2.5 increased dramatically, with a slope of 1.512. The interannual variations in ROI-a were shown in Figure 4a. The correlation analysis showed a significantly positive correlation between LST and SINDVI (r = 0.661, p < 0.01) and a significantly negative correlation between LST and PM_2.5 (r = −0.746, p < 0.01), which means that both SINDVI and PM_2.5 contributed to the changes in LST. However, the partial correlation analyses indicated lower confidence levels (p > 0.05,

Table 3. Statistical analysis of relation between LST with SINDVI and PM_2.5 in the regions of interest (ROIs).

	Simple Correlation Coefficient				Partial Correlation Coefficient
	r LST SINDVI	p	r LST PM2.5	p	r LST SINDVI (PM2.5)	p	r LST PM2.5 (SINDVI)	p
ROI-a	0.661 **	0.005	−0.746 **	0.001	0.249	0.371	−0.511	0.051
ROI-b	−0.877 **	0.000	−0.630 **	0.009	−0.789 **	0.000	−0.106	0.708
ROI-c	0.336	0.203	−0.416	0.109	0.569 *	0.027	−0.608 *	0.016
ROI-d	−0.849 **	0.000	0.743 **	0.001	−0.736 **	0.002	0.513	0.051

* indicates a statistically significant p-value (p < 0.05); **indicates p < 0.01.

). Therefore, the effects of vegetation cover cannot be ignored in this area.

3.4.2. India

LST experienced a marked decrease as SINDVI and PM_2.5 increased in India (Figure 1). The interannual variations in ROI-b were shown in Figure 4b. The variation in vegetation cover played the leading roles in affecting LST changes during the period from 2001 to 2016 (Figure 3). The correlation analysis showed that LSTs were significantly correlated to both SINDVI and PM_2.5, with Pearson correlation coefficients of −0.877 (p < 0.01) and −0.630 (p < 0.01), respectively. However, the partial correlation analysis indicated that the correlation between LST and SINDVI (r _{LST SINDVI (PM2.5)} = −0.789, p < 0.01) was more significant than that between LST and PM_2.5 (r _{LST PM2.5 (SINDVI)} = −0.106).

3.4.3. China

LST showed an increasing trend in West China, but a reverse trend in East China, particularly in Northeast China (Figure 1a). The interannual variations in ROI-c were shown in Figure 4c. Interestingly, the partial correlation analyses indicated a significantly positive correlation between LST and SINDVI and a significantly negative correlation between LST and PM_2.5, while the simple correlation coefficients were smaller than the partial correlation coefficients (Table 3). Therefore, the interaction between SINDVI and PM_2.5 weakened their respective effects on LST in the Northeast China.

3.4.4. Australia

LST in Australia showed a decreasing trend, with a noticeable increase in SINDVI and a slight decrease in PM_2.5 (Figure 1). The interannual variations in ROI-d were shown in Figure 4d. In 2010, the value of LST and PM_2.5 were low, but SINDVI reached its peak (Figure 4d). The simple correlation analysis showed significant relationship between LST with SINDVI and PM_2.5, but the partial correlation analysis indicated that SINDVI had a larger impact on LST than PM_2.5 (Table 3).

4. Discussion

4.1. Potential Causes for Variations in LST, SINDIV, and PM_2.5

LST is very sensitive to the global climate changes. It had a prominent increasing trend in the Arctic, but a decreasing trend over the Antarctic. This is consistent with the conclusion of Overland et al. [64] that the North and South poles experienced opposite climate changes. Authors identified anthropogenic influence as the primary cause for sea ice reduction in the Arctic [65]. Under the background of global warming, researchers have paid more attention to the warming in the Arctic due to the importance of global sea level rise, especially regarding the Greenland ice sheet [66]. LST could also be used to monitor glacial melt zones based on remote sensing technology. LSTs increased in the Arctic, and the surface temperature of North Greenland had risen by approximately 1 K over the past 16 years (Figure 1a). The main indicators of Arctic temperature increase include reducing sea ice and permafrost. In addition, LST could indicate the burn severity of forest fires [67]. The increase in LST (0.073–0.325 K) has deteriorated forest fires in Siberian boreal forests. In addition to direct warming by fires, lower albedo with darker surfaces and lower evapotranspiration due to biomass burning can increase energy absorption by land surface and decrease the evaporative cooling, with the consequence of higher LST [68].

Although both the Arctic and Antarctic were subjected to similar solar radiation and increasing greenhouse gas concentrations, the two regions had undergone different variations [64]. Studies showed that decreasing LSTs in the Antarctic may be related to the ozone hole [69]. Turner et al. [70] suggested that the ozone hole would shield global warming in the Antarctic. Ozone depletion in the stratosphere had intensified polar vortices, and winds over the center of the Antarctic continent had strengthened the change of weather patterns over the mainland. Strong winds have effectively separated the Antarctic from global warming, with a result of a decrease in LST and an increase in snowfall [71].

SINDVI value was used to denote interannual variations of vegetation productivity [58]. It integrated NDVI over the growth season and weaken the details of NDVI variation. Vegetation cover measured by SINDVI showed a greening trend at the global scale from 2001 to 2016 (Figure 1b). Specifically, vegetation greening trends increased in midlatitude areas of the Northern Hemisphere. Other studies demonstrated that snow cover decreased dramatically in those regions [57]. As a result of global warming, the melting of snow and ice produced a large amount of water that could be used for irrigating plant. In addition, longer growing season caused by more heat can also led to the increases in SINDVI. A dramatic browning trend was found in central Asia (Figure 1b). In this area, a warming climate in the summer worsened the aridity [72]. The decreases in precipitation and aridity were the main reasons for the decreases in SINDVI [72].

The sources of air pollutants were diversified, such as biomass burning, fossil fuel combustion, and industrialization [73,74]. PM_2.5 in the atmosphere could be emitted directly by combustion or mechanical processes and it can also be generated by chemical reactions of inorganic or organic molecules in the atmosphere [73]. The time series analyses of the interannual mean PM_2.5 concentration demonstrated an increasing trend at a global scale (Figure 1c and Table 1). Dramatic decreases in PM_2.5 were observed in the Eastern USA and Europe because of effective pollution control measures, but an opposite trend was found in Sub-Saharan Africa due to biomass burning and higher concentration of mineral dust [60]. New research showed that the emissions of sulfur dioxide had increased by 50% in India from 2007 to 2017 [75]. The sharp increase in PM_2.5 in India was influenced by both anthropogenic and natural factors such as chemical processes [73,75] and mineral dust [60]. The regions where PM_2.5 was the dominate factor for the change in LST were mainly distributed in central Russia and Kazakhstan where the concentration of PM_2.5 increased because of the biomass burning [76], and Saudi Arabia where a lot of mineral dust contributed to the increase in PM_2.5.

4.2. Influence of SINDVI and PM_2.5 on LST

This study analyzed the effects of the vegetation cover measured by SINDVI and air pollution measured by PM_2.5 on the variation in LST. Previous studies showed the lower LSTs with higher vegetation cover [18,77]. Similarly, the current study found a significant negative relationship between SINDVI and LST in the Southern Hemisphere and low latitude areas of the Northern Hemisphere. On the contrary, a positive relationship between LST and SINDVI was found in high latitude areas. Different from the low latitude of Northern Hemisphere, high latitude areas are characterized by lesser solar radiation and lower temperature. Plant photosynthesis is limited particularly by lower temperature and increase in temperature including LST will promote plant growth and NDVI values in high latitude. Therefore, a positive relationship between LST and SINDVI was found in high latitude. A study found that although most regions experienced overall greenness, the increase rate decreased and browning trend expanded in the Northern Hemisphere especially in northern mid-low latitudes [26]. It could be the reason for that the negative relationship between SINDVI and LST in the Northern Hemisphere was weaker than that in the Southern Hemisphere (Figure 3). Furthermore, compared with past time, the increasing temperature caused a longer plant growth season. It partially caused the increase in SINDVI [78]. Therefore, the browning trend was hidden in the greenness in the Northern Hemisphere. In addition, drought can cause higher LST, but the legacies of drought were different among forest, shrubs, and grass in the North Hemisphere [79].

SUHI and the air pollution have become two major problems along with urbanization without proper environmental protection [43]. The positive correlation between LST and PM_2.5 was found in most areas from 2001 to 2016 (Figure 2b). However, negative correlation was also identified in Saudi Arabia, northeast China and south of Sahara Desert. It might be caused by a different rate of the variation in PM_2.5 concentrations. Cooling effect of aerosol was found in the urban areas [44]. Radiative effects of aerosols, reducing solar shortwave radiation and increasing atmospheric longwave radiation, have a marked cooling effect on air temperature, especially downwind of the city cluster [42]. The cooling effects of aerosol can explain the negative relation between LST and PM_2.5 in Saudi Arabia, northeast China and south of Sahara Desert where PM_2.5 increased markedly. The increase of PM_2.5 can reduce solar radiation reaching land by scattering more radiation back to space, but it could also reduce the surface longwave flux and latent heat flux which may increase LST [80,81].

In addition, Hajiloo et al. [37] used a geographically weighted regression model to assess the relationship between environmental parameters and PM_2.5 concentrations. They found a negative correlation between LST with PM_2.5 in western Tehran and a positive correlation in eastern Tehran which resulted from asymmetry distribution of vegetation. However, different trends of SINDVI did not cause opposite correlation in ROI-a and ROI-b (Figure 4a,b). Therefore, other factors may play a more important part in the negative correlation between LST and PM_2.5 in the two areas.

4.3. Limitations and Future Research

Same as many studies, there are limitations in the current research. Due to the spatial heterogeneity, LST can be affected by many factors, such as elevation, albedo, meteorological parameters, human activities, and others. [9,18,57,82]. However, because of the limitation of high spatial resolution remote sensing images and statistical data at a global scale, other factors rather than SINDVI and PM_2.5 have not yet been considered in this study. In addition, the enhanced vegetation index (EVI) was more sensitive than NDVI in some regions, such as tropical rainforests with high biomass [83]. The saturation of NDVI might underestimate the correlation between LST and vegetation cover. In future studies, other vegetation indexes can be used to analyze vegetation changes in the regions with saturation of NDVI. In addition, it is important to develop models including more climatic, topographical, and socio-economic factors such as precipitation, wind speed, elevation, SUHI, population, and others for better estimate the changes of LST.

5. Conclusions

In this study, the spatiotemporal trends of global land surface temperature, vegetation cover, and air quality were analyzed during the period of 2001–2016. The research evaluated the influence of vegetation cover and PM_2.5 on land surface temperature. A total of 62.4% of the world’s land surface temperature exhibited a increasing trend from 2001 to 2016. A total of 58.2% of regions had an increasing trend of SINDVI, indicating that vegetation coverage increased from a global aspect. The concentrations of PM_2.5 showed an increasing pattern mainly in developing countries, for example China and India, and a decreasing pattern mostly in developed countries, such as the USA. A significant negative correlation between vegetation and land surface temperature indicated vegetation was the main driver for land surface temperature changes in most areas on the earth. The effect of PM_2.5 on surface temperature was not uniform. A significant increase in PM_2.5 concentrations caused a decrease in surface temperature, but a slight increase in PM_2.5 concentrations increased surface temperature. In future studies, models including more parameters should be developed to more accurately analyze the variation in land surface temperature.