Impacts of terrain on land surface phenology derived from Harmonized Landsat 8 and Sentinel-2 in the Tianshan Mountains, China

ABSTRACT

Knowledge of terrain impacts on land surface phenology (LSP) is crucial for understanding the responses of mountainous ecosystems to environmental changes. While the effects of terrain factors on LSP spatial patterns have been observed to vary across regions due to their different climate and terrain conditions, the specific effects across different elevations are still largely unclear, especially in regions with diverse hydrothermal conditions, such as the Tianshan Mountains located in arid and semiarid region. Here, we investigated the spatial relationships between LSP metrics and terrain factors (i.e. elevation and aspect) in the Tianshan Mountains in Xinjiang, China. Our analysis utilized surface reflectance at a 30 m spatial resolution from the Harmonized Landsat 8 and Sentinel-2 dataset for 2021 and 2022. We focused on two LSP metrics, vegetation greenup (GU20) and maturity (GU90), which were estimated using 20% and 90% thresholds of seasonal amplitude of the enhanced vegetation index (EVI) time series, respectively. We modeled the spatial relationships using ordinary least square (OLS) linear regression for the entire study region and then applied geographically weighted regression (GWR) with a 2.5 km bandwidth to explore local variations. Our results suggest that, at a large scale, elevation played a primary role in controlling the spatial variations in both LSP metrics, overshadowing the role of aspect. However, when examined at a local scale using GWR, aspect emerged as an important factor, with south-facing aspects associated with earlier dates of GU20 and GU90 for most regions. Furthermore, we found that the influences of terrain on the LSP metrics varied across elevations. The explanatory power of terrain was stronger at middle elevations (approximately 2000–3000 m) than at lower (<2000 m) and higher (>3000 m) elevations. In addition, the sensitivities of the LSP metrics to elevation and aspect demonstrated varying patterns at elevations above 2000 m. Our findings highlight the diverse environmental controls on LSP across elevations, with a particular emphasis on the phenological sensitivities to aspect-induced local climatic differences.

KEYWORDS:

• Vegetation phenology
• drylands
• Harmonized Landsat and Sentinel-2
• elevation
• aspect
• Tianshan Mountains

1 Introduction

Mountainous regions are experiencing rapid climate change (Pepin et al. 2022). Vegetation phenology has emerged as an essential variable for studying mountain ecosystem functions under the influences of climate change (Hwang et al. 2011; Kraus, Zang, and Menzel 2016; Vitasse, Signarbieux, and Fu 2018). For example, mountainous regions are typically characterized by high biodiversity (Perrigo, Hoorn, and Antonelli 2020), and shift in vegetation phenology is a critical factor affecting the suitability of habitats for certain species (Ettinger, Chamberlain, and Wolkovich 2022; Viña et al. 2016).

In mountainous regions, the large elevation gradient and complex terrain lead to a variety of sunlight, temperature, and moisture combinations, which consequently affect the spatial patterns of vegetation phenology (Adams et al. 2021; Dunn and de Beurs 2011; Hua, Ohlemüller, and Sirguey 2022; An et al. 2018). Knowledge of terrain controls (e.g. elevation, aspect, and slope) on vegetation phenology is critical to understanding the responses of vegetation phenology to climate variability and mountainous ecosystem functions (Hwang et al. 2011; Vitasse, Signarbieux, and Fu 2018). Satellite remote sensing serves as an indispensable tool for observing vegetation phenology in mountainous regions (Dunn and de Beurs 2011; Shahgedanova et al. 2021). Additionally, earth observation satellites generally measure the seasonality of a landscape or ecosystem rather than individual plants (Guan et al. 2014; Reed, Schwartz, and Xiao 2009); hence, satellite-derived vegetation phenology is also referred to as land surface phenology (LSP, de Beurs and Henebry 2004; Henebry and de Beurs 2013).

