Decoupling of Ecological and Hydrological Drought Conditions in the Limpopo River Basin Inferred from Groundwater Storage and NDVI Anomalies

: Droughts are projected to increase in intensity and frequency with the rise of global mean temperatures. However, not all drought indices equally capture the variety of inﬂuences that each hydrologic component has on the duration and magnitude of a period of water deﬁcit. While such indices often agree with one another due to precipitation being the major input, heterogeneous responses caused by groundwater recharge, soil moisture memory, and vegetation dynamics may lead to a decoupling of identiﬁable drought conditions. As a semi-arid basin, the Limpopo River Basin (LRB) is a severely water-stressed region associated with unique climate patterns that regularly affect hydrological extremes. In this study, we ﬁnd that vegetation indices show no signiﬁcant long-term trends (S-statistic 9; p -value 0.779), opposing that of the modeled groundwater anomalies (S-statistic -57; p -value 0.05) in the growing season for a period of 18 years (2004–2022). Although the Mann-Kendall time series statistics for NDVI and drought indices are non-signiﬁcant when basin-averaged, spatial heterogeneity further reveals that such a decoupling trend between vegetation and groundwater anomalies is indeed signiﬁcant ( p -value < 0.05) in colluvial, low-land aquifers to the southeast, while they remain more coupled in the central-west LRB, where more bedrock aquifers dominate. The conclusions of this study highlight the importance of ecological conditions with respect to water availability and suggest that water management must be informed by local vegetation species, especially in the face of depleting groundwater resources.


Introduction 1.Background
Global warming will lead to unprecedented hydrological extremes this century due to the increasing water holding capacity of the troposphere [1,2].This is especially concerning in warm, dry climates, where evaporative demand is naturally high and home to vulnerable populations as defined by international socioeconomic standards [3].The Limpopo River Basin (LRB) is one such region as a critical transboundary watershed for more than 18 million people across four countries: Botswana, Mozambique, South Africa, and Zimbabwe.Characterized by a predominantly semi-arid climate, the LRB receives about 530 mm of annual rainfall, which must be distributed across pre-existing water demand between its agricultural, mining, and urban sectors [4][5][6].Drought conditions attributed to climate teleconnection events such as El Niño-Southern Oscillation and the Indian Ocean Dipole regularly exacerbate the water deficit in this closed basin as well [7,8].Although mesoscale weather patterns such as cloud bands and tropical low-pressure systems typical of January through March provide seasonal relief, associated storms also run the risk of causing devastating floods [9,10].Given this dynamic variance in water quantity, additional studies on this basin's hydrologic anomalies are required to better understand the uncertainty of and the watershed's response to wet or dry spells, for the benefit of current and future stakeholders.
A drought is generally defined by the duration and magnitude of a regional deficit in available water and may be additionally categorized as one of the following: meteorological, agricultural, hydrological, or, more recently, ecological [11,12].Depending on the hydrologic component in question, different indices may be calculated to predict and contextualize the drought condition.For example, the Standardized Precipitation Index (SPI) uses a historical precipitation dataset for reference, while the Palmer Drought Severity Index (PDSI) includes a temperature component to inform soil moisture and evapotranspiration conditions, though both indices are used to assess meteorological droughts [11,13,14].With the advent of satellite-based Earth observations, remotely sensed products derived from the Moderate Resolution Imaging Spectroradiometers (MODIS) have also become useful for monitoring droughts by way of determining vegetation and temperature conditions based on the Normalized Difference Vegetation Index (NDVI) and land surface temperature retrievals, respectively-two variables that may be further combined to generate the Vegetation Health Index (VHI), interpreted to be a measure of ecological drought in this study [15][16][17].In recent years, the Standardized Precipitation Evapotranspiration Index (SPEI) has become a well-known meteorological drought index that considers the effects of both available water (precipitation, P) and atmospheric demand (potential evapotranspiration, PET) [18,19].Precipitation plays such a significant role in determining drought conditions because of its temporal relationship to a basin's water balance components.For example, drought-stricken watersheds with flashier hydrographs may still be vulnerable even after intense rainfall events due to low infiltration capacities and high runoff volumes.The basin's response to precipitation is critical in determining hydrological droughts, which may not be captured by rainfall and temperature conditions alone.Subsurface interactions between soil moisture, groundwater recharge, and plant water uptake are therefore important to consider for more comprehensive drought assessments.
Fortunately, the National Atmospheric and Space Administration's (NASA) Gravity and Recovery and Climate Experiment (GRACE) and its successor, GRACE Follow-On (FO), satellite missions have become instrumental in assessing total water storage anomalies (TWSA) around the globe.By inferring the temporal variations in Earth's gravity field as mass changes caused by the movement and storage of water, GRACE observations have been discovered to detect trends in groundwater depletions on the regional scale [20,21].Despite its coarse spatial resolution (>150,000 km 2 ), GRACE has been used to generate TWSA products based on regional mass concentration (mascon) functions at a higher gridded resolution of 0.5 degree.The GRACE mascon data facilitates the assimilation of data into land surface models for the purpose of simulating a seamless terrestrial water budget for drought monitoring and analysis [22].The most recent efforts by Li et al. (2019) [23] incorporate GRACE TWSA into NASA's Catchment Land Surface Model (CLSM) using mascon solutions to develop a global groundwater storage data product.

