Prospects for Long-Term Agriculture in Southern Africa: Emergent Dynamics of Savannah Ecosystems from Remote Sensing Observations

: Hydro-climatic resilience is an essential element of food security. The miombo ecosystem in Southern Africa supports varied land uses for a growing population. Albedo, Leaf Area Index (LAI), Fractional Vegetation Cover (FVC), Solar-Induced chlorophyll Fluorescence (SIF), and precipitation remote-sensing data for current climate were jointly analyzed to explore vegetation dynamics and water availability feedbacks. Changes in the surface energy balance tied to vegetation status were examined in the light of an hourly albedo product with improved atmospheric correction derived for this study. Phase-space analysis shows that the albedo’s seasonality tracks the landscape-scale functional stability of miombo and woody savanna with respect to precipitation variations. Miombo exhibits the best adaptive traits to water stress which highlights synergies among root-system water uptake capacity, vegetation architecture, and landscape hydro-geomorphology. This explains why efforts to conserve the spatial structure of the miombo forest in sustainable farming of seasonal wetlands have led to signiﬁcant crop yield increases. Grass savanna’s high vulnerability to water stress is illustrative of potential run-away impacts of miombo deforestation. This study suggests that phase-space analysis of albedo, SIF, and FVC can be used as operational diagnostics of ecosystem health.


