The GLTFCA is located in southeastern Africa and spans across three countries—South Africa (ZA), Mozambique (MOZ), and Zimbabwe (ZIM), and five national parks—Banhine (MOZ), Gonarezhou (ZIM), Kruger (ZA), Limpopo (MOZ), and Zinave (MOZ) (Figure 1) (Home - Peace Parks Foundation, 2022). The national parks alone within the GLTFCA are approximately 45,000 km² in area (Table 1), but the entirety of the GLTFCA, which connects and surrounds the national parks, acting as functional buffer zones, is almost 100,000 km² in area (Home - Peace Parks Foundation, 2022). The GLTFCA is home to more than 850 animal species as well as over 2,000 plant species, all of which are now carefully managed for conservation (Home - Peace Parks Foundation, 2022).

This area is one of the earliest established and most successful examples worldwide of the conservation value of peace parks. Peace parks are defined as conservation areas that are managed cooperatively across geopolitical boundaries, and in this case specifically, are coordinated by a neutral third party, the Peace Parks Foundation (PPF) (Home - Peace Parks Foundation, 2022). The long-term goal of the PPF is to preserve large and functional ecosystems, such as the GLTFCA, so that the natural resources that humans depend on, such as food, medicine, clean water, and clean air are preserved for current and future generations (Home - Peace Parks Foundation, 2022). “Our dream is to reconnect Africa’s wild spaces to create a future for man in harmony with nature” (Home - Peace Parks Foundation, 2022).

Each of these ecosystems within the GLTFCA can be categorized as some type of savanna system, which range from grassland dominated to woodland dominated. One of the most important factors that influences the type of savanna is precipitation, and this whole area is categorized as dryland (Sankaran et al., 2005; Andela et al., 2013; Campo-Bescós et al., 2013b; Lehmann et al., 2014). Precipitation varies in these parks from Mean Annual Precipitation (MAP) of 439 mm in Bahine to 602 mm in Zinave (Figure 2), but as all these MAP values are below the well-published 650 mm threshold, the vegetation types are grass dominated savanna systems (Scholes and Archer, 1997; Sankaran et al., 2005; Staver et al., 2011; Campo-Bescos et al., 2013a). Southern Africa has two primary seasons, a rainy season from November to April and a dry season from May to October. Therefore, there is a high degree of seasonal fluctuation of precipitation across the year in these national parks (Figure 2).

Precipitation data was a key part of this analysis, given its important relationship with savanna vegetation, and healthy vegetation is the basis of conservation for all other species in this study area (Sankaran et al., 2005; Andela et al., 2013; Campo-Bescós et al., 2013b; Lehmann et al., 2014). The precipitation data was acquired from the Climate Hazards Group InfraRed Precipitation with Station Data (CHIRPS) website (Funk et al., 2015). This gridded dataset is available at 0.05° spatial resolution beginning in 1981 to present (Funk et al., 2015). To account for the discrepancy in spatial resolution of the precipitation data and the size of the parks, a mean value for each park over time and by seasons was extracted using the park shapefile subset over the CHIRPS grids.

In southern Africa, the dry season occurs from approximately May to October, and the wet season occurs from approximately November to April. Vegetation green-up lags 1 month behind precipitation occurrence (Zhu and Southworth, 2013). Triads were created to represent the seasons (beginning wet, end wet, beginning dry, end dry). Therefore, precipitation/vegetation triads were defined respectively as November-December-January (NDJ)/December-January-February (DJF); February-March-April (FMA)/March-April-May (MAM); May-June-July (MJJ)/June-July-August (JJA); August-September-October (ASO)/September-October-November (SON).

The two time periods evaluated are 2000–2020 to match the beginning of the available satellite imagery, and 1981–2020 which is the longest available climate data to establish longer term trends for this region through CHIRPS. Evaluating longer-term climate trends using a continuous time period is a standard practice in savanna system research (Scholes and Archer, 1997; Sankaran et al., 2005; Staver et al., 2011; Herrero et al., 2016; Herrero et al., 2020b). CHIRPS data was downloaded at the monthly time scale and accumulated in the following ways. This includes: a) total annual precipitation (defined as February 2000 to January 2021 to match the start of satellite vegetation data in March 2000), b) total precipitation by each season from 1981 to 2020 for each park (where mean seasonal precipitation was also calculated), and c) similarly, total precipitation by each season over time from 1999 to 2020 by park. To determine any significant changes in precipitation over time, Z-scores were calculated to identify anomalies, and linear trendlines were creat ed (Herrero et al., 2020b). This was done in R and represents 40 statistical tests: the five national parks, each tested for significant increases or decreases in precipitation for each of the four triads at both time periods a) 1981–2020 and b) 1999–2020.