Motivation
In the context of Southern Africa, only a handful of studies exist that leverage GRACE data products and have primarily focused on total water storage trends, precipitation anomalies, and climate teleconnections [42][43][44][45][46]. Without the inclusion of other hydrologic components, a knowledge gap on long-term interactions between total water storage and its constituents remains for the Limpopo River Basin.Water available to vegetation, or green water, is especially important in comparable climates, where soil moisture may be insufficient during dry spells and root systems may tap into deeper sources of water [47][48][49].The widespread phenomenon of greening in the form of shrub encroachment across the greater region raises further questions about its relationship with water availability in this region [50][51][52].
Although NDVI trends were considered alongside TWSA by Scanlon et al. (2022) [46] for the major aquifers on the African continent, the LRB was not included.The NDVI of the Karoo basin, which is located southwest of the LRB, was determined to have minimal correlation (r = 0.3) with TWSA, while midlatitude African aquifers had higher values with coefficients ranging between 0.6 and 0.9.One overlooked component of their study is aquifer type.Because the LRB is underlain by both alluvial sand (downstream river valley) and fractured crystalline (upstream tributary catchments) aquifers, the recharge potential will vary with both human abstraction and seasonal precipitation [53].Saveca et al. (2021) [54] and Walker et al. (2018) [55] investigated sand river deposits along the Limpopo River in Mozambique and the Molototsi sand river catchment in South Africa, respectively, and found that shallow, alluvial aquifers have significant potential for groundwater abstraction, as they have the added benefits of frequent recharge potential with both aquifer thicknesses ranging between 0 and 6 m and sustained lateral groundwater flow occurring in both the wet and dry seasons through surface runoff.Such aquifers can become fully saturated even after consecutive periods of droughts but require robust management practices to prevent the reduction of perennial flow conditions [54,55].On the other hand, stable isotope analyses by Abiye et al. (2020) [56] show that short-term intense rainfall events (non-seasonal) are crucial recharge events for bedrock aquifers due to the infiltrated water's long residence time, especially with groundwater recharge being a smaller component of the water balance for deep aquifers (2-5% of the annual rainfall) [57].
The primary objective of this study is therefore to determine drought conditions in the LRB by vegetation health, groundwater storage anomalies, and meteorological components and assess their potential interactions.First, a series of the drought indices are calculated using the following datasets: the GRACE-assimilated CLSM outputs, MODIS-based optical and thermal products, and precipitation estimates based on retrievals from the Global Precipitation Measurement (GPM) mission, which are calibrated and processed under the Integrated Multi-satellitE Retrievals for GPM (IMERG) algorithm (Table 1) [23,58,59].Then, both growing season and monthly anomalies are assessed with time series analyses and correlation tests to observe potential trends and relationships.Lastly, the spatial heterogeneity of the indices and variables are compared for the LRB for the growing season based on a reference period of 2004-2022.The variety and long-term record of the selected datasets offer a unique study on droughts and their impact on remotely sensed vegetation in the Limpopo River Basin.

