The Loess Plateau (LP) region is in the middle part of the Yellow River Basin (Figure 1). It is adjacent to the Qin Mountains in the south, the Yinshan Mountains in the north, the Taihang Mountains in the east, and the Wuqiaoling and Riyueshan mountains in the west. This topographically complex region covers 287 counties across Shaanxi, Shanxi, Inner Mongolia, Ningxia, Qinghai, Gansu, and Henan. Its territory comprises 620,000 km². The climate is arid in the northwest and semiarid in the southeast. The annual precipitation in this region ranges from 50 mm in the arid area to 700 mm in the semiarid area (Figure 1). Most rainfall occurred during June and September (Chen et al., 2007; Fu et al., 2017). The Loess Plateau is covered by thick loess and is well-known as the largest region with loess deposition, which is suitable for agriculture (Li et al., 2015). Gullies developed well in this topographically fragmented plateau (Cai, 2001; Li et al., 2003). In addition, the soil has a weak structure and is highly erodible, which makes the Loess Plateau a global hotspot of soil and water loss (He et al., 2004; Boardman, 2006). Revegetation was a dominant countermeasure adopted to prevent severe soil erosion on the Loess Plateau. Since the 1990s, driven by the national ecological restoration actions, the fast-growing economy, and the change in industrial composition, the Loess Plateau, particularly the gully and hilly areas, has entered an era of extensive and rapid revegetation through afforestation and grazing management (Fu et al., 2017; Wang et al., 2017).

The heterogeneity index of coupling (HIC) is proposed based on the reality of vegetation adaptation to environmental conditions adjusted by topographical variation. Under natural succession, topography with a great spatial heterogeneity shapes heterogeneous local settings, such as soil properties, available water, and solar radiation (Aguilar et al., 2010; Hoylman et al., 2019; Chen et al., 2024), for plant colonization and growth. While topography with low variation shapes homogeneous local settings and thus uniform vegetation cover. Consequently, the spatial heterogeneity of vegetation (H _v) tends to coincide with the spatial heterogeneity of topography (H _T). However, human interventions, such as afforestation in mountainous areas and diversity land use in plains, and natural disturbances such as fire and wind scouring may change this relationship and lead to the decoupling of H _v and H _T.

In arid and semiarid environments, the vegetation is water-limited (Wang et al., 2011; Wang et al., 2013), and topography plays a critical role in shaping the heterogeneity of soil water at the local scale, which leads to the variation in vegetation status among aspects and slope positions (Liu and Fu, 2013). Without anthropogenic intervention, the heterogeneity of vegetation growth status ( H v ) tends to couple with the heterogeneity of topographical features ( H T ), which can be quantified as the standard deviation of elevation within the unit area. Thus, the vegetation index in areas with small topographical variation shows low spatial heterogeneity, while a great topographical variation leads to a great spatial heterogeneity of the vegetation index. As illustrated by Figure 2, the external disturbances and interventions, such as planting trees in the Loess Plateau (Xiao, 2014; Zhou et al., 2019), may change the status of coupling between H _v and H _T. The afforestation usually shapes a vegetation cover with simplified species compositions of similar age and may lead to spatial homogeneity in the vegetation growth status (Figure 3). In addition, the diverse utilization of the land may create highly heterogeneous vegetation cover and vegetation growth status but seldom change the topography. Thus, the heterogeneity of vegetation cover decouples with the heterogeneity of topographical heterogeneity when there were natural disturbances such as fire and human activities such as afforestation and land utilization. The natural disturbances also shape the heterogeneity of vegetation cover and topography. Therefore, the consistency between H v and H T may reflect human interventions and various disturbances to the landscape.

In this study, the HIC was developed to quantify the coupling between the spatial heterogeneity of the vegetation index and spatial variation in topography. The coupling and decoupling between the spatial heterogeneity of topography and that of vegetation growth status can be observed by the coincidence and deviation between H _v and H _T. H _T and H _v is indicated by the spatial variation in elevation and the vegetation index, respectively. HIC is quantified as the difference between H v and H T , as shown by Eq. 1. H I C = H V − H T (1)

