Moisture and temperature influences on nonlinear vegetation trends in Serengeti National Park

While long-term vegetation greening trends have appeared across large land areas over the late 20th century, uncertainty remains in identifying and attributing finer-scale vegetation changes and trends, particularly across protected areas. Serengeti National Park (SNP) is a critical East African protected area, where seasonal vegetation cycles support vast populations of grazing herbivores and a host of ecosystem dynamics. Previous work has shown how non-climate drivers (e.g. land use) shape the SNP ecosystem, but it is still unclear to what extent changing climate conditions influence SNP vegetation, particularly at finer spatial and temporal scales. We fill this research gap by evaluating long-term (1982–2016) changes in SNP leaf area index (LAI) in relation to both temperature and moisture availability using Ensemble Empirical Mode Decomposition and Principal Component Analysis with regression techniques. We find that SNP LAI trends are nonlinear, display high sub-seasonal variation, and are influenced by lagged changes in both moisture and temperature variables and their interactions. LAI during the long rains (e.g. March) exhibits a greening-to-browning trend reversal starting in the early 2000s, partly due to antecedent precipitation declines. In contrast, LAI during the short rains (e.g. November, December) displays browning-to-greening alongside increasing moisture availability. Rising temperature trends also have important, secondary interactions with moisture variables to shape these SNP vegetation trends. Our findings show complex vegetation-climate interactions occurring at important temporal and spatial scales of the SNP, and our rigorous statistical approaches detect these complex climate-vegetation trends and interactions, while guarding against spurious vegetation signals.


