A Twenty-Year Assessment of Spatiotemporal Variation of Surface Temperature in the Yangtze River Delta, China

Abstract

A good understanding of the processes of land surface temperature (LST) change is important for assessing regional climate change. In the present study, we obtained the MODIS MOD11A2 LST products over the Yangtze River Delta (YRD) from 2001 to 2020. In order to comprehensively assess the spatial and temporal variability of LST in the YRD region over the past two decades, the Theil–Sen Median trend analysis and Mann–Kendall test, BFAST01 trend decomposition, and landscape pattern analysis were used in this study. We show that the rate of linear change in LST in the YRD ranges from −0.019 °C/month to 0.046 °C/month. The BFAST01 trend decomposition identifies more details of LST change and monotonic increases, reversal increase, and interruption increase are the main warming trends. The distribution of the different trend types shows strong aggregation with high spatial heterogeneity. The LST breakpoints are mainly located in the northern and southern YRD, which frequently occurred during 2010–2013. Of the various land types, breakpoints occur most frequently in cropland and high NDVI (0.5–0.7) areas, and the intensity of most of them is within 2 °C. In addition, much stronger warming occurs in urban areas than in other land types. Our study provides a better understanding of the dynamics of LST in the YRD region over the past 20 years and highlights that breakpoints cannot be circumvented in regional temperature assessment.

1. Introduction

Climate warming has been a hot issue of public and academic concern in recent decades and is gradually being confirmed by numerous studies and observations [1,2]. In response, many international organizations and countries have taken measures to mitigate global temperature rise. Land surface temperature (LST), a key indicator of global warming, is an important characteristic parameter for the physical property of the earth, which has been used in surface heat balance and ground-air energy interaction [3,4,5,6]. Tracking changes in temperature trends based on historical LST observations has been already an important component of global warming assessment [7,8].

Traditionally, temperature was often identified using in situ measurements or mobile measurements in the urban canopy layer [9,10]. With the development of remote sensing (RS) technology, RS inversion of land surface temperature (LST) has been increasingly used to examine the surface thermal environment [6]. LST based on satellite RS, e.g., Terra, Aqua, Landsat, and Sentinel, can provide continuous observational data, reflecting spatial and temporal evolutions of LST over a long term and on a large scale [11,12,13]. This provides an important pathway to assess regional and even global warming. There are often significant differences in temperature variation across regions and periods for temperature can be highly disturbed by local environments. When it comes to “temperature variation”, the fluctuating variability of temperature has been confirmed by many historical observations [14,15,16]. However, existing studies have generally based their judgments on the general trend of temperature change, and detailed information on the process of temperature change (abrupt changes, non-linear changes, etc.) is often ignored [17,18]. Compared to gradual warming/cooling trends, abrupt changes in temperature result in much stronger temperature fluctuations, which is important for the assessment of regional temperature change.

The trend analysis of RS images is a general judgment of the long-term change process of LST, which can be generally divided into linear trend analysis and non-linear trend analysis [19,20]. The common linear trend analysis includes univariate linear regression analysis and Theil–Sen Median trend analysis [21,22]. Linear trend analysis methods can evaluate the overall direction and change rate of a trend, which has been widely used in long-time series analysis. For example, Yu et al. [22] examined the trends of annual mean LST and annual amplitude in China from 2003–2018 using the Mann–Kendall test and the Theil–Sen estimator, which showed an increasing trend at a rate of 0.02 K/year. Fang et al. [23] analyzed the trend changes of the maximum, minimum, and average temperature in China from 1979 to 2018 and found that the number of cold nights and warm days gradually decreased and increased, respectively. Zhao et al. [24] studied the linear trend of surface air temperature for 16 cities in eastern China from 1909 to 2010 and reported a warming trend of 1.52 °C (100 yr)⁻¹ in the regional mean surface air temperature. El Kenawy et al. [25] assessed the long-term change and variability of temperature over northeastern Spain (1920–2006) and a significant increase in temperatures was found, but the diurnal temperature range showed less coherent variability with both positive and negative trends. Linear trend analysis can only analyze the overall direction and speed of change during long time series changes, but detailed information on trend changes is ignored. To comprehensively analyze the dynamic characteristics of features during long-term changes, new detection methods have been proposed, such as the breakpoint detection methods including LandTrendr [26] and DBEST [27], the trend decomposition methods including BFAST [28] and STL [29], and the trajectory fitting methods including Rescaled Range Analysis (R/S) [30], Empirical Orthogonal Function (EOF) [31], and Empirical Model Decomposition (EEMD) [32]. Most of these algorithms were initially applied in the dynamic monitoring of vegetation with widespread applications [29,33,34]. However, studies using these methods in the assessment of spatial and temporal variations of LST have rarely been reported.

The past decades have witnessed very rapid urbanization in the Yangtze River Delta (YRD), China, where human activities have significantly reshaped the land surface and the environment [35,36,37]. These changes perturb the LST strongly and contribute to the alteration in the dynamic trends of LST [38,39,40]. Existing studies have confirmed the warming trend of LST in the YRD, especially in urban agglomerations [40,41,42]. However, these studies, from an overall trend change perspective of LST, do not adequately reflect the dynamic processes of the regional thermal environment [43,44]. Since the detailed information, such as segmented linear transformations and abrupt changes, is ignored in the overall LST trend change analysis, these studies do not provide an in-depth assessment of LST change processes. An abrupt change is offset by an integral of a long-term linear trend spanning several years or even decades [45,46,47]. This plays an essential role in the assessment of regional and global warming but has been often underappreciated in existing studies [48].

In this study, we obtained the MODIS MOD11A2 LST products over the YRD region from 2001 to 2020 and performed a comprehensive assessment of the spatial and temporal variability of the LST through an analysis and decomposition of the trend. LST breakpoints, abrupt changes in temperature where LST exceeds previous averages, were of focus in this study, and the effects of land cover type (LCT) on the type for LST trend changes and breakpoints were analyzed. In addition, the inconsistent warming over different LCTs was also explored. Section 2 introduces the data and methodology. Section 3 presents the results and discussion and Section 4 draws the main conclusions.

2. Date and Methods

2.1. Study Area

The YRD is one of the most economically developed regions in China, consisting of Shanghai City, Jiangsu province, Zhejiang province and Anhui province. It is located in eastern China, bordering the East China Sea. The climate is subtropical monsoon and average annual temperatures range from 13–18 °C. The northern of the YRD is mainly plain, while the southern is hilly. As shown in Figure 1, there are obvious spatial differences among the LCT in the YRD.

2.2. LST, NDVI and LCT Data

Daytime MODIS MOD11A2 LST data over the YRD from 1 January 2001 to 31 December 2020 were retrieved on the Google Earth Engine (GEE). A total of 240 images of the LST data were processed to obtain monthly daytime LST data using the maximum synthesis method [49,50,51]. To avoid noise and perturbations caused by missing values, the processed data were then filtered using the Savitsky–Golay method [52,53], for best performance, the optimal filter window is set to 2 and the number of smoothing polynomials is set to 4 after testing.

The annual NDVI data and LCT data in 2010, the middle year of the study period (2001–2020), were also involved in the present study [53]. The raw NDVI data were from the 2010 MODIS MOD13A2 product and the raw LCT data were from the 2010 MODIS MCD12Q1 product. NDVI was ranked equally to five levels of <0, 0–0.25, 0.25–0.5, 0.5–0.75, and 0.75–1, and The LCT reclassification is based on the nature of the subsurface as woodland, grassland, cropland, built-up area, and others, as shown in Table 1.

2.3. Theil–Sen Median Trend Analysis and Mann–Kendall Test