Numerous studies have investigated the impacts of elevation on the spatial patterns of LSP metrics (e.g. Ding et al. 2022; Misra, Asam, and Menzel 2021; Shen et al. 2014; Thompson and Paull 2017). In most cases, spring vegetation greenup is delayed along the elevation due to elevation-induced temperature differences (Shen et al. 2014; Vitasse, Signarbieux, and Fu 2018). In addition, the influences of elevation on LSP have been found to differ across elevations (Misra, Asam, and Menzel 2021). On the other hand, some regions also exhibited negative elevation gradients of spring LSP (Elmore et al. 2012; Fisher, Mustard, and Vadeboncoeur 2006). For instance, later greenup dates of deciduous forests at lower elevations were observed in a small valley in the USA (Fisher, Mustard, and Vadeboncoeur 2006). This phenomenon was attributed to a unique microclimate condition in which cold air drainage results in a cooler valley bottom (Fisher, Mustard, and Vadeboncoeur 2006). Moreover, Elmore et al. (2012) reported a conversion of elevation gradients of spring LSP with a conversion point at 275 m. In addition to cold air drainage, Elmore et al. (2012) suggest that the negative gradient below 275 m may be related to winter chilling and atmospheric boundary conditions.

Solar radiation intensity and photoperiod exhibit variations with aspect, which can lead to diverse near surface temperature and soil moisture conditions (Newman et al. 2006; Zapata-Rios et al. 2016). These aspect-induced local climatic differences reportedly affect LSP patterns in some mountains (e.g. Dunn and de Beurs 2011; Hua, Ohlemüller, and Sirguey 2022; Hwang et al. 2011; An et al. 2018). Specifically, spring LSP occurred earlier on south-facing slopes in some regions (Dunn and de Beurs 2011; Hwang et al. 2011; An et al. 2018). Hwang et al. (2011) inferred that disparities in photoperiod may explain the variations in growing season length of forests caused by aspect in North Carolina, USA. Furthermore, different influences of aspect on LSP of subtypes of grasslands at different elevations were observed (Hua, Ohlemüller, and Sirguey 2022; An et al. 2018). The diverse roles of aspect in shaping LSP patterns are tied to the local environmental conditions (Dunn and de Beurs 2011; An et al. 2018). Although terrain slope can affect solar geometry and soil moisture (Beven and Kirkby 1979; Newman et al. 2006), its impacts on LSP spatial patterns were found to be minor for some regions (Adams et al. 2021; Dunn and de Beurs 2011; Reaves et al. 2018). Overall, the effects of aspect and slope on LSP were generally weaker than that of elevation in many studies (e.g. Elmore et al. 2012; Hwang et al. 2011; Misra et al. 2018).

Some uncertainties and problems persist in explaining the impacts of terrain on LSP spatial patterns. Firstly, the presence of mixed pixels, particularly in coarse spatial resolution satellite images, can lead to ambiguous results when modeling the influences of terrain on LSP. For example, many negative elevation gradients of the start of growing season in the Nanling Mountains of southern China were caused by mixed pixels at lower elevations (Ding et al. 2022). Secondly, a mountain slope may be represented by several coarse resolution pixels, or even a single pixel in complex terrain regions. This is likely to impede the analysis of the effects of aspect on LSP patterns (Hwang et al. 2011). Thirdly, the relative importance of terrain factors in shaping LSP patterns could be linked to the spatial extent of the study area (Elmore et al. 2012). Finally, it is essential to characterize how the terrain controls on LSP differ across environmental gradients (Misra et al. 2018). Accordingly, a multiscale spatial analysis that uses LSP metrics derived from high or median spatial resolution satellite measurements is needed (Elmore et al. 2012; An et al. 2018).

Here, we aim to investigate the impacts of terrain factors (elevation and aspect) on the spatial patterns of LSP metrics in the Tianshan Mountains in Xinjiang, China, by employing satellite observations with a median spatial resolution of 30 m. Terrain slope was excluded from this study due to its correlation with elevation over the study area. The Tianshan Mountains show a large environmental gradient with diverse hydrothermal conditions (Hu 2004). Hence, terrain controls on LSP patterns are expected to vary with hydrothermal conditions across elevations. To this end, this study seeks to answer two main questions: (1) What are the impacts of elevation and aspect on LSP metrics? (2) How do these impacts vary across elevations?