Study Area
The Limpopo River Basin spans a drainage area of about 412,000 km 2 with an elevation band ranging from sea level to above 2000 m.The region may be subdivided into four climate regions: the lowveld and central river valley (arid), middleveld (semi-arid), escarpment (humid subtropical), and coastal plains (wet tropical).A general north-south and west-east rainfall gradient exists following the orography, and ranges between about 200 mm of annual rainfall in the dry and 1500 mm in the wetter areas.Potential evapotranspiration for the basin over open water ranges between 800 and 2400 mm per year, with the highest rates occurring in the river valley and during rainy seasons in the austral summer (October through March), defined as this region's growing season.Mean daily temperatures for the basin range between 0 • C in the winter and 36 • C in the summer [4,6,64,65].
The dominant land cover classes are savannah-grass and shrublands, which include woody trees (deep-rooted, ≥4 m) and shrub species (shallow-rooted) [48].Croplands across the LRB are populated by primarily rainfed maize, sorghum, and millet, as well as some irrigated fields for cultivating tropical fruits such as mangos and bananas year-round."Blue water" irrigation is made possible by both groundwater and surface water abstraction, which is supplemented by reservoirs formed by dams, one of which is the Massingir Dam in the Oliphants sub-watershed [4,[66][67][68][69][70]. Forested regions are made up of both native forests and commercial plantations of exotic species (e.g., pine, eucalyptus) in the escarpment and mountainous region of South Africa [4,66].Studies have also identified a trend in decreasing natural vegetation in exchange for increasing cropland cover in the past decades in the Limpopo [71,72].
Figure 1 shows a land cover map based on 2017 imagery from Sentinel-2.Although no formal land cover accuracy assessment was conducted, a qualitative comparison between it and South Africa's custom 2014 and 2020 maps revealed sufficient overlap [73].

In Situ Datasets
A groundwater level dataset from the South African Department of Water and Sanitation (SA DWS) containing more than 1000 borehole sites matching the South African extent of the Limpopo River Basin was selected to assess the groundwater storage anomalies derived from the GRACE-assimilated CLSM groundwater data-further validated with remotely sensed GRACE TWSA mascon solutions (Figure 2)-and build confidence in using the CLSM outputs for further drought index analyses.The borehole sites were filtered to 172 records based on location, limited to records north of greater Pretoria to remove unnatural groundwater level fluctuations due to irrigation and urban abstraction.The water level data were then linearly interpolated for temporal continuity and standardized for comparison with the basin-averaged CLSM groundwater storage anomalies, as shown in Equation ( 1): where z equals the water level/storage anomaly index, x is the water level/storage for a given month, x is the average water level/storage for the entire record, and σ is the standard deviation of the record for that borehole or the LRB.After standardization, all records were averaged to create a mean and standard deviation to characterize uncertainty.Growing season anomalies were calculated by filtering for observations in October through March by water year.

In Situ Datasets
A groundwater level dataset from the South African Department of Water and Sanitation (SA DWS) containing more than 1000 borehole sites matching the South African extent of the Limpopo River Basin was selected to assess the groundwater storage anomalies derived from the GRACE-assimilated CLSM groundwater data-further validated with remotely sensed GRACE TWSA mascon solutions (Figure 2)-and build confidence in using the CLSM outputs for further drought index analyses.The borehole sites were filtered to 172 records based on location, limited to records north of greater Pretoria to remove unnatural groundwater level fluctuations due to irrigation and urban abstraction.The water level data were then linearly interpolated for temporal continuity and standardized for comparison with the basin-averaged CLSM groundwater storage anomalies, as shown in Equation ( 1): where z equals the water level/storage anomaly index,  is the water level/storage for a given month, ̅ is the average water level/storage for the entire record, and  is the standard deviation of the record for that borehole or the LRB.After standardization, all records were averaged to create a mean and standard deviation to characterize uncertainty.Growing season anomalies were calculated by filtering for observations in October through March by water year.
Daily discharge rates from the Global Runoff Data Centre (GRDC) were selected from 36 gauges across the LRB to compare with the CLSM surface runoff dataset.Each record was multiplied by upstream areas, summed by month, log-transformed for greater interpretability, and then standardized like the groundwater level dataset.Fifteen daily pre-    Daily discharge rates from the Global Runoff Data Centre (GRDC) were selected from 36 gauges across the LRB to compare with the CLSM surface runoff dataset.Each record was multiplied by upstream areas, summed by month, log-transformed for greater interpretability, and then standardized like the groundwater level dataset.Fifteen daily precipitation records, also from the SA DWS, were used to assess the accuracy of GPM IMERG for the basin.A monthly validation with both its given unit (mm) and a logtransformed unit were evaluated for context (Figure 3).Table 2 shows the in situ datasets used for this study.One limitation of these in situ datasets is that they only represent the South African region of the LRB.Thus, they are only used as partial references to compare against the modeled and remotely sensed products, which cover the entire basin.

Time Series of Additional Remotely Sensed and Modeled Datasets
Monthly Land Surface Temperature (LST) estimates were calculated using the averages of daily MODIS Terra (MOD11) and Aqua (MYD11) products, which are generated using the generalized split-window algorithm and brightness temperature values observed in the sensors' thermal infrared bands.These products are thoroughly validated for various land covers and approach an error of ±1 K [60]).NDVI estimates are based on the averages of monthly MOD and MYD products, which take advantage of MODIS' near infrared and red bands to produce a standardized reflectance ratio (−1 to 1), minimizing error [58].
Monthly precipitation accumulation estimates from GPM IMERG are derived from a constellation of passive and active microwave-sensing satellites with a core observatory that houses both the GPM Microwave Imager and Dual Frequency Precipitation Radar [57].The available daily accumulation final run products were summed on a monthly timestep for October 2003 through September 2021, after which the late run products were required.Potential evapotranspiration products are calculated using the Penman-Monteith equation, which incorporates daily meteorological reanalysis data generated by the National Atmospheric and Space Administration's Global Modeling and Assimilation Office and 8-day products from MODIS, including land surface temperature, land cover, and albedo [59].Monthly estimates were aggregated using a weighted average algorithm to appropriately partition the given 8-day estimates that overlap on a monthly timestep.
Soil moisture estimates at both the surface and root zone layer from CLSM were investigated.Due to the lack of in situ soil moisture data, downscaled 1 km volumetric water content based on remotely sensed brightness temperatures captured with the passive L-band radiometer on board the Soil Moisture Active Passion (SMAP) mission were assessed for similar drying and wetting signals [76].The selected soil moisture product has a spatial resolution of 1 km at near daily resolution, with a reported unbiased root mean squared error of 0.063 m 3 m −3 .This soil moisture product is engineered using the thermal inertia principle and requires both MODIS LST and NDVI products [62].Iterations of this product have been generated using a combination of various sensors, including the Advanced Microwave Scanning Radiometer (AMSR) for Earth Observing System [76], AMSR-2 [77], and aircraft and field experiments [78,79].Both the surface soil moisture estimates from CLSM and SMAP estimate the water content within the top 5 cm, which must be multiplied to the volumetric soil moisture to get a liquid water equivalent length.All variables mentioned so far are plotted in Figure 4 to provide context.

Drought Indices (Standardized Anomalies, VHI, SPEI)
The time series datasets for this study period were standardized appropriately (at both monthly and growing season time steps) to compare across variables at their respective frequencies, using the following equation, Equation (2): where z equals the variable anomaly, x i is the estimate for month/season i, x i is the monthly estimate for that given month/season, and σ i is the standard deviation of the monthly/seasonal estimate.Notice that this equation slightly differs from Equation (1).
Monthly LST and NDVI values were used to calculate the Temperature Condition Index (TCI) (Equation ( 3)) and Vegetation Condition Index (VCI) (Equation ( 4)), respectively.The former is an indicator of temperature stress, while the latter is an indicator of moisture content.Lower values in either signal drought-like conditions during the window of observation.Their composite index, VHI (Equation ( 5)), can then be calculated using a weighted average of the TCI and VCI.Typically, an alpha value of 0.5 is used for studies lacking more nuanced information regarding the relative importance of either temperature or moisture conditions in the basin [15,16].
The SPEI is a climatic drought index developed by Vincente-Serrano et al. (2010) and is informative in that it includes the effects of temperature variability based on the aggregated deficit or surplus of water availability at a selected timescale (e.g., 3, 6, 12, 48 months) [18] (Equation ( 6)).
where D is the monthly difference between precipitation, P, and potential evapotranspiration, PET, for a given month i.The differences are then normalized based on a distribution (e.g., Gamma, Pearson III, Log-Logistic) calculated with parameters generated from the L-moment procedure.Over a long-enough record, the SPEI is able successfully capture the intensity and duration of drought events over a region [19].For this study, SPEIs based on the Gamma distribution at timescales of 12 months were calculated using an open-source python library called climate-indices developed at the National Ocean and Atmospheric Administration [80].The 12-month scale was selected for its possible relationship to groundwater anomalies [41].

Statistical Analyses
To identify trends in the time series, a Mann-Kendall test was performed for each calculated index and anomaly.It is a widely used non-parametric test that can determine the presence of a negative or positive monotonic trend based on the S-statistic (Equation ( 7)).The null hypothesis, which assumes the existence of a non-monotonic trend in a series without autocorrelation, must be rejected.A seasonal test function (to remove serial correlation) under the pyMannKendall package was used for the monthly datasets [81].
where n is the number of total observations, i is ranked from 1 to n − 1, and j is ranked from i to n.A Spearman's Rank correlation (or rho) test, another non-parametric test, was performed for a more robust time series analysis.The Spearman's Rank test's null hypothesis assumes the rho coefficient to equal 0, meaning that no trend exists, such as the Mann-Kendall test, and must be rejected to confirm any monotonic trend.Both statistical tests have become well documented since their popularization [82].Correlation coefficients were also calculated between the standardized time series values.

Results
Figure 5a shows that the gridded GRACE-assimilated CLSM groundwater storage anomalies capture the decreasing trend in groundwater levels, matching the annual trend by growing season, of the distribution of in situ groundwater level anomalies to within one standard deviation (highlighted in a lighter shade of blue).The CLSM anomalies exhibit a relatively greater annual magnitude in the growing season, oscillating between −1.5 and +1.9.The in situ anomalies hover above 0 between 2006 and 2016, after which there is a steady decreasing trend to about −0.5 by 2021.However, the CLSM anomalies further show a recovery beyond 2020, trending positive to a peak anomaly of about +0.8 in 2021; the in situ anomalies do not reach a baseline of 0. Interestingly, the trough in 2016 is not captured with the in situ data, and the rebound in groundwater levels is not as dramatic between 2020 and 2022.Below, (b) shows the transformed and standardized seasonal in situ discharge and CLSM surface runoff estimates.Compared to the groundwater anomalies, the in situ discharge observations show more resilience, returning to a baseline of 0 after the 2016 drought.The CLSM surface runoff also exhibits a similar pattern as the monthly discharge totals, except the peak in 2008 and trough in 2010-2012.The naturally low flows of this basin likely cause the minimal variation in discharge or runoff anomalies, in contrast to the visibly more negative groundwater anomalies of 2016 and 2020, based on Figure 5a.
Table 3 shows all Mann-Kendall and Spearman's Rank statistics, which confirm a couple of the trends already depicted in Figure 5.While no obvious trend exists at the annual frequency (based on the growing season), groundwater and both soil moisture anomalies appear to have a significantly decreasing monotonic trend, with both non-parametric tests.Additionally, even in the growing season, the significance of the groundwater storage anomalies straddles a p-value of 0.05.On the other hand, NDVI and evapotranspiration monthly anomalies show strong signs of a monotonically increasing trend, though not as significant based solely on the growing season.Figure 6a shows the trends visually with monthly anomalies in NDVI and components from CLSM.Although the soil moisture component varies across a wider distribution of z-scores and trends toward lower z-scores, the NDVI appears relatively consistent, even increasing in the latter months.The rising and falling signals in tandem suggest that NDVI and soil moisture are well correlated; however, that is not consistently the case.Figure 7 shows the plot of drought indices calculated for the LRB.The SPEI-12 and VHI are somewhat correlated, but the groundwater storage anomalies appear to show more intense wet and dry periods, approaching a maximum z-score of 2.2 in early 2008 and a minimum z-score of −1.7 in late 2019.The VHI has greater variance due to its seasonal cycle, and at times do not agree with either the SPEI-12 or groundwater anomalies.For example, in late 2012, healthier vegetation conditions (0.5) were observed in the growing season months despite lower groundwater (−1.2) and SPEI-12 z-scores (−1).A similar example is shown during the rebound in water storage after the 2020 drought, in which the VHI approaches positive z-scores sooner than either the groundwater anomalies or SPEI-12.Given that the VHI is influenced by NDVI, the plot in Figure 6b reveals a similar decoupling as Figures 6a and 7, in which the NDVI anomalies are relatively higher during some droughts (z = −0.6 for 2020) but lower for the rainier growing seasons in 2006 (z = 0) and 2008 (z = 0.75) compared to the CLSM soil moisture, groundwater, and runoff anomalies.Figure 7 shows the plot of drought indices calculated for the LRB.The SPEI-12 and VHI are somewhat correlated, but the groundwater storage anomalies appear to show more intense wet and dry periods, approaching a maximum z-score of 2.2 in early 2008 and a minimum z-score of −1.7 in late 2019.The VHI has greater variance due to its seasonal cycle, and at times do not agree with either the SPEI-12 or groundwater anomalies.For example, in late 2012, healthier vegetation conditions (0.5) were observed in the growing season months despite lower groundwater (−1.2) and SPEI-12 z-scores (−1).A similar example is shown during the rebound in water storage after the 2020 drought, in which the VHI approaches positive z-scores sooner than either the groundwater anomalies or SPEI-12.Given that the VHI is influenced by NDVI, the plot in Figure 6b reveals a similar decoupling as Figures 6a and 7, in which the NDVI anomalies are relatively higher during some droughts (z = −0.6 for 2020) but lower for the rainier growing seasons in 2006 (z = 0) The spatial distribution of the variables in Figures 8 and 9 reveals some distinct patterns relative to the anomalies' time series.For example, for both 2016 and 2020, despite experiencing drought conditions, the NDVI remained high in the mountainous region of the escarpment in central-east LRB, as shown in Figure 8a,b.Precipitation patterns are also well correlated with vegetation health, with higher percentages in VHI showing spatial association with greater monthly precipitation accumulation.Groundwater and soil moisture estimates appear to be especially low to the east in Mozambique, but again, relatively high in the central-east, likely due to orographic influence.Figure 9a,b depicts the opposite conditions with much higher total monthly rainfalls.However, some discrepancy exists between the vegetation and groundwater and soil moisture distributions in Figure 9b.Despite showing relatively good vegetation conditions in the east, groundwater storage is noticeably lower, compared to Figure 9a.In both anomalously wet and dry growing seasons, the spatial distributions of rainfall, groundwater, and soil moisture are only moderately consistent.Table 4 lists the correlation coefficients across all variables as a matrix, which all exhibit statistically significant relationships with at least 95% confidence.Here, surprising results include the lower correlation between precipitation and NDVI (r = 0.510) and high correlation between surface runoff and root zone soil moisture (r = 0.913) but not NDVI (r = 0.684).Lower correlations were expected between groundwater and precipitation due to the lagging influence of infiltration (r = 0.475), but a moderate relationship exists between it and NDVI anomalies (r = 0.705).While all correlation coefficients were considered significant (p-value < 0.05), precipitation and groundwater anomalies showed the least significance (p-value = 0.04).Table 5 is the correlation matrix of monthly anomalies and indices for all variables.Here, we see that the strongest relationships (besides VHI-NDVI and the soil moisture components, which are dependent by default) exist between the surface soil moisture and VHI (r = 0.774), and the least between precipitation and groundwater (r = 0.064).The relationship between groundwater and NDVI anomalies is moderate (r = 0.671).Weak linear relationships also appear between NDVI and both precipitation and runoff (r < 0.4).Precipitation is mostly correlated with runoff (r = 0.724) and moderately with the surface soil moisture (r = 0.509).All relationships except for P and GW (p-value > 0.3) had a p-value of < 0.001, removing the need for a sub-table.
Figure 10a,b shows the spatial growing season anomalies for the groundwater storage and NDVI, two variables that were determined to have significant monotonic trends based on their monthly anomalies.Although Table 3 had shown non-significant results for the growing season, the anomalies' spatial distributions reveal a different story.The decreasing trend in groundwater storage appears to be strongest in the eastern LRB, within Mozambique's boundaries.Despite the wetter growing seasons of 2021 and 2022, it appears that the groundwater reserves have not quite recovered, likely matching the in situ borehole anomalies from South Africa (Figure 3a).However, a bimodal trend exists for the NDVI, in which different parts of the basin are experiencing greening or browning.Figure 10c,d shows the calculated p-value for each pixel and trend, confirming a significant trend in decreasing groundwater storage, but more spatial heterogeneity with the NDVI.Table 4 lists the correlation coefficients across all variables as a matrix, which all exhibit statistically significant relationships with at least 95% confidence.Here, surprising results include the lower correlation between precipitation and NDVI (r = 0.510) and high correlation between surface runoff and root zone soil moisture (r = 0.913) but not NDVI (r = 0.684).Lower correlations were expected between groundwater and precipitation due

Discussion
The results show unexpected relationships across the analyzed hydrologic variables.Although precipitation, surface runoff, and evapotranspiration did not exhibit any significant trends, NDVI, groundwater storage, and soil moisture conditions did (Table 3).This was confirmed with the SPEI-12, which did not reveal a specific trend, but was informative for both drought duration and magnitude, tracking similarly to the groundwater storage anomalies.On the other hand, NDVI showed a significant increasing trend, while LST showed no such trend during this reference period.Surprisingly, the VHI also showed no significant trends; nevertheless, its peaks and troughs still informed drought conditions such as the SPEI-12, albeit at a higher frequency (Figure 7).Assuming all datasets analyzed in this study are reliable enough to exhibit the appropriate signals for drought conditions-as verified with the in situ datasets for rainfall, groundwater levels, and discharge, and remotely sensed datasets for TWSA and soil moisture (Figures 2, 3, 4a and 5)-these results point to a potential decoupling of hydrological and ecological droughts.This is further strengthened by the comparison of growing season and monthly analyses.As shown in Figure 6, the corresponding offset between NDVI and groundwater anomalies are visible at both timescales.Because Table 3 reveals that the monthly anomalies exhibit strong trends, the lack of obvious trends in the growing season anomalies suggests that interannual interactions between subsurface processes (i.e., plant water uptake, recharge, and human abstraction) may be maintaining some sense of equilibrium through times of high and low annual rainfall.
In the first half of the time series in Figure 6a, the NDVI peak anomalies were relatively lower than the soil moisture components, averaging a z-score of about 1 compared to 1.5, respectively.As the time series approached 2017, this relationship inversed, in which the NDVI peak anomalies were higher than the soil moisture anomalies, suggesting that, on the basin-scale, vegetation in the LRB can withstand these more long-term drought events.However, that may not be the case for more short-term, intense droughts, as shown in January 2016, where the NDVI anomaly was in fact marginally lower than either soil moisture anomaly.This pattern becomes more obvious for the growing season anomalies in Figure 6b, in which the NDVI z-score reaches a global minimum of about −1.3 in 2016 but does not again in the following droughts.However, the groundwater and soil moisture anomalies drop again for the second time during the 2019 drought.Based on the spatial distribution, Figure 8b shows how precipitation may have been the cause for this rebound.However, this is not consistent for Figure 9b in the western LRB, in which a relatively high NDVI is associated with comparably drier surface soil moisture and groundwater storage conditions, with patterns such as the distribution of values depicted in Figure 8 for the drier areas.
A lag of approximately one month exists between precipitation and NDVI, while the runoff peaks align well with the monthly precipitation accumulations (Figure 4).The correlation coefficients in Table 5 confirm these observations and reveal further subsurface flow dynamics.For example, the groundwater anomalies appear to be uncorrelated with precipitation anomalies, which confirms that seasonal precipitation is insufficient for groundwater recharge [56].However, the soil moisture and vegetation components appear to retain some water post rainfall, as demonstrated by their correlations, albeit weak (r ≤ 0.5).
Perhaps the most interesting finding of this study is the confirmation of this decoupled relationship between groundwater and NDVI anomalies.Although the basin-averaged growing season anomalies for groundwater storage and NDVI showed no significant trends, this was proven to be false at the distributed level, as shown in Figure 10.Decreasing groundwater storage anomalies were especially significant in the eastern basin, just north of Limpopo River's mouth at the Mozambiquan coast.However, the NDVI appeared to be increasing during the growing season, despite depleting groundwater storage.The decreasing trend in NDVI in the central LRB was more consistent with what is expected of decreased groundwater storage [83,84].This pattern could reflect mean annual rainfall.In the wetter portions of the study region to the southeast, vegetation water subsidies from groundwater would likely be less important and the ample rainfall maintains the greening of Krueger National Park at the border between South Africa and Mozambique, and Banhine National Park near the eastern boundary of the LRB in Mozambique.In the dryer portion of the basin to the northwest, deeper rooted plants (≥40 m) may be more dependent on groundwater storage [85].
A more threatening explanation may be that the increasing trend in NDVI is attributed to groundwater loss due to plants exploiting limited resources as an adaptation strategy.Key processes regarding the vegetation-groundwater dynamic include (1) partitioning of water by plants during precipitation and evapotranspiration, (2) extraction of groundwater in either the saturated or unsaturated zone, and (3) hydraulic redistribution via the root system [47,86].Past studies have noted that woody invasive species in South Africa, previously introduced in past centuries to address timber shortages, can lead to reduced recharge and exploit groundwater at deeper root depths [50][51][52].Therefore, this explanation may offer an insight into how or why depleting water resources is not correlated with browning in some areas of the LRB.Another point to consider is the difference in land cover, assuming the Sentinel-2-based map is accurate in Figure 1.Given the more densely forested region of the eastern LRB compared to the grassy-rangeland cover of the central LRB, future studies that relate equilibrium between woody vegetation and grasses may be necessary to better characterize the relationship between groundwater (or precipitation) deficit and vegetation species [87].
The middle-ground hypothesis is that plant species in such regions can leverage the more sustainable shallow, unconfined aquifers, which are recharged with every growing season, given normal precipitation conditions, while more intense precipitation events are required to offset groundwater loss [54,55,88].Given the surface lithology of the LRB, as shown in Figure 11, this may be a plausible explanation, because the majority of the eastern LRB is composed of sand dune alluvium (deposited river sediment), and is primarily a phreatic aquifer [85,88].Upstream of the Limpopo River, metaigneous rocks, which have hard crystalline structures, may be associated with fractured, deeper aquifers where root systems depend on groundwater resources, leading to a browning in each growing season that does not receive enough intense precipitation, and thus limiting groundwater recharge [56,85,89].Because the statistics for precipitation anomalies remain inconclusive with negative trends but high p-values, decreasing annual precipitation cannot be accepted as the sole reason for these groundwater anomalies.In fact, human abstraction must also be accounted for, as the semi-arid nature of the LRB requires using groundwater for agricultural irrigation as an alternative to surface water, which is relatively low throughout most of the year [57,85].Because seasonal precipitation in the growing season is only able to sustain the vegetation so much and cannot adequately recharge groundwater storage (Table 4), this browning phenomenon appears to be more significant in the arid climate of the central river valley (Figure 10).This is further supported by the non-significant, decorrelated relationship between monthly precipitation and groundwater anomalies (Table 5).
Therefore, the relationships and statistical metrics for NDVI and CLSM-based anomalies suggest that ecological drought conditions with respect to vegetation health may not inherently follow hydrological drought conditions, at least according to the available record of modeled groundwater storage anomalies in the LRB.One concerning conclusion is that groundwater storage overall has depleted, despite the temporary relief in 2021 (Figures 3a and 10).Because the record-breaking drought of 2016 was particularly intense due to a strong El Niño event, water management practices must become even more precautious in preparation for such future events [46,90,91].Therefore, vegetation health in the eastern LRB cannot be used to assess groundwater availability, though it may be in the more arid climate of the central west region.
ord of modeled groundwater storage anomalies in the LRB.One concerning conclusion is that groundwater storage overall has depleted, despite the temporary relief in 2021 (Figures 3a and 10).Because the record-breaking drought of 2016 was particularly intense due to a strong El Niño event, water management practices must become even more precautious in preparation for such future events [46,90,91].Therefore, vegetation health in the eastern LRB cannot be used to assess groundwater availability, though it may be in the more arid climate of the central west region.

Conclusions
This study analyzed a collection of monthly and annual growing season anomalies related to meteorological, hydrological, and ecological droughts, including land surface temperature, NDVI, precipitation, evapotranspiration, soil moisture, surface runoff, groundwater storage, and two drought indices, the VHI and SPEI-12, for the Limpopo River Basin from a water year of 2004 through 2022.We highlight the importance of vegetation dynamics with respect to groundwater anomalies and recommend that future groundwater studies be investigated in the context of vegetation health to understand the

Conclusions
This study analyzed a collection of monthly and annual growing season anomalies related to meteorological, hydrological, and ecological droughts, including land surface temperature, NDVI, precipitation, evapotranspiration, soil moisture, surface runoff, groundwater storage, and two drought indices, the VHI and SPEI-12, for the Limpopo River Basin from a water year of 2004 through 2022.We highlight the importance of vegetation dynamics with respect to groundwater anomalies and recommend that future groundwater studies be investigated in the context of vegetation health to understand the impact of both humans and plants on terrestrial water storage trends.Significant trends in both vegetation greening and browning appear through the study period's growing seasons, which exhibit signs of both negative and positive correlation with the other significant trend in groundwater depletion based on aquifer type and surficial lithology.Although this research is limited by both the lack of in situ data in Mozambique, Botswana, and Zimbabwe and the unaccounted influence of human impact by either abstractions or dam control, basin-wide in situ data from South Africa's DWS and the GRDC, and the agreement in drought signals across modeled and remotely sensed data products increase the validity of these research conclusions.Therefore, a method to isolate surface discharge is necessary to characterize its role in linking the hydrologic cycle with not only human impact but also vegetation health.Assuming all excess precipitation is converted to either runoff or river flow, reanalyzing the studied factors at shorter timescales (hourly to daily) is required to understand the intermediate responses of the system.Other limitations may include the validity of the land cover map and the surface lithology map provided by the USGS; however, these maps are not essential to the conclusions made in this study.As one of the few reports in the LRB that consider vegetation health in the context of groundwater, this study calls for renewed attention to the heterogeneous response of vegetation to groundwater availability and seasonal precipitation.Thus, droughts must be researched with the greater ecosystem to better assess present-day and future water management practices in a water-stressed, semi-arid basin such as the Limpopo under current climate change projections.

Figure 1 .
Figure 1.A land cover map of the Limpopo based on Sentinel-2 from ESRI, overlaid with in situ sites and the main channel of the Limpopo River.

Figure 1 .
Figure 1.A land cover map of the Limpopo based on Sentinel-2 from ESRI, overlaid with in situ sites and the main channel of the Limpopo River.

Figure 2 .
Figure 2. Time series of water storage anomalies from GRACE and CLSM.Figure 2. Time series of water storage anomalies from GRACE and CLSM.

Figure 2 .
Figure 2. Time series of water storage anomalies from GRACE and CLSM.Figure 2. Time series of water storage anomalies from GRACE and CLSM.

Figure 4 .
Figure 4. Basin-averaged monthly time series of hydrologic components: (a) Monthly precipitation accumulation and average monthly soil moisture in liquid water equivalent (mm) from both CLSM and SMAP are plotted for comparison.(b) Evapotranspiration, NDVI, and land surface temperature estimates from MODIS.(c) Groundwater storage and surface runoff from CLSM.

Figure 4 .
Figure 4. Basin-averaged monthly time series of hydrologic components: (a) Mo accumulation and average monthly soil moisture in liquid water equivalent (mm and SMAP are plotted for comparison.(b) Evapotranspiration, NDVI, and land s estimates from MODIS.(c) Groundwater storage and surface runoff from CLSM

Figure 4 .
Figure 4. Basin-averaged monthly time series of hydrologic components: (a) Monthly precipitation accumulation and average monthly soil moisture in liquid water equivalent (mm) from both CLSM and SMAP are plotted for comparison.(b) Evapotranspiration, NDVI, and land surface temperature estimates from MODIS.(c) Groundwater storage and surface runoff from CLSM.

Figure 5 .
Figure 5. Normalized seasonal anomalies: (a) groundwater storage and level anomalies are plotted above, while (b) log-transformed runoff and discharge anomalies are plotted below.

Figure 5 .
Figure 5. Normalized seasonal anomalies: (a) groundwater storage and level anomalies are plotted above, while (b) log-transformed runoff and discharge anomalies are plotted below.

Figure 6 .
Figure 6.NDVI and precipitation anomalies are plotted with components from CLSM plotted (a) for all months and (b) for the growing season.

Figure 6 .
Figure 6.NDVI and precipitation anomalies are plotted with components from CLSM plotted (a) for all months and (b) for the growing season.

Figure 5 .
Figure 5. Normalized seasonal anomalies: (a) groundwater storage and level anomalies are plotted above, while (b) log-transformed runoff and discharge anomalies are plotted below.

Table 3 .
Time series analysis of anomalies and indices with significant trends bolded.

Table 5 .
Monthly Anomalies and Drought Indices Correlation Matrix.