Introduction
Recently, a great deal of work has been carried out to identify and attribute late 20th century vegetation changes, namely an inter-annual 'greening' trend at global and regional spatial scales. While the full extent and causes of these trends are still under investigation (Cortés et al 2021, Walker et al 2020, CO 2 -fertilization, land cover change, and agricultural intensification may be important regional drivers (Forkel et al 2016, Zhu et al 2016, Chen et al 2019, Walker et al 2020. At smaller spatial scales, however, particularly across global protected areas such as IUCN Wilderness and National Parks (Stolton et al 2013), challenges remain in detecting and attributing vegetation changes and their drivers (Mondal et al 2020). Protected areas provide critical ecosystem services and sources of human livelihood, and quantifying their vegetation changes, particularly in response to anthropogenic forcings, are important for biodiversity goals.
Prior work has identified increasing land use as one major driver of ecosystem change in protected areas (Venter et al 2016, Jones et al 2018, particularly across Africa. In East Africa, land use change is critical to explain recent ecosystem shifts, e.g. to woodier vegetation in savannas, and vegetation productivity declines across protected areas (Landmann and Dubovyk 2014, Midgley and Bond 2015, Probert et al 2019, Kija et al 2020. Yet, additional work is also needed to identify how climate change may contribute to vegetation trends in East African protected areas (Midgley andBond 2015, Adole et al 2018). These ecosystems' responses to intensifying climate change are uncertain due in part to an incomplete understanding of climate-vegetation interactions, difficulty in detecting robust climate change signals (e.g. precipitation trends), differences in the scales and appropriateness of ecosystem indicators, and limitations on widespread, frequent ground observations (Pettorelli et al 2012, Rannow et al 2014, Mondal et al 2020.
However, advances in satellite and other remotesensing measures are enabling improved quantification of historical vegetation change and variability across Africa (Pettorelli et al 2012, Adole et al 2016, Hawinkel et al 2016. These products have been used to evaluate vegetation sensitivities at large spatial scales and identify complex vegetation-climate interactions. While many African ecosystems are waterlimited (Nemani et al 2003, Zhang et al 2006, Ugbaje and Bishop 2020, recent studies have explored vegetation dependencies on other environment and climate variables and interactions between them (Ryan et al 2017, Adole et al 2018. For example, precipitation responses of African vegetation, even in water-limited ecosystems, may be strongly mediated through interactions with thermal plant-growth thresholds and temperature effects (Adole et al 2016(Adole et al , 2019, which warrant further investigation under changing climate conditions. While these water-temperature interactions are generally important across Africa, in East Africa, precipitation variability remains a key driver of vegetation change (Nicholson et al 2017, Wei et al 2018. The El Niño-southern oscillation (ENSO) and Indian Ocean surface temperatures drive precipitation variability that explains much of the historic normalized difference vegetation index (NDVI) variability in East Africa (Hawinkel et al 2016). Additionally, during the late 20th century, intra-annual vegetation variability may have increased across East Africa (Ritchie et al 2008), particularly in protected areas in response to precipitation trends (Pettorelli et al 2012). A recent study controlling for land use change shows that East African woodland and savannas areas exhibit a delayed start to the growing season, thereby contributing to enhanced seasonal vegetation variability (Adole et al 2018). Embedded in this increasing variability may be long-term vegetation trends and even nonlinear responses to changing climate variables, in which temporal trends display 'reversals' (e.g. from greening to browning) (Pettorelli et al 2012, Wei et al 2018, Kalisa et al 2019.
Furthermore, seasonal timescales, the short rains (November-December) and long rains (March-May), are also important when evaluating East African vegetation responses to climate change, particularly precipitation. Nicholson et al (2017) highlighted key differences in East African atmospheric circulation between these seasons, which are subject to different modes of externally-forced variability (e.g. pronounced ENSO impacts on short rains) with distinct responses to climate change. Even within the long rains, March, April, and May exhibit substantial differences in the governing modes of rainfall variability, relatively low spatial coherence, and weak correlations in precipitation across the three months (Camberlin and Philippon 2002, Camberlin and Okoola 2003, Camberlin et al 2009, Nicholson 2015, Nicholson et al 2017. Thus, vegetation responses may also be distinct on monthly timescales.
We herein investigate vegetation changes within Serengeti National Park (SNP), a highly-protected ecosystem where land use activities are relatively restricted per IUCN Category II (Stolton et al 2013, Kija et al 2020, allowing us to better identify potential climate-driven trends. The SNP supports one of the world's 'great migrations' of herbivores following seasonal vegetation availability . Increasing climate change and variability, in addition to other anthropogenic pressures, may disrupt key ecosystem processes related to vegetation growth and phenology, with subsequent impacts on both wildlife and human populations (Ritchie et al 2008, Sinclair et al 2014, Byrom et al 2015, Dublin and Ogutu 2015, Kilungu et al 2017, Hunninck et al 2020. There is uncertainty, however, regarding the presence and timing of SNP vegetation trends. Some studies indicate annual NDVI declines since 1982 across the Serengeti-Mara ecosystem (Pettorelli et al 2012, Kalisa et al 2019, while others suggest nonlinear trends, with annual browning since the early 2000 s (Wei et al 2018). In addition, intra-annual variability of SNP-averaged net primary production may be increasing (Pettorelli et al 2012), possibly as a result of increased rainfall variability across the greater Serengeti-Mara ecosystem surrounding the SNP (Ogutu et al 2008, Bartzke et al 2018. These sitebased studies also show increased dry season rainfall alongside annual rainfall declines since the early 2000s, possibly a response to rising Indian Ocean sea surface temperatures (Williams and Funk 2011). Notably, grazing pressure at the edges of the Serengei-Mara ecosystem may indirectly reduce maximum leaf area within the park boundaries, potentially changing vegetation-precipitation interactions (Estes et al 2012, Probert et al 2019, Veldhuis et al 2019, Kija et al 2020. We therefore acknowledge that the SNP ecosystem is not fully isolated from land use effects. Nevertheless, there remain outstanding research gaps to characterize long-term vegetation trends across the SNP at seasonal and monthly scales, and to identify how changing, interactive climate variables contribute to these trends (Adole et al 2016). Identifying these vegetation trends requires rigorous statistical methods, for prior work using simplified techniques (Chen et al 2019) has likely overestimated the prevalence of global greening (Cortés et al 2021). It is also unclear if there are nonlinear vegetation responses (i.e. trend reversals) in the SNP, similar to those appearing across greater East Africa. Here, we fill these research gaps with a novel investigation of spatially-explicit, time-varying, and sub-seasonal SNP vegetation trends in response to changing precipitation, soil moisture, and temperature. Our methods provide a robust, rigorous approach to evaluate vegetation trends, and account for the lagged interactions between moisture and temperature variables in their contributions to these trends. More specifically, we apply ensemble empirical mode decomposition (EEMD) along with principal component analysis (PCA) and regression techniques to a widely-used satellite-derived leaf area index (LAI) product to answer the following research questions: (1) what types of vegetation trends (monotonic greening/browning or nonlinear) are occurring at annual, seasonal, and monthly scales across the SNP? and (2) how do climate variables and their interactions contribute to these SNP LAI trends?

Study area: SNP
The SNP is situated on the border of Kenya and Tanzania (figure 1), covering 14 763 km 2 at altitudes of 1000-1800 m . A west-toeast precipitation gradient partly shapes transitions in vegetation classes between acacia-dominated woodlands in the west and south, grasslands to the southeast, and grass and shrublands in the east and north (figure 1). Consequently, LAI (figure A1(a)), the ratio of green leaf area to unit ground area (m 2 /m 2 ), ranges from 0.75 for grass and shrubland regions to 2.5 in savanna areas with a mix of grass and trees. Figure A1 also displays annual average (2003-2016) spatial patterns of soil moisture, precipitation, and surface temperature to further contextualize regional LAI.
The SNP's annual precipitation cycle is characterized by the shifting Intertropical Convergence Zone, which results in two general rainy seasons : the short rains (November-December) and the long rains (March-May) (figure 2(c)). Climatic influences on precipitation can also be distinct on monthly timescales (Nicholson et al 2017), however, and thus we evaluate SNP vegetation and climate trends for each month. The seasonal soil moisture cycle (figure 2(b)) is consistent with rainfall, with a slight decrease in January-February and a peak around April. Regional SNP LAI increases at the beginning of short rains and stays relatively high until May, at which point it declines into the dry season (figure 2(a)), although we note there is much interannual variability in the timing and spatial distribution of green-up and maximum SNP LAI (figure C3).

Data
We evaluate vegetation change in the SNP using satellite-derived estimates of LAI, as it provides a more direct physiological interpretation of vegetation growth than the NDVI, which is a proxy for greenness (Cook and Pau 2013). LAI is positively and mostly linearly related to net primary productivity (NPP) (Goward andN 1985, Li et al 2017), although this linearity diminishes for dense canopies. We use monthly mean LAI from the global LAI3g product (Zhu et al 2013) (table 1), which spans July 1981 to December 2016 and is derived from the third generation GIMMS NDVI3g data set at 1/12 degree resolution, 15-day composites. A set of neural networks is used to predict LAI from NDVI, supervised by an improved version of Terra Moderate Resolution Imaging Spectroradiometer LAI (MODIS BNU, Yuan et al (2011)) for the overlapping period 2000-2009 period. The trained neural networks are then used to infer LAI for the remaining years. LAI3g compares well with ground LAI measurements and other satellite-based estimates (Zhu et al 2013) and has been used widely in many previous ecosystem studies (Cook and Pau 2013, Zhu et al 2016, Chen et al 2019. We also conduct spatial and temporal correlation analyses using the MODIS LAI as a 'benchmark' to assess the suitability of LAI3g for this study and its similarity to other LAI products (table 1).

Evaluating time-varying changes to seasonal LAI
To identify the magnitude and direction of any intrinsic vegetation trend, we seek to remove the periodic components, e.g. due to climate variability such as ENSO, from the overall LAI signal. Moreover, we do not limit our trend identification to linear trends only, as nonlinear trends may capture important vegetation-climate interactions. Thus, we decompose the underlying vegetation trend using an established method: EEMD (Wu and Huang 2009). EEMD has been employed in previous studies to identify key modes of vegetation variability, including in East Africa, particularly in relation to precipitation variability (e.g. Hawinkel et al 2016, Zhu et al 2016, Wei et al 2018. EEMD has been shown to be more effective in trend identification than Mann-Kendall tests, especially if the trend magnitude is small compared to other components of the signal (Sang et al 2014). Its core algorithm, empirical mode decomposition (EMD) (Huang et al 1998), decomposes the LAI time   series into a finite number of oscillatory, intrinsic mode functions (IMFs), and a smooth residual term: We hereafter use 'Trend' to refer to the decomposed Trend component from our EEMD specifically ( figure 3(a)). To obtain a more robust decomposition, EEMD ensembles multiple trials of EMD, each of which decomposes the original time series perturbed with white noise. We perform EEMD on SNP LAI for each month at each pixel (196 grid cells, 35 years) using 100 trials with a noise amplitude 0.05 of the absolute signal amplitude. We decompose at most three IMF modes and, similar to Wei et al (2018), categorize the remaining Trend into: strictly greening or strictly browning (i.e. monotonic Trend); greeningto-browning or browning-to-greening (i.e. nonlinear or inflected Trend); and undetermined (i.e. more than one extrema). To determine if strictly greening or browning montonic Trends are significant, we assess whether they are statistically different from random noise that preserves the original LAI distribution and autocorrelation structures (i.e. surrogates), measured by a test statistic. We use iteratively refined amplitude adjusted Fourier transform method (IAAFT, appendix B) (Schreiber and Schmitz 1996) to generate these surrogates for monthly LAI series. We choose the time reversibility (TR) as our test statistic: where X t is the value at time t, and n is the number of time steps. If a significant monotonic Trend is present, the TR magnitude will be large since reversing the data in time will differ substantially from the original trend. We performed a one-sided test at the 5% significance level by generating 19 surrogates from the original monthly LAI and extracting their Trends with EEMD. The LAI Trend is significant when its TR statistic magnitude is greater than the TR statistics from all the other surrogates' trends. By testing on the LAI Trend alone, this procedure avoids false positive trends found in other studies of global greening (Cortés et al 2021), which may result from longer-term modes of ecosystem variability. Details on the surrogate generation and significance test can be found in appendix B. Nonlinear Trends may constitute another interesting response in which the LAI Trend has 'reversed' over time. Such  Lastly, we evaluate how the Trend contributes to the overall LAI signal, i.e. how much of the total LAI variance the Trend explains, using the variance contribution rate (VCR), a standard metric used with EEMD approaches in previous climate studies (e.g. Bai et al 2015, Guo et al 2016, Liu et al 2020: where σ 2 represents the variance for the IMF or the trend component.

Evaluating the contribution of changing climate variables to SNP LAI change
We further investigate the contributions to the SNP LAI Trends by key regional climate variables (figure 2, table 1): root zone soil moisture from the global land evaporation Amsterdam model (GLEAM) (Martens et al 2017); precipitation from the climate hazards group infrared precipitation with stations (CHIRPS) (Funk et al 2015); and surface temperature from CRU Surface Temperature Version 4 (CRU TS) (Harris et al 2020).
We first conduct the above EEMD analysis on each of the climate variables (figure 3(a)). Due to the higher variability present in precipitation and soil moisture, we use a maximum of four IMF modes for all the climate variables. We then perform the same Trend categorization to identify systematic patterns that may relate to the LAI Trends.
Next, we identify months with spatially-uniform Trends (e.g. domain-wide browning) from EEMD analysis. For these month(s), we then take the SNP area-average for the time-varying (1982-2016) climate variables and LAI, with intention of regressing the LAI timeseries on those of the climate variables. However, multicollinearities likely exist between soil moisture, precipitation, and temperature, and their lags, which complicate interpreting their effects as independent variables in a regression. To mitigate these effects, we first apply a PCA to these climate variables, inclusive of their values up to two months prior to capture potential memory effects (e.g. in moisture) and lagged ecosystem responses. We then use the first two principal components (PCs), also interpreted as explanatory climate components of regional climate variability and change (figure 3(b)), in our LAI regression. The PCA preprocessing allows us to: (1) identify independent climate components for our regressions, and (2) reveal how these components are comprised of the interacting climate variables. The latter is important because recent work (Adole et al 2016(Adole et al , 2019 highlight how interactions between moisture and temperature mediate vegetation responses across Africa, even in water-limited domains. We then regress LAI on these PC climate variables to investigate their contributions to the LAI Trend and the overall LAI signal (equation (4)).
Specifically, to assess the importance of the respective PC1 and PC2 climate Trends to the overall LAI timeseries, we further fit the PC regression by using EEMD to remove the Trend component from PC1 while keeping PC2 fully intact (equation (5)) and vice versa (equation (6)): Trend categorization: Brown = monotonic browning; GtoB = reversal from greening to browning; BtoG = reversal from browning to greening; Green = monotonic greening; NaN denotes the remaining Trend has more than one extrema and cannot be classified in any of these categories; overplotted on some grid cells is black hatch lines, denoting that our surrogate test indicates the monotonic greening or browning Trend is significant at 0.05 level.
Finally, we evaluate how Trends from these climate variables affect the LAI Trend using the standard multiple EEMD-Regression (Yang et al 2011, Masselot et al 2018:

Time-varying changes to seasonal LAI
From 1982 to 2016, SNP LAI exhibits slight greeningto-browning Trend reversals and monotonic greening Trends at the annual and wet season temporal scales (figure E2), respectively. However, climate change and its interaction with modes of climate variability may have distinct impacts at monthly timescales (Nicholson et al 2017), which could further impact vegetation responses. Therefore, we further examine time-varying SNP LAI Trends for each month using the EEMD approach described in section 2.3 (see figure 3(a)). Firstly, there exists large spatial variation in the LAI Trend categorization across months, particularly during the dry season (June-October) (figure 4). Moreover, there is a lack of widespread significant monotonic Trends (indicated by black hatch lines) across the domain and months. Nevertheless, browning-to-greening and greening-to-browning Trend reversals are spatially widespread during November-December in the short-rains (figures 4(k) and (l)) and in March-April of the long-rains (figures 4(c) and (d)), respectively. These responses differ markedly from those indicated by a simple Mann-Kendall test (figure C1), as has been used in prior work (e.g. Chen et al 2019, Cortés et al 2021), which show more widespread greening trends (especially during the wet season from January to April) and significant gridcells.
Furthermore, while there are few significant monotonic Trends in monthly mean SNP LAI, the intra-annual LAI variability, i.e. the relative range defined as RREL = (max-min)/mean (Pettorelli et al 2012), does show widespread and significant monotonic increases (black hatching, per a Mann-Kendall test) over the 1982-2016 period. This measure captures the LAI amplitude during the annual and wet season cycles. This enhanced intra-annual variability is particularly concentrated in the southern and eastern areas (figure C2) of the domain which are broadly dominated by herbaceous vegetation and shrubs. Monthly plots of GLEAM soil moisture Trend categorization: Dry = monotonic drying; WtoD = reversal from wetting to drying; DtoW = reversal from drying to wetting; Wet = monotonic wetting; NaN denotes the remaining Trend has more than one extrema and cannot be classified in any of these categories.
Yet, the timing of LAI maximum and minimum still displays large interannual variability without clear trends (figure C3). We provide context for this response in section 4, although a quantified attribution is beyond the scope of our analysis.
We also examine the timing of the prevailing LAI Trend reversals (figures C4 and C5): Greening-tobrowning Trend reversals in the long rains season emerge around 2001-2002 (darker brown grid cells on March and April maps in figure C6 and the reversal points in the area-aggregated figures C5(g) and (h)). In contrast, the browning-to-greening that characterizes much of the short rains occurs in the late 1980s/early 1990s (figures C5(k) and (l)). We note that for the short rains, November and December, the VCRs (or the amount of total LAI variance the Trend captures) are relatively small (2.84%, 0.91%). March, in particular, exhibits the most uniform greeningto-browning across the domain with a comparatively higher 4.28% VCR.

The contribution of changing climate variables to SNP LAI change
The EEMD Trend categorization for soil moisture (figure 5) and precipitation (figure 6) reveal distinct responses during the individual months in the wet season, similar to the LAI results.
In November and December, the soil moisture and precipitation Trends (figures 5(k) and (l)) and 6(k) and (l)) show a drying-to-wetting and monotonic wetting Trend, respectively, which is consistent with prior work (Nicholson et al 2017, Bartzke et al 2018. These moisture Trends are broadly consistent with browning-to-greening during these months (figures 4(k) and (l)). We also note increased October precipitation (figure 6(j)) and drying-to-wetting Trends in soil moisture (figure 5(j)), preceding the November/December LAI responses. The short rains Trend reversal for soil moisture occurs near the beginning of the time series (figure C7), similar to the timing of the browning-to-greening Trends in LAI ( figure C6). Qualitatively, this suggests that the LAI increases during short rains (figures C5(g) and (h)) may be related to the largely co-occurring increases in moisture availability that begin early in the 1982-2016 period. We note however that the soil moisture dataset is more coarsely spatially resolved compared to LAI3g dataset, which limits comparative interpretations of reversal timings between these variables.
The potential relationships between LAI and the moisture variables appear more complicated during the long rains. There are widespread wetting-todrying Trends in March for soil moisture ( figure 5(c)) that are consistent with the greening-to-browning Figure 6. Monthly plots of CHIRPS precipitation Trend categorization: Dry = monotonic drying; WtoD = reversal from wetting to drying; DtoW = reversal from drying to wetting; Wet = monotonic wetting; NaN inside the SNP polygon denotes the remaining Trend has more than one extrema and cannot be classified in any of these categories.
LAI Trends (figure 4(c)). However, the wet season precipitation Trends are less temporally coherent with these LAI and soil moisture Trends, particularly from January to April. March (figure 6(c)) largely shows monotonic increases in precipitation throughout much of the domain, while January and February (figures 6(a) and (b)) show monotonic declines and wetting-to-drying reversals in precipitation, respectively. In contrast, soil moisture largely shows increasing Trends over nearly the entire SNP in January and much of the domain in February (figures 5(a) and (b)). Qualitatively, this may indicate a potentially lagged relationship between precipitation and LAI Trends in the long rains, as mediated through soil moisture.
Given this response, we evaluate more explicitly the possible (lagged) relationships between the climate variables and LAI using PC-Regression and EEMD-Regression techniques (section 2.4, figure 3(b)). We focus our regression analyses on better characterizing the March LAI response during the long rains, motivated by its spatially widespread greening-to-browning Trend reversals and relatively high VCR (figure 4(c)).
Our PC analysis ( figure 3(b)) shows that PC1 loads heavily on soil moisture and precipitation while PC2 loads more strongly on temperature (table D1). We therefore interpret these as the 'moisture' and 'temperature' components (of regional climate variability and change), respectively. Furthermore, moisture variables in January and February display relatively strong correlations with March LAI (figure D3), supporting a lagged relationship between moisture availability and early long rains green-up, and are also important components of the PCs (table D1). PCregression results (tables 2 and 4) indicate both PC1 and PC2, i.e. the moisture and temperature components combined, are significant and explain March LAI moderately well (adj. R 2 = 0.6). Sensitivity analyses (figure 7(a)) show that the detrended moisture component (table 2, equation (5)) yields a worse model than the detrended temperature component (table 2, equation (6)), evidenced by lower adjusted R 2 and worse AIC/BIC scoring criteria. This suggests that the moisture Trend is more important than the temperature Trend for predicting the overall March LAI. Nevertheless, Trends belonging to both the moisture and temperature components are highly important to explain the greening-to-browning Trend reversal in March LAI, evidenced by the good fit in EEMD Trend Regression (adj. R 2 = 0.975) and their statistical significance (table 2, figure 7(b)).

Discussion
Our results show no linear declines in 1982-2016 annual SNP LAI ( figure E1(a) specifically, a browning-to-greening in the short rains (November-December) starting in the late 1980s/early 1990s and a greening-to-browning in March-April in 2001-2002. Furthermore, when measured from their inflection points, the approximately linear browning trend in March is around 0.1 m 2 /m 2 /year (figure C5(g)) and the November/ December greening (figures C5(k)/(l)) is 0.08-0.1 m 2 /m 2 /year. While we did not find similar monthly analyses against which to compare these Trends, they are generally comparable in magnitude to the global, growing-season trends reported by other studies, such as approximately 0.07 m 2 /m 2 /year reported by Zhu et al (2016).
These findings demonstrate the utility and importance of high-resolution, spatially-explicit, sub-seasonal and monthly analyses for assessing climate change-ecosystem interactions, particularly in protected areas where land use change is, ostensibly, minimal (Pelkey et al 2003, Nicholson et al 2017. Our EEMD and surrogate testing approach captures nonlinear Trends without the restrictive assumptions of the Mann-Kendall Test (e.g. serial-independence), which has been previously used to assess vegetation trends in Africa and globally (Pettorelli et al 2012, Chen et al 2019. Importantly, our approach is a robust procedure that minimizes false positive trends, which have resulted from simpler approaches to evaluate vegetation greening effects (Cortés et al 2021).
To identify climate influences on mean SNP LAI Trends, our combined PCA and EEMD regressions move beyond simply correlating LAI and climate variables, and are unique in explicitly capturing interactions between climate variables and their lagged SNP LAI impacts. A novel finding of this approach is the importance of both SNP moisture and temperature Trends to the SNP March LAI Trend, as well as the interaction between these climate variables (and their lags) evidenced by their loading in both PCs. Increasing water limitation in March across the SNP, i.e. soil moisture declines that lag declining January and February precipitation Trends, and rising temperature Trends across these months all contribute to March browning. To our knowledge, this is the first analysis of how interacting climate variables, both in their monthly lagged effects and Trends, contribute to subseasonal vegetation variability and Trends in SNP.
These findings align with previous work suggesting that interactions between moisture, temperature, and other physical variables shape intra-seasonal vegetation phenology across Africa writ large (Ryan et al 2017, Adole et al 2018. The browningto-greening during the short rains likely results from increasing precipitation that enhances soil moisture, which we also show (figures 5 and 6) and has previously been reported around the SNP (Ogutu et al 2008, Bartzke et al 2018, although these studies did not focus on subsequent impacts to vegetation. Notwithstanding important moisture-temperature interactions, our findings that moisture availability and memory are strong determinants of SNP vegetation responses is generally consistent with our understanding of water-limited African ecosystem dynamics (Zhang et al 2005, 2006, Wei et al 2018, Adole et al 2019.
In addition, our findings of increasing intraannual LAI variability from 1982 to 2016 is consistent with NPP analyses by Pettorelli et al (2012), averaged for the SNP domain. Our spatially-explicit analysis further reveals that this increasing variability occurs in predominantly herbaceous vegetationand shrub-dominated eastern and southern SNP. As this study primarily focused on identifying mean LAI Trends and their drivers, attributing these long-term changes in intra-annual LAI variability is beyond our scope. We nevertheless suggest that increasing precipitation variability is likely to play a role. For example, prior studies suggested that regional precipitation may be displaying more intra-and interannual variability (Bartzke et al 2018), partly relating to changing Indo-Pacific sea surface temperatures and subsequent impacts to regional atmospheric circulation and moisture transport (e.g. Williams andFunk 2011, Nicholson et al 2017). Future work using process-based Earth system and/or ecosystem models (e.g. Liu et al 2015) that include (African) vegetation sensitivities to a range of climate variables (Adole et al 2018(Adole et al , 2019 could better isolate natural variability versus anthropogenic versus forcings to attribute change in multiple SNP vegetation characteristics. Global and larger-scale regional studies of greening and browning may obscure or even overestimate important vegetation trends in protected areas, where climate change poses an increasing threat (Mondal et al 2020). Our results augment the increasing body of work on East African ecosystems to show that SNP vegetation Trends are nonlinear, and thus require regular monitoring and rigorous evaluation at fine spatial and temporal scales to tease out the many sources of variability, inclusive of climate change. Indeed the long-term monitoring and evaluation of vegetation (and ecosystem) health may be just as important for biodiversity and conservation goals, e.g. Sustainable Development Goal 15: Life on Land, as the amount of area protected , UN 2015, Elsen et al 2020. While the monthly LAI Trends we identify account for a relatively small portion of the overall inter-annual vegetation variability, our results nevertheless reveal contrasting Trends between the short and long rains. Along with enhanced intra-annual LAI variability, these differences may influence ecosystem phenology and dynamics (Scheffers et al 2016, Elsen et al 2020 as 'slow changes and trends' . For example, monthly vegetation Trends and their differences between the short and long rains may incite shifts in the timing and locations of herbivore migration, contributing to changing food web dynamics (Ritchie et al 2008).
We further note limitations to this study that complicate a comprehensive attribution of these LAI Trends and climate variables. Firstly, prior work comparing different LAI products demonstrates that different Trends and responses may emerge depending upon the product used (Zhu et al 2016). Therefore, evaluating agreement across products is critical to assessing vegetation changes, particularly in protected ecosystems with relatively small area footprints (∼15%) and where widespread, frequent ground monitoring is not always available.
Secondly, we note that our two-step approach of PCA followed by EEMD-Regression does not exclude possible multicollinearities between the IMFs of each PC, which could limit the statistical test power on the covariates in EEMD-Regression. One might consider a more flexible model that does not impose assumptions on independent covariates and better capture their interactions. Regardless, advanced statistical methods to identify climate change impacts on SNP vegetation should focus on interpretability while also capturing important interactions between climate variables.
Lastly, previous studies have documented changes across a host of other important abiotic and biotic determinants of Serengeti ecosystem dynamics. Among the most impactful is increasing human land use around the greater Serengeti-Mara ecosystem, which includes the SNP, that may threaten SNP ecosystem dynamics and functionality (Sinclair et al 2008, Li et

Conclusions
SNP is a critical hotspot of biodiversity, characterized by complex ecosystem dynamics and intensifying human pressures, including land use and climate change. We conducted rigorous analyses using EEMD to identify sub-seasonal, spatially-explicit Trends in widely-used remote sensing datasets of leaf area, precipitation, temperature, and model-produced soil moisture in the SNP from 1982 to 2016. We further employed novel PCA and regression techniques to determine how temperature and moisture variables, and their lagged interactions, contribute to these trends. We found long-term SNP Trends to be nonlinear, with spatially-uniform greening-tobrowning during March and browning-to-greening during November-December. While the SNP ecosystem is water-limited, both moisture and temperature Trends, and their values from January-March, were critical to explain March greening-to-browning vegetation Trends. Our findings demonstrate the importance of evaluating long-term ecosystem trends in the SNP at monthly timescales, when vegetation responses can differ substantially, and of including multiple, interactive climate variables in attributing these vegetation changes.

Data availability statement
No new data were created or analysed in this study.

Appendix A. Main manuscript section 2.1: Serengeti National Park (SNP)
Average annual precipitation in the SNP is highest in the northern and western regions ( figure A1(c)). The spatial pattern of average annual soil moisture is consistent with precipitation ( figure A1(b)), highest in the northern SNP corresponding well with the regional maximum in LAI. Annual average temperatures are less variable over the domain ( figure A1(d)). The highest temperatures occur in the central and western areas, while slightly cooler conditions prevail in the east and north.
For completeness, we compare LAI3g with MCD15A3H Version 6 Moderate Resolution Imaging Spectroradiometer Level 4 LAI (MODIS LAI) (Myneni et al 2015), GLASS LAI From Time-Series MODIS Surface Reflectance (GLASS) (Liang et al 2013), and NOAA Climate Data Record of AVHRR LAI (CDR) (Claverie et al 2016). MODIS LAI is a 4-day composite data set with 500-meter pixel size. It chooses the best pixel available from all the acquisitions of both MODIS sensors located on NASA's Terra and Aqua satellites from within the 4-day period. GLASS is an 8-day composite data with 0.05-degree resolution. It generates one-year LAI estimates by fitting general regression neural networks on one-year MODIS Surface Reflectance data, which improves quality by considering the temporal pattern (Xiao et al 2017). CDR is a daily composite data with 0.05degree resolution, derived from the NOAA AVHRR Surface Reflectance product. CDR estimates LAI values globally over land surfaces, but not over bare or very sparsely vegetated areas, which results in several grid cells of missing data within our defined SNP domain. Brief descriptions and references for all four LAI products considered for our analyses are given in table 1. MODIS LAI and CDR datasets are queried from Google Earth Engine while LAI3g and GLASS are downloaded from source webpages.
To evaluate agreement between LAI products, we first compare their temporal patterns shown in figure A2. The mean values for each data product time series are calculated by averaging all the grid cells with centroids overlapping the SNP polygon at the respective time interval. Although MODIS LAI only has data available from 2002 and such short  time frame is not suitable for our long-term trend study, it is a high quality dataset with the finest spatial resolution (of those products considered here) and relatively less measurement uncertainty (Fang et al 2013). Therefore, we conduct a set of spatial and temporal correlation analyses using the MODIS LAI as a 'benchmark' to assess how similar these LAI products are both to MODIS LAI and to each other, during their overlapping periods. We employ the Spearman rank correlation test here as it is suitable for Figure A3. 1000 sampling points from SNP, used to calculate the spatial correlation among LAI products. nonlinear data, and this analysis enables us to identify which LAI product may be best suited for our analyses to detect long-term trends changes.
To standardize our correlation analysis in the temporal dimension, we first compute the monthly mean for each product (16 years × 12 months) spatially aggregated over the domain for each point in time. Given the variation in geospatial resolutions of the LAI datasets used, we adopt the following sampling strategy for assessing the products' spatial  figure A3) uniformly from the SNP domain using 'rejection sampling' . Specifically, we define a rectangular boundary by the minimum and maximum of x-axis and y-axis in SNP polygon, and then generate a point within this rectangle selecting those that specifically fall within the SNP polygon. Points are generated until the total number of points equals N. The spatial correlation is then computed using the N sampling points for the annual or seasonal mean. We find that N = 1000 is a satisfactory number of sampling points, after which the improvements in correlation cease and we risk duplicating sampling points in coarser products (such as LAI3g). The spatial distribution of these 1000 sampling points, shown in figure A3, does well to represent the SNP domain. Figure A4 and table A1 present the spatial and temporal correlations results. Temporally, MODIS LAI is most closely correlated to CDR and LAI3g, but less so with GLASS. However, we note that CDR has several missing values over the SNP domain, which may complicate its use in this analysis. Spatially, however, MODIS LAI is more correlated to LAI3g, which is expected given that algorithms used to create the LAI3g dataset were calibrated and trained with respect to MODIS LAI (Zhu et al 2013). Given these considerations and correlative findings, we use LAI3g throughout our analysis to answer our research questions.

Appendix B. Main manuscript section 2.3: evaluating time-varying changes to seasonal LAI
We describe our surrogate hypothesis testing method for the nonlinear, monotonic EEMD Trend using iteratively refined amplitude adjusted Fourier transform method (IAAFT, (Schreiber and Schmitz 1996)).
Surrogate hypothesis testing (Theiler et al 1992) is a statistical approach to evaluate nonlinearity in timeseries data. Time series data is typically composed with long-term nonlinear determinism, shortterm linear correlations, chaos and noise. To test for the significance of nonlinearity, the method first specifies some linear stochastic process as a null hypothesis, then generates surrogate data sets based on this null hypothesis, and finally computes a discriminating test statistic for the original and for each of the surrogate data sets. If the value computed for the original data is significantly different than those computed for the surrogate data, then the null hypothesis is rejected and significant nonlinearity is detected.
To test on nonlinear, monotonic LAI Trends, we need to specify our null model (which may be due to other sources of climate variability or other ecosystem dynamics). We assume the time series noise is generated from a linear stochastic process, and measured through a potentially nonlinear function. The null hypothesis is that the underlying process is from linear noise, while the alternative being a nonlinear deterministic trend.
We then apply IAAFT to generate surrogates using Fourier transform. Recall the key relationship between the time domain and the frequency domain: The linear correlations in the time domain are represented by the power spectrum amplitudes in the frequency domain, while the nonlinear correlations are captured in the phase angles. In light of this, the algorithm first computes the Fourier transform of the signal. Next, it randomizes the phase angles of the Fourier transformed signal so as to destroy any original nonlinearity, while preserving the amplitudes that capture original linear signals. It then reconstructs the surrogates through Inverse Fourier transform. Lastly, it corrects the spectrum and the distribution of the surrogate to align with the original signal, so as to avoid mixing potential nonlinear effects from the measurement function. In the end, IAAFT generates a surrogate that has same linear property as the original signal while randomizes the nonlinearity.
After generating the surrogates, we test the significance of the nonlinearity of the LAI Trend versus its surrogates by computing the test statistic. For a one-sided test, if the test statistic of LAI Trend is the maximum or minimum out of all K surrogates, then it yields a p-value of 1/(1 + K). Here, we fix our p-value as 0.05 and thus generate K = 19 surrogates (Schreiber and Schmitz 1996).

Appendix C. Main manuscript section 3.1: changes and trends in interannual SNP LAI
For comparison, we perform Mann-Kendall trend test on SNP LAI. Figure C1 shows the monthly trend slope from 1982 to 2016 per grid cell over SNP. Contrary to the EEMD Trends, the Mann-Kendall trends display quite different patterns, with a predominantly greening trend in the wet season and more significant grid cells.
To identify the spatially-explicit changes in intraannual LAI variability, we measure annual and wet season LAI variability using relative range: RREL = [MAX − MIN]/MEAN, where MAX, MIN, and MEAN are the annual or wet season maximum, minimum and mean LAI for each year during 1982-2016 (Pettorelli et al 2012). We then test for significant monotonic trend in the annual/seasonal LAI RREL using a Mann-Kendall trend test (MK test) (Kendall 1948, Pettorelli et al 2012. We compute RREL over the 1982-2016 period for the annual and wet season (November-May) (figure C2). The LAI RREL has been mostly increasing and a number of grid cells show significant monotonic trends. Earlier studies have also suggested domain-averaged increases SNP vegetation variability, (Ritchie et al 2008, Pettorelli et al 2012. Our results augment these previous findings to show that the increase in variability is predominantly in the eastern, central and south, parts of the domain, which are mainly characterized by grasses and shrubs. To further evaluate SNP vegetaion variability, we compute the empirical distributions of the maximum Figure C2. Intra-annual variability of LAI3g measured by Mann-Kendall trend test on relative range (RREL: (max-min)/mean in the given season). Green grid cells indicate positive slope from Mann-Kendall trend test, while brown grid cells indicate negative slope. Black hatch lines mark out grid cells with significant trend at 5% level. For each year, we calculate the month that achieves peak LAI per pixel, and then aggregate them to compute the empirical distribution. Each grid cell represents the empirical probability of this month being the maximum out of all 12 months.
LAI timing from 1982 to 2016. As shown in figure C3, the timing of peak LAI in SNP varies drastically in different years.
To evaluate the validity of monthly time-scale of LAI response, we use bi-weekly LAI3g data to perform EEMD and Trend categorization during the wet season. As shown in figure C4, the Trend categorization remains largely the same between the first half and the second half of each month during the wet season, except May. These bi-weekly patterns are Figure C4. Bi-weekly plots of LAI3g Trend categorization: Brown = monotonic browning; GtoB = reversal from greening to browning; BtoG = reversal from browning to greening; Green = monotonic greening; NaN denotes the remaining Trend has more than one extrema and cannot be classified in any of these categories. Top (bottom) row shows the decomposed Trend from the first (last) two weeks of each month.   highly consistent with the monthly patterns shown in figure 4. Apart from the pure data-driven evidence, monthly response carries more distinct climate signals while finer time-scales could introduce more short-term climate noises. Hence, we rely on monthly responses in our analysis. Consistent with figure 4, the trends on the spatial average of LAI over SNP domain show the similar reversal patterns. The VCR analysis suggests that the intrinsic trends are relatively subtle compared to the overall variability of LAI, which are largely driven by external forces such as climate factors.
Monthly maps (figures C6-C8) below indicate the timing of the reversal point revealed for each of the variables, LAI, soil moisture, and precipitation, using our EEMD analysis 7 . Gridcells shown are only for inflected trends, not for monotonic trends. Browning-to-greening reversals are denoted by positive values and greening-to-browning reversals are denoted by negative values. The absolute value of the pixel indicates the year of reversal in a numeric sequence, while the darker the color, the more recent      Figure D4. Projection of the climate variables data matrix (35 × 9) to the two-dimensional subspace spanned by the first two PCs. Points colored by LAI_Mar values.  (6)).
Co-variate Coefficients P-value PC1 0.160 0.000 PC2_IMF −0.066 0.041 Figure D5. EEMD Decomposition of LAI_March (a), its PC1 (b) and PC2 (c) from the driver matrix. Interestingly, we can see PC1_Trend and PC2_Trend are very similar. loads more strongly on temperature. We also note that Precip_Feb contributes strongly and negatively to PC2. We further illustrate the nine driver variables with 35 data points (i.e. 1982-2016) projected to the first two PCs, color-coded by March LAI values (figure D4). March LAI variation is mainly captured in the PC1 direction. Tables D2-D4 show the results of PC-Regression (original), PC-R with detrended PC1, and PC-R with detrended PC2 respectively. Figure D5 shows the EEMD for March LAI, its PC1 and PC2 from the driver matrix. Table D5 shows the EEMD-Regression on Trend component: Trends from both PCs are highly significant. However, the high VIF factors indicates the multicollinearity issue between the covariates, and thus their coefficient values are biased due to confounding effects (Graham 2003).

Appendix E. Other supplemental figures
Below shows the full ensemble empirical mode decompositions for LAI in the SNP (using LAI3g), for annual/seasonal values (figure E1) and for each month (figure E2).