2 Materials and methods

2.1 Study area

This study was performed in the western of the Tianshan Mountains in Xinjiang, China (Figure 1). This region encompasses an impressive elevation span exceeding 3500 m (Figure 1(b)). The region below 2000 m a.s.l. presents a flat topography. However, the landscape transforms dramatically beyond 2000 m, exhibiting a steep elevation gradient (Figure A1). The region predominantly features western and eastern terrain aspects. Located in a temperate arid and semiarid region, the vegetation is dominated by grasslands. The species composition shows noticeable variations along the elevation gradient (Hu 2004). At low-lying elevations, the major vegetation type is temperate desert steppe (Hu 2004; Xu et al. 2011). As the elevation increases, a transition occurs to mountain steppe, followed by forests, alpine meadows, and ultimately, alpine tundra (Hu 2004; Xu et al. 2011).

Figure 1. An overview of the study area. (a) Location of the study area, (b) elevation derived from the NASADEM product, (c) transformed aspect (TA), and (d) land cover map at 30 m spatial resolution compiled from the ESA WorldCover 2021 land cover map V2.

2.2 Datasets

2.2.1 Surface reflectance time series

We used the Harmonized Landsat 8 and Sentinel-2 (HLS) surface reflectance product with a spatial resolution of 30 m generated by NASA (Claverie et al. 2018) to estimate LSP metrics. The HLS product has a revisit interval of approximately 2–3 days, a significant improvement over Landsat 8 and Sentinel-2A/2B. This enhancement is due to the integration of available images from all three satellites (Claverie et al. 2018). The shortened revisit interval is advantageous for LSP studies that require high frequency observations (Bolton et al. 2020; Li and Roy 2017; Zhang et al. 2020). The HLS product already incorporates atmospheric correction and bidirectional reflectance distribution function (BRDF) reflectance normalization. In addition, the HLS product adjusts the surface reflectance values of Sentinel-2 visible and infrared bands to align with the corresponding bandpass of the Landsat 8 OLI sensor using empirical functions (Claverie et al. 2018).