HIC > 0 implies the greater spatial heterogeneity of the vegetation index than that of elevation. While HIC < 0 indicates a lower spatial heterogeneity of the vegetation index than that of elevation. HIC equals 0 means a total coupling between the spatial heterogeneity of vegetation index and the spatial heterogeneity of elevation.

The Block Statistics tool of ArcGIS 10.3 (Environmental Systems Research Institute Inc., California, United States) was applied to calculate HIC. The Block Statistics tool performs a neighborhood operation that calculates a statistic for input grid cells within a fixed set of non-overlapping windows. The statistic (e.g., standard deviation, maximum, average, or sum) is calculated for all input cells contained within each window. The resulting value for an individual block is assigned to all cells located in the minimum bounding rectangle window. Since the windows do not overlap, any cell is included in the calculations for one block only. The shape of a window was rectangle in this study. Possible statistics that can be calculated within a window are mean, majority, maximum, median, minimum, minority, range, standard deviation, sum, and variety. In this study, the MODIS NDVI with a spatial resolution of 250 m and the grid DEM with a cell size of 90 m were used to calculate the HIC of rectangle windows with a size of 2.5 km × 2.5 km. Therefore, each window covers 100 grid cells of NDVI and approximately 771 grid cells of DEM. H v and H T of each window are expressed as the normalized standard deviation of NDVI (STD _v) and elevation (STD _T), respectively, followed by Eqs 2, 3. H v = S T D v − S T D v _ m i n S T D v _ m a x − S T D v _ m i n (2) H T = S T D T − S T D T _ m i n S T D T _ m a x − S T D T _ m i n (3) where S T D v and S T D T stand for the original standard deviation of the vegetation index (STD _v) and elevation (STD _T) within a rectangle block, respectively; S T D v _ m i n and S T D T _ m i n denote the minimum and maximum of STD _v and STD _T, respectively; and S T D v _ m a x and S T D T _ m a x represent the minimum and maximum of the standard deviation of the vegetation index and elevation of the study area, respectively.

The Loess Plateau is an ideal region to test the effectiveness of HIC for several reasons. First, the topography in this area varied from south to north and from west to east (Figure 1). Second, the Loess Plateau covered areas experienced rapid revegetation after 1999 (Fu et al., 2017) and areas covered by forest experienced long-term succession (Zheng, 2006; Liu et al., 2020). Third, this region included high plains, thick loess-covered gully and hilly mountains, sandy areas, rocky mountainous areas, and basins with large plains. In addition, the remote area with low population density and spot residential areas and the areas with high population density and sprawl built-up areas concurred in this region. So the Loess Plateau provided an ideal test scenario for HIC.