Introduction
Global warming due to anthropogenic causes such as greenhouse gas emissions impacts the interconnected climate, hydrological, environmental/ecological, and social systems, that is the Earth system [1,2]. Land-use and land cover dynamics are crucial for modeling climate change impacts, particularly in areas with intensive agricultural practices and population growth. Agricultural lands cover approximately 40% of the land surface of the Earth, and population increases worldwide present challenges for sustainable food production and water resource management [3,4]. One of the main issues hindering African development is drought, making up fewer than 20% of natural disasters in Sub-Saharan Africa but more than 80% of the population affected by natural disasters in this region [5,6]. For example, drought can intensify in dry areas resulting from increased evaporation and enhanced atmospheric moisture holding capacity caused by warming [1]. Projected increases in water scarcity in combination with restricted access to groundwater can in turn lead to increased risk of heat waves and wildfires with devastating impacts on African communities dependent on rainwater as a primary water source [1,2,[7][8][9].
In addition to estimating albedo to improve radiative uncertainty and examine surface energy changes, vegetation parameters were retrieved to map plant dynamics in response to different hydrometeorological conditions from year to year. Data such as the Leaf Area Index (LAI) and Fractional Vegetation Cover (FVC) facilitate the tracking and evaluation of changes in the biophysical properties and photosynthetic rates of canopy and other plant cover [50,51]. LAI and FVC are based on remotely sensed spectral signatures, and Solar-Induced chlorophyll Fluorescence (SIF) provides a more direct approach for predicting changes in photosynthetic efficiency that are not easily detected by changes in reflectance [51]. Depending on canopy condition that is measured by LAI, increases in vegetation cover increase shade effectively leading to lower temperatures near the ground. Phase-space analyses were conducted with albedo, SIF, LAI, and FVC to understand the vegetation life cycle and diurnal land-atmosphere dynamics of the miombo ecosystem. Joint, spatial analyses of these multi-sensor data with precipitation at different nominal spatial resolutions enabled intercomparisons to examine the consistency of the representation of spatial heterogeneity across the different sensors.
The main objectives of this study were: (1) to illustrate the utility of remote sensing data for capturing feedbacks between precipitation variations, vegetation dynamics, and water availability; (2) to provide a quantitative understanding of natural ecosystem resilience in the region to water stress and changes in climate and its potential agricultural applications through remotely-sensed land surface monitoring. Section 2 covers the materials and methods, namely the study area's topography and climatology and data processing methodology for precipitation data, vegetation parameters, land cover mapping, albedo, SIF analysis, and joint data analyses. Section 3 covers the main results and discussion regarding the spatiotemporal variability of rainfall, spatial variability of vegetation data, SIF analysis, spatiotemporal variability of albedo, phase-space analysis, seasonal and interannual variability of albedo and vegetation, data anomalies, uncertainty, and future directions. Section 4 contains the main conclusions regarding the interactions among rainfall variability, soil moisture availability, and root water uptake demonstrated by analysis of remotely sensed land-surface properties and the functional behaviors of different land cover types as illustrated by phase-space analysis. The interplay between the hydrologic behavior of dambos (e.g., seasonal wetlands in the headwaters) and the miombo trees is further discussed in the light of potential implications for long-term, sustainable agriculture.

Study Area
The Upper Zambezi River Basin's (UZRB) ecosystem consists of miombo, tree savanna, grasslands, and wetlands in an arid to semi-arid transition zone between the Kalahari Desert to the south and the tropical rainforest of the Congo to the north [52]. Figure 1 shows a digital elevation map of the UZRB. In arid ecosystems, productivity is driven primarily by soil water content, and the seasonal cycle of productivity is driven by soil moisture in larger part than by vapor pressure deficit [53]. The climate of this area is characterized by complex rainfall gradients driven by the movement of the Inter-Tropical Convergence Zone (ITCZ), a narrow, clearly defined band of clouds that stretches around the globe parallel to the equator [54,55]. The ITCZ moves northward during the cooler, dry season (May to October), generating an average accumulated precipitation of less than 100 mm. During the warmer, rainy season (November to April), the ITCZ moves southward, resulting in an average accumulated precipitation between 700-1000 mm [52] (see Table S1 in the Supplementary Materials for average accumulated seasonal rainfall values).  Table 1 summarizes the data products used for this project. The UZRB spans two tiles in sinusoidal projection used in the canonical MODIS products. All MODIS data used in this study were stitched and re-projected to UTM 34S projection using the HDF-EOS to GeoTIFF (HEG) Conversion Tool [56]. Since the ITCZ over this region often results in cloudcontaminated pixels, a software package called TIMESAT performed adaptive temporal filtering using the Savitzky-Golay method to smooth out anomalies and discontinuities in the data impacted by cloud contamination [57]. This is accomplished by using time series data, quality control information, and moving windows to fit a quadratic polynomial function. All data products were re-scaled and re-projected to 1 km, UTM 34S projection with various methods depending on the underlying physics and original spatial resolution of each product described in the following sections.   1 The MODerate resolution Imaging Spectroradiometer (MODIS) was launched on two satellites: Terra, which descends across the equator in the morning local time [71], and Aqua, which ascends across the equator in the afternoon local time [72]. MCD refers to observations collected by the combined MODIS instruments from both Terra and Aqua satellites. MOD refers to data retrieved by the Terra satellite only, and MYD refers to data retrieved by the Aqua satellite only. 2

Precipitation Data
Integrated Multi-satellitE Retrievals for Global Precipitation Measurement (IMERG) precipitation intensity data at 0.1 • (∼10 km), monthly resolution were obtained from the National Aeronautics and Space Administration (NASA) Precipitation Processing System's (PPS) Science Team On-Line Request Module (STORM) [70]. The original data were bilinearly interpolated to 1 km spatial resolution. Accumulated precipitation data were calculated by multiplying the monthly mean precipitation intensities by the time period of interest (i.e., converted from mm/h to mm/month, mm/season, mm/period, or mm/year).

Vegetation Data
MODIS MCD15A3H Version 6 Leaf Area Index (LAI) products at 4-day, 500 m resolution were re-projected with the HEG tool, cloud-corrected in TIMESAT (Figure 2), and aggregated to 1 km spatial resolution by taking the average of every four pixels with a final 4-day, 1 km resolution [65]. LAI is representative of canopy cover, which is linked to surface energy balance and modulates interactions between vegetation, hydrologic processes, and the atmosphere, such as evapotranspiration, rainfall interception, throughfall, and stemflow [73,74].
Global 'Orbiting Carbon Observatory-2 (OCO-2)' solar-induced chlorophyll fluorescence (GOSIF) data were obtained from the Global Ecology Group's Data Repository at 8-day, 0.05 • resolution and OCO-2 SIF data were obtained from the Goddard Earth Sciences Data and Information Services Center (GES DISC) at 16-day, 2.25 km by 1.29 km resolution [68,69]. For the purposes of this manuscript, SIF and GOSIF are used interchangeably unless otherwise specified (i.e., OCO-2 SIF used in SIF analysis is a separate data product from GOSIF). This fluorescence is emitted directly from the core of photosynthetic machinery and is thus indicative of vegetative photosynthetic status, unlike other vegetation indices that describe "greenness" [68]. SIF can serve as a proxy for ecosystem gross primary productivity (GPP), with strong correlations between SIF and GPP climatologies on the monthly scale for deciduous broadleaf, mixed forest, evergreen needleleaf forests, and croplands [53]. Significant correlations with higher variability have been observed for land cover types including savanna, evergreen broadleaf, and shrublands [53]. SIF data were re-projected in the Quantum Geographic Information System (QGIS) to UTM 34S projection, and then interpolated using nearest-neighbor to 1 km spatial resolution. OCO-2 SIF are kept at the original resolution, and comparisons between the two data are made by selecting GOSIF values that are closest geographically to each OCO-2 orbit overpass. Fractional Vegetation Cover (FVC) was calculated using a semi-empirical method (Equation (1)) based on the Normalized Difference Vegetation Index (NDVI) (Equation (2)) [75]. NDVI is a measure of "greenness" that exploits chlorophyll reflectance of near-infrared light and absorption of visible red light [76]. Any FVC values over 1 were set to 1, and any pixels whose LAI values were 3 or higher were also assumed to have 100% FVC.
where NDV I 0 is the NDVI corresponding to bare soil (NDV I 0 = 0.3) and NDV I s is the NDVI corresponding to full vegetation cover (NDV I s = 0.8), which were approximated after systematic inspection of MODIS NDVI time series graphs. NDV I is the time-varying NDVI value calculated from TIMESAT-corrected MODIS MCD43A4 Version 6 Band 1 (red) and Band 2 (N IR 1 ) reflectance data that was aggregated by taking the average of every four 500 m pixels to 1 km resolution with daily time steps [76], given by:

Land Cover Re-Classification
The MODIS MCD12Q1 Version 6 product's International Geosphere-Biosphere Programme (IGBP) land cover classification system was used to characterize land cover in the UZRB [58]. The data were downloaded at annual, 500 m spatial resolution and ag-gregated by taking the mode of every four pixels to 1 km resolution, then re-classified into six main land cover types: miombo, tree savanna, grass savanna, wetland, water, and barren/built-up. The re-classification is presented in Table 2. Larger trees and shrublands in the IGBP classification were assumed to be representative of miombo because the understory in miombo woodlands is scattered and composed of suppressed tree saplings and broadleaved shrubs [18,77]. Sedges, forbs, and grasses sparsely but continuously cover the ground layer [18,77]. The dominant tree species in miombo are of the genera Brachystegia, Julbernardia, and/or Isoberlinia, which are unique to miombo ecosystems [18,[77][78][79]. Other important species include Burkea africana, Pseudolachnostylis maprouneifolia, Diplorhynchus condilocarpon, and more [18,77]. The large broadleaf trees are characteristic of wet miombo, where canopy cover is dense and trees are over 15 m in height [18,77]. In the IGBP classification, open and closed shrublands contain woody perennials 1-2 m in height, which is characteristic of Eriosema, Sphenostylis, Kotschya, Dolichos, and Indigofera perennial plants that dominate the shrub composition of dry miombo [18]. Since the IGBP classification provides woody savanna, savanna, grassland, and cropland land cover classes that are distinct from shrubland land cover, the IGBP classes were assumed to cover all woody savanna and grass savanna land cover types. The same re-classification was used by Lowman et al. (2018) [52] to represent Southern Africa using MCD12Q1 Version 5 data instead of MCD12Q1 Version 6 data.
Using the same methodology as Lowman et al. (2018) [52], dynamic wetland maps were created using five wetland predictor variables derived from MODIS MCD43A4 Version 6 Nadir BRDF-adjusted reflectance data after aggregation to 1 km resolution with daily temporal resolution [59,60]. These indices aid in visualizing seasonal changes in wetland area, and stationary land cover data for three consecutive wet years (2010-2012) were used as training data for wetland probability mapping (Appendix A) considering six wetland indices (Figure 3c-h). Two distinct peaks in the histograms generally signifies a "good" wetland predictor. As such, NDVI, which has the greatest amount of overlap between wetland and non-wetland values and the most irregular peaks, was not used as a wetland predictor variable. As shown in Figure 3b, threshold area curves were generated for 2009 to 2018, indicating that, at scales above 150 km, there is less uncertainty with wetland mapping. At finer resolutions, uncertainty increases greatly and varies according to interannual variability in climate and weather. For example, a 60% threshold yields uncertainty on the order of 100%. Based on the black line in Figure 3b marking the permanent wetland area derived from stationary land cover between 2010 and 2012 ( Figure 3a), a threshold value of 99% was selected. Compared to previous work on wetland probability mapping [52] that used MCD12Q1 Version 5 data processed with the nowretired MODIS Reprojection Tool (MRT), the MCD12Q1 Version 6 data processed with the HDF-EOS to GeoTIFF Conversion Tool refines the permanent wetland class in the stationary land cover and thus highlights the interannual variability in wetland extent tied to precipitation. These dynamic wetland area maps were overlaid onto the land cover maps to complete the reclassification, as shown in Figure 4. Figure 5 displays varying vegetation densities of the three land cover types of interest in this study: grass savanna, miombo, and woody savanna. Miombo trees are generally taller and more densely packed, with extensive root systems that can reach up to 27 m in lateral length and over 5 m in vertical depth [17]. On the other hand, woody savanna trees are shorter and more spread out, with smaller root systems ranging 2-12 m laterally and 0.7-1.5 m vertically [80]. Grass savanna vegetation are the least dense with shallow root systems, often containing seasonally waterlogged wetland areas and shallow depressions called dambos that intersperse the landscape [17]. Figure 6 illustrates relative sizes and root zone distributions of woody savanna, grass savanna, and miombo.

Albedo Algorithm
Given the lack of hourly albedo data available for this region, surface albedo data were derived from MODIS MAIAC MCD19A2 aerosol optical depth (AOD) and MCD19A3 bi-directional reflectance distribution function (BRDF) model parameter data using existing methodology [61,62,82]. Figure 7 depicts the surface albedo data processing and calculations flowchart with final spatial and temporal resolutions specified for each data product.
Black-Sky Albedo (BSA), or directional-hemispherical reflectance, is the albedo of a surface illuminated in one direction. Equation (3) describes how BSA is calculated by a simplified polynomial expression to represent the integration of the BRDF over the exitance hemisphere for a single irradiance direction [83]. White-Sky Albedo (WSA), or bi-hemispherical reflectance, is the surface albedo generated by diffuse illumination. Equation (4) describes how WSA is calculated by a simplified polynomial expression to represent the integration of the BRDF over all viewing and irradiance directions [83]. Famine Early Warning Systems Network (FEWS NET) Land Data Assimilation System (FLDAS) surface pressure, Global 30 Arc-Second Elevation (GTOPO30) digital elevation, and MODIS MOD11A1/MYD11A1 land surface temperature data were inputted into the National Renewable Energy Laboratory's Solar Position Algorithm (NREL SPA) to calculate the hourly solar zenith angles (SZA) for each pixel [63,64,66,67,84].
To obtain fraction of diffuse skylight (SKYL) values, a look-up table retrieved from the MODIS User Tools section of Professor Crystal Schaaf's Lab site was used to bilinearly interpolate SKYL values based on AOD and SZA values [86]. These SKYL values act as weights for calculating blue-sky albedo, or albedo corresponding to actual atmospheric conditions, from an interpolation between BSA and WSA, as shown in Equation (5) [83]: where α(θ, λ) is narrowband blue-sky albedo, S(θ, τ(λ)) is the fraction of diffuse skylight, and τ(λ) is the aerosol optical depth. The seven resulting narrowband blue-sky albedos were then combined to calculate shortwave broadband albedo using Equation (6) [83]: where A(θ) is the broadband albedo and c i is the conversion coefficient corresponding to MODIS reflectance band i presented in Table 4. Table 4. Narrowband spectral-to-broadband albedo conversion coefficients from [48,85]. Since GOSIF data were derived from MODIS data, Modern-Era Retrospective Analysis for Research and Applications, Version 2 (MERRA-2) meteorological reanalysis data, and OCO2-SIF data, both GOSIF and OCO-2 SIF datasets were compared in an orbit-based analysis. GOSIF data were corrected with a 24 h daily correction factor, while OCO-2 SIF data were corrected with the following daytime-only daily correction factor (any SZAs over 90 • were set to 90 • ) in Equation (7) [87]:

MODIS Band
whereF s is the average daily SIF value, F s is the instantaneous OCO-2 SIF value at 757 nm, t 1 is the beginning time (sunrise for daytime-only), t 2 is the end time (sunset for daytimeonly), SZA(t 0 ) is the instantaneous solar zenith angle at the time the observation was taken, t 0 , and SZA(t) is the solar zenith angle at time t. To match the 24 h daily correction factor calculations provided in the OCO-2 User Guide as closely as possible, the solar zenith angles were calculated at 10-minute intervals using the NREL SPA, and the integral was approximated via the trapezoidal method for daytime hours only [84,87]. Only OCO-2 SIF observations collected under nadir and for clear skies conditions were used, since the Cubist regression tree model used to produce GOSIF from OCO-2 soundings averaged OCO-2 SIF data collected under similar conditions [68].

Data Analysis
For each year of interest (e.g., 2010-2011 and 2015-2016), the data were compared in phase-space plots by masking the data for each land cover type and calculating the mean values and 95% confidence intervals over time. Ellipses were then overlaid onto the graphs to approximate the shape of the limit-cycles with arrows to indicate the direction of time. The time period for each graph was not based off of the calendar year so that the data points would represent two consecutive rainy and dry seasons.
To explain anomalous seasonality in the albedo data, the slope of the terrain was examined using GTOPO30 digital elevation data in QGIS to determine where flat land or water may impact the data values. When other anomalous behavior was identified in specific locations in the landscape, phase-space graphs for albedo and fractional vegetation cover were plotted as well.

Spatial and Temporal Variability of Precipitation
Spatial and bar graph comparisons of mean monthly accumulated precipitation are displayed in Figures 8 and 9, respectively.  Table S1 in the Supplementary Materials). Based on the spatial interannual comparisons (Figure 8) and land cover comparisons (Figure 9), the northern areas of the basin where miombo and woody savanna dominate receive more accumulated rainfall on average than the southern portion where grass savanna dominates.