The Theil–Sen Median trend analysis method can be combined with the Mann–Kendall test, reflecting the trend changes of each pixel in the time series. This combination does not require the data to obey a certain distribution. The Theil–Sen median trend analysis is a robust statistical method for the extraction of trend, and the formula is as follows:

(1)$S_{LST}=\mathrm{Median}\left(\frac{LS{T}_j-LS{T}_i}{j-i}\right),$

where $i$ and $j$ represent different months, January 2001–December 2020 in the present study. A positive (negative) $S_{LST}$ represents a warming (cooling) trend.

The Mann–Kendall test is used to evaluate the significance of a trend. It is a non-parametric statistical test and is free from the interference of outliers. This method has the advantage that samples do not need to obey certain distributions. The calculation of the Mann–Kendall test follows:

(2)$Z=\left\{\begin{array}{c}\frac{S-1}{\sqrt{s\left(S\right)}},S>0\\ 0,S=0\\ \frac{S+1}{\sqrt{s\left(S\right)}},S<0\end{array}\right.,$

where, $S=\sum _{j=1}^{n-1}\sum _{i=j-1}^n\mathrm{sgn}\left(LS{T}_j-LS{T}_i\right)$, $s\left(S\right)=\frac{n\left(n-1\right)\left(2n+5\right)}{18}$, and

(3)$\mathrm{sgn}\left(LS{T}_j-LS{T}_i\right)=\left\{\begin{array}{ccc}1& ,\hfill & LS{T}_j-LS{T}_i>0\hfill \\ 0& ,\hfill & LS{T}_j-LS{T}_i=0\hfill \\ -1& ,\hfill & LS{T}_j-LS{T}_i<0\hfill \end{array}\right.,$

In the equation, $i$ and $j$ indicate different months, $n$ denotes the length of the time series, $\mathrm{sgn}$ is the sign function, and the statistic $Z$ has a range of (−∞, +∞). At a given significance level α, when $\left|Z\right|>u_{1-\alpha /2}$, indicates a significant change in the time series at the $\mathsf{\alpha }$ level. In this study, we choose $\alpha =0.05$, meaning that we measure the significance of the LST trend from 2001 to 2020 on a pixel scale at a confidence level of 95%. A significant trend is defined with $\left|Z\right|>1.96$, otherwise a slight trend is obtained.

2.4. BFAST Algorithm

To detect and characterize changes in LST variations, the BFAST algorithm was applied to analyze the LST series. The BFAST algorithm was initially used to capture the dynamic changes of vegetation. This algorithm iteratively decomposes the original data series into the trend, seasonal, and residual components, which can detect abrupt changes in the trend and seasonal components. The process of implementing the BFAST algorithm is as follows:

• (1). An additive model is used to decompose the original time series into a trend component, a seasonal component, and a residual component. The algorithm is formulated as follows:

  (4)$Y_t=T_t+S_t+e_t,\hspace{1em}t=1,\dots ,n,$

  where $Y_t$ is the observed value at time $t$, $T_t$ is the trend component, $S_t$ is the seasonal component, and $e_t$ is the residual component.

• (2). A segmented linear fit is used to fit the trend component $T_t$. For each trend segment $t_i^{*}<t\le t_i+1^{*}$ after defining $T_0=0$, the linear model algorithm is as follows:

  (5)$T_t={\alpha }_i+{\beta }_{i}t$

In the above equation, $i=1\dots m$, $m$ is the number of mutation points, ${\alpha }_i$ represents the intercept of the linear model after segment fitting, and ${\beta }_i$ represents its slope.

• (3). For the seasonal component $S_t$, the harmonic model is fitted due to the obvious periodic variation of LST. For each trend segment $t_i^{*}<t\le t_i+1^{*}$, after defining $T_0=0$, the harmonic model can be expressed as:

  (6)$S_t=\sum _{k=1}^k{a}_{j,k}\sin \left(\frac{2\pi kt}{f}+{\delta }_{j,k}\right)$

In the above equation, $j$ is the position of the breakpoint $\left(j=1,\dots ,q\right)$. $q$ is the number of breakpoints in the seasonal component, $k$ is the number of harmonic terms in the harmonic model, and the amplitude $a_{j,k}$ and phase ${\delta }_{j,k}$ need to be obtained by a linear regression model. $f$ is the known observed frequency in the input time-series image.

More specifically, we use the BFAST01 algorithm, which is a variant of the BFAST algorithm. BFAST01 is implemented following the same procedure as BFAST, which detects the most significant breakpoints in long time series [54]. The algorithm not only detects abrupt changes in the LST time series but also distinguishes different trend change types [55]. The formula argument of BFAST01 was set to “response ~ trend + harmon”, because we used monthly time series with seasonality [55], and the h value is set to 1/9 [33,53,54]. The trend test was performed using four structural tests: Bayesian Information Criterion (BIC), Lagrange Multiplier (SupLM), Least Squares Moving Summation (OLS-MOSUM), and F-test (supF). A breakpoint is detected if any single test identifies a significant break (p < 0.05) [53,56]. In addition, LST pixels where the change in trend was not significant either before or after the breakpoint are defined as insignificant trends. In this study, we classified the types of LST dynamic trends detected by BFAST01 into six types. Monotonic changes indicate that the LST trend has always maintained the same trend. Interruption changes suggest an abrupt change in the opposite direction of the trend but with the same overall direction of the trend change. Reversal changes mean a shift in the direction of trend change during the study period. See Table 2 for more details.

2.5. Landscape Pattern Analysis

Landscape pattern is the spatial arrangement of different patterns among a series of landscape elements of different sizes and shapes, including the number, type and spatial distribution and configuration of landscape components [57].

By calculating the landscape indices, the basic spatial characteristics of landscape patterns can be understood. In this study, the grid method of Fragstats 4.2 [58] was used to calculate the overall Landscape pattern based on the types derived from the BFAST01 trend decomposition, while the Class metrics are calculated directly using Fragstats 4.2.

More specifically, at the overall landscape metrics, the four landscape indices SHEI, SHDI, SPLIT and CONTAG were selected and the overall landscape distribution pattern was spatially mapped using the grid method. The selected landscape indices and corresponding landscape ecological implications are shown in Table 3.

3. Results and Discussion

3.1. Linear LST Trends Based on Theil–Sen Median Trend Analysis and the Mann–Kendall

The long-term variation of LST on each pixel can be effectively traced using the Theil–Sen Median trend analysis and the Mann–Kendall test. The trend analysis shows that the range of LST variation in the YRD varied between −0.019 °C/month and 0.046 °C/month. According to the meaning of LST in reality and previous related studies [59], the types of LST trend changes are classified into four types (Table 4).

We conducted statistics on the monthly daytime LST in Shanghai, Hangzhou, Nanjing and Hefei from 2001 to 2020 and found that the LST in all four cities showed a slowly increasing trend (Figure 2). To gain further insight into the overall temperature change in the YRD region. Figure 3 illustrates the spatial distribution of the LST trend from 2001 to 2020 in the YRD. Significant warming (1.83%) occurred in the Suzhou–Wuxi–Changzhou urban agglomeration and the Hefei urban agglomeration, which are highly urbanized areas. Most of these areas have experienced rapid urbanization over the past decades, with urban expansion, population growth, and energy consumption recognized as the main causes of rapid warming [60,61,62,63]. Non-significant cooling warming (71.46%) mainly occurred in the northern YRD with large areas of cropland. The areas characterized by non-significant cooling (26.70%) were mainly concentrated in southwest Zhejiang and western Anhui. These areas are generally higher in elevation and the ground cover is dominated by vegetation. This result suggests that vegetation probably has a cooling effect on the LST [62]. Recent numerous studies also proposed that there is a warming trend for the LST over most of the YRD region [44,47,64].

3.2. LST Variations Based on BFAST01 Decomposition

3.2.1. LST Trends Based on BFAST01 Decomposition

The BFAST01 algorithm was used to decompose the variations of LST in the YRD region over the period 2001–2020. The different trend types obtained from BFAST01 trend decomposition can be used to investigate abrupt changes in trend direction and rate during LST variations, to monitor more detailed information about the LST change process and to identify different patterns of regional warming. The spatial distribution of the type derived by BFAST01 trend decomposition is shown in Figure 4. The most remarkable finding is that LST is highly dynamic (although the warming is very slow), with 63.4 % of pixels having a statistically significant trend. Non-significant trends accounted for 36.6% of all trend types. The type derived by BFAST01 trend decomposition is further divided into warming trends (monotonic increases, interruption increase, reversal increase) and cooling trends (monotonic decreases, interruption decrease, reversal decrease) based on the overall temperature change. Monotonic increases (27.3%), reversal increase (19.3%), and interruption increase (10.64%) constitute the majority of significant trends in the warming trend. As a comparison, the cooling trend (6.16%) accounts for little. In the spatial distribution, monotonic increases are concentrated in the YRD estuary urban agglomeration and central Anhui Province, which have higher levels of urbanization. Reversal increase, suggesting a reversal in the direction of the trend change, is mainly found in the northern YRD, where the land cover type is predominantly cropland. Interruption increase, reversal increase and interruption decrease mosaics are distributed in central Zhejiang and southern Anhui, indicating that external disturbance events are more frequent in these areas.

Compared to a linear trend change, BFAST01 decomposition can reveal underlying changes that are masked in linear trends and provide a more comprehensive view of the dynamic changes of LST [55,65,66,67]. In most of the warming/cooling regions, the BFAST01 method monitored more types of trend changes. Both the linear trend and the BFAST01 decomposition show a warming trend in the YRD estuary urban agglomeration and in the northern region, consisting mainly of monotonic increases, and reversal increases. The trend types of non-significant types, interrupted decreases and reversed decreases are musically distributed in the southern part of the YRD. However, the linear trend is dominated here by regions of slight cooling. To further demonstrate the differences in linear trend/BFAST01 trend decomposition, two trend-type conversion tables were calculated (Table 5). There are significant differences in LCTs of the YRD region. The southern part is mainly mountainous and hilly, with extensive woodland, while the eastern part is dominated by dense urban agglomerations and the northern part by large areas of agricultural land. The LCT plays a key role in the dynamic evolution of the LST, more details of which are given in Section 3.3.

3.2.2. Landscape Pattern Analysis

Landscape pattern is a high level of condensation of information on spatial distribution [57,68]. To further reveal the distribution pattern of LST dynamic changes in the YRD region, for the type derived by BFAST01 trend decomposition, we have further analyzed its landscape distribution.

Regarding the overall landscape metrics, Figure 5 shows the landscape pattern map of the LST trend types in the YRD after trend decomposition. SHEI and SHDI are relatively similar in their spatial distribution patterns, with higher values in central Zhejiang Province, north-central Anhui Province and northern Jiangsu Province. This indicates that the LST trend types in these regions are very diverse, but the dominance of the trend types is not high. The value of SPLIT reaches its highest value in the south-central region of Zhejiang, suggesting that the fragmentation of LST trend types in this region is very high, with a mosaic distribution of different trend types. CONTAG values are higher in southern Anhui, the YRD urban agglomeration, and northern Anhui. The LST trend types in these regions are generally relatively homogeneous, with strong connectivity of trend types. These indicators suggest that external disturbance activity on LST is somewhat spatially concentrated but with extensive spatial heterogeneity.

For the different trend types, Table 6 shows the landscape index values for the different trend types at class metrics. the monotonic increasing and non-significant trend patches are very similar in all indicators, with the AREA_MN and the LPI being significantly higher than the other trend types. The two trend types, both of which have relatively large and highly connected patches, are predominantly located in the east and southwest part of the YRD. This is mainly related to the highly centralized urban agglomerations in the east area and the forests in the southwest part. The interrupted increasing type has the highest number of patches, but the smaller AREA_MN and lower degree of AI indicate a very dispersed distribution of the trend type. The lower PD for all indices further suggests a mosaic distribution between patches of each trend type. Overall, the warming trend type is higher than the non-significant trend and the cooling trend in terms of the mean patch area, maximum patch index and degree of aggregation, respectively, indicating that the overall warming effect is dominant in the YRD, followed by the non-significant trend and finally the cooling trend.

3.2.3. Breakpoint Strength, Occurrence Times and Spatial Distribution

Anomalous increases or decreases in LST during long time series temperature changes often reflect strong external disturbances to the local environment. To detect temperature anomalies during LST variations, further analysis was performed on the timing and intensity of temperature breakpoints (breakpoints in monotonic trends were also counted) during the type derived by BFAST01 trend decomposition in the LST. Figure 6 illustrates the spatial distribution of breakpoint strength. Breakpoints are widely distributed in the YRD, concentrated in the northwestern and southern regions. The magnitude of the temperature change at the breakpoints tends to reflect the strength of the region’s exposure to perturbations from the external environment. The LST breakpoint intensity varies between −9 °C and 4.1 °C but does not exceed 2 °C in most areas (91.54%). With respect to the rate of the LST trend (Section 3.1, it takes about 3.62 to 8.77 years to compensate for the temperature decrease or increase at a breakpoint. This finding reveals the range of variation in the intensity of LST breakpoints and emphasizes the important role of the abrupt LST change in local temperature assessments. As a comparison, there are few areas with breakpoint strength above 2 °C, with a sporadic distribution. Further investigation was conducted into the timing of breakpoints in the type derived by BFAST01 trend decomposition (breakpoints in monotonic trends were also counted). The breakpoint occurrence times are distributed almost every year (Figure 7), by further aggregating the breakpoint occurrence times, one can find that the breakpoint occurrence times are highly concentrated in the periods 2010–2011 and 2012–2013. More specifically, the breakpoints in 2010–2011 were more concentrated in Anhui Province and northern Jiangsu Province. In contrast, the breakpoints in 2012–2013 are more concentrated in central Zhejiang.

3.3. Attributions of LST Trends

As a key parameter of the physical and chemical properties of the surface, the variability of LST is closely related to the local environment. In previous studies, LCT and NDVI have been reported as the most important influences of LST, and they play an important role in long-term changes in LST [69,70,71]. Therefore, we have analyzed the influences of these two factors on LST breakpoints and trend types.

3.3.1. NDVI, LCT and the LST Breakpoints

We further analyzed the frequency of breakpoints on different NDVI levels and LCTs (Figure 8). For NDVI, the occurrence of breakpoints is highly concentrated in the range of NDVI of 0.5–0.75. The occurrence of breakpoints is closely related to NDVI, with higher NDVI corresponding to a higher probability of breakpoints occurrence. Higher vegetation cover is more likely to induce dynamic changes in LST because human activities, pest and disease hazards, and hill fires are more likely to erupt in these areas [33,72,73]. It has been confirmed that the higher variation of vegetation variations in the southern and northern parts of the YRD region [74]. This is thought to be a potential cause of LST breakpoints. Breakpoints statistics for LUCC show that the distribution of breakpoints is highly concentrated in cropland (50.4%), followed by grassland (27.7%), woodland (10.7%), built-up areas (7.1%) and others (4.1%). Processes of LCT change, such as rapid urbanization, afforestation and reclamation of wasteland, can have a greater impact on the nature of the land surface, which changes the biophysical state of the land and causes feedback in the LST. Numerous studies have confirmed that the YRD region has experienced a large reduction in arable land over the past decades, which plays a key role in the occurrence of LST breakpoints [75,76].

3.3.2. NDVI, LCT and the Type Derived by BFAST01 Trend Decomposition

Further analysis of the types derived from the BFAST01 trend decomposition at different NDVI levels and LCTs (Figure 9). The statistics on the type of temperature change can help to understand the details of LST change on different LCTs. Overall, the proportion of type derived by BFAST01 trend decomposition was non-significant change (36.64%), monotonic increase (27.27%), reversal increase (19.30%), interrupted increase (10.58%), interrupted decrease (5.95%), monotonic decrease (0.21%) and reversal decrease (0.05%). The warming-type trend dominates the temperature change. However, there are significant differences in type derived by BFAST01 trend decomposition at different LCTs and NDVI levels. Figure 9 shows the frequency occurrence of the LST trend types decomposed by BFAST01 at different LCTs. The built-up area is dominated by monotonic increases (54.56%) and interrupted increases (19.53%). Monotonic increases and interrupted increases mean that there are no negative disturbances to LST during the temperature rise. This is because increases in population and emissions are continuous and slow. While external disturbances such as urban expansion, changes in the subsurface are often more rapid and intense, to which temperature changes respond more strongly. The main types of trend changes on cropland are monotonic increases (33.61%), decreases to increases (28.85%) and non-significant changes (26.05%). Compared to other LCTs, the frequency in the type of decreases to increases is the highest in cropland. Numerous studies have shown that the main target of urban expansion in the YRD is cropland [77,78]. The replacement of natural vegetation by urban buildings in urban expansions causes a significant change in temperature variations in terms of trend direction and rate. The non-significant change was the dominant type over grasslands and woodlands. As a comparison, for NDVI, the proportion of non-significant change generally increases as the NDVI value gets higher. In the background of global warming, regions with higher NDVI values show non-significant responses to warming, and monotonic increases are predominant in the range of NDVI less than 0.5.

3.4. The Inconsistent Warming of Different LCTs

To reveal the effect of different LCTs on regional warming, we compared the difference between the average monthly 2001 LST and the average monthly 2020 LST in the YRD region and further analyzed the LST changes in woodland, grassland, cropland, built-up area for warming-type trends (Figure 10). There were significant differences in the mean LST changes under different LCTs, with the largest changes occurring in built-up areas (3.14 °C), followed by grassland (2.01 °C), cropland (1.94 °C) and woodland (0.91 °C). There is also some variation in the warming contribution of the different warming-type trends within LCTs. For built-up areas and cropland, the highest warming contribution of the interrupted warming type was 3.46 °C and 2.09 °C, respectively. This suggests that external disturbances such as urban expansion and subsurface changes have a greater impact on the warming in built-up areas and cropland [79]. However, the warming contribution of monotonic warming is highest in grassland and woodland. Uneven warming is a phenomenon that cannot be ignored in global warming and has been widely reported, especially in urban areas. Numerous studies have confirmed that warming in urban contexts is significantly higher than in rural contexts, with urban expansion contributing to the major warming in China’s cities [5]. Regional warming due to rapid urban warming is the main cause of global warming. Among the different warming trends, monotonic warming and interrupted warming are more pronounced and contribute to the main warming during temperature change.

4. Conclusions

The analysis of the Spatiotemporal variations of the LST has been widely analyzed and discussed, but the variability of the LST is often complex and susceptible to external disturbances, especially at long timescale. In this study, In this study, we use the BFAST01 algorithm to decompose LST trends and monitor break-points, compare the spatial and temporal patterns of LST variability in the YRD region from different perspectives, and further analyze them using a landscape pattern approach for non-linear patterns of LST variability. Our study not only deepens the understanding of the spatiotemporal variability of LST in the region but also highlights that LST breakpoints caused by external disturbances cannot be ignored in temperature assessment. The main conclusions obtained from this paper are as follows:

• (1). The linear rate of change of LST in the YRD ranged from −0.019 °C/month to 0.046 °C/month, with a more pronounced warming trend in the north and near urban agglomerations. However, within the warming trend, it is mainly composed of monotonic increases (27.3%), reversal increases (19.3%) and interruption increases (10.64%). The landscape index shows a strong aggregation of the type derived by BFAST01 trend decomposition, but low connectivity and high spatial heterogeneity. Monotonic increases and non-significant trends are more dominant.

• (2). The breakpoints are widely distributed in the YRD but are more concentrated in the southern and northern regions. The intensity of the breakpoints is mostly within 2 °C, with reference to the linear trend rate of change, which typically takes 3.62–8.77 years to offset an abrupt change. The breakpoints are highly concentrated in the period 2010–2013, suggesting stronger external disturbances in this period. Breakpoints occurred more frequently over cropland and the NDVI range of 0.5–0.7, indicating more disturbances over these areas.

• (3). The types of LST trends varied considerably for different NDVI levels and LCTs. In general, the proportion of non-significant trends generally increases gradually as the NDVI level increases. Within a global warming background, this suggests a suppressive effect of vegetation on LST warming. The warming in the built-up area is significantly higher than in the other LCTS, with monotonic warming and interrupted warming contributing more to warming.

Author Contributions

Data curation, Q.Z., M.W.; Formal analysis, Q.Z., G.Y. and W.S.; Funding acquisition, T.F. and H.L.; Investigation, T.F. and H.L.; Methodology, Q.Z. and T.F.; Resources, W.S.; Supervision, G.Y.; Validation, G.Y. and W.S.; Writing—original draft, Q.Z.; Writing—review and editing, T.F., M.W. and H.L. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

In this study, the data were acquired and processed via GEE. The raw data used are publicly available from the sources provided in following links. The processed data generated from this study are available upon request to the corresponding author. MOD11A2 data: https://lpdaac.usgs.gov/products/mod11a2v006/ (accessed on 7 July 2022). MOD13A2 data: https://lpdaac.usgs.gov/products/mod13a2v006/ (accessed on 12 August 2022). MCD12Q1 data: https://lpdaac.usgs.gov/products/mcd12q1v006/ (accessed on 12 August 2022).

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Land use map of the study area.

Figure 2. Temporal variations of monthly daytime LST in Shanghai, Hangzhou, Nanjing and Hefei in the YRD from 2001 to 2020.

Figure 3. Distributions of the LST linear trend from the Theil–Sen Median trend analysis and the Mann-Kendall test in the YRD during 2001–2020. Islands are not included in the analysis.

Figure 4. Trend types of the LST derived from the BFAST01 trend decomposition in the YRD during 2001–2020. Islands are not included in the analysis.

Figure 5. Distributions of landscape pattern of LST trend types in the YRD from the trend decomposition.

Figure 6. Distributions of the LST breakpoint intensity over the YRD region during 2001–2020. Islands are not included in the analysis.

Figure 7. Distributions of the LST breakpoint time over the YRD region during 2001–2020. Islands are not included in the analysis.

Figure 8. Percentages of different NDVI levels (a) and LCTs (b) in the LST breakpoints during 2001–2020 in the YRD region.

Figure 9. Percentages of the LST trend types derived from the BFAST01 trend decomposition in different LCTs and NDVI levels during 2001–2020 in the YRD region.

Figure 10. LST changes of various trend types derived from the BFAST01 trend decomposition over different LCTs.

Table 1

MCD12Q1 reclassification scheme.

Reclassification	MCD12Q1 IGBP Classification
Woodland	Evergreen Needleleaf Forests, Evergreen Broadleaf Forests, Deciduous Needleleaf Forests, Deciduous Broadleaf Forests, Mixed Forests
Grassland	Closed Shrublands, Open Shrublands, Woody Savannas, Savannas, Grasslands
Cropland	Croplands, Cropland/Natural Vegetation Mosaics
Built-up area	Urban and Built-up Lands
Others	Barren, Permanent Wetlands, Permanent Snow and Ice, Water Bodies

Table 2

The different types of LST change detected by BFAST01.

Type of Change	Example	Description
Monotonic increase	[Image omitted. Please see PDF.]	A significant increase with one significant break or none
Monotonic decrease	[Image omitted. Please see PDF.]	A significant decrease with one significant break or none
Interruption increase	[Image omitted. Please see PDF.]	An increasing trend with a negative breakpoint
Interruption decrease	[Image omitted. Please see PDF.]	A decreasing trend with a positive breakpoint
Reversal decrease	[Image omitted. Please see PDF.]	An increasing trend disturbed by a breakpoint and followed by a decrease trend
Reversal increase	[Image omitted. Please see PDF.]	A decreasing trend disturbed by a breakpoint and followed by an increasing trend

Table 3

Landscape indices and the meaning in landscape ecology.

	Landscape Indices	Value	Meaning
Class metrics	NP	≥1	The number of patches in the landscape
	AREA_MN	≥0	Average area of patches
	LPI	0 ≤ LPI ≤ 100	The percentage of the landscape comprised by The largest patch
	PD	≥0	Patch density
	LSI	≥0	Complexity of patch shape
	AI	0 ≤ AI ≤ 100	Degree of aggregation of patches
Landscape metrics	NP	≥1	The number of patches in the landscape
	SPLIT	0 ≤ SPLIT ≤ NP2	Higher values indicate greater landscape fragmentation
	CONTAG	0 ≤ CONTAG ≤ 100	Higher values indicate greater landscape connectivity
	SHDI	≥0	Higher values indicate more landscape types
	SHEI	0 ≤ SHEI ≤ 1	Higher values indicate lower landscape dominance

Table 4

LST linear trend variation statistics.

$\mathit{S}\mathit{l}\mathit{o}\mathit{p}{\mathit{e}}_{\mathit{L}\mathit{S}\mathit{T}}$	Z	LST Trend	Area Percentage (%)
−0.019–0	≥1.96	Significant warming	1.83%
−0.019–0	−1.96–1.96	Non-significant warming	71.46%
0–0.046	≥1.96	Significant cooling	0.01%
0–0.046	−1.96–1.96	Non-significant cooling	26.70%

Table 5

Linear and the type derived by BFAST01 trend decomposition conversion table.

		Non-Linear Trends
		Monotonic Increases	Monotonic Decreases	Interruption Increase	Interruption Decrease	Reversal Increase	ReversalDecrease	Non-Significant Change
Linear trends	Significant cooling	/	34.29	2.86	/	22.86	/	40
	Non-significant cooling	1.52	0.75	6.97	13.54	17.10	0.01	60.10
	Significant warming	60.73	/	28.01	0.12	4.37	0.24	6.53
	Non-significant warming	36.09	0.01	11.54	3.23	20.49	0.06	28.58

Table 6

Landscape index values for the different trend types at class metrics.

Class Metrics	Monotonic Increases	Monotonic Decreases	Interruption Increase	Interruption Decrease	Reversal Increase	ReversalDecrease	Non-Significant
NP	2357	170	3969	1523	3445	95	3886
AREA_MN	3993.7537	419.1853	924.5099	1335.1613	1932.0885	182.5037	3247.9073
LPI	16.7232	0.0209	0.5316	1.0246	10.7034	0.0045	23.9916
PD	0.0068	0.0005	0.0115	0.0044	0.01	0.0003	0.0113
LSI	62.8771	14.9683	77.0671	45.4685	63.9435	10.3548	71.7793
AI	82.6631	53.1167	65.8089	73.0886	79.0188	34.2404	82.8819

References

1. Diffenbaugh, N.S.; Burke, M. Global warming has increased global economic inequality. Proc. Natl. Acad. Sci. USA; 2019; 116, pp. 9808-9813. [DOI: https://dx.doi.org/10.1073/pnas.1816020116] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31010922]

2. King, A.D.; Karoly, D.J.; Henley, B.J. Australian climate extremes at 1.5 degrees C and 2 degrees C of global warming. Nat. Clim. Chang.; 2017; 12, 114031. [DOI: https://dx.doi.org/10.1038/NCLIMATE3296]

3. Hansen, J.; Ruedy, R.; Sato, M.; Lo, K. Global Surface Temperature Change. Rev. Geophys.; 2010; 48, [DOI: https://dx.doi.org/10.1029/2010RG000345]

4. Li, Z.L.; Tang, B.H.; Wu, H.; Ren, H.Z.; Yan, G.J.; Wan, Z.M.; Trigo, I.F.; Sobrino, J.A. Satellite-derived land surface temperature: Current status and perspectives. Remote Sens. Environ.; 2013; 131, pp. 14-37. [DOI: https://dx.doi.org/10.1016/j.rse.2012.12.008]

5. Liu, Z.; Zhan, W.; Bechtel, B.; Voogt, J.; Lai, J.; Chakraborty, T.; Wang, Z.-H.; Li, M.; Huang, F.; Lee, X. Surface warming in global cities is substantially more rapid than in rural background areas. Commun. Earth Environ.; 2022; 3, 219. [DOI: https://dx.doi.org/10.1038/s43247-022-00539-x]

6. Zhan, W.; Chen, Y.; Zhou, J.; Wang, J.; Liu, W.; Voogt, J.; Zhu, X.; Quan, J.; Li, J. Disaggregation of remotely sensed land surface temperature: Literature survey, taxonomy, issues, and caveats. Remote Sens. Environ.; 2013; 131, pp. 119-139. [DOI: https://dx.doi.org/10.1016/j.rse.2012.12.014]

7. Aguilar-Lome, J.; Espinoza-Villar, R.; Espinoza, J.-C.; Rojas-Acuña, J.; Willems, B.L.; Leyva-Molina, W.-M. Elevation-dependent warming of land surface temperatures in the Andes assessed using MODIS LST time series (2000–2017). Int. J. Appl. Earth Obs. Geoinf.; 2019; 77, pp. 119-128. [DOI: https://dx.doi.org/10.1016/j.jag.2018.12.013]

8. Zhou, D.; Xiao, J.; Frolking, S.; Liu, S.; Zhang, L.; Cui, Y.; Zhou, G. Croplands intensify regional and global warming according to satellite observations. Remote Sens. Environ.; 2021; 264, 112585. [DOI: https://dx.doi.org/10.1016/j.rse.2021.112585]

9. Camuffo, D.; Bertolin, C. The earliest temperature observations in the world: The Medici Network (1654–1670). Clim. Chang.; 2012; 111, pp. 335-363. [DOI: https://dx.doi.org/10.1007/s10584-011-0142-5]

10. Stahl, K.; Moore, R.D.; Floyer, J.A.; Asplin, M.G.; McKendry, I.G. Comparison of approaches for spatial interpolation of daily air temperature in a large region with complex topography and highly variable station density. Agric. For. Meteorol.; 2006; 139, pp. 224-236. [DOI: https://dx.doi.org/10.1016/j.agrformet.2006.07.004]

11. Eleftheriou, D.; Kiachidis, K.; Kalmintzis, G.; Kalea, A.; Bantasis, C.; Koumadoraki, P.; Spathara, M.E.; Tsolaki, A.; Tzampazidou, M.I.; Gemitzi, A. Determination of annual and seasonal daytime and nighttime trends of MODIS LST over Greece—Climate change implications. Sci. Total Environ.; 2018; 616, pp. 937-947. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2017.10.226]

12. Laraby, K.G.; Schott, J.R. Uncertainty estimation method and Landsat 7 global validation for the Landsat surface temperature product. Remote Sens. Environ.; 2018; 216, pp. 472-481. [DOI: https://dx.doi.org/10.1016/j.rse.2018.06.026]

13. Malakar, N.K.; Hulley, G.C.; Hook, S.J.; Laraby, K.; Cook, M.; Schott, J.R. An Operational Land Surface Temperature Product for Landsat Thermal Data: Methodology and Validation. IEEE Trans. Geosci. Remote Sens.; 2018; 56, pp. 5717-5735. [DOI: https://dx.doi.org/10.1109/TGRS.2018.2824828]

14. Ding, Y.; Zhang, S.; Zhao, L.; Li, Z.; Kang, S. Global warming weakening the inherent stability of glaciers and permafrost. Sci. Bull.; 2019; 64, pp. 245-253. [DOI: https://dx.doi.org/10.1016/j.scib.2018.12.028]

15. Du, Q.; Zhang, M.; Wang, S.; Che, C.; Ma, R.; Ma, Z. Changes in air temperature over China in response to the recent global warming hiatus. J. Geogr. Sci.; 2019; 29, pp. 496-516. [DOI: https://dx.doi.org/10.1007/s11442-019-1612-3]

16. Ji, F.; Wu, Z.; Huang, J.; Chassignet, E.P. Evolution of land surface air temperature trend. Nat. Clim. Chang.; 2014; 4, pp. 462-466. [DOI: https://dx.doi.org/10.1038/nclimate2223]

17. Croitoru, A.-E.; Holobaca, I.-H.; Lazar, C.; Moldovan, F.; Imbroane, A. Air temperature trend and the impact on winter wheat phenology in Romania. Clim. Chang.; 2012; 111, pp. 393-410. [DOI: https://dx.doi.org/10.1007/s10584-011-0133-6]

18. Sonali, P.; Nagesh Kumar, D. Review of trend detection methods and their application to detect temperature changes in India. J. Hydrol.; 2013; 476, pp. 212-227. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2012.10.034]

19. Li, Y.; Yao, N.; Chau, H.W. Influences of removing linear and nonlinear trends from climatic variables on temporal variations of annual reference crop evapotranspiration in Xinjiang, China. Sci. Total Environ.; 2017; 592, pp. 680-692. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2017.02.196]

20. Wu, Z.; Huang, N.E.; Long, S.R.; Peng, C.-K. On the trend, detrending, and variability of nonlinear and nonstationary time series. Proc. Natl. Acad. Sci. USA; 2007; 104, pp. 14889-14894. [DOI: https://dx.doi.org/10.1073/pnas.0701020104]

21. Panwar, M.; Agarwal, A.; Devadas, V. Analyzing land surface temperature trends using non-parametric approach: A case of Delhi, India. Urban Clim.; 2018; 24, pp. 19-25. [DOI: https://dx.doi.org/10.1016/j.uclim.2018.01.003]

22. Yu, Y.; Duan, S.-B.; Li, Z.-L.; Chang, S.; Xing, Z.; Leng, P.; Gao, M. Interannual Spatiotemporal Variations of Land Surface Temperature in China From 2003 to 2018. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens.; 2021; 14, pp. 1783-1795. [DOI: https://dx.doi.org/10.1109/JSTARS.2020.3048823]

23. Fang, S.; Mao, K.; Xia, X.; Wang, P.; Shi, J.; Bateni, S.M.; Xu, T.; Cao, M.; Heggy, E.; Qin, Z. Dataset of daily near-surface air temperature in China from 1979 to 2018. Earth Syst. Sci. Data; 2022; 14, pp. 1413-1432. [DOI: https://dx.doi.org/10.5194/essd-14-1413-2022]

24. Zhao, P.; Jones, P.; Cao, L.; Yan, Z.; Zha, S.; Zhu, Y.; Yu, Y.; Tang, G. Trend of Surface Air Temperature in Eastern China and Associated Large-Scale Climate Variability over the Last 100 Years. J. Clim.; 2014; 27, pp. 4693-4703. [DOI: https://dx.doi.org/10.1175/JCLI-D-13-00397.1]

25. El Kenawy, A.; López-Moreno, J.I.; Vicente-Serrano, S.M. Trend and variability of surface air temperature in northeastern Spain (1920–2006): Linkage to atmospheric circulation. Atmos. Res.; 2012; 106, pp. 159-180. [DOI: https://dx.doi.org/10.1016/j.atmosres.2011.12.006]

26. Kennedy, R.E.; Yang, Z.; Cohen, W.B. Detecting trends in forest disturbance and recovery using yearly Landsat time series: 1. LandTrendr—Temporal segmentation algorithms. Remote Sens. Environ.; 2010; 114, pp. 2897-2910. [DOI: https://dx.doi.org/10.1016/j.rse.2010.07.008]

27. Jamali, S.; Jönsson, P.; Eklundh, L.; Ardö, J.; Seaquist, J. Detecting changes in vegetation trends using time series segmentation. Remote Sens. Environ.; 2015; 156, pp. 182-195. [DOI: https://dx.doi.org/10.1016/j.rse.2014.09.010]

28. Verbesselt, J.; Hyndman, R.; Newnham, G.; Culvenor, D. Detecting trend and seasonal changes in satellite image time series. Remote Sens. Environ.; 2010; 114, pp. 106-115. [DOI: https://dx.doi.org/10.1016/j.rse.2009.08.014]

29. Jacquin, A.; Sheeren, D.; Lacombe, J.-P. Vegetation cover degradation assessment in Madagascar savanna based on trend analysis of MODIS NDVI time series. Int. J. Appl. Earth Obs. Geoinf.; 2010; 12, pp. S3-S10. [DOI: https://dx.doi.org/10.1016/j.jag.2009.11.004]

30. Bassingthwaighte, J.B.; Raymond, G.M. Evaluating rescaled range analysis for time series. Ann. Biomed. Eng.; 1994; 22, pp. 432-444. [DOI: https://dx.doi.org/10.1007/BF02368250]

31. Iida, T.; Saitoh, S.-I. Temporal and spatial variability of chlorophyll concentrations in the Bering Sea using empirical orthogonal function (EOF) analysis of remote sensing data. Deep. Sea Res. Part II Top. Stud. Oceanogr.; 2007; 54, pp. 2657-2671. [DOI: https://dx.doi.org/10.1016/j.dsr2.2007.07.031]

32. Hawinkel, P.; Swinnen, E.; Lhermitte, S.; Verbist, B.; Van Orshoven, J.; Muys, B. A time series processing tool to extract climate-driven interannual vegetation dynamics using Ensemble Empirical Mode Decomposition (EEMD). Remote Sens. Environ.; 2015; 169, pp. 375-389. [DOI: https://dx.doi.org/10.1016/j.rse.2015.08.024]

33. Fang, X.; Zhu, Q.; Ren, L.; Chen, H.; Wang, K.; Peng, C. Large-scale detection of vegetation dynamics and their potential drivers using MODIS images and BFAST: A case study in Quebec, Canada. Remote Sens. Environ.; 2018; 206, pp. 391-402. [DOI: https://dx.doi.org/10.1016/j.rse.2017.11.017]

34. Shen, X.; An, R.; Feng, L.; Ye, N.; Zhu, L.; Li, M. Vegetation changes in the Three-River Headwaters Region of the Tibetan Plateau of China. Ecol. Indic.; 2018; 93, pp. 804-812. [DOI: https://dx.doi.org/10.1016/j.ecolind.2018.05.065]

35. Xie, M.; Liao, J.B.; Wang, T.J.; Zhu, K.G.; Zhuang, B.L.; Han, Y.; Li, M.M.; Li, S. Modeling of the anthropogenic heat flux and its effect on regional meteorology and air quality over the Yangtze River Delta region, China. Atmos. Chem. Phys.; 2016; 16, pp. 6071-6089. [DOI: https://dx.doi.org/10.5194/acp-16-6071-2016]

36. Yang, X.C.; Leung, L.R.; Zhao, N.Z.; Zhao, C.; Qian, Y.; Hu, K.J.; Liu, X.P.; Chen, B.D. Contribution of urbanization to the increase of extreme heat events in an urban agglomeration in east China. Geophys. Res. Lett.; 2017; 44, pp. 6940-6950. [DOI: https://dx.doi.org/10.1002/2017GL074084]

37. Yao, R.; Hu, Y.Q.; Sun, P.; Bian, Y.J.; Liu, R.L.; Zhang, S.L. Effects of urbanization on heat waves based on the wet-bulb temperature in the Yangtze River Delta urban agglomeration, China. Urban Clim.; 2022; 41, 101067. [DOI: https://dx.doi.org/10.1016/j.uclim.2021.101067]

38. Qian, L.X.; Cui, H.S.; Jie, C. Impacts of land use and cover change on land surface temperature in the Zhujiang Delta. Pedosphere; 2006; 16, pp. 681-689. [DOI: https://dx.doi.org/10.1016/S1002-0160(06)60103-3]

39. Song, Z.J.; Li, R.H.; Qiu, R.Y.; Liu, S.Y.; Tan, C.; Li, Q.P.; Ge, W.; Han, X.J.; Tang, X.G.; Shi, W.Y. et al. Global Land Surface Temperature Influenced by Vegetation Cover and PM2.5 from 2001 to 2016. Remote Sens.; 2018; 10, 2034. [DOI: https://dx.doi.org/10.3390/rs10122034]

40. Yang, B.; Yang, X.C.; Leung, L.R.; Zhong, S.Y.; Qian, Y.; Zhao, C.; Chen, F.; Zhang, Y.C.; Qi, J.G. Modeling the Impacts of Urbanization on Summer Thermal Comfort: The Role of Urban Land Use and Anthropogenic Heat. J. Geophys. Res. Atmos.; 2019; 124, pp. 6681-6697. [DOI: https://dx.doi.org/10.1029/2018JD029829]

41. Lin, S.; Feng, J.; Wang, J.; Hu, Y. Modeling the contribution of long-term urbanization to temperature increase in three extensive urban agglomerations in China. J. Geophys. Res. Atmos.; 2016; 121, pp. 1683-1697. [DOI: https://dx.doi.org/10.1002/2015JD024227]

42. Luo, M.; Lau, N.-C. Increasing Heat Stress in Urban Areas of Eastern China: Acceleration by Urbanization. Geophys. Res. Lett.; 2018; 45, pp. 13060-13069. [DOI: https://dx.doi.org/10.1029/2018GL080306]

43. Peng, X.; She, Q.; Long, L.; Liu, M.; Xu, Q.; Zhang, J.; Xiang, W. Long-term trend in ground-based air temperature and its responses to atmospheric circulation and anthropogenic activity in the Yangtze River Delta, China. Atmos. Res.; 2017; 195, pp. 20-30. [DOI: https://dx.doi.org/10.1016/j.atmosres.2017.05.013]

44. Sang, Y.-F. Spatial and temporal variability of daily temperature in the Yangtze River Delta, China. Atmos. Res.; 2012; 112, pp. 12-24. [DOI: https://dx.doi.org/10.1016/j.atmosres.2012.04.006]

45. Vincent, L.A.; Wang, X.L.; Milewska, E.J.; Wan, H.; Yang, F.; Swail, V. A second generation of homogenized Canadian monthly surface air temperature for climate trend analysis. J. Geophys. Res. Atmos.; 2012; 117, [DOI: https://dx.doi.org/10.1029/2012JD017859]

46. Yan, Z.; Ding, Y.; Zhai, P.; Song, L.; Cao, L.; Li, Z. Re-Assessing Climatic Warming in China since 1900. J. Meteorol. Res.; 2020; 34, pp. 243-251. [DOI: https://dx.doi.org/10.1007/s13351-020-9839-6]

47. Zhao, B.; Mao, K.; Cai, Y.; Shi, J.; Li, Z.; Qin, Z.; Meng, X.; Shen, X.; Guo, Z. A combined Terra and Aqua MODIS land surface temperature and meteorological station data product for China from 2003 to 2017. Earth Syst. Sci. Data; 2020; 12, pp. 2555-2577. [DOI: https://dx.doi.org/10.5194/essd-12-2555-2020]

48. Li, L.; Zhang, Y.; Liu, Q.; Ding, M.; Mondal, P.P. Regional differences in shifts of temperature trends across China between 1980 and 2017. Int. J. Climatol.; 2019; 39, pp. 1157-1165. [DOI: https://dx.doi.org/10.1002/joc.5868]

49. Bai, X.; Zhang, L.; He, C.; Zhu, Y. Estimating Regional Soil Moisture Distribution Based on NDVI and Land Surface Temperature Time Series Data in the Upstream of the Heihe River Watershed, Northwest China. Remote Sens.; 2020; 12, 2414. [DOI: https://dx.doi.org/10.3390/rs12152414]

50. Han, G.; Xu, J. Land Surface Phenology and Land Surface Temperature Changes Along an Urban–Rural Gradient in Yangtze River Delta, China. Environ. Manag.; 2013; 52, pp. 234-249. [DOI: https://dx.doi.org/10.1007/s00267-013-0097-6]

51. Peng, J.; Ma, J.; Liu, Q.; Liu, Y.; Hu, Y.; Li, Y.; Yue, Y. Spatial-temporal change of land surface temperature across 285 cities in China: An urban-rural contrast perspective. Sci. Total Environ.; 2018; 635, pp. 487-497. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2018.04.105]

52. Deng, Y.; Wang, M.; Yousefpour, R.; Hanewinkel, M. Abiotic disturbances affect forest short-term vegetation cover and phenology in Southwest China. Ecol. Indic.; 2021; 124, 107393. [DOI: https://dx.doi.org/10.1016/j.ecolind.2021.107393]

53. Higginbottom, T.P.; Symeonakis, E. Identifying Ecosystem Function Shifts in Africa Using Breakpoint Analysis of Long-Term NDVI and RUE Data. Remote Sens.; 2020; 12, 1894. [DOI: https://dx.doi.org/10.3390/rs12111894]

54. De Jong, R.; Verbesselt, J.; Zeileis, A.; Schaepman, M.E. Shifts in Global Vegetation Activity Trends. Remote Sens.; 2013; 5, pp. 1117-1133. [DOI: https://dx.doi.org/10.3390/rs5031117]

55. Bernardino, P.N.; De Keersmaecker, W.; Fensholt, R.; Verbesselt, J.; Somers, B.; Horion, S.; Silva, T. Global-scale characterization of turning points in arid and semi-arid ecosystem functioning. Glob. Ecol. Biogeogr.; 2020; 29, pp. 1230-1245. [DOI: https://dx.doi.org/10.1111/geb.13099]

56. Zeileis, A. A Unified Approach to Structural Change Tests Based on ML Scores, F Statistics, and OLS Residuals. Econom. Rev.; 2005; 24, pp. 445-466. [DOI: https://dx.doi.org/10.1080/07474930500406053]

57. Riitters, K.H.; O’Neill, R.V.; Hunsaker, C.T.; Wickham, J.D.; Yankee, D.H.; Timmins, S.P.; Jones, K.B.; Jackson, B.L. A factor analysis of landscape pattern and structure metrics. Landsc. Ecol.; 1995; 10, pp. 23-39. [DOI: https://dx.doi.org/10.1007/BF00158551]

58. McGarigal, K.; Cushman, S.A.; Ene, E. FRAGSTATS v4: Spatial Pattern Analysis Program for Categorical and Continuous Maps. Computer Software Program Produced by the Authors at the University of Massachusetts, Amherst, MA, USA. 2012; Available online: http://www.umass.edu/landeco/research/fragstats/fragstats.html (accessed on 12 November 2022).

59. Jiang, W.; Yuan, L.; Wang, W.; Cao, R.; Zhang, Y.; Shen, W. Spatio-temporal analysis of vegetation variation in the Yellow River Basin. Ecol. Indic.; 2015; 51, pp. 117-126. [DOI: https://dx.doi.org/10.1016/j.ecolind.2014.07.031]

60. Bokaie, M.; Zarkesh, M.K.; Arasteh, P.D.; Hosseini, A. Assessment of Urban Heat Island based on the relationship between land surface temperature and Land Use/Land Cover in Tehran. Sustain. Cities Soc.; 2016; 23, pp. 94-104. [DOI: https://dx.doi.org/10.1016/j.scs.2016.03.009]

61. Dewan, A.; Kiselev, G.; Botje, D.; Mahmud, G.I.; Bhuian, M.H.; Hassan, Q.K. Surface urban heat island intensity in five major cities of Bangladesh: Patterns, drivers and trends. Sustain. Cities Soc.; 2021; 71, 102926. [DOI: https://dx.doi.org/10.1016/j.scs.2021.102926]

62. Peng, S.S.; Piao, S.L.; Zeng, Z.Z.; Ciais, P.; Zhou, L.M.; Li, L.Z.X.; Myneni, R.B.; Yin, Y.; Zeng, H. Afforestation in China cools local land surface temperature. Proc. Natl. Acad. Sci. USA; 2014; 111, pp. 2915-2919. [DOI: https://dx.doi.org/10.1073/pnas.1315126111] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24516135]

63. Tran, D.X.; Pla, F.; Latorre-Carmona, P.; Myint, S.W.; Gaetano, M.; Kieu, H.V. Characterizing the relationship between land use land cover change and land surface temperature. ISPRS J. Photogramm. Remote Sens.; 2017; 124, pp. 119-132. [DOI: https://dx.doi.org/10.1016/j.isprsjprs.2017.01.001]

64. Wang, J.; Yan, Z.; Quan, X.-W.; Feng, J. Urban warming in the 2013 summer heat wave in eastern China. Clim. Dyn.; 2016; 48, pp. 3015-3033. [DOI: https://dx.doi.org/10.1007/s00382-016-3248-7]

65. Ghorbanian, A.; Mohammadzadeh, A.; Jamali, S. Linear and Non-Linear Vegetation Trend Analysis throughout Iran Using Two Decades of MODIS NDVI Imagery. Remote Sens.; 2022; 14, 3683. [DOI: https://dx.doi.org/10.3390/rs14153683]

66. Holtvoeth, J.; Vogel, H.; Valsecchi, V.; Lindhorst, K.; Schouten, S.; Wagner, B.; Wolff, G.A. Linear and non-linear responses of vegetation and soils to glacial-interglacial climate change in a Mediterranean refuge. Sci. Rep.; 2017; 7, 8121. [DOI: https://dx.doi.org/10.1038/s41598-017-08101-y]

67. Jenkins, L.K.; Barry, T.; Bosse, K.R.; Currie, W.S.; Christensen, T.; Longan, S.; Shuchman, R.A.; Tanzer, D.; Taylor, J.J. Satellite-based decadal change assessments of pan-Arctic environments. Ambio; 2020; 49, pp. 820-832. [DOI: https://dx.doi.org/10.1007/s13280-019-01249-z]

68. Weng, Y.-C. Spatiotemporal changes of landscape pattern in response to urbanization. Landsc. Urban Plan.; 2007; 81, pp. 341-353. [DOI: https://dx.doi.org/10.1016/j.landurbplan.2007.01.009]

69. Marzban, F.; Sodoudi, S.; Preusker, R. The influence of land-cover type on the relationship between NDVI–LST and LST-Tair. Int. J. Remote Sens.; 2018; 39, pp. 1377-1398. [DOI: https://dx.doi.org/10.1080/01431161.2017.1402386]

70. Du, H.; Wang, D.; Wang, Y.; Zhao, X.; Qin, F.; Jiang, H.; Cai, Y. Influences of land cover types, meteorological conditions, anthropogenic heat and urban area on surface urban heat island in the Yangtze River Delta Urban Agglomeration. Sci. Total Environ.; 2016; 571, pp. 461-470. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2016.07.012]

71. Sun, D.; Kafatos, M. Note on the NDVI-LST relationship and the use of temperature-related drought indices over North America. Geophys. Res. Lett.; 2007; 34, [DOI: https://dx.doi.org/10.1029/2007GL031485]

72. Smith, V.; Portillo-Quintero, C.; Sanchez-Azofeifa, A.; Hernandez-Stefanoni, J.L. Assessing the accuracy of detected breaks in Landsat time series as predictors of small scale deforestation in tropical dry forests of Mexico and Costa Rica. Remote Sens. Environ.; 2019; 221, pp. 707-721. [DOI: https://dx.doi.org/10.1016/j.rse.2018.12.020]

73. Yan, J.; He, H.; Wang, L.; Zhang, H.; Liang, D.; Zhang, J. Inter-Comparison of Four Models for Detecting Forest Fire Disturbance from MOD13A2 Time Series. Remote Sens.; 2022; 14, 1446. [DOI: https://dx.doi.org/10.3390/rs14061446]

74. Yao, R.; Wang, L.; Huang, X.; Chen, X.; Liu, Z. Increased spatial heterogeneity in vegetation greenness due to vegetation greening in mainland China. Ecol. Indic.; 2019; 99, pp. 240-250. [DOI: https://dx.doi.org/10.1016/j.ecolind.2018.12.039]

75. Ning, J.; Liu, J.; Kuang, W.; Xu, X.; Zhang, S.; Yan, C.; Li, R.; Wu, S.; Hu, Y.; Du, G. et al. Spatiotemporal patterns and characteristics of land-use change in China during 2010–2015. J. Geogr. Sci.; 2018; 28, pp. 547-562. [DOI: https://dx.doi.org/10.1007/s11442-018-1490-0]

76. Wang, X.; Xin, L.; Tan, M.; Li, X.; Wang, J. Impact of spatiotemporal change of cultivated land on food-water relations in China during 1990-2015. Sci. Total Environ.; 2020; 716, 137119. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2020.137119]

77. He, C.; Liu, Z.; Xu, M.; Ma, Q.; Dou, Y. Urban expansion brought stress to food security in China: Evidence from decreased cropland net primary productivity. Sci. Total Environ.; 2017; 576, pp. 660-670. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2016.10.107]

78. Li, J.; Wang, Z.; Lai, C.; Wu, X.; Zeng, Z.; Chen, X.; Lian, Y. Response of net primary production to land use and land cover change in mainland China since the late 1980s. Sci. Total Environ.; 2018; 639, pp. 237-247. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2018.05.155]

79. Sun, Y.; Zhang, X.; Ren, G.; Zwiers, F.W.; Hu, T. Contribution of urbanization to warming in China. Nat. Clim. Chang.; 2016; 6, pp. 706-709. [DOI: https://dx.doi.org/10.1038/nclimate2956]