Version 2.0 of the HLS surface reflectance time series product contains two sub-datasets: the HLS L30 (derived from Landsat 8, Masek et al. 2021a) and S30 (derived from Sentinel-2A/2B, Masek et al. 2021b) data. We collected the HLS data for the years 2021 and 2022 from the Earthdata website (https://search.earthdata.nasa.gov/). It should be noted that the revisit interval of the HLS time series varies across locations. For the majority of pixels within our study area, there are more than 100 observations each year.

2.2.2 Digital elevation model (DEM) data

We obtained the DEM data provided by the NASADEM product (NASA JPL 2020) from https://search.earthdata.nasa.gov/search. This product is a refined version of the Shuttle Radar Topography Mission (SRTM) DEM data, delivered at a spatial resolution of 30′ under the WGS 84 geographical coordinate system. We adjusted the projection of the DEM data to match that of the HLS data.

First, we calculated the aspect in degrees, oriented clockwise from north. We then transformed the aspect to a linear format, with values ranging from 0 (corresponding to a northeast slope) to 1 (representing a southwest slope) through Eq. 1 (Neitlich and McCune 1997).

(1) $TA=0.5\ast \left(1-cos\left(A-45\right)\right)$(1)

where TA represents the transformed aspect, and A represents the original aspect. TA is also referred to as the heat load index since the southwest slope tends to receive stronger solar radiation than other aspects (Neitlich and McCune 1997). Figure 1(c) displays the spatial distribution of TA across the study area. Additionally, the terrain slope is presented in Figure A1. Notably, a strong correlation exists between slope and elevation in this region. To reduce the effects of multicollinearity, we opted to exclude terrain slope when modeling the spatial relationships between LSP metrics and terrain factors (Hwang et al. 2011).

2.2.3 Land cover map

The ESA WorldCover 2021 land cover map V2 with a 10 m spatial resolution (Zanaga et al. 2022) was downloaded from Zenodo (https://zenodo.org/record/7254221). This land cover product uses Sentinel-2 and Sentinel-1 time series images as primary inputs. We rescaled the land cover map to a 30 m resolution using a majority resampling rule (Figure 1(d)). It is noteworthy that forests predominantly exist at mid-range elevations. Here, we focus on herbaceous vegetation that spans a broad elevation spectrum and generally shows gradual changes in species composition along the elevation gradient. Therefore, we confined our LSP analysis to grasslands and moss and lichen.

2.3 Processing of HLS time series data and estimation of LSP metrics

We adopted the enhanced vegetation index (EVI) to estimate LSP metrics due to its advantage of reducing the effects from soil background and residual aerosols (Huete et al. 2002). In addition, we calculated the normalized difference snow index (NDSI, Hall, Riggs, and Salomonson 1995) to identify snow-affected EVI values. The equations for calculating the EVI and NDSI are denoted as follows.

(2) $EVI=2.5\ast \frac{{\rho }_{NIR}-{\rho }_{Red}}{{\rho }_{NIR}+6\ast {\rho }_{Red}-7.5\ast {\rho }_{Blue}+1}$(2)

(3) $NDSI=\frac{{\rho }_{Green}-{\rho }_{SWIR}}{{\rho }_{Green}+{\rho }_{SWIR}}$(3)

Here, ${\rho }_{Blue}$, ${\rho }_{Green}$, ${\rho }_{Red}$, ${\rho }_{NIR}$, and ${\rho }_{SWIR}$ represent the surface reflectance values of the spectral bands blue (0.45–0.51 μm), green (0.53–0.59 μm), red (0.64–0.67 μm), near infrared (0.85–0.88 μm), and shortwave infrared (1.57–1.65 μm), respectively.

For each pixel, observations with NDSI > 0.1 during the day of year (DOY) periods 1–151 and 274–365 were considered as snow-affected values (Riggs, Hall, and Román 2016). We used the cloud mask provided in the quality assessment (QA) layer of the HLS product and the surface reflectance of the blue band to identify cloudy observations. This operation was taken to minimize the instances of misidentifying clear pixels as clouds, which can occur when solely using the QA layer. The criterion for identifying clouds is shown in Eq. 4.

(4) $Cloud=\left(cloudmask=true\right)\text{amp}({\rho }_{Blue}>0.15)$(4)

The threshold of the reflectance of the blue band was set following Oreopoulos et al. (2011). Cloud shadows were identified exclusively using the QA layer. In addition, it is important to note that in this mountainous region, there is a high prevalence of mountain shadows (i.e. dark pixels) in the HLS images, especially during the winter season. These dark pixels are generally associated with abnormal EVI values. We identified dark pixels using the surface reflectance values of the five HLS bands necessary for calculating the EVI and NDSI. Pixels with a surface reflectance value smaller than 0.02 in any of these five bands were regarded as dark pixels, and their corresponding EVI values were subsequently removed from the time series. Snow-affected EVI values were filled by the 5^th percentile of the remaining clear EVI values (Gray, Sulla-Menashe, and Friedl 2019). Nonetheless, noisy observations may still exist due to the effects of undetected clouds and residual aerosols (Alvarez-Mendoza, Teodoro, and Ramirez-Cando 2019; Zhang et al. 2020). To mitigate the remaining noises in the EVI time series, we conducted a 3-day maximum value composite (Zhang et al. 2020).

We employed the asymmetric Gaussian function fitting in the TIMESAT 3.3 software to further smooth the EVI time series (Eklundh and Jönsson 2017; Jönsson and Eklundh 2002, 2004). The median filter with a spike value of 1.0 was implemented in TIMESAT to eliminate the remaining outliers in the EVI time series. The iteration times of fitting the envelope of the EVI time series were set as 3, with an adaptation strength of 3. We utilized the relative threshold method to estimate LSP metrics. Specifically, the timings of the EVI values corresponding to the 20% and 90% of EVI amplitude plus the minimum EVI value during the vegetation green-up period were identified as GU20 and GU90, respectively (Guan et al. 2014; Xie et al. 2022). The GU20 metric has been widely regarded as indicating the onset of vegetation greenup in spring in many studies (e.g. O’Connor et al. 2012; Xie et al. 2022). The GU90 metric represents the onset of maturity, which is close to the timing of peak greenness during the growing season (Gray, Sulla-Menashe, and Friedl 2019; Guan et al. 2014).

2.4 Statistical analysis

We analyzed the spatial relationships between the two LSP metrics and terrain factors (elevation and TA) at two spatial scales, herein referred to as spatial extents. Ordinary least square (OLS) linear regression between the LSP metrics and terrain factors was performed for the entire study region. We also used geographically weighted regression (GWR), a local form of regression model, to quantify spatially nonstationary local relationships (Brunsdon, Fotheringham, and Charlton 1996; Fotheringham, Brunsdon, and Charlton 2002). GWR executes a local linear regression for each data point. For the local regression of a target point, an input data point is assigned a weight (W) based on its distance (d) to the target point. Here we used a bisquare kernel function with a bandwidth of 2.5 km (b = 2.5 km) to quantify W (Eq. 5, Fotheringham, Brunsdon, and Charlton 2002)

(5) $W=\left\{\begin{array}{c}{\left(1-{\left(\frac{d}{b}\right)}^2\right)}^2,d\le b\\ 0,d>b\end{array}\right\}$(5)

Using this kernel function, only the points located within a 2.5 km buffer of the target point were included in the local regression. Although GWR was performed for the entire region, we used the local regression coefficients and the local determination coefficient (R²) within each spatial window with a radius of 2.5 km. Hence, the GWR outputs for each window can be regarded as local scale results.

We used the boxplot method to eliminate outliers of GU20 and GU90 (Duarte et al. 2018). Then, we performed OLS regression and GWR via a spatially random selection of 20% of the pixels (accounting for 709,773 pixels). The regression slopes derived from OLS and GWR can be interpreted as the sensitivities of the LSP metrics to terrain factors.

3 Results

3.1 Spatial patterns of the LSP metrics

The spatial distributions of GU20 and GU90 in 2021 and 2022 are shown in Figure 2. The spatial patterns of each LSP metric remained relatively consistent between the two years. Nevertheless, considerable spatial variabilities were observed in both LSP metrics. For example, the interquartile range of GU20 in 2021 exceeded 30 days across the region (Figure 2(e)). The mean dates of GU20 over the region were DOY 124 and 125 in 2021 and 2022, respectively. Moreover, the mean dates of GU90 were delayed by 36 and 35 days in comparison to GU20 in 2021 and 2022, respectively.

Figure 2. Spatial patterns of GU20 and GU90. (a) GU20 in 2021, (b) GU20 in 2022, (c) GU90 in 2021, (d) GU90 in 2022, and (e) box chart of the LSP metrics. The black square in the box chart represents the mean value.

3.2 Spatial variations in the LSP metrics modeled by OLS

Across the entire study region, the spatial variations in GU20 were highly influenced by elevation for both years, as indicated by the results of the OLS linear regression (Table 1). For example, elevation alone explained 71% of the spatial variability in GU20 in 2021. In contrast, TA exhibited minor effects on the spatial patterns of GU20 at the scale of the entire study region, as evidenced by the very small R² values of OLS regression (Table 1).

Table 1. Statistics of the OLS linear regression models between the two LSP metrics and terrain factors. All OLS regression slopes were statistically significant (p < 0.05). ELEV represents elevation.

LSP metric	Regression model	Year	R^2	Slope of ELEV (days/m)	Slope of TA
GU20	GU20=a*ELEV+b	2021	0.71	0.030	-
		2022	0.73	0.024	-
	GU20=a*TA+b	2021	0.01	-	−7.4
		2022	0.01	-	−4.2
	GU20=a*ELEV+b*TA+c	2021	0.71	0.030	−3.8
		2022	0.73	0.024	−1.3
GU90	GU90=a*ELEV+b	2021	0.37	0.015	-
		2022	0.31	0.012	-
	GU90=a*TA+b	2021	0.01	-	−4.1
		2022	0.01	-	−3.2
	GU90=a*ELEV+b*TA+c	2021	0.37	0.015	−2.3
		2022	0.32	0.012	−1.8

The spatial variations in GU90 were also significantly influenced by elevation, rather than TA, across the entire region (Table 1). However, the explanatory ability of elevation decreased with vegetation growth, as seen in the considerably lower R² values of OLS regressions for GU90. For example, elevation explained only 31% of the variation in GU90 in 2022. Scatter plots illustrating the relationships between the LSP metrics and elevation are presented in Figure 3. We found numerous instances of later GU90 dates at lower elevations (<2000 m) than at elevations of 2000–2500 m (Figure 3). The spatial belts of later GU90 at lower elevations can also be clearly observed in Figure 2(c-d).

Figure 3. Scatter plots between the LSP metrics and elevation for the entire region. (a) GU20 in 2021, (b) GU20 in 2022, (c) GU90 in 2021, and (d) GU90 in 2022. The gray lines represent the OLS linear regression lines.

3.3 Spatial variations in the LSP metrics modeled by GWR

Statistics of the local R² and regression slopes of the GWR models are presented in Figure 4. Unlike the OLS models, the GWR models exhibited comparable local R² values between the regression models, with elevation and TA each functioning as a single explanatory variable (Figure 4(a-b)). Over 10% of the local R² values for TA exceed 0.3 for both LSP metrics (Figure 4(b)), indicating that TA can substantially affect GU20 and GU90 at a local scale (2.5 km radius in this study) in some regions.

Figure 4. Box charts for the GWR model outputs. Local R² values for the GWR models with (a) elevation, (b) TA, and (c) elevation and TA as explanatory variables, and the local regression slopes of (d) elevation and (e) TA. The gray and dark cyan boxes represent data in 2021 and 2022, respectively. The black square in the box chart indicates the mean value.

Consistent with the OLS models, most regression slopes of elevation were positive for the GWR models (Figure 4(d)). This suggests that in most local regions, elevation increases can lead to later GU20 and GU90 dates. The regression slopes of TA for GU20 and GU90 were mostly negative, indicating that grasslands on south-facing slopes generally experience greenup and maturity earlier than those on north-facing slopes (Figure 4(e)). Specifically, more than 10% of the samples in 2021 showed TA regression slopes for GU20 steeper than −10. This result indicates that GU20 can occur over 10 days earlier on the southwest slope than on the northeast slope in some regions. On the other hand, the proportions of positive slopes for TA regarding both LSP metrics were nonnegligible.

The spatial patterns of the local R² for the GWR models with both elevation and TA as explanatory variables in 2021 are displayed in Figure 5. The relationships between terrain factors (elevation and TA) and the two LSP metrics (GU20 and GU90) were dependent on elevation. More specifically, stronger correlations (indicated by local R² values > 0.4) were mostly found at elevations ranging from approximately 2000 m to 3000 m. The fact that these phenomena were consistent in 2022 (as shown in Figure A2) provides additional evidence supporting this finding.

Figure 5. Spatial patterns of the local R² for the GWR models with elevation and TA as explanatory variables in 2021. Subfigures (a) and (b) depict GU20, while (c) and (d) illustrate GU90. The brown circle and error bar in (b) and (d) represent the mean value and one standard deviation of the local R² values within the 200 m elevation bin, respectively.

In analyzing the sensitivities of the LSP metrics to elevation, as shown through the GWR regression slopes, we found that negative values primarily occurred in regions with lower elevations, specifically those less than 2000 m (Figures 6 and A3). For areas with elevations above 2000 m, a noticeable pattern emerged in the regression slopes of GU90, showing a trend of increasing and then decreasing (resembling a parabolic shape) along with the elevation. The specific elevation bins that exhibited the largest mean regression slopes showed differences between GU20 and GU90.

Figure 6. Spatial patterns of the local regression slopes of elevation for the GWR models with elevation and TA as explanatory variables in 2021. Subfigures (a) and (b) depict GU20, while (c) and (d) illustrate GU90. The brown circle and error bar in (b) and (d) represent the mean value and one standard deviation of the regression slope within the 200 m elevation bin, respectively.

The regression slopes of TA for GU20 and GU90 were also elevation-dependent (Figures 7 and A4). Positive regression slopes of TA were primarily found at lower elevations (<2000 m). The sensitivities to TA demonstrated an increasing trend with elevations between 2000 and 3000 m for both LSP metrics. For areas above 3000 m, the sensitivities of GU90 to TA slightly decreased with elevation. In addition, GU20 typically exhibited stronger sensitivities to TA compared to GU90 for areas above 2000 m.

Figure 7. Spatial patterns of the local regression slopes of TA for the GWR models with elevation and TA as explanatory variables in 2021. Subfigures (a) and (b) depict GU20, while (c) and (d) illustrate GU90. The brown circle and error bar in (c) and (d) represent the mean value and one standard deviation of the regression slope within the 200 m elevation bin, respectively.

4 Discussion

Spatial variations in GU20 and GU90 were observed to be predominantly influenced by elevation across the whole study area. This is in agreement with previous findings that have highlighted elevation as a principal determinant of LSP spatial variations, compared to other terrain factors (Dunn and de Beurs 2011; Elmore et al. 2012; Hua, Ohlemüller, and Sirguey 2022; Hwang et al. 2011; Misra et al. 2018). Furthermore, the degree to which elevation explained GU90 was considerably less than its influence on GU20 (Table 1). This was largely due to a higher incidence of later GU90 dates at lower elevations (<2000 m) as opposed to those observed between 2000 m and 2500 m elevations (Figure 3). There are two potential explanations for this finding. Firstly, grasslands at lower elevations within this region are generally drier compared to their counterparts at higher elevations (Hu 2004). Such water limitation may hinder vegetation growth and thereby lead to later maturity dates (Liu et al. 2022). Secondly, human activities, such as livestock grazing and forage planting, may also influence phenological patterns at lower elevations (Han et al. 2015; Zhao, Li, and Qi 2007). Evidence of human influences can be inferred from the LSP metrics images, specifically from the sharp and erratic phenological transitions (Figure 2).

With respect to TA, its explanatory powers for GU20 and GU90 were very limited over the whole region, as indicated by the R² values of the global OLS regression (Table 1). In contrast, the local R² values of elevation and TA from the GWR models with a bandwidth of 2.5 km were close (Figure 4a-b), which highlights the importance of TA in shaping the spatial patterns of GU20 and GU90 at a local scale. The ability to perform this local analysis in complex terrain was facilitated by the 30 m spatial resolution of the HLS dataset. The GWR results reveal that both GU20 and GU90 generally occur earlier on south-facing slopes (Figures 4 and 7), a pattern that echoes findings from previous studies in the Rocky Mountains in North America (Dunn and de Beurs 2011) and the meadow area on the Qinghai-Tibet Plateau (An et al. 2018). This impact could be partially attributed to the fact that south-facing slopes, which receive more solar radiation, tend to be warmer than north-facing slopes (An et al. 2018).

The influences of terrain on both LSP metrics were found to be elevation dependent (Figures 5-7). In lower elevation regions, the local regression slopes of elevation and TA showed a heterogeneous mix of both negative and positive values. As mentioned above, water availability and human activities might be the primary drivers of the spatial variations in LSP for lower elevation regions (<2000 m). In addition, the large percentage of positive regression coefficients for TA in lower elevation regions may reflect the impacts of water availability on LSP spatial patterns, as sun-exposed slopes are generally drier than shaded slopes. The mechanisms differ from those reported in some small valleys, where cold air drainage may lead to the varying controls of elevation on spring LSP (Elmore et al. 2012; Fisher, Mustard, and Vadeboncoeur 2006). It is noteworthy that the sensitivities of both LSP metrics to elevation and TA demonstrated nonlinear changes along elevation for regions above 2000 m. These patterns may stem from two interconnected factors. Firstly, different responses of vegetation growth to terrain factors may be caused by differences in species composition (Shen et al. 2022; Xie et al. 2015) and soil conditions (Liu et al. 2022) across elevations. For example, soil moisture has been reported to regulate the response of spring phenology to air temperature in grasslands (Liu et al. 2022). Secondly, the changing impacts of terrain on climate conditions, such as local temperature and snow cover durations, across elevations may also partially explain the varying sensitivities (Misra, Asam, and Menzel 2021; Shen et al. 2016; Tomaszewska, Nguyen, and Henebry 2020; Xie et al. 2021). These findings suggest that LSP at the highest elevations may not be the most sensitive to terrain-induced climatic differences. Further studies are required to elucidate the underlying mechanisms driving the varying responses of LSP to terrain across elevations.

5 Conclusions

In this study, we investigated the impacts of elevation and TA on the spatial patterns of two LSP metrics during the vegetation green-up period (GU20 and GU90) in the Tianshan Mountains, China using the HLS dataset with a 30 m spatial resolution. The results showed that the impacts of elevation and TA on the spatial variations in GU20 and GU90 were dependent on the spatial scale of analysis. At the scale of the entire study area, spatial variations in both LSP metrics were primarily controlled by elevation. However, when the GWR model with a 2.5 km bandwidth was applied, TA emerged as a significant factor affecting the two LSP metrics at a more local scale. For most areas, a larger TA was associated with an earlier date of GU20 and GU90. Additionally, the effects of terrain on GU20 and GU90 varied across elevations. Stronger explanatory powers of elevation and TA were observed at middle elevations (approximately 2000–3000 m) than at lower and higher elevations. The sensitivities of GU20 and GU90 to elevation and TA exhibited nonlinear but different change patterns along elevation. Our results highlight the diverse responses of the LSP metrics to terrain-induced local climatic differences in the arid and semiarid Tianshan Mountains. Specifically, LSP at the highest elevations may not be most sensitive to terrain-induced climatic differences.

Acknowledgments

We are grateful to the anonymous reviewers for their constructive comments.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data that support the findings of this study were obtained from publicly available datasets, and the websites providing access to these datasets are included in the article.

Additional information

Funding

This work was supported by the National Natural Science Foundation of China [Grant No. 42201348, 42201349, and 32201349].

Appendices

Figure A1. Terrain slope of the study area.

Figure A2. Spatial patterns of the local R² values for the GWR models with elevation and TA as explanatory variables in 2022. Subfigures (a) and (b) depict GU20, while (c) and (d) illustrate GU90. The brown circle and error bar in (b) and (d) represent the mean value and one standard deviation of the local R² values within the 200 m elevation bin, respectively.

Figure A3. Spatial patterns of the local regression slopes of elevation for the GWR models with elevation and TA as explanatory variables in 2022. Subfigures (a) and (b) depict GU20, while (c) and (d) illustrate GU90. The brown circle and error bar in (b) and (d) represent the mean value and one standard deviation of the regression slope within the 200 m elevation bin, respectively..

Figure A4. Spatial patterns of the local regression slopes of TA for the GWR models with elevation and TA as explanatory variables in 2022. Subfigures (a) and (b) depict GU20, while (c) and (d) illustrate GU90. The brown circle and error bar in (b) and (d) represent the mean value and one standard deviation of the regression slope within the 200 m elevation bin, respectively..