The normalized difference vegetation index (NDVI) is a measure of vegetation abundance and greenness and is often used in the context of evaluating the health of conservation areas (Herrero et al., 2016; Herrero et al., 2019; Herrero et al., 2020b; Herrero et al., 2020a; Blentlinger and Herrero, 2020). NDVI uses a ratio between the red band (R, absorptive in healthy vegetation) and the near infrared band (NIR, reflective in healthy vegetation) from multispectral satellite imagery to determine a pixel value between −1 and +1, where higher values indicate greener/denser vegetation and values less than 0 indicate no vegetation (Measuring Vegetation, 2000). N D V I = N I R − R / N I R + R

The satellite imagery used for this analysis was an NDVI composite product from the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor. This product, MOD13Q1.006 Terra Vegetation Index, has a 250 m spatial resolution and is made from a 16-day composite of the maximum NDVI value per pixel. This data from Google Earth Engine comes atmospherically corrected as well as being masked for clouds, cloud shadows, heavy aerosols, and water (Google Earth Engine Explorer, 2021). This MODIS data was available from 2000 to 2020. The NDVI data was then aggregated into triads by mean value, hereafter referred to as seasons, to make it compatible with the precipitation regimes. The wet season is represented by DJF and MAM and the dry season is represented by JJA and SON, the spatial results of which are below in Figure 6.

The mean NDVI value for each season was calculated by pixel and then averaged by park. This was plotted as a continuous time series by park from 2000 to 2020 to demonstrate annual phenology in these protected areas. The mean NDVI pixel value was also calculated by seasons and park for the periods 2000–2004, 2005–2008, 2009–2012, 2013–2016, and 2017–2020 to further identify differences over finer time periods (multiple years to account for any potential anomalies). These values were also presented spatially through maps.

Greening studies have utilized remotely sensed time-series of vegetation indices, which included seasonality and serial auto-correlation. These studies have attempted to correct for trends using multiple statistical methods such as harmonic regression, time series using calendar days, linear models with non-parametric components for seasonality, etc. However, in a study by De Jong et al. (2012) they determined that greening or browning results varied significantly, depending on the method used in the analysis, and not necessarily the actual change on the landscape (De Jong et al., 2012). In response to these difficulties, the creation of a straightforward, repeatable, and statistically valid method was warranted (Southworth et al., 2023). NDVI time series can be used to study global vegetation change in several ways. D is a per pixel continuous measure of the direction of change in NDVI over time and identifies significant positive or negative cumulative change, as compared to a baseline condition (Waylen et al., 2014; Southworth et al., 2016; Herrero et al., 2020b, 2020a). The NDVI product (MOD13Q1.006) used for calculating D is described above in Section 2.3. This novel technique allows for both the space and time components to be addressed in and around these protected areas at a pixel scale over large landscapes and for extended periods of time. This technique provides a spatially broad and temporally long metric to determine patterns in vegetation stability or change. The application of this technique for understanding dynamics of these critical conservation areas, such as the GLTFCA, can assist in their management and protection (Herrero et al., 2020b; Herrero et al., 2020a).

To calculate the D statistic, a multi-year seasonal mean (2000–2004) is used as a baseline of comparison. Including multiple years in the baseline accounts for any anomalous years (Waylen et al., 2014). Following this, each ensuing season/year is compared to the baseline, creating a random walk statistic for each pixel (Waylen et al., 2014). This indicates that there is an equal chance that a pixel will be above or below the mean. If there is a “success” (NDVI is higher than the baseline), that pixel is assigned a value of +1. If there is a “failure” (NDVI is lower than the baseline), that pixel is assigned a value of −1. These pixel values are summed through time producing a net D value; in this study comparing 2005–2020 to the baseline creates N = 15. A pixel value of +15 would mean that every subsequent year was higher than the baseline (a very strong positive trend in vegetation health of the protected areas over time), and a pixel value of −15 would mean that every subsequent year was lower than the baseline (a very strong negative trend in vegetation health of the protected areas over time).

Where the value of statistical significance lies is based on a hypergeometric distribution (Waylen et al., 2014). Given the nature of the application of this study and investigating how the vegetation health of these conservation areas have changed over time and how these changes could make a meaningful difference in the parks, the lowest conventionally accepted statistical significance of 0.10 was used to allow for examining more change. Therefore, any pixels with a D value of ±7 are deemed significant.

The general relationship between precipitation and savanna vegetation using NDVI as a proxy has been established (Tucker, 1979; Carlson and Ripley, 1997; Lambin and Ehrlich, 1997; Scholes and Archer, 1997; Sankaran et al., 2005; Wang et al., 2005; Jiang et al., 2006; Bégué et al., 2011; Staver et al., 2011; Campo-Bescós et al., 2013a; Bunting et al., 2016; Southworth et al., 2016). However, in this study the specific question is determining if the savanna vegetation in these five national parks have a similar relationship with changing precipitation over time. Total seasonal precipitation was quantitatively compared to the mean NDVI over time in each park by regressing mean NDVI against precipitation and statistically interpreting the results (Tucker, 1979; Carlson and Ripley, 1997; Lambin and Ehrlich, 1997; Scholes and Archer, 1997; Sankaran et al., 2005; Wang et al., 2005; Jiang et al., 2006; Bégué et al., 2011; Staver et al., 2011; Campo-Bescós et al., 2013a; Southworth et al., 2016).