The 16-day composited MODIS NDVI dataset (MOD13Q1) with a resolution of 250 m during 2000–2015 was downloaded from the website of NASA (https://earthdata.nasa.gov/). This NDVI dataset was produced based on a surface reflectance algorithm applying a set of filters based on quality, cloud, and viewing geometry (Didan, 2021). The production flow of this dataset also included minimizing the BRDF (bidirectional reflectance distribution function) effects, the impact of cloud/smoke contamination, and the impact of aerosols. The annual sums of the 16-day composited NDVI were calculated for the Block Statistics in ArcGIS 10.3 (Environmental Systems Research Institute Inc., California, United States). The annual maximum NDVI was calculated based on the 16-day composited MODIS NDVI data. Geomorphological zonation data were downloaded from the Resource and Environment Data Cloud Platform (http://www.resdc.cn/). The land cover data on the Loess Plateau in 2000 and 2015 were produced by using the Landsat imageries (Zhang et al., 2014; Zhu et al., 2017). Specifically, the Landsat TM imageries were adopted to produce the land cover in 2000. For producing the land cover of 2015, Landsat 8 OLI imageries were used. The object-based classification (Zhang et al., 2014) was used to extract land covers, including forest land, shrubland, grassland, farmland, built-up area, water bodies, and bare area. The 90-m SRTM DEM was downloaded from the geospatial data cloud (http://www.gscloud.cn). The zonal statistics tool in ArcGIS 10.3 (Environmental Systems Research Institute Inc., California, United States) was used to summarize the HIC for geomorphological zones, land cover categories, and land cover change types. The consistence of HIC was evaluated by comparing HIC calculated by using NDVI, and EVI. The impact of spatial resolution was analyzed by adopting 250-m MODIS NDVI data and 10-m Sentinel-2A NDVI data in a sampled area in the central area of the Loess Plateau. Excel 2010 was employed in statistical analysis.

The land cover in the Loess Plateau region changed in large magnitude and with great spatial heterogeneity from 2000 to 2015. In this period, land cover change featured with expanded vegetated area and built-up area (Figures 4A, C, D), shrunk farmland and declined bare area (Figure 4B, C, D). Several urban agglomerations formed in the northern high plain and plains in the southern and eastern parts. In the gully and hilly areas that extend from the northeastern to southwestern part of the Loess Plateau, farmland shrank dramatically, accompanied by the vast expansion of grassland, shrubland and forest land (Figure 4). The bare area (mostly sandy area and desert) also retreated obviously in the northwest of the Loess Plateau. Forest cover expanded rapidly in the hilly and gully regions in the central area and rock mountains in the eastern part (Figures 4A, B), which was driven by extensive afforestation (Lü et al., 2012; Fu et al., 2017).

The spatial pattern of the HIC is consistent among datasets with different spatial resolutions. The heterogeneity of vegetation measured ( H v ) by using the EVI and NDVI of MODIS data and Sentinel-2 data was compared in a sampled area in the Loess Plateau. As shown in Figures 5, 6, the spatial pattern of H v was similar between EVI and NDVI, either based on MODIS data or Sentinel-2 data. However, the range of vegetation heterogeneity was greater when calculated using Sentinel-2 data than that derived from MODIS data (Figures 5, 6). Correspondingly, the spatial pattern and the range of HIC showed similar pattern with the spatial pattern of H v (Figures 5, 6E–H).

The HIC varies among land cover categories (Figure 7A). In 2000, the value of HIC ordered as bare area > forestland > shrubland > farmland > grassland. In 2015, the order changed into bare area > farmland > shrubland > grassland > forest land. In 2000, the HIC of all land cover types was greater than zero (Figure 7A). However, in 2015, the HIC of forest land, shrubland, and grassland was in minus, which implied that the spatial heterogeneity of NDVI was lower than that of topography. The HIC of farmland and bare area was greater than zero both in 2000 and 2015, which implied that the spatial heterogeneity of the NDVI was greater than the spatial variation of topography. The HIC of forest land, shrubland, grassland, and bare area in 2000 was greater than that in 2015. This indicates that the vegetation growth status in forest land, shrubland, grassland, and bare area became more homogeneous. However, the HIC of farmland has increased from 2000 to 2015, which may be attributed to the more diverse utilization of farmland.

Land cover change included shifts among the land cover categories and changes in properties such as configuration and coverage while maintaining the same land cover category. Both changes reshaped the coupling of spatial heterogeneity between vegetation indices and topography. As shown in Figure 7B, HIC decreased and deviated more from zero in the areas with the land cover shifting from farmland to forest, shrub, and grassland compared with areas with unchanged land cover types. However, the HIC of bare areas transferred from farmland was lower than that of unchanged bare areas, which can be attributed to the great spatial heterogeneity of vegetation in sandy bare land in the northwest part of this region.

The temporal dynamic of HIC in areas with unchanged land cover types varied among land cover types (Figure 8A). Generally, the HIC of bare areas showed apparently different temporal trajectory compared with other land cover types. It was above zero all the time and increased from 2001 to 2012, and then declined. HIC in forest land and shrubland generally declined and was gradually close to zero in this period. HIC in grassland was the lowest among land cover types in the studied period and close to zero. Farmland HIC was greater than zero all the time and declined before 2004, and then gradually moved into a steady stage (Figure 8A).

The temporal dynamics of HIC from 2001 to 2015 differed between land use change modes. Retiring farmland for revegetation and desertification of farmland are two avenues of farmland shrinking in the Loess Plateau. It led to a distinguished trajectory of HIC evolution (Figure 8B). In areas where farmland shifted into bare areas, HIC fluctuated greatly and was greater than zero during 2000–2015. However, in the areas where farmlands were retired for revegetation, HIC remained below zero and in a fluctuating pattern with a smaller magnitude than those areas transformed into bare areas (Figure 8B). After 2009, it declined (Figure 8B).

The spatial heterogeneity of topography and the vegetation index showed inconsistent patterns in the Loess Plateau region. In the northwestern and the southeastern parts of this region, the plain is dominant (Figure 9A). In the gully and hilly zones extended from southwest to northeast part of this region, the topography has great spatial variation (Figure 9A). The spatial heterogeneity of mean NDVI during 2001–2015 was great in the rocky mountainous area in the eastern part of this plateau (Figure 9B). However, in the gully and hilly zones, the spatial heterogeneity of NDVI was low (Figure 9B) and inconsistent with the highly heterogeneous topography. Consequently, the HIC was low in the hilly and gully areas. In areas covered by forests in the gully and hilly regions that experienced revegetation over 60 years (Zheng, 2006), HIC was apparently greater than in other areas with similar topographical patterns (Figure 9C). In plains such as those around Xi’an and Taiyuan, where farmland and residential landscapes dominated, the HIC also was rather high (Figure 9C).

As shown by the HIC, the temporal dynamic of the coupling between H v and H T varied significantly among geomorphic zones in the Loess Plateau region (Figure 9D). There are five geomorphic categories in this region, including the sandy area, hilly area, rocky mountainous area, Loess tableland, and plain (Figure 9). As illustrated by Figure 9D, during 2001–2015, the HIC in plain areas was the highest and ranged from 0.018 to 0.027, followed by the HIC of the sandy area (ranged from 0.009 to 0.018), rocky mountainous area (ranged from −0.03 to 0.004), hilly area (ranged from −0.011 to −0.006), and loess tableland (ranged from −0.026 to −0.017). The order of the five geomorphic zones remained the same during the study period (Figure 9). Except in the hilly area, the temporal trajectories of HIC in the geomorphic zones declined (Figure 9D). In the hilly area and the loess tableland, the HIC was below zero all the time in the study period. It means that the spatial heterogeneity of NDVI was lower than that of elevation all the time. In plain and sandy areas in the northern part of the Loess Plateau, HIC was above zero all the time during the study period, which implies that the spatial heterogeneity of NDVI was greater than that of topography all the time.

The HIC measuring the coupling between the spatial heterogeneity of vegetation status and the spatial variation in topography is an effective indicator to detect the human intervention of landscape in large extension. The spatial pattern of HIC was consistent with the regional pattern of revegetation. As revealed in this study, afforestation led to a lower NDVI spatial heterogeneity than the topographical heterogeneity (Figures 8, 9). There were afforestation in the abandoned farmlands in the hilly and gully areas and the sandy areas (Figures 3A, B, D, E, F), nursing the natural vegetation (Figure 3C) and rapid land use change in the past decades (Feng et al., 2016a; Fu et al., 2017). As shown in this study, during 2001–2015, the HIC of land use types consistently declined except of sandy areas. This implied a continuous declined spatial heterogeneity of NDVI caused by the extensive farmland retirement for revegetation. The declined HIC in areas with farmland abandonment for revegetation revealed that the revegetation of the farmlands reduced the spatial heterogeneity of vegetation. The Loess tableland area, hilly area, and rocky mountainous area were prior areas to implement the revegetation projects by planting trees (Figures 3A, B, D, E, F) (Cao et al., 2009; Chen et al., 2015; Zhou et al., 2019). The vegetation status was rather homogeneous, as revealed by the NDVI. Consequently, the HIC was relatively low, even though the topography was of great spatial heterogeneity. In the plains in the southern and northern parts, the land use types are diverse (Figure 4). Since most of the areas in plains are agricultural lands and residential areas (Figure 4), the land use intensity was high. Consequently, the spatial variation in vegetation was greater than that of terrain and was greatest among all the geomorphological zones. Due to the long history of cultivation, the HIC of plains was high and rather stable during the study period. In the plains, the sprawl of urban areas and built-up areas, and the diversification of the agricultural plots promoted the spatial heterogeneity of the vegetation status (Yan et al., 2009; Kou et al., 2021) and led to the greater NDVI heterogeneity than topographical variation (Figures 8, 9). Consequently, HIC was very high in these areas.

Resistance to disturbances sourced from physical processes or anthropogenic activities and intensity of disturbances are of critical importance for the coupling between the spatial heterogeneity of vegetation status and spatial variation of topography. As revealed by the spatial and temporal patterns of HIC among geomorphological zones, in the area where vegetation is vulnerable to natural disturbances or intensively altered by human interventions, the heterogeneity of vegetation status was greater than the topographical heterogeneity (Figure 9). The vegetation in the sandy area has a great inter-annual variation and spatial heterogeneity due to the intensive disturbance of sand movement (Richards et al., 2022). Therefore, the HIC of the bare area was the greatest among land cover types during the studied period (Figure 9).

As the ecosystem succession progresses, the homogeneous status of planted vegetation will change into heterogeneous patterns (Liu et al., 2020) due to the spatial variation in available water, radiation, and soil nutrients, as well the interaction of plants and various disturbances. The human intervention-shaped homogeneous land will change into the heterogeneous land due to adaptation to the local environment. Therefore, the decoupling between the heterogeneity of vegetation status and the topographical variation in revegetated areas will be replaced by the coupling as the ecosystem succession progresses (Figure 3C). The age of vegetation established by revegetation in most areas of the Loess Plateau was shorter than 30 years. From the perspective of ecosystem succession, the climate climaxes of species composition and coupling of vegetation cover to topography have not been reached. Therefore, during the mapped period, there was a lower heterogeneity of NDVI than the heterogeneity of topography in these areas.

Vegetation–topography coupling is an issue at various spatial scales (Yousefi Lalimi et al., 2017; Xiong and Wang, 2022; Brocard et al., 2023). As reported by various literatures (Radeloff et al., 2019; Shen et al., 2020; Shahzaman et al., 2021), different vegetation indices reflect different aspects of vegetation status. Therefore, the spatial heterogeneity of vegetation growth status may differ because of using different vegetation indices. As shown in Figure 5, the spatial patterns of H _v and HIC estimated by using the enhanced vegetation index (EVI) are well-consistent with those estimated by using the NDVI. The statistical comparative analysis of H _v (r = 0.89 for MODIS data and r = 0.93 for Sentinel-2 data) and HIC (r = 0.87 for MODIS data and r = 0.92 for Sentinel-2 data) also showed a good correlation between NDVI and EVI. However, the pattern of HIC may apparently be different when using the vegetation index with a different resolution (Figures 5, 6B). For H _v calculated by using NDVI and EVI, the correlation coefficient between MODIS data and Sentinel-2 data was 0.72 and 0.67, respectively. For HIC, they were 0.69 and 0.64, respectively. Generally, the consistency of HIC between the two data sources is better when using NDVI rather than using EVI. Since this study was conducted at regional scale, the resolution of the analysis is coarse. The further classification of forests was not adopted. However, the spatial pattern of land cover categories should be included in studies at the local scale because the high resolution is needed. For instance, more detailed forest classifications (Ivanova et al., 2022) related to specific local settings (Pfister and Arno, 1980) should be adopted.