Spatial Variability of Vegetation
Seasonal spatial comparisons of LAI, SIF, and FVC are displayed in Figure 10. The use of an NDVI-based fractional vegetation cover calculation appears to highlight the flush of grasses during the rainy season and their subsequent senescence during the dry season characteristic of semi-arid regions in Africa [88]. Woody savanna and miombo retain some fractional vegetation cover during the dry season, which could be attributed to frequent accidental and deliberate fires as people prepare land for harvest or for livestock grazing. The miombo ecosystem likely contains the largest burned area globally estimated at approximately 1 million km 2 yr −1 [17]. These fires generally decimate grasses, shrubs, and fine woody biomass, but trees and other woody plants often persist, re-sprouting quickly from any remaining root systems and stems [17]. LAI and SIF data reflect similar seasonal contrasts, although the tree savanna LAI data in the northwest portion of the UZRB appear to be less sensitive across both seasons. This could be attributed to the smoothing effect of the TIMESAT cloud correction since many pixels were cloud-contaminated during the peak of the rainy season demonstrated by the large fluctuations in original LAI values, particularly for the woody savanna pixel shown in Figure 2. In addition, the seasonally inundated grass savanna/wetland areas in the center of the basin consistently have low values of LAI, GOSIF, and FVC that are distinct from the seasonal greening and senescence observed in other grass savanna areas in the southern portion of the basin. The wetland areas mapped in Figure 4 are based on a 99% threshold from the wetland probability maps determined using stationary MODIS MCD12Q1 permanent wetland pixels (Figure 3b), but the actual wetland areal extent may cover a larger, sparsely vegetated area marked by the low (red) values in Figure 10 that the MODIS (LAI and FVC) and OCO-2/MERRA-2 instruments (GOSIF) are able to detect.