In Figure 3, total seasonal precipitation is plotted by park over the total time period of 1981–2020. The long term trends in precipitation were important to study in order to establish the historical conditions and climatological norms of each park, as conservation efforts are heavily influenced by climate. This was further broken down in Figure 3 into 1999–2020 to compare seasonally specific precipitation trends with satellite-based vegetation data in the form of NDVI. Figure 3 shows that there are strong negative trends during the dry seasons (MJJ and especially ASO) over time across all parks. During the wet season, and especially the onset of the rainy season (NDJ) there was a slightly stronger positive trend over time. This is particularly true in Zinave, which has the greatest mean annual and seasonal precipitation, Table 2. In FMA over the longer time period there was a slightly positive trend, however, in FMA from 2000 to 2020, there was a strong negative trend, especially in Zinave, indicating that there may be more variability in this park. These stronger trends at the beginning of the wet season and end of the dry season again support that precipitation is changing most likely in terms of temporal distribution, for example, early or late rains, etc., (Herrero et al., 2020b).

Results presented in Table 3 show that there are statistically significant negative trends in both the longer time series of 1981–2020 and from 1999 to 2020. The parks with the greatest (Zinave) and least (Bahine) mean annual and seasonal precipitation, had the most significantly negative trends. In every park there was a significant negative trend in precipitation during the later part of the dry season (ASO), occurring across the whole time period. Similarly, in the beginning of the dry season (MJJ) there was a significant downward trend in either the period 1981–2020 or 1999–2020 in every park, except Kruger, which only exhibited significant change in precipitation during ASO 1981–2020). In only one park, Zinave, during the beginning of the wet season (NDJ) over the whole time period of 1981–2020 was there any significant positive change in precipitation, which supports results in Figure 3. This further supports the theory that changes in precipitation most prominently occur at the edges of the seasons. In total, of the 40 statistical scenarios tested, there were 14 that were statistically significant 13 negative and 1 positive change in precipitation over time.

When approaching the Normalized Difference Vegetation Index (NDVI) as a time series, Figure 4 shows that there is a great deal of intra-annual variability in NDVI, which follows the rains, i.e., there is greater abundance of vegetation during the rainy season and a lower abundance of healthy vegetation during the dry season. When looking at these seasonal trends over time, each park, buffer, and the total landscape has a negative trend in the abundance of healthy vegetation (Figure 4). However, Zinave had the least negative trend, which aligns with the results of wet season precipitation above.

Figure 5 shows seasonal NDVI by season and park and depicts the expected variability between parks. During MAM and JJA, consistently across parks, the first and last time periods had the highest mean NDVI, with the lowest mean NDVI occurring during the second period. During SON, the first time period had the highest mean NDVI. In DJF, the third and fourth periods had the highest mean NDVI.

Figure 6 represents the spatialized context of NDVI across southern Africa by season and time period; there is a clear trend of minimum NDVI values in the JJA and SON seasons. The latter part of the wet season (MAM) represents the highest values of the year, and the latter part of the dry season (SON), represents the lowest values of the year.

There is a strong negative mean value of D (2000–2020, 2000–2004 baseline and 2005–2020 comparison to benchmark) across all five of the national parks and the buffer zone in all seasons (except for DJF) ranging from −4 to −12 (Figure 7). From the maps in Figure 8, it can be identified that the strongest negative trend is proximate to our study area, where South Africa, Mozambique, and Zimbabwe meet. Furthermore, their trends of negative persistence are far more extensive than positive persistence. Table 4 highlights that there are large portions of the National Parks and buffer zone in the GLTFCA with significant negative persistence across all seasons, whereas there are only very small portions of the National Parks in the GLTFCA with significant positive persistence across the seasons, again, except for DJF, where about 25% of the landscape has a significant positive persistence value. Spatialized D shows a very significant negative pattern in vegetation health across the study area in all seasons except those in the onset of the rainy season (DJF), where there are more neutral values (Figure 8).

Figure 9 demonstrates that there is a positive relationship between NDVI and precipitation within the five national parks from 2000 to 2020. Table 5 indicates that there is a significant positive relationship between precipitation and NDVI in each of the parks (p-values = 0.00). Notably, there are varying degrees of explanatory power of precipitation on NDVI, from 55% [Bahine (a), relatively smallest slope] to 79% [Kruger (c), relatively steepest slope], which is strong for a real-world relationship.