The directional increase of woody plants, in terms of biomass, stem density, and canopy coverage, has long been observed across open areas of continents and biomes (Stevens et al., 2017). This phenomenon is known as woody plant encroachment and causes social and ecological degradation in various ways (Archer et al., 1995(Archer et al., , 2017)). Gray and Bond (2013) showed in a conservation area of Africa that increase in tree density negatively impacted visitors' ability to see megafauna and undermined the tourism economy. Abreu et al. (2017) reported acute declines in both plant and animal species in the Cerrado of South America due to woody plant encroachment. Espunyes et al. (2019) suggested that woody plant expansion in the alpine grassland of Spanish Pyrenees disturbs the feeding efficiency of free-ranging livestock.
In the United States, the southern Great Plains (SGP) is receiving increased attention with regard to woody plant encroachment (Twidwell et al., 2013;Ge and Zou, 2013;Scholtz et al., 2018). It is mainly because this region has been largely altered by the encroachment, at an overall rate over five times higher than other regions of the US (Barger et al., 2011). The SGP spans three states of Kansas, Oklahoma, and Texas from north to south. In Kansas, 56% of tallgrass prairie in the Flint Hills are at the risk of converting to juniper (Juniperus) woodland due to low fire frequency (less than every 3 years) (Ratajczak et al., 2016). In Oklahoma, Wang et al. (2018) found that juniper forest (above 5 m in height) encroached into the subhumid and semiarid prairies at a regular pace of ~40 km 2 /year during 1984 and 2010. This encroachment was shown to affect land surface temperature and evapotranspiration (Wang et al., 2021). In Texas, the Edwards Plateau has mostly transitioned to juniper and oak (Quercus) woodland due to encroachment (Diamond and True, 2008). The transition has led to habitat fragmentation, native herbaceous species loss, and invasive grass establishment (Alofs andFowler, 2010, 2013).
The above studies in the SGP are small scale and focused on areas with aggressive encroachment. However, sub-regions in the SGP have distinct encroachment trajectories, owing to different land use and management history, precipitation level, and soil type (Archer et al., 2011;Wilcox et al., 2018). At state level, for instance, Kansas and Oklahoma have a much shorter history of woody plant encroachment than Texas. This is because a lot of areas in the former two states were largely cultivated by Europeans since their settlement in the late 1800s. Extensive encroachment did not occur until those areas were returned to grassland later when cultivation proved unsustainable (Wilcox et al., 2018). At county level, Wang et al. (2018) found substantial variation of juniper encroachment in Oklahoma. Within the Texas drylands, Asner et al. (2003) revealed that the woody cover increase between 1937 and 1999 varied widely according to management history.
Moreover, the existing studies usually span different time periods, making it even harder to extrapolate the revealed local encroachment trend to the whole SGP. Consequently, it is unclear whether or not the overall encroachment in the SGP is intensifying in response to climate change, as reported or predicted in other parts of the world (Donohue et al., 2013;Caracciolo et al., 2016;García Criado et al., 2020). Accordingly, the first objective of this study is to quantify the encroachment rate across the SGP over the recent decades.
Knowledge of the key factors and their role in encroachment is needed for targeted management plans and for prediction of future encroachment trends (Yang and Crews, 2020). Rosan et al. (2019) concluded that fire suppression and land use abandonment are the most important factors facilitating woody plant expansion in the Brazilian savanna in the last 15 years. In the SGP, the decreased fire occurrence due to the disruption of the fuel-fire feedback has been identified as a major cause of the encroachment (Fuhlendorf et al., 2008;Wilcox et al., 2018). Through long-term fire experiments in Kansas, Ratajczak et al. (2016) found that the tallgrass prairie with rare fire occurrence (less than every 10 years) can become woodland in 30-50 years. Collins et al. (2021) suggests that annual burning can slow the encroachment in the Flint Hills ecoregion.
In addition, precipitation affects woody plant encroachment. Yang and Crews (2020) revealed that precipitation exerts a significant positive effect on the juniper encroachment rate in the semiarid part of the Texas savanna. Woody plant encroachment in sub-Saharan Africa has also been correlated with higher precipitation (Brandt et al., 2017;Venter et al., 2018). Over time, Weber-Grullon et al. (2022) found that increased precipitation facilitated the germination and survival of honey mesquite in southern New Mexico, USA. According to precipitation projection, Yang et al. (2020b) suggests that the drying trend of Texas savanna in the 21st century will lower the potential woody cover, and therefore suppress the encroachment to some extent.
Initial forest area, the primary seed source, is a third factor influencing the encroachment rate (Yang and Crews, 2020). Kepfer-Rojas et al. (2014) showed that tree and shrub densities were lower in areas farther away from seed sources in a heathland of southwest Denmark. By model simulation, Caracciolo et al. (2016) demonstrated that initial shrub cover condition has a strong impact on the encroachment rate in the northern Chihuahuan desert. Venter et al. (2018) found in sub-Saharan Africa that regions with moderate initial woody cover (e.g. 30-60%) experienced the most rapid encroachment. Woods et al. (2019) demonstrated that woody plant encroachment rate has strong correlation with seed dispersal in a coastal grassland.
Despite the recognized importance of the above three factors, their role in the encroachment in the SGP have been rarely evaluated comprehensively and quantitatively. As such, the second objective of this study is to ascertain the role of fire frequency, precipitation level, and initial forest area in relation to encroachment rate across the SGP by quantitative analysis. For several reasons, this study will focus on the encroachment of evergreen forest in the SGP. Firstly, and most importantly, evergreen forest (primarily junipers) has been the main encroaching component (Engle et al., 2008;Barger et al., 2011;Twidwell et al., 2016;Zou et al., 2016). Secondly, the focus on evergreen forest can minimize the confounding effect of deciduous species (e.g. length of growing season, nitrogen-fixing ability) on the encroachment trend (Rogers et al., 2009;Buitenwerf et al., 2015). Thirdly, the unique phenology of evergreen forest (e.g. green foliage in winter season) exhibits strong spectral imprint (e.g. high NDVI value in winter season) (Wang et al., 2017(Wang et al., , 2018)). The availability of long history and large scale optical remote sensing dataset (e.g. Landsat) enables us to trace evergreen forest encroachment over the past several decades across the SGP. Overall, this study will examine the spatiotemporal variability of evergreen forest encroachment across the SGP in recent decades, and disentangle the role of fire, precipitation, and initial evergreen forest area in the encroachment.
The SGP (Kansas, Oklahoma, Texas) lies between 26 • and 40 • latitudes, and between -106 • and -94 • longitudes (Fig. 1). This vast expanse comprises arid, semiarid, subhumid, and humid climates (Wang et al., 2018;Yang et al., 2021). Historically, the SGP was primarily grassland and open savanna, maintained by the endogenous processes of fire and herbivory (e.g. Bison; Bond and Parr, 2010). Since the late 1800s, the European settlement brought changes to the landscape. The most important was active fire suppression across the whole region (Walker and Janssen, 2002). In addition, while Texas was mostly used as rangeland for livestock grazing, Kansas and Oklahoma were largely cultivated as cropland (Wilcox et al., 2018). The combination of fire suppression and overgrazing in Texas initiated woody plant encroachment in the state since then. But evident encroachment did not occur in Kansas or Oklahoma until the lands were returned to grassland later for sustainability purposes (Engle et al., 2008).
As aforementioned, junipers (Juniperus) have been the major encroaching component in the SGP. Three juniper species are prevalent enough, namely eastern red cedar (Juniperus virginiana), ashe juniper (Juniperus ashei), and redberry juniper (Juniperus pinchotii) (Lyons et al., 2009). These juniper species differ slightly in physical and growth characteristics. Besides junipers, southern live oak (Quercus virginiana) is an important encroaching evergreen species in this region (Starr et al., 2019).
The SGP contains a lot of US Level III ecoregions, which are distinguished by pattern and composition of climate, land use, vegetation, soil, hydrology, or wildlife, among others (Omernik, 2004). Pronounced differences in the encroachment stage and rate also exist among the ecoregions. For instance, woody plant cover in much of the Edwards Plateau in south Texas has reached its maximum according to the mean annual precipitation (MAP; Yang et al., 2020a). In contrast, the High Plains in west Texas has little encroachment so far (Yang et al., 2021). In Oklahoma, junipers have been expanding mainly in Cross Timbers and Central Great Plains (DeSantis et al., 2011;Williams et al., 2013). The Flint Hills region in Kansas has had little encroachment and maintains the largest tallgrass prairie remnant (Twidwell et al., 2016).
The ecoregions (US Level III) in southeast Oklahoma and southeast Texas were excluded from this study. This is because many of them are under extensive human cultivation (e.g. evergreen pine tree plantation and harvest; Fagan et al., 2018;Shephard et al., 2021). Their exclusion minimizes the bias that human cultivation may induce to estimates of the true encroachment rate. As for the other excluded ecoregions (e.g. Edwards Plateau), where forest or woodland has been widely established, restoration back to the original status is extremely difficult (Ratajczak et al., 2016). These ecoregions are of less urgency in terms of conservation. The remaining 35 ecoregions under study encompass broad physiographic gradients. Their MAP ranges between 241 and 1286 mm, while surface soil moisture varies from 2.9 to 19.7 mm. Surface temperature ranges between 287 and 312 K, while elevation increases from 141 to 2463 m (Fig. S1).
Yang et al. ( 2021) generated a 30 m resolution forest map (above 2 m in height) of the SGP for the time period of 2015-2017, with integration of radar and optical remote sensing data. Then they extracted evergreen forest from the forest map, by a threshold (0.3) of seasonal NDVI change (annual maximum NDVI minus winter season mean NDVI). Surface reflectance data (2015-2017) from Landsat 7/8 was utilized. The resulting evergreen forest map is displayed in Fig. 1. The inclusiveness (in terms of height) of this validated evergreen forest map (referred to as evergreen_base) makes it more suitable for studying woody plant encroachment, as compared to other evergreen forest maps (e.g. from NLCD) that target at trees above 5 m in height. This is because much of the encroaching evergreen forest (e.g. junipers) in the SGP is shorter than 5 m due to the limitation imposed by climate, soil type, tree species, and life stage (Scholtz et al., 2018).
The orthorectified and atmospherically corrected surface reflectance data from the Thematic Mapper (TM) of Landsat 4 and 5, and the Enhanced Thematic Mapper (ETM+) of Landsat 7 was used in this study. They share the same spatial resolution (30 m) and spectral resolution (e. g. 0.63-0.69 μm for red band, 0.77-0.90 μm for near infrared band). The surface reflectance data is available from 1982 to 1993 for Landsat 4, from 1984 to 2012 for Landsat 5, and from 1999 to 2021 for Landsat 7. The long history of these Landsat surface reflectance data provides us a great opportunity to trace evergreen forest encroachment across the SGP in recent decades. They were processed in Google Earth Engine.
Ecoregions are a spatial framework to assess and monitor the quality and quantity of environmental resources. They been widely applied in research and management (Fusco et al., 2019;Lacher et al., 2019). In this study, the US Level III ecoregions, subdivision of the coarser levels I and II, were adopted (Omernik and Griffith, 2014). The US Level III ecoregions shapefile with state boundaries was downloaded from the website of United States Environmental Protection Agency (https ://www.epa.gov/eco-research/level-iii-and-iv-ecoregions-continental-united-states). In this shapefile, the level III ecoregions crossing states were subdivided by state boundaries. This shapefile was clipped to our study area (Fig. 1).
Validation data was prepared for assessing the accuracy of the two composite evergreen forest maps of 2005-2009 and 2015-2017 at subregion scale. The validation data were collected in reference to timeseries (winter season) very high spatial resolution imagery in Google Earth Pro. During each of the two time periods, random plots of representative land cover types (according to NLCD: evergreen forest, deciduous forest, grassland, cropland, built-up, water body, and shrubland) were selected across each of the five sub-regions separately (Tables S1 andS2). The land cover types of all the validation plots were
The PRISM Climate Group's product of monthly total precipitation (rain + melted snow) for recent years (since 1981) was utilized (https:// prism.oregonstate.edu/recent/). This product is developed with the approach of climatologically-aided interpolation (CAI), which takes into account long-term climatic spatial patterns (Daly et al., 2008). It has a spatial resolution of 2.5 min (~4 km). In this study, monthly precipitation data from the three time periods 1990-1999, 2000-2009, and 2010-2017 was aggregated to respective MAP (Fig. 2). With each MAP layer, average MAP was calculated at ecoregion and sub-region scales respectively.
The MODIS Fire_cci Burned Area Pixel product (version 5.1, Fire-CCI51) was accessed in Google Earth Engine. It identifies burned area across the globe at ~250 m spatial resolution on a monthly basis since January 2001 (Chuvieco et al., 2018). The applied burned area algorithm takes into account the varying land surface cover and uses an adaptive thresholding approach (Lizundia-Loiola et al., 2020). This product has been validated and utilized in different scenarios (Hall et al., 2021;Oto´n et al., 2021). In this study, the FireCCI51 product within the two time periods 2001 to 2009, and 2010 to 2017 was composited into two burned area maps, which indicate the number of burned years for each pixel (Fig. 3). From each burned area map, a mean annual burned area (MABA) was calculated at ecoregion and sub-region scales separately.
This study aims to characterize the spatiotemporal variability of evergreen forest encroachment across the SGP, and ascertain the role of MAP, MABA, and initial evergreen forest area (Fig. 4). Due to the availability of historical Landsat surface reflectance data, this study will focus on the time period of 1985-2017. The encroachment will be monitored at decadal-scale time resolution (Axelsson and Hanan, 2018;Yang and Crews, 2020). First, on the basis of evergreen_base, this study will develop a single evergreen forest map for each of the four time periods 1985 to 1989 (evergreen_8589), 1995 to 1999 (evergreen_9599), 2005 to 2009 (evergreen_0509), and 2015 to 2017 (evergreen_1517), with the same set of region-specific NDVIwinter thresholds for evergreen forest. Second, this study will quantify the mean annual encroachment rate of evergreen forest at each of the three temporal stages 1990-1999, 2000-2009, and 2010-2017, at both ecoregion and sub-region scales. Third, it will relate the encroachment rate at each spatial scale to corresponding initial evergreen forest area, MAP, and MABA.
2.4.1. Region-specific NDVIwinter thresholds for evergreen forest Wang et al. (2017Wang et al. ( , 2018) ) suggest that evergreen forest exhibits a certain minimum level of NDVIwinter. As such, this study applied NDVIwinter threshold to identify the presence/absence of each evergreen forest pixel of evergreen_base in the past three time periods (1985-1989, 1995-1999, 2005-2009). According to our observation, herbaceous vegetation (e.g. grasses) in the SGP, especially in the southern part, is still vigorous and green in December. To minimize the effect of background (e.g. understory grasses) on NDVIwinter, this study defines winter season as January and February of each year. In addition, it has been observed that evergreen forest across the SGP varies a lot in NDVIwinter, which might be attributed to the encompassed broad physiographic gradients (e.g. precipitation, land surface temperature). To obtain accurate historical evergreen forest maps, the study area was divided into 5 sub-regions (Kansas, Oklahoma, west Texas, north Texas, southwest (SW) Texas) by Level III ecoregions, mainly according to their physiographic characteristics (Fig. 1). Then a NDVIwinter threshold was derived for each sub-region by the following steps.
First, sample evergreen forest plots (as of 2017) were collected for each sub-region. The number of plots/Landsat pixels is 67/3990 for Kansas, 54/5360 for Oklahoma, 51/4553 for west Texas, 41/4838 for north Texas, and 21/2231 for southwest Texas (Fig. S4). Second, all the high quality Landsat 7 observations in January and February of 2015-2017 were used to calculate NDVIwinter for the sample evergreen forest pixels. Third, an NDVIwinter density graph was developed for each of the five sub-regions (Fig. 5). Last, since we have high confidence in the sample data, a confidence interval of 99.7% was applied to obtain NDVIwinter thresholds. That is, the 0.3% percentile NDVIwinter value in each density graph was chosen as respective NDVIwinter threshold. The resulting NDVIwinter threshold is 0.33 for Kansas, 0.37 for Oklahoma, 0.31 for west Texas, 0.34 for north Texas, and 0.29 for southwest Texas.
In this study, historical evergreen forest maps refer to the presence/ absence of evergreen_base's evergreen forest pixels (during 2015-2017) in each of the three past time periods (1985-1989, 1995 to 1999, 2005 to 2009). First, for each year within the three time periods, a NDVIwinter layer was developed for the whole study area with available good quality observations from Landsat 4, 5 and 7 in January and February of the year. The number of images from each sensor across the years is summarized in Table S3. Second, for each of those years, an evergreen forest map was created. That is, for each evergreen forest pixel in ever-green_base, if its NDVIwinter in a given year is above corresponding region-specific threshold, it is identified as evergreen forest in that year. All other pixels will be marked as 'others' in that year.
Third, for each historical time period, a composite evergreen forest map was developed. For a given time period, a pixel will be classified as evergreen forest as long as it is identified as such in one or more years within the time period. This approach can minimize the effect of possible flash drought (drought intensification at sub-seasonal scale) that may result in abnormally lower NDVIwinter value (Christian et al., 2021). It can also reduce the effect of missing good quality Landsat observations in some winter seasons (Fig. S5). For easy description, the resulting evergreen forest maps are hereinafter referred to as ever-green_8589 for the time period 1985-1989, as evergreen_9599 for 1995-1999, and as evergreen_0509 for 2005-2009.
It is worth noting that although evergreen_base is for the time period 2015 to 2017, it was developed by applying a threshold of seasonal NDVI change (rather than NDVIwinter threshold) to a validated forest map (Yang et al., 2021). To ensure comparability among the time series evergreen maps, a new evergreen forest map for the time period 2015 to 2017 was developed by repeating the above three steps. It is referred to as evergreen_1517 and will be used in the quantification of evergreen forest encroachment.
The above four resulting evergreen forest maps (evergreen_8589, evergreen_9599, evergreen_0509, evergreen_1517) were used to quantify evergreen forest encroachment. The encroachment was identified between each two consecutive time periods, in reverse chronological order (Fig. 6). For each evergreen forest pixel in evergreen_1517, if it was not identified as evergreen forest in evergreen_0509, it will be marked as encroachment between 2010 and 2017; If it was identified as such in evergreen_0509 but not in evergreen_9599, it will be marked as encroachment between 2000 and 2009; If it was identified as such in evergreen_0509 and evergreen_9599, but not in evergreen_8589, it will be marked as encroachment between 1990 and 1999; If it was identified as such in evergreen_0509, evergreen_9599, and evergreen_8589, it will be marked as existing evergreen forest during 1985 and 1989. The maps of evergreen forest encroachment at the three temporal stages (1990-1999, 2000-2009, 2010-2017) can be found in Fig. S6. The mean annual encroachment rate (km 2 /year) at the three temporal stages were summarized at ecoregion scale and sub-region scale separately, with the function Zonal Statistics as Table in ArcMap.
An initial evergreen forest map was developed for each of the three temporal stages, namely initial_1990 for 1990-1999, initial_2000 for 2000-2009, and initial_2010 for 2010-2017. While initial_1990 equals the existing evergreen forest map during 1985 and 1989, initial_2000 is the sum of initial_1990 and the encroachment during 1990 and 1999, initial _2010 is the sum of initial_2000 and the encroachment during 2000 and 2009. The area of initial evergreen forest for each temporal stage is summarized at ecoregion and sub-region scales.
To understand the spatiotemporal variability of the evergreen forest encroachment in the SGP, both quantitative and qualitative analyses were performed. First, over each of the three temporal stages, the mean annual encroachment rate was related to the three factors of initial evergreen forest area, average MAP, and MABA at ecoregion-scale, by linear regression. These quantitative analyses were aimed to ascertain the determining factors of the encroachment rate. Second, pairwise comparison was conducted between the changing trend of encroachment rate and that of average MAP and MABA across the three temporal stages, at ecoregion and sub-region scales. These qualitative analyses were designed to elucidate the role of MAP and MABA in the encroachment. All these analyses were carried out in RStudio.
The resulting evergreen forest maps of the four time periods under the uniform criterion are displayed in red in Fig. 7. The area of evergreen forest is 5654 km 2 in evergreen_8589 (Fig. 7a), 6755 km 2 in ever-green_9599 (Fig. 7b), 8467 km 2 in evergreen_0509 (Fig. 7c), 9123 km 2 in evergreen_1517 (Fig. 7d). Among the three states, it is evident that Texas (west Texas, north Texas, and southwest Texas combined) has the most evergreen forest area in all the four time periods, while Kansas has the least amount. For a clearer picture, zoom-in views of a sample site (blue polygon in Fig. 7a) are displayed in Fig. 7e-h. This sample site is centered at (34.9312 • , -101.127 • ) and has a size of 18 km by 18 km.
Compared to evergreen_base (2015-2017), evergreen_1517 is more parsimonious. Specifically, the evergreen forest area in evergreen_1517 accounts for 83-89% of that in evergreen_base, depending on specific subregions (Table S4). The accuracy of evergreen_1517 and evergreen_0509 was assessed at sub-region scale, with the prepared validation plots of representative land cover types. As shown in Table 1, both ever-green_1517 and evergreen_0509 have very high accuracy across all the five sub-regions. The three initial evergreen forest maps initial_1990, initial_2000, and initial_2010 are displayed in Fig. S7.
The encroachment rate of evergreen forest in the SGP exhibits both spatial and temporal variability at sub-region scale (Fig. 8). Among the five sub-regions, west Texas experienced the most rapid encroachment at all the three temporal stages (66.2 km 2 /year, 92.0 km 2 /year, 53.8 km 2 /year), while southwest Texas saw the slowest encroachment (12.4 km 2 /year, 10.9 km 2 /year, 7.4 km 2 /year). The largest fluctuation in annual encroachment rate occurred in west Texas, while the least variation was observed in southwest Texas. Across the three temporal stages, all the five sub-regions exhibit an overall declining trend of evergreen forest encroachment. For instance, the mean annual encroachment rate in the sub-region Oklahoma varied from 46.7 km 2 /year during 1990 and 1999, to 38.3 km 2 /year during 2000 and 2009, and to 36.1 km 2 /year during 2010 and 2017. In Kansas, the respective encroachment rates are 14.5 km 2 /year, 14.3 km 2 /year, and 7.5 km 2 /year. Nevertheless, there are two local and temporary exceptions. The encroachment accelerated in west Texas at temporal stage 2 and in north Texas at stage 3. Among the 35 ecoregions under study, 17 have no evergreen forest throughout the study period. The spatiotemporal variability of evergreen forest encroachment rate in the remaining 18 ecoregions can be found in Figs. S8, S9, and S10.
Linear regression analysis suggests that the mean annual encroachment rate at ecoregion scale is strongly correlated with initial evergreen forest area at all the three temporal stages (Fig. 9). However, it shows no significant relationship with average MAP or MABA (Figs. S11 andS12). The initial evergreen forest area explains 68% of variation among the 18 ecoregions in mean annual encroachment rate of 1990-1999. It accounts for 93% and 70% of the ecoregion-scale variation over 2000-2009 and 2010-2017, respectively. The pairwise comparison of mean annual encroachment rate to concurrent average MAP at sub-region scale is displayed in Fig. 10. It is evident that the trend of encroachment rate across the three temporal stages is consistent with that of average MAP in all the sub-regions except southwest Texas (Fig. 10e). In southwest Texas, while the average MAP increased from 319 mm at stage 1-334 mm at stage 2, and to 338 mm at stage 3, the mean annual encroachment rate decreased from 12.4 km 2 /year to 10.9 km 2 /year, and to 7.4 km 2 /year. The most pronounced fluctuation in average MAP was observed in Oklahoma (Fig. 10b; from 870 mm to 823 mm, and to 787 mm), while the least variation was seen in north Texas (Fig. 10c; from 687 mm to 682 mm, and to 685 mm). It is interesting that the relatively wet sub-regions (Kansas, Oklahoma, north Texas) experienced an overall drying trend since 1990, while the dry sub-regions (west Texas, southwest Texas) saw a wetting trend.
The pairwise comparison of the mean annual encroachment rate to the concurrent MABA at sub-region scale over the latter two temporal stages is displayed in Fig. 11. The trend of the encroachment rate (over the latter two temporal stages) is opposite to that of MABA. In southwest Texas, the increase of MABA from 2 km 2 /year to 6 km 2 /year (Fig. 11e) might explain the decreased encroachment rate, despite its slight increase in average MAP from 334 mm to 338 mm (Fig. 10e). The largest fluctuation in MABA was seen in west Texas (from 7 km 2 /year to 12 km 2 /year). It is notable that the burned area increased over the latter two temporal stages in all the five sub-regions except north Texas, where MABA stayed at 2 km 2 /year (Fig. 11c). Kansas had significantly higher MABA than the other sub-regions. It has to be noted that the MABA of 2001-2009 was used as the proxy of 2000-2009 in this study, due to the lack of MODIS burned area product (FireCCI51) in 2000. Pairwise comparison at ecoregion scale also suggests that the trend of the encroachment rate tends to be consistent with that of average MAP, but opposite to that of MABA (Figs. S8, S9, S10).
The existing evergreen forest during 2015 and 2017 (evergreen_1517) is the net result of encroachment, climate impact, and human intervention. The formation process (pace) of these existing evergreen forest across the study area over the past several decades afford a unique lens to examine the spatiotemporal variability of evergreen forest encroachment in the SGP. It also provides an opportunity to assess the role of climate and human intervention in the encroachment. In this study, climate is represented by mean annual precipitation, while human intervention is reflected in mean annual burned area.
As shown in Fig. 8, all the 5 sub-regions of SGP exhibit an overall slowing trend of evergreen forest encroachment. The slowing trend is also exemplified by the zoom-in views (Fig. 7e-h) of the historical evergreen forest maps at the sample site (blue polygon in Fig. 7a). This finding is consistent with the conclusion by Deng et al. (2021) that global woody plant encroachment slowed down after 2010. However, it differs from the report and prediction in other studies. By synthesizing existing literature on woody plant encroachment, García Criado et al. (2020) concluded that the encroachment is intensifying in tundra and savanna biomes in recent decades due to increased precipitation. This discrepancy could be attributed to the prominent drying trend in the vast majority of the SGP (Figs. 2 and 10). In particular, the drying trend of wet sub-regions and the wetting trend of arid sub-regions in the SGP are contrary to the reported wetting trend of wet areas and drying trend of arid areas in other regions (Dore, 2005;Trenberth, 2011;Feng and Zhang, 2015).
The spatial variability of the encroachment rate among ecoregions in the SGP is largely associated with initial evergreen forest area (Fig. 9), but not related to MAP or MABA (Figs. S11 and S12). This is because initial evergreen forest is the major source of seeds, which are required for recruitment (Santos, 2010;Chaneton et al., 2012). This result agrees with the finding by Vitali et al. (2017) that the upward treeline shift in mountain landscapes is determined by distance to reproductive age trees (acting as a seed source) and microsite topography, rather than by climate. It is also consistent with the conclusion by Hanan et al. (2010) that the presence of adult trees is the most important driver of tree recruitment in the prairies of Florida.
The temporal variability of the encroachment rate in the SGP, however, is under the influence of MAP and MABA (Figs. 10 and 11), as reported in African savannas (Axelsson and Hanan, 2018). Specifically, it accelerates with MAP, but slows down with MABA. This finding agrees well with existing literature in the field. García Criado et al. (2020) revealed that woody plant encroachment in savannas is positively associated with precipitation. The encroachment rate in Kansas is slightly higher than that in southwest Texas at all the three temporal stages, although the latter has about four times higher evergreen forest area in 1990 (480 km 2 vs. 128 km 2 according to initial_1990). This is mainly because the average MAP in Kansas is above 720 mm, while the average MAP in southwest Texas is below 340 mm. With regard to fire, Miller et al. (2017) analyzed woody cover change in glade grasslands of Missouri, USA, in response to reintroduction of prescribed fire. It was found that woody cover increased in unburned glades but stayed the same in burned glades.
In comparison, fire plays a more dominant role than precipitation in the encroachment rate (Brunsell et al., 2017), as exemplified in southwest Texas. That is, although both initial evergreen forest area and average MAP increased from temporal stage 2-3, the encroachment rate declined by 32% at temporal stage 3. This is because the MABA expanded from 2 km 2 /year to 6 km 2 /year over the two temporal stages. The dominant role of fire can also explain the low encroachment rate in Kansas (MABA >72 km 2 /year), despite its high average MAP (>720 mm). The remaining intact tallgrass prairie in the Flint Hills should be attributed to its high frequency of prescribed fire (Fig. 3; Twidwell et al., 2016). These results suggest that prescribed fire could be an effective means to combat encroachment in the SGP before high level woody plant cover is established (Ratajczak et al., 2014). It is also applicable since social constraints that curb prescribed fire in this region are being overcome (Twidwell et al., 2013).
The role of other factors in the encroachment rate should not be overlooked. Firstly, different land cover types (e.g. barren land, grassland) exist in the SGP (Yang et al., 2021). Some land cover types may be more favorable than the others for evergreen forest encroachment. Secondly, discrepancy in response to fire exists among the evergreen species. For instance, after the tops are killed by fire, ashe juniper does not sprout while redberry juniper can sprout from a bud zone (Lyons et al., 2009). In contrast, low-intensity fire in dormant season can facilitate regeneration of oaks (DeSantis and Hallgren, 2011).
The four time-series evergreen forest maps of the uniform criterion (Fig. 7) are the linchpin of this study. Across the study area, the NDVIwinter of evergreen forest pixels varies significantly (Fig. 5). It reflects the effect of physiographic characteristics such as precipitation level, soil type, and topography on NDVIwinter (Onema and Taigbenu, 2009;Svoray and Karnieli, 2011;Meng et al., 2020). The region-specific NDVIwinter thresholds were proposed to minimize these effects and are key to the high accuracy of the four time-series evergreen forest maps. Additionally, their 100% user's accuracy (Table 1) should be credited to the double thresholds of seasonal NDVI change and NDVIwinter. That is, while evergreen_base was derived from a validated forest map by the threshold of seasonal NDVI change (NDVImax -NDVIwinter < 0.3; Yang et al., 2021), the four time-series evergreen forest maps (Fig. 7) were developed by further applying region-specific NDVIwinter thresholds to the evergreen_base for different time periods.
Compared to evergreen_base, evergreen_1517 is a refined evergreen forest map for the time period 2015 to 2017 and has slightly lower evergreen forest area (Table S4). The discrepant evergreen forest pixels are those that meet the seasonal NDVI change threshold, but their NDVIwinter does not reach the region-specific thresholds. While some of the discrepant evergreen forest pixels could be commission error of evergreen_base, the others could represent real evergreen forest. Nevertheless, evergreen_1517 is under the same standard as the other three historical evergreen forest maps, which ensures its validity for tracing evergreen forest encroachment. Although the seasonal NDVI change threshold approach can effectively extract evergreen forest out of forest map, it falls short of judging the historical status of current evergreen forest pixels. This is because other land cover types such as barren land and water body also have low seasonal NDVI change. This comparison reflects the strengths and shortcomings of different remote sensing approaches.
Limitations exist in this study. First, the four historical evergreen forest maps could be even more accurate if NDVIwinter thresholds were derived for more sub-regions. Second, among the four time periods (1985-1989, 1995-1999, 2005-2009, and 2015-2017) and across the study area, the number of winter seasons having good quality Landsat observations varies (Fig. S5). More winter seasons with good quality observation means less susceptibility of classification accuracy to shortterm abnormal conditions (e.g. flash drought). Third, an accuracy assessment was not performed for evergreen_8589 and evergreen_9599 due to the scarcity of high resolution (winter season) imagery in Google Earth over the earlier time periods. Although it is reasonable to assume that evergreen_8589 and evergreen_9599 have as good accuracy as ever-green_0509 and evergreen_1517 that developed with the same approach.
In conclusion, this study reveals the spatiotemporal variability of evergreen forest encroachment across the SGP. It suggests that the encroachment rate is strongly correlated with initial evergreen forest area that acts as seed source. The overall declining trend of evergreen forest encroachment since 1990 is consistent with the concurrent trend of decreasing precipitation and increasing burned area. With regard to the dominant role that fire plays in the encroachment rate, prescribed fire could be restored in the SGP to combat the encroachment there. These results also improve our ability to forecast future woody plant encroachment in other regions, under the context of global climate change. Lastly, this study embodies the usefulness of remote sensing approaches to tracing and understanding ecological processes of broad scale.