Orbit-Based SIF Analysis
The results of the orbit-based SIF analysis in Figure 11a,c display GOSIF values that are more densely concentrated between 0 and 0.4 W m −2 µm −1 sr −1 , while the OCO-2 SIF values are more variable over a wider range of SIF values. Despite the many advantages of GOSIF in its global, continuous coverage and longer record (2000-2020), these results suggest that the application of the daily correction factor over 24 h may be underestimating the SIF signal because nighttime values were taken into account. On the other hand, the GOSIF data appear to capture seasonal trends more clearly (Figure 11a (Figure 14c) for reference. The locations of the three pixels within the basin used for the interannual monthly mean diurnal cycle comparison (Figure 14a) for each land cover type are marked in Figure 14c with a black plus sign for miombo, a black cross for tree savanna, and a black circle for grass savanna. Figure 14d   Based on these results, the hourly shortwave broadband albedo exhibits a diurnal cycle expected of actual albedo because of the data's dependence on instantaneous solar zenith angles. The spatial graphs ( Figure 13) suggest that albedo values are generally under 0.12 for inundated areas in the wetlands during the rainy season around noon local time (see Table S2 for mean annual albedo values corresponding to permanent wetlands). Surface albedo was also consistent with land cover differences during peak sunlight hours (i.e., noon local time), which makes sense given that the solar zenith angle and subsequently the fraction of diffuse skylight at those times is smaller. Assuming albedo estimates at noon local time are representative of actual albedo, increases in solar zenith angle (i.e., after sunrise and before sunset in Figures 12 and 13) appear to overestimate albedo because the fraction of diffuse skylight is higher. In Figure 12, spectral albedo for near-infrared wavelengths (MODIS Band 2) is also generally lower during the dry season. The decreased albedo leads to increases in the net radiation available for latent and sensible heat fluxes, particularly for sensible heat fluxes because the vegetation is dry. These results suggest a positive feedback between dynamic radiative properties of the landscape and rainfall that captures seasonal variability.  A grassy portion of the temporary wetlands in the UZRB outlined by the gray rectangle in Figure 15 had higher albedo values than would be expected during the dry season. The results from the slope analysis suggest that the anomalous albedo values may be explained by increased turbidity in flat water [89]. The slopes also provide a potential explanation for why the monthly mean diurnal cycles for the grass savanna pixel in Figure 14 near the start of the rainy season (October-December) appear to flip upside down. The grass savanna pixel (marked by a black circle) is in a flat area close to the natural outlet in the southeastern portion of the UZRB, and it turns into a wetland area during the rainy season. This inundated, grassy, flat terrain could slow down waters carrying sediments as they drain out of the basin, increasing the turbidity of the water and causing surface albedo to increase.

Phase-Space Analysis
Phase-space maps of albedo, FVC, and SIF for miombo, grass savanna, and woody savanna across different years are shown in Figures 16 and 17. Each point represents the mean value for all pixels of that land cover type, with colors representing the corresponding day of the year. The 95% confidence intervals were plotted as horizontal and vertical error bars, although many of the confidence interval values were so close to the mean that the short error bars were obscured by the markers. Dark blue shapes show the approximate form of each set of pixels, with arrows overlaid on the shapes to mark the direction of time. The offset between the two branches of the seasonal limit-cycle ("loop") indicates the degree of hysteresis in the relationship between the two variables. The land cover classes used here incorporate dynamic wetland mapping since each point represents a 4-day (for graphs with LAI or FVC but without GOSIF) or 8-day (for graphs with GOSIF) time step, depending on which dataset has the coarser temporal resolution. This enables a more realistic comparison of the interactions between albedo, LAI, FVC, and SIF and rainfall variability for each land cover type as wetland areas become inundated.

Seasonal Variability of Albedo
Across all land cover types, inspection of the smooth limit-cycles in Figure 16a indicates large hysteresis (offset of the trajectory in phase-space) consistent with strong seasonal patterns in the wet year and decreasing albedo in the rainy season (tied to increasing FVC and increased LAI and SIF in Figure 17a) with increasing albedo ranges from woody to grass savanna and miombo in-between. In the dry year, however, the limit-cycles collapse (Figure 16b) with flipped behavior for grass and woody savanna land cover types. For grass savanna, albedo decreases more dramatically in the wet season when the precipitation deficit is at a maximum during the dry year (see Table S2 in the Supplementary Materials) and darker, greener grass flushes the landscape. For woody savanna, the limit-cycles exhibit smaller hysteresis in dry conditions (Figure 16b). Thus, drought depresses the seasonality of albedo. Nevertheless, note that, while the albedo offset is reduced, the miombo limit-cycles are stable and exhibit seasonal hysteresis.
The collapse of the radiative hysteresis (albedo offset) under dry conditions during 2015-2016, particularly for grass savanna pixels, is indicative of the high vulnerability of this ecosystem to water stress. Indeed, Figure 18 compares the mean limit-cycle for grassland under wet and dry conditions with the limit-cycles for grassland in the areas that receive daily mid-day rainfall as per Lowman et al. (2018) [52]. Under dry conditions, the local afternoon rainfall is critical for grass sustainability. The decreased hysteresis offsets during 2015-2016 for woody savanna and miombo reflects the prominent decrease in precipitation in the northeastern portion of the UZRB during 2015-2016, where miombo and woody savanna dominate. As a result, reduced rainfall and soil moisture differences between seasons decrease the difference in albedo during transitional periods. The drop in albedo values for grass savanna in the southern portion of the basin could be due to the relative increase in the contribution of flooded dambos to the landscapescale albedo when grass savanna is stressed under dry conditions in the rainy season of 2015-2016. These results suggest that these data capture both seasonal and interannual variability in response to the amount of precipitation in a given year.

Seasonal Variability of Vegetation
LAI, SIF, and FVC values were highest during the peak of the rainy season and lowest towards the end of the dry season. Both FVC and SIF values dropped close to 0 for grass savanna during the dry season, which can be explained by seasonal senescence characteristic of herbaceous vegetation in semi-arid regions [88]. Based on comparisons of the vegetation data depicted in Figure 17, the degree of hysteresis for each graph varied depending on the time of year and the land cover type. In general, LAI and SIF comparisons display a positive linear dependence without hysteresis for both miombo and woody savanna in the wet year. This behavior highlights increases in LAI due to increases in soil moisture that result in increased photosynthesis and SIF activity. Note however that, due to cloud cover in the rainy season, hysteresis emerges for grass savanna with the lower branch (lower LAI, lower SIF) in the wet season, highlighting seasonal increases in soil moisture in response to precipitation.   [52]. The 95% confidence interval values were plotted as horizontal and vertical error bars, with bar lengths representing the distance of the interval endpoints to the mean. (b) On the right is a conceptual diagram that synthesizes the impact of afternoon convection (Anomaly) produced by local circulations organized by land-atmosphere interactions and groundwater convergence [52] on the seasonal cycle of albedo. This shows the importance of landscape organization of water fluxes in producing hot-spots of enhanced ecosystem resilience.
FVC saturates in the peak of the rainy season (cross-over loop). LAI and FVC comparisons indicate strong hysteresis for most of the year. The phase-space trajectories of (LAI, FVC) and (LAI, SIF) do not change significantly between the wet and dry year suggesting that they capture functional behavior of the three vegetation types, which could be used to parameterize regional phenology in models.
Woody savanna and miombo land cover exhibit similar LAI and SIF values, which can be explained by their similar root system structures. Notably, miombo pixels display enhanced hysteresis offsets compared to woody savanna pixels, emphasizing its more extensive lateral and deep root system [17]. Lateral roots near the surface effectively function as shallow grass roots because water does not have to infiltrate as deep into the soil to reach tap roots. Grass savanna has the strongest hysteresis during the rainy season, which is expected because the shallow root systems are more sensitive to changes in soil moisture. Near the surface, soil moisture is governed by rainfall variability because the topsoil is more susceptible to water loss through evaporation compared to deeper soils. In addition, the development of tap roots in woody savanna and miombo vegetation allow for root uptake of deep-water storage, resulting in a decreased sensitivity to changes in precipitation. Overall, hysteresis decreases between the wet and dry year, with greater impacts on grass savanna because of the higher soil moisture sensitivity due to the shallow root system characteristic of herbaceous vegetation [88].

Sources of Uncertainty
Although all data were interpolated or aggregated to 1 km spatial resolution for data analysis, Figure 19 illustrates how comparisons between coarser resolution (i.e., SIF in Figure 19a) and finer resolution (i.e., land cover in Figure 19b) data may decrease the accuracy of remotely-sensed land surface properties. Detailed, dynamic physical features such as wetland areal extent are lost as spatial resolution decreases, and smaller features such as dambos can be obscured even at 1 km resolution. In addition, the shapes in Figure 17 resemble a cone and an ellipse for LAI vs. SIF and LAI vs. FVC, respectively, which appear as though they are folded over. The sharp drops in LAI during the peak of the rainy season indicate that the interpolation of widely variable LAI values in the TIMESAT cloud correction algorithm for each year could be skewed by outliers at the transition between calendar years (see the first data points of Figure 2b for woody savanna), resulting in LAI that is not representative of the landscape during the rainy season when cloud contamination issues are most severe. Figure 2a shows a decrease in LAI values to the east of the wetlands which may be a result of TIMESAT filtering at short time scales and cloud contamination. Increased values to the west of and within the wetlands after TIMESAT correction are a product of averaging in TIMESAT, suggesting a loss of information. Of the re-classified miombo pixels, mixed forest pixels (based on the original IGBP classification) are more influenced by averaging issues, while evergreen broadleaf pixel values (based on the original IGBP classification) are preserved (see Figure 3a of the stationary land cover for the original IGBP land cover classes).

Future Work
Future studies could further examine the changes in surface energy fluxes discussed in this paper by, for example, analyzing in situ data from flux-towers co-located with the three main ecosystems across the UZRB to examine the phase-space of radiative properties and vegetation phenology. Lowman et al. (2018) [52] used satellite data to specify vegetation parameters including radiative properties and compared their simulations using a hydrologic model with dynamic vegetation [90] against flux-tower data from Mongu, Zambia with good results (see their Figure 7). Lowman and Barros (2019) [91] added a predictive phenology module to the same model to simulate the feedbacks among vegetation dynamics, climate variability (temperature, humidity, and rainfall), and the surface energy budget to be compared against satellite and flux-tower observations. Furthermore, it would be interesting to examine the strength of local land-atmosphere coupling and to characterize the impact of changes in the spatial organization of albedo on the spatial organization of clouds and rainfall. Because of the large areal extent of savanna ecosystems and because seasonal changes in albedo imply seasonal changes in Bowen ratio, and there-fore in low-level relative humidity and moist instability, there is also potential for regional impacts on convective activity tied to vegetation condition (e.g., [92,93]). Whereas previous studies have generally focused on humid climatic regions, the UZRB is located in a region of transition between the Congo River Basin and the Kalahari Desert.
Future work could also involve analyzing finer resolution (i.e., 30 m) Landsat imagery in addition to the land surface data examined in this manuscript to develop more dynamic land cover maps and improve comparisons across different land cover types. Agricultural areas are of great interest since local populations in the UZRB rely on subsistence agriculture. The land use and land cover changes associated with small-scale agriculture can have profound impacts on the water and surface energy balance in the region as population increases and climate change add additional stressors to food security [10,11,17,94,95]. A better understanding of the land-atmosphere feedbacks and natural system resiliency could inform climate-smart agricultural practices aligned with the United Nations Sustainable Development Goals to address climate mitigation, adaptation, and poverty reduction.

Conclusions
A joint, multi-sensor analysis of albedo, leaf area index, fractional vegetation cover, solar-induced chlorophyll fluorescence, and precipitation was conducted to understand how land surface properties are interrelated in maintaining the surface energy balance and water budget. The estimated albedo exhibit realistic diurnal cycles, and the results from the analysis indicate that strong seasonality and inter-annual variability in albedo are directly tied to landscape greening, and in particular soil moisture availability in the root zone and root water uptake following precipitation events. This study establishes regionalscale relationships that illustrate the functional behaviors of miombo, grass savanna, and woody savanna land cover types. Deep-rooted vegetation is less vulnerable to drought. Miombo exhibits optimal adaptive characteristics with robust seasonality for both dry and wet conditions in the region of study (see Figure 6). Operational phase-space analysis of albedo-FVC and albedo-SIF limit cycles can be used to monitor ecosystem stability and to detect changes. For example, grass savanna's high vulnerability to water stress is indicative of potential run-away feedbacks of miombo deforestation.
Dambos (seasonal headwater wetlands) provide localized surface storage for runoff in this region. The hydrologic functionality of dambos is similar to flood irrigation with distributed runoff control, that is the natural system counterpart to micro-dams in smallholder agriculture [96]. Infiltration of ponded water recharges the water table and soil moisture in the rainy season; in the dry season, miombo's deep roots physiologically pump water from deep soil layers upward that is then laterally redistributed at shallow soil depths. During the rainy season, miombo root systems anchor the soil catenas reducing erosion, and thus preserve the topsoil and landform. Dambos and miombo trees play therefore complementary roles in regional hydrology that increase the overall landscape resilience to water stress throughout the year.
The importance of dambos in food production in Southern Africa, and in particular Zambia, Zimbabwe, and Malawi, has long been studied (e.g., [22,[97][98][99][100][101][102]). The synergy between miombo and dambos provides a framework for sustainable agriculture to emulate the resiliency of natural systems in this region. Indeed, recent data show significant increases in crop yield (30-60%) from cultivated dambos when miombo trees are preserved in the landscape as opposed to those that have been cleared [103].
In addition to increasing surface water storage through cultivated dambos, accessing blue water resources, or water withdrawn from surface water and aquifers without depleting environmental flows or groundwater, in a sustainable fashion is crucial for long-term agriculture. Economic water scarcity corresponds to conditions in which blue water resources are physically available, but society's use of that water is hindered by the lack of institutional and economic capacity required to access such water resources [9,[104][105][106]. Agricultural economically water scarce croplands are rainfed agricultural areas where sustainable irrigation expansion is suitable [9]. In Zimbabwe, irrigation by recharge from thousands of human-made earth micro-dams has long been recognized as key to high agricultural productivity [107][108][109]. Dambos play in natural systems an irrigation role not unlike the small earth dams in Zimbabwe farms. Abandonment and inadequate maintenance of these micro-dams in the last twenty years has been in large part associated with the severe decline in agricultural productivity (e.g., [107,108]). Therefore, prospects for long-term subsistence agriculture in the region encompass the preservation of dambos' ecohydrological functionality.
As shown by Rosa et al. (2020) [9], croplands in the savannah ecosystems to the north and northeast of the Upper Zambezi River Basin experience agricultural economic water scarcity, particularly in the dry season. Irrigation expansion for Sub-Saharan croplands facing economically water scarce conditions along with sustainable deficit irrigation practices (i.e., irrigation that reduces blue water supply to below maximum levels and enables the growth of mildly water-stressed crops with minimized impact on crop yield) from areas not currently facing water scarcity such as the UZRB can generate enough food for 189-235 million more people with a roughly 24-96% increase in irrigation water consumption [9]. Adopting sustainable deficit irrigation practices in conjunction with cultivated shallow dambos and miombo trees can enhance the resilience of the ecosystem and surrounding communities, improving food and water security as climate hazards intensify in the coming years [1,9,103,110].

Acknowledgments:
The authors would like to thank the reviewers for their helpful comments and constructive criticisms of this manuscript. Lauren E.L. Lowman, as well as other past and present members of the Barros Research Group (Malarvizhi Arulraj, Masih Eghdami, Steven Chavez, Yueqian Cao, Mochi Liao, and others), provided invaluable guidance to T.M.W. in the early stages of this work.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

UZRB
Upper Zambezi River Basin WSA White-Sky Albedo

Appendix A. Wetland Probability Mapping Algorithm
Using seven bands of TIMESAT-corrected MODIS MCD43A4 Version 6 Nadir BRDFadjusted reflectance data (Table A1), several indices sensitive to the greenness and water content of the surface were tested to determine how well they predicted the dynamic wetland areal extent using a similar methodology as Lowman et al. (2018) [52]. Stationary MODIS MCD12Q1 Version 6 land cover data from consecutive years with the largest permanent wetland area (2010-2012) were used to create training data of each index, and, based on how well the peaks of wetland and non-wetland pixels were separated in the histograms in Figure 3, five wetland predictor variables were selected: Near-Infrared Reflectance (NIR 1 , or MODIS Band 2), the Normalized Difference Water Index (NDWI), Tasseled Cap (TC) Brightness, TC Greenness, and TC Wetness. These indices, described in Equations (A1) and (A2), were incorporated into a logistic regression model (Equations (A3) and (A4)) to produce seasonally-varying wetland maps.
where NDW I is the normalized difference water index, green is MODIS reflectance Band 4 representing green light, and MIR 1 is MODIS reflectance Band 6 representing mid-infrared light [111]. TC x = w x1 (red) + w x2 (N IR 1 ) + w x3 (blue) + w x4 (green) + w x5 (N IR 2 ) + w x6 (MIR 1 ) where TC x is the tasseled cap index value for index x (x = brightness, greenness, wetness) and w xi is the tasseled cap coefficient for MODIS band i. MODIS instrument specifications and TC index coefficients are presented in Table A1.
Log e P(X) 1 − P(X) = ω 0 + ω 1 x 1 + · · · + ω n x n = ω T X, P(X) ∈ (0, 1) , (A3) P(X) = exp(ω T X) (1 + exp(ω T X)) , where P(X) is the probability that wetland area is present within a given pixel based on predictor variables X and ω n is the regression coefficient for predictor variable x representing the variable's weight toward the resulting wetland presence probability [52]. Based on the intersection of the stationary wetland area line with the threshold-area curves for the three wet years (2010-2012) in Figure 3b, a 99% probability threshold was selected to determine where wetland areas were located.