Evaluation of Selected Sentinel‐2 Remotely Sensed Vegetation Indices and MODIS GPP in Representing Productivity in Semi‐Arid South African Ecosystems

The ability to validate satellite observations with ground‐based data sets is vital for the spatiotemporal assessment of productivity trends in semi‐arid ecosystems. Modeling ecosystem scale parameters such as gross primary production (GPP) with the combination of satellite and ground‐based data however requires a comprehensive understanding of the associated drivers of how the carbon balance of these ecosystems is impacted under climate change. We used GPP estimates from the partitioning of net ecosystem measurements (net ecosystem exchange) from three Eddy Covariance (EC) flux tower sites and applied linear regressions to evaluate the ability of Sentinel‐2 vegetation indices (VIs) retrieved from Google Earth Engine to estimate GPP in semi‐arid ecosystems. The Sentinel‐2 normalized difference vegetation index (NDVI), enhanced vegetation index (EVI) and the land surface water index (LSWI) were each assessed separately, and also in combination with selected meteorological variables (incoming radiation, soil water content, air temperature, vapor pressure deficit) using a bi‐directional stepwise linear regression to test whether this can improve GPP estimates. The performance of the MOD17AH2 8‐day GPP was also tested across the sites. NDVI, EVI and LSWI were able to track the phase and amplitude patterns of EC estimated gross primary production (GPPEC) across all sites, albeit with phase delays observed especially at the Benfontein Savanna site (Ben_Sav). In all cases, the VI estimates improved with the addition of meteorological variables except for LSWI at Middleburg Karoo (Mid_Kar). The least improvement in R2 was observed in all EVI‐based estimates—indicating the suitability of EVI as a single VI to estimate GPP. Our results suggest that while productivity assessments using a single VI may be more favorable, the inclusion of meteorological variables can be applied to improve single VIs estimates to accurately detect and characterize changes in GPP. In addition, we found that standard MODIS products better represent the phase than amplitude of productivity in semi‐arid ecosystems, explaining between 68% and 83% of GPP variability.


Introduction
The amount of carbon captured from the atmosphere by terrestrial plants during the process of photosynthesis is referred to as terrestrial gross primary production (Chapin et al., 2011).This is the largest global carbon flux and it drives important ecosystem processes such as respiration and growth (Beer et al., 2010).Importantly, the direct contribution of GPP to human welfare is critical since it is the basis of food, fiber and wood production (Beer et al., 2010;Chapin et al., 2011).In addition, by altering the amount of carbon dioxide in the atmosphere, the role of terrestrial plants in photosynthesis and respiration continue to play a crucial part in the carbon cycle, with around 30% of anthropogenic emissions being absorbed by terrestrial ecosystems each decade (Friedlingstein et al., 2022).Dryland ecosystems occupy about 40% of the global land surface and are expected to expand in the twenty-first century (Feng & Fu, 2013;Huang et al., 2017).These systems are responsible for about 50% of the trend in global terrestrial carbon sink with at least 40% of global terrestrial carbon inter-annual variability also attributed to these ecosystems (Ahlström et al., 2015).
Importantly, the opportunistic growth tactics of desert shrubs and ephemeral plants, including grasses and weeds, which are acclimated to highly fluctuating water availability, account for a significant portion of the contribution of dryland ecosystems to the inter-annual variability of the global carbon flux.For example, dryland vegetation in the southwestern region of the United States of America show high sensitivity to precipitation, where wet periods are associated with a tendency of these systems becoming CO 2 sinks, while becoming CO 2 sources during dry periods (Scott et al., 2015).With climate projections indicating that climate change is likely to deepen drought periods and intensify heatwaves in dryland ecosystems, this is a major threat to these water-limited environments (C.J. Engelbrecht & Engelbrecht, 2016;Engelbrecht et al., 2015;Huang et al., 2016).For instance, surface temperatures are projected to increase in southern Africa with the region also projected to become generally drier under low-mitigation climate change futures (Engelbrecht et al., 2009(Engelbrecht et al., , 2015)).Therefore, the global ability to understand the behavior, vulnerability and resilience of terrestrial ecosystems and their role as carbon sinks is a major focus of earth system science especially in the context of the responses of carbon sinks to the anthropogenic pressures such as agriculture, urban expansion and climate change (Laurance et al., 2014).
Eddy Covariance (EC) CO 2 flux measurements have offered a continuous and direct method to study the global carbon cycle.EC flux towers generally measure continuous exchanges of CO 2 , latent heat and sensible heat between ecosystems and the atmosphere, with global data archive lengths ranging from hours to decades, thus enabling the evaluation of seasonal and interannual variability in these exchanges while also elucidating their climatic controls (Baldocchi, 2003).A number of studies have applied the EC technique to study the temporal and spatial characteristics of the carbon cycle in African ecosystems (Archibald et al., 2009;Kutsch et al., 2008;Merbold et al., 2009;Räsänen et al., 2017;Scanlon & Albertson, 2004;Veenendaal et al., 2004).While these studies have focused on savanna woodlands across differing precipitation gradients, there is still a great need for representation of other biomes, especially in semi-arid, water-limited ecosystems.
The use of satellite observations in combination with climate data along EC data has been minimally explored in African dryland ecosystems (Abdi et al., 2017(Abdi et al., , 2019;;Jin et al., 2013;T. Kato et al., 2013;Sjöström et al., 2009Sjöström et al., , 2011Sjöström et al., , 2013)).Modeling GPP empirically by associating vegetation greenness indices like the Enhanced Vegetation Indices (EVI) and the Normalized Difference Vegetation Index (NDVI), as well as water-based indices such as the Land Surface Water Index (LSWI) to EC estimated GPP is a widely used approach (Sjöström et al., 2011).These greenness-based VIs have been effective in monitoring vegetation parameters such as phenology in semi-arid ecosystems (Archibald & Scholes, 2007;Higgins et al., 2011;A. Kato et al., 2021), as well as estimating productivity (Sjöström et al., 2009(Sjöström et al., , 2011(Sjöström et al., , 2013)).MODIS NDVI and EVI were found to be highly correlated with 8-day sums of EC GPP in a sparse savanna in central Sudan with R 2 values of 0.90 and 0.93, respectively, while the applicability of LSWI as a water-stress indicator was limited due to low fractional vegetation cover-making the retrieval of accurate information on canopy water content a challenge (Sjöström et al., 2009).
The relationship between GPP and EVI was particularly strong on a site-by-site basis and was further improved by combining EVI with tower-measured photosynthetically active radiation (PAR) and Evaporative Fraction (EF) (a measure of water sufficiency), with the relationship between GPP and EVI × PAR × EF at Skukuza, in South Africa observed with R 2 = 0.60 and being the most variable amongst the sites (Sjöström et al., 2011).Further studies looking at estimating productivity using remote-sensing driven models in dryland ecosystems in China continue to indicate that while remote sensing products can successfully capture temporal GPP for grasslands and shrublands, there are still large overestimation biases in amplitude present (Li et al., 2021).Similarly for a dryland site in North America (Yan et al., 2019), the limitations of greenness-GPP relationships especially over short timescales were observed as plant physiological functions change rapidly in response to biophysical drivers compared to greenness which changes more slowly-highlighting the potential and limitations of modeling GPP using remotely sensed greenness indices over dryland ecosystems.
The MOD17A2 GPP 1,000 m product was assessed for 12 African eddy sites where GPP seasonality was well captured over wet and dry sites but its amplitude was underestimated for many of the dry sites specifically in the Sahel (Sjöström et al., 2013).This underestimation persisted even in the updated MOD17A2H GPP (collection 6) product featuring an increased resolution to 500 m with two data input substitutes for meteorological data and fraction of photosynthetically active radiation (FPAR) products (Tagesson et al., 2017).This has been attributed to the maximum light use efficiency being set too low for semi-arid ecosystems in the MODIS algorithm (Tagesson et al., 2017).In addition, the performance of MOD17A2H was assessed on a global scale by comparisons with EC-derived GPP measurements and it continued to underestimate EC-derived GPP at most sites with poor performances in estimating both annual (R 2 = 0.62) and 8-day GPP (R 2 = 0.52) for sites featuring evergreen needleleaf forests, evergreen broadleaf forests, deciduous broad forests, mixed forests, grasslands and croplands (L.Wang et al., 2017)-necessitating the need for further updates and improvements to MODIS standard products (Zhu et al., 2018).
The contribution of dryland ecosystems to the global carbon cycle continues to be important, especially under a changing climate.With much of South Africa dominated by semi-arid ecosystems, which are expected to face warmer and drier periods under low mitigation efforts, it is vital to gain an improved spatial and temporal understanding of productivity.The strong relationship between water availability and the productivity of dryland ecosystems means that any likely changes in phenological events in these ecosystems such as shifts in the growing season length under the backdrop of these climate projections, there are likely impacts on the seasonal exchange of CO 2 , energy and water between terrestrial ecosystems and the atmosphere (Richardson et al., 2013).Being able to test the possibility to upscale site-specific EC measurements to regional scales using satellite-based data will offer such an understanding and provide further progress on satellite driven primary production modeling for semi-arid ecosystems as well as improve on data gap-filling approaches during periods of instrument failure.
The study intends to (a) compute site-based multilinear regression models to test the applicability of Sentinel-2 remotely sensed VIs (NDVI, EVI and LSWI) to track phase and amplitude patterns of EC estimated GPP and, (b) to test whether the combination of VIs (NDVI, EVI, LSWI) with in situ meteorological variables, namely: vapor pressure deficit (VPD), incoming solar radiation, soil water content (SWC) and air temperature can improve the prediction of GPP patterns and; (c) to assess the performance of a Light Use Efficiency (LUE)-based model (MOD17AH2) to simulate GPP of Savanna and Nama-Karoo sites differing in vegetation composition and annual precipitation in South Africa.

Site Characteristics
The study uses three (3) EC flux sites in South Africa where 2, namely: Benfontein Savanna (Ben_Sav) (28°5 3′26.2″S,24°51′40.0″E)and Nama-Karoo (Ben_Kar) (28°51′23.0″S,24°50′23.3″E)are located in the transition zone between the Savanna and the Nama-Karoo biomes, while the Middleburg site is located in the heart of the Nama-Karoo biome.These sites are illustrated in Figure 1 along others biomes in South Africa (Mucina & Rutherford, 2006).The Benfontein Savanna and Nama-Karoo towers were established in January 2020 and are located about 4 km from each other at the Benfontein Nature Reserve (BNR) along the border of the Free State and Northern Cape provinces of South Africa, about 15 km from the city of Kimberley.This area predominately receives summer rainfall with a mean annual rainfall of 419 ± 134 mm (Kamler et al., 2012), with a mean annual temperature of 18.5°C.Thunderstorms are common in summer and frost on occasion during winter.The savanna footprint represents Kimberley thornveld vegetation, with the vegetation being open camel thorn, Vachellia erioloba, and mix of mainly two grass species, namely Schmidtia pappophoroides and Stipagrostis uniplumis on deep Kalahari sands.The Nama-Karoo footprint is a fairly flat and homogeneous area dominated by Pentzia globosa.The periods of interest were data selected from the 4 February 2020 to the 10 March 2023 for the Ben_Sav and Ben_Kar sites.
The leniently grazed Middleburg site (Mid_Kar) (31°25′20.97″S,25°1′46.38″E) is located in the Eastern Upper Karoo vegetation type of the Nama-Karoo biome and was installed in October 2015.The site receives an average of 374 mm per year mainly falling during mid-to late-summer, with March usually being the wettest month, with a mean annual temperature of 15°C.The footprint of this leniently grazed site is dominated by Digitaria eriantha and Pentzia globosa over loamy soils.Data from the 4 February 2019 to the 10 March 2022 from Mid_Kar were used in this study and is referred to as the period of interest.

Instrumentation and Measurements
The EC is a micrometeorological technique that relies on high frequency measurements of the fluctuating components of vertical wind and a gas of interest (CO 2 ) in the constant flux region of the surface boundary layer to directly measure the net ecosystem exchange (NEE) of CO 2 (Baldocchi, 2003;Foken & Wichura, 1996;Madsen et al., 2009).NEE was directly measured using a Campbell Scientific IRGASON® Integrated CO 2 /H 2 O Open-Path Gas Analyzer and three-dimensional Sonic Anemometer (Campbell Scientific Inc., Logan, Utah, USA) at 10.5 and 3.5 m above ground at the Benfontein Savanna (Ben_Sav) and Nama-Karoo (Ben_Kar) sites, respectively.The data were acquired with a CR6 Campbell Scientific data logger at a frequency of 20 Hz and postprocessed into 30-min averages.The BNR sites were paired with identical instrumentation.Air temperature was measured with a HygroVUE™ 10 (Campbell Scientific, Logan, Utah, United States) probe at 2.0 m, while incoming solar radiation was measured using a CNR4 net radiometer (Kipp & Zonen, Delft, Netherlands).Soil moisture content measurements at 2.5 cm were averaged from two CS616 reflectometers (Campbell Scientific, Logan, Utah, United States) in two separate soil pits.Rainfall was measured using a TE525 tipping bucket rain gauge (Campbell Scientific, Logan, Utah, United States).
A Li-7200 enclosed path fast-response Infra-Red Gas Analyzer was coupled with a CSAT3 three-dimensional sonic anemometer (Campbell Scientific Inc., Logan, Utah, USA) to measure NEE at the Middleburg Karoo (Mid_Kar) site at a height of 2.5 m and a frequency of 20 Hz.Air temperature was recorded using a Vaisala HMP155 probe (Vaisala, Helsinki, Finland) while precipitation was measured using a TR 525 tipping bucket rain gauge (Texas Electronics, Texas, USA).Incoming solar radiation was measured using a CNR4 net radiometer (Kipp & Zonen, Delft, Netherlands).Soil moisture at 10 cm was measured using ML3x soil moisture probes (Delta T, EcoTech, Bonn, Germany).(Mucina & Rutherford, 2006).Panel 1 indicates the relative positions of the Benfontein Savanna and Nama-Karoo sites at the Benfontein Nature Reserve.Panel 2 shows the positioning of the Middleburg Karoo site.Panels (a-c) show the footprint extents of each site as well as the positioning of the polygons used to extract satellite data in relation to footprint extent.

Data Processing
The datalogger programs at the Ben_Sav and Ben_Kar sites were set to perform a double rotation coordinate system to transform the vertical and horizontal velocity components to zero and adjust for the misalignment of the sonic anemometer (Wilczak et al., 2001).Block averaging was applied to detrend turbulent fluctuations by removing the mean value from the time series (Kaimal et al., 1989;Massman, 2000).In addition, fluxes were corrected for the spectral losses in the low and high frequencies due likely to sensor responses, data processing choices as well as system characteristics (Moncrieff et al., 1997;Moore, 1986;Van Dijk, 2002).The Foken et al. (2012) quality control and flagging methodology was applied to all fluxes.These flux correction steps resulted in half-hourly values of CO 2 fluxes or NEE.The data post-processing steps for Mid_Kar are described in Rybchak et al. (2023).The half-hourly flux data across the sites were filtered according to the Foken et al. (2012) for Ben_Sav and Ben_Kar, while the Mauder and Foken (2011) was applied to discard bad quality data at Mid Kar.Data flagged between "0-4" were used for the Ben_Sav and Ben_Kar sites, while data flagged between "0-1" system were used for Mid_Kar-the remaining bad quality data were discarded.

Eddy Covariance Gross Primary Production Data
The measured NEE can be partitioned into fluxes of interest, namely: gross primary production and ecosystem respiration (ER).The REddyProc flux partitioning tool was used to partition NEE by applying the Lasslop et al. (2010) daytime-based method which uses the common rectangular hyperbolic light-response curve to estimate GPP (Falge et al., 2001;Lasslop et al., 2010) for all three sites (https://www.bgcjena.mpg.de/REddyProc/ui/REddyProc.php).The resulting half-hourly GPP estimates were aggregated to 5-day averages, while halfhourly meteorological variables, namely incoming solar radiation, SWC, air temperature and VPD, were also aggregated to 5-day averages to temporally match with Sentinel-2 VIs.The half-hourly GPP data were further aggregated to 8-day sums to temporally match the MOD17AH2 8-day GPP product.R statistical software version 4.2.3 (R core Team, 2023) was used to align GPP aggregations with the date sequences of Sentinel-2 and MOD17AH2 missions over each site-ensuring the data had the same temporal resolution.Following quality control of removing values out of range on the satellite data, these were replaced with the mean of the previous and the subsequent value in the each VI time series, with an average of six Sentinel-2 observations each month.

Ecosystem Footprint Calculation
Mean ecosystem footprints were predicted over each site's period of interest using the simple two-dimensional parameterization flux footprint model based on novel scaling approach for the crosswind distribution of the flux footprint (Kljun et al., 2015).The following variables were required for the footprint calculation: (zm) measurement height above ground in meters; (d) displacement height; (z0) roughness length; (u_mean) mean wind speed at measuring height; (L) Obukhov length; (σ v ) standard deviation of lateral velocity fluctuations after rotation; (u * ) frictional velocity; and wind direction.Polygons for the extraction of pixels for satellite data were created on the mean footprint dimensions to create regions of interest.The mean footprints for each site were in the dimensions of approximately 400 × 400 m (∼16 ha) for the Savanna site, 200 × 200 m (∼4 ha) for the Nama-Karoo site, and 200 × 100 m (∼2 ha) for the Middleburg Karoo site.

Satellite Data
The Copernicus Sentinel-2 Multispectral Instrument (MSI) Level-2A data were used for the computing of VIs (Copernicus, 2023).The Sentinel-2 MSI has a revisit time of 5 days while monitoring variability in land surface conditions.The Level-2A orthorectified atmospherically corrected surface reflectance data were retrieved from Google Earth Engine using square polygons with a 20 m negative buffer located at the center of the mean flux footprint over each study site.The negative 20 m buffer was set to minimize on the uncertainties of footprint and pixel mismatch wherein the mean value of the pixels covering the extent of the tower footprint were retrieved at each site.
A cloud filter was applied to the data for data with less than 20% cloud cover for the duration of the period of interest.A further filtering function was applied within the polygons to mask out non-vegetation and non-soil pixels.To remove the influence of cloud cover, the function "maskS2clouds" on Google Earth Engine enables cloudy and cloud-free pixels to be identified for both dense and cirrus clouds with an indicator specifying the cloud type.The method identifies cloud pixels based on high reflectance in the blue spectral region (information on this cloud mask can be found at https://developers.google.com/earthengine/datasets/catalog/COPERNICUS_S2_SR#description).The VIs detailed below were then retrieved by applying band calculations.

Normalized Difference Vegetation Index
VI are considered robust and empirical measurements of vegetation activity at the land surface, making use of measured spectral responses through the combination of multiple wavebands usually in the red (0.6-0.7 μm) and near-infrared (NIR) radiation wavelengths (0.7-1.1 μm) regions (Didan et al., 2015).One of the most common VI is the normalized difference vegetation index (NDVI): where NIR is reflectance in the near-infrared (NIR) band (band 8) and Red is the reflectance in the red band (band 4).NIR radiance is reflected by leaf cells since absorption of these wavelengths would result in overheating of the plant, while red radiance is absorbed by chlorophyll.The index is then normalized to improve the vegetation signal as well as to reduce atmospheric errors, solar zenith angles and sensor viewing geometry.A welldocumented problem of NDVI is the saturation at high biomass due to absorption of red light at ∼670 nm peaking at higher biomass loads.To counter this, the Enhanced Vegetation Index (EVI) is designed to maintain sensitivity in high biomass regions.

Enhanced Vegetation Index
In dealing with some of the atmospheric attenuation by aerosols, EVI minimizes the atmospheric effect by estimating the atmospheric influence level by the difference in blue and red reflectances.When aerosol concentrations are high, the difference between the two bands is larger and this influence is used to stabilize the index values against variations in aerosol concentration levels.By reducing atmospheric and soil background influences, EVI increases the vegetation signal and maintains sensitivity in high biomass regions (Huete et al., 2002).
where Blue is the reflectance in the blue band (band 2).L is the canopy background adjustment for correcting the nonlinear, differential NIR and red radiant transfer through a canopy; C 1 and C 2 are the coefficients of the aerosol resistance term and G is a gain or scaling factor.The coefficients adopted for the EVI algorithm are L = 1, C 1 = 6, C 2 = 7.5, and G = 2.5.

Land Surface Water Index
Leaf water content influences the measured leaf reflectance in several regions (0.4-2.5 μm) of the electromagnetic spectrum therefore shortwave infrared reflectance (SWIR) is negatively related to leaf water content.The region between 1.3 and 2.5 μm is the largest in the SWIR interval where increased reflectance in these wavelengths indicates a response to plant stress in general, including water stress.The LSWI uses the shortwave SWIR (band 12) and the NIR regions of the electromagnetic spectrum and is known to be sensitive to the total amount of liquid water contained in vegetation and the soil background (Chandrasekar et al., 2010), and is computed as the following: Since Sentinel-2 provides the NIR and SWIR bands at different spatial resolutions (10 and 20 m), the NIR band was resampled to 20 m on Google Earth Engine.Canopy water content is not the only parameter responsible for variations in reflectance in the SWIR channels since variations in leaf structure and dry matter content influence SWIR reflectance-effectively making the retrieval of vegetation water content using SWIR reflectance alone not ideal.Therefore, to isolate the contribution of these parameters to SWIR reflectance, the use of NIR information Journal of Geophysical Research: Biogeosciences which is affected by leaf structure and dry matter content but not by water is more useful for the combination of these signals to improve the retrieval of vegetation water content (Fensholt & Sandholt, 2003).

MOD17AH2 GPP Product
The MOD17AH2 Version 6 Gross Primary Productivity (GPP) product is provided by NASA LP DAAC at the USGS EROS Center and was retrieved using the Google Earth Engine (Running & Zhao, 2015).This collection 6 version features updated Biome Property look-up tables (BPLUT) as well as an updated version of the daily Global Modeling and Assimilation Office meteorological data.The product is based on the LUE concept (Monteith, 1972), which uses the relation between incident photosynthetically active radiation (PAR), the fraction of photosynthetically active radiation absorbed by plants (FPAR) and the actual light use efficiency term of the vegetation in the following model: where APAR, a product of FPAR and PAR, is the absorbed photosynthetically active radiation, where an 8-day estimate of FPAR is retrieved from MOD15 and daily estimates of PAR values are retrieved from the Global Modeling Assimilation Office reanalysis for each pixel (Running & Zhao, 2015).The PAR conversion efficiency, ε, values are derived from the attenuation from its maximum value, ε max , which is constrained by two environmental stresses, namely, (a) minimum temperature, T min , which can truncate GPP on days when air temperature is below 0°C and (b) VPD where high VPD can reduce stomatal conductance.Ultimately, the MOD17AH2 GPP product offers a cumulative 8-day composite of GPP values per 500-m pixel from the MODIS sensor onboard the Terra satellite.The polygons were applied here without the buffer since the extent of the footprints fell within one pixel of the product across all sites.A timeseries of the GPP pixel values in kg C m 2 were extracted at each site and converted to g C m 2 for consistent units between EC GPP and the MODIS GPP product.

Statistical Analysis
The EC estimated GPP, meteorological data and satellite-based VIs from each site were partitioned into training (70%) and (30%) testing sets.Simple linear regressions were computed for individual VIs with GPP EC to build initial VI-based models across each site on the training data set.To improve these models, a bi-directional (forward and backward) stepwise linear regression was used on the training set to build vegetation-index coupled with meteorological data from individual VIs (NDVI, EVI, LSWI) in combination with a varied selection of meteorological data (incoming solar radiation, air temperature, SWC, VPD (not available for Mid_-Kar)).To rank the applicability of these experimental models, the Akaike information criterion (AIC) was applied to select models with the lowest AIC score at each site during the model selection process at each site, for each VIbased model.The models with the lowest AIC score were selected, with their respective β coefficients for each VI and the selected supporting meteorological variables from each site.
The simulated data from both model experiments and observed data were assessed using the coefficient of determination (R 2 ) and residual standard error (RSE) to describe their performance against EC estimated gross primary productivity (GPP EC ) for both training and testing data sets.To evaluate the performance of all computed VI-based models, Taylor diagrams were computed to assess the correlation, the root-mean-square differences and the ratio of the variance for the test data set.The same performance test protocol was applied between MOD17AH2 against GPP EC evaluations.All analyses were performed using R statistical software version 4.2.3 (R core Team, 2023).Highest NEE, GPP and ER values were observed during the wet season, with NEE however remaining above 5 g C m 2 day 1 (negative values indicate more carbon uptake by the ecosystems, while positive values indicate carbon release into the atmosphere) across all sites, except for the 2021 wet season at the Ben_Sav site.GPP and ER were highest during the wet season at Ben_Sav, followed by Ben_Kar and Mid_Kar, reaching values above ∼6 g C m 2 day 1 across Ben_Sav and Ben_Kar.GPP and ER values were the lowest at Mid_Kar (below 4 g C m 2 day 1 ), with the 2020 wet season being the most productive compared to the other years.Air temperature was closely matched at the Benfontein sites and remained within range of the Mid_Kar site.Ben_Sav and Ben_Kar were observed to receive more rainfall compared to Mid_Kar for year-on-year comparisons.SWC was variable across the Benfontein sites, due likely to soil properties, while SWC at 10 cm at Mid_Kar was considerably higher during dry periods compared to the top layers (2.5 cm) measured at the Benfontein sites.NDVI, EVI and LSWI estimates also followed the seasonality observed for ecosystem productivity, rainfall, SWC and air temperature across the sites.

Correlation Matrix Between GPP EC Versus Vegetation Indices and Meteorological Variables
Pearson correlations were computed between GPP EC and VIs and meteorological variables for the training portion of the data set (70%) as shown in Table 1.GPP EC was positively correlated with all variables, except for VPD at Ben_Kar, with the strongest correlations observed between GPP EC and the greenness-based VIs; where r ranged between 0.85 (Ben_Sav) and 0.86 (Ben_Kar) for GPP EC and NDVI and, ranged between 0.88 (Ben_Sav) and 0.91 (Ben_Kar) for EVI.
Strong positive correlations were observed between GPP EC and LSWI across all sites, especially for the Nama-Karoo type vegetation sites (0.84 for Ben_Kar and 0.83 for Mid_Kar) compared to the savanna-type site (0.72).SWC correlations with GPP EC were however not as strong as with LSWI at Ben_Sav (0.46) compared to Ben_Kar (0.67) and Mid_Kar (0.76), with no major redundancy observed between LSWI and SWC across all sitescorrelations ranged between 0.34 (Ben_Sav) and 0.65 (Mid_Kar).GPP EC correlations with T air and R g followed similar patterns across the sites, where the highest correlations were observed at Ben_Sav for T air and R g (0.70 and 0.60), followed by 0.56 and 0.40 at Ben_Kar while the lowest were observed at Mid_Kar (0.36 and 0.19) for T air and R g , respectively.On all occasions, the strong positive correlations were significant (p < 0.001).

Model Training for VIs
The applicability of NDVI, EVI and LSWI as well as the inclusion of meteorological variables to improve the estimates was tested.Notably, all the VIs were able to track the seasonal patterns of GPP across the sites.However, there was a temporal delay and periods of overestimation in GPP NDVI and GPP EVI estimates particularly in periods approaching the peak of the wet season as well as the periods leading up to the dry season especially at Ben_Sav and Ben_Kar.LSWI estimates of GPP performed poorly at Ben_Sav and Ben_Karoverestimating GPP during the dry season and lagging behind wet season increases especially at Ben_Sav.Nevertheless, LSWI was able to track the fire period during September 2021 at both Ben_Sav and Ben_Kar where sharp decreases in GPP LSWI were observed-indicating substantial losses in vegetation, especially at Ben_Sav where the grass layer was completely burnt from the fire.
The model selection results are shown in Table 2, where the included meteorological variables for each improved model are indicated across all the sites with their respective β coefficients.The improved VI-based GPP estimates with selected meteorological variables are referred to as GPP NDVI , GPP EVI and GPP LSWI (an addition of an underscore from the previous estimates).The structure of the applied equation for each VI and selected meteorological variables for each site is indicated as follows: where β 0 indicates the intercept of the equation, β VI indicates the selected VI, β Tair indicates the coefficient for temperature, β Rg indicates the incoming solar radiation coefficient, β VPD indicates the VPD coefficient (except for Mid_Kar), and β SWC indicates the coefficient for SWC for each site as reported in Table 2.
The simulated GPP NDVI , GPP EVI estimates indicated considerable improvement particularly in phase than in amplitude agreement with GPP EC at Ben_Sav as shown in Figure 3.While both phase and amplitude improvements were observed for GPP LSWI at Ben_Sav, amplitude disagreement between GPP EC peaks during the third wet season and GPP EVI and GPP LSWI was a still notable limitation at Ben_Sav.The simulated estimates at Ben_Kar and Mid_Kar matched in phase and in amplitude with GPP EC during the training process.Linear regressions were computed for NDVI, EVI and LSWI as well as GPP NDVI , GPP EVI and GPP LSWI with GPP EC at each site and the results for estimates with and without meteorological variables for the complete period and between the wet and dry seasons are shown in Table 3.More than 80% of GPP EC variability was explained by the improved models across all the sites.

Model Evaluation
All the models were evaluated on the remaining 30% and were tested against GPP EC across the sites and the time series of the results are shown in Figure 4. Considerable phase agreement was observed in all model estimates across the sites.Notably, underestimation of GPP EC amplitude was observed for GPP NDVI at Ben_Sav (R 2 = 0.77;

Table 3 Coefficient of Determination and Residual Standard Error Between GPP EC and GPP Estimates From Vegetation Index (VI)-Based Estimates (GPP NDVI , GPP EVI , and GPP LSWI ) and Model Estimates Improved With Meteorological Variables (GPP NDVI , GPP EVI and GPP LSWI ) for Complete Periods and Between Wet and Dry Seasons During the Training Period
Site Period RSE = 0.54 g C m 2 5 day 1 ) and Ben_Kar (R 2 = 0.87; RSE = 0.74 g C m 2 5 day 1 ).GPP EVI better estimated GPP EC across the sites compared to GPP NDVI in terms of amplitude, however initial periods leading up to the wet season were overestimated at Ben_Sav.GPP LSWI was within phase at Ben_Sav but overestimated during the dry season while underestimating during the wet season.While all the models without the inclusion of meteorological variables were in phase agreement at Mid_Kar, they underestimated wet season GPP EC .
Considerable improvements were observed in the performance of all models with the inclusion of meteorological variables at Ben_Sav in terms of amplitude agreement with GPP EC however GPP NDVI underestimation continued still from the start of the wet season.GPP EVI performed well overall, explaining 94% (RSE = 0.61 g C m 2 5 day 1 ) of the variability in GPP EC at Ben_Sav, while GPP NDVI and GPP LSWI both explained 91% (RSE = 0.51 g C m 2 5 day 1 ) and 90% (RSE = 0.64 g C m 2 5 day 1 ), respectively.GPP NDVI showed amplitude improvements too at Ben_Kar, however wet season GPP EC was still underestimated (R 2 = 0.79; 0.62 g C m 2 5 day 1 ).Negligible improvements were observed for GPP EVI and GPP LSWI at Ben_Kar.
The performance of GPP NDVI at Mid_Kar was the weakest compared to the other sites-explaining only 54% (RSE = 0.58 g C m 2 5 day 1 ) of variability in GPP EC .The performance of GPP EVI was however better than GPP NDVI , explaining 68% (RSE = 0.57 g C m 2 5 day 1 ) of GPP EC variability, while GPP LSWI explained only 62% (RSE = 0.63 g C m 2 5 day 1 ) of variability in GPP EC -a reduction in R 2 from GPP LSWI .During periods of low GPP during the dry season, all the improved models estimated values below zero at Mid_Kar.Coefficient of determination and RSE results for the evaluated estimates with and without meteorological variables for the complete period and between the wet and dry seasons are shown in Table 4.
Taylor Diagrams were computed to illustrate how three complementary model performance statistics: correlation coefficient, the standard deviation and the centered root-mean-square error vary simultaneously across the sites as shown in Figure 5.The GPP NDVI , GPP EVI GPP LSWI as well as GPP NDVI , GPP EVI and GPP LSWI predictions were compared against GPP EC observations.At Ben_Sav, the closest simulated data to the observed values were from GPP NDVI in terms of standard deviation, a high correlation and a CRMSE close to zero compared to all other models.While GPP LSWI , GPP EVI and GPP EVI were also in close proximity, GPP LSWI and GPP EVI were lesser correlated and GPP EVI had a higher standard deviation.
GPP EVI and GPP EVI showed the highest standard deviation compared to observed GPP EC values, while GPP NDVI ; GPP NDVI and GPP LSWI estimates were clustered, indicating less similarity with observed values in terms of their

GPP EC Versus LUE-Based MODIS MOD17AH2 GPP 8-Day Cumulative Product
MODIS MOD17AH2 GPP 8-day product (GPP MODIS ) versus EC derived GPP 8-day sums (GPP EC ) for February 2020 to January 2022 were compared for Ben_Sav, Ben_Kar and Mid_Kar as shown in Figure 6.MODIS_GPP was observed to follow the seasonal cycles of GPP EC throughout the measuring period.However, moments where GPP EC receded during the fire at the end of September 2021 were not tracked by GPP MODIS at Ben_Sav and Ben_Kar-indicating some limitations in GPP MODIS capturing key ecosystem disturbances especially over a short period.Overall the relation between the two data sets indicated that variability in GPP EC was captured by GPP MODIS across all sites, however there were periods of disagreement during both the wet and dry seasons.
Coefficient of determination values (R 2 ) were 0.68 at Ben_Sav, 0.70 at Ben_Kar-indicating moderate fit between the GPP MODIS product and GPP EC 8-day cumulative sums.An R 2 value of 0.83 was observed at Mid_Kar, with a better fit between the GPP MODIS product and 8-day sums of GPP EC as shown on Figure 7. Differences in correlation, standard deviation and RSE error between the cumulative GPP EC and GPP MODIS values were simultaneously compared using Taylor diagrams for each site as shown in Figure 8.The data were normalized for the wet and dry season comparison by dividing both the RMS difference and the standard deviation of GPP MODIS by the standard deviation of GPP EC .While good correlations between GPP MODIS and GPP EC were observed across the sites overall (complete), there were differences between GPP EC and GPP MODIS and were attributed to standard deviation differences on a seasonal basis.The largest standard deviation differences were mainly observed during the wet season for GPP MODIS and GPP EC values at Ben_Sav, while dry season GPP EC and GPP MODIS values were moderately related.This can be attributed to the substantial underestimation of GPP EC peaks by GPP MODIS during the wet season leading to increased variability between the two data sets, albeit they are temporally in good agreement.At Ben_Kar, MODIS_GPP  The 8-day GPP MODIS and GPP EC sums were cumulated for the whole study period at each site, where GPP MODIS had a cumulative sum of 1,756 g C m 2 versus 3,244 g C m 2 by GPP EC between 2 February 2020 and 10 February 2023 at Ben_Sav (difference of 1,488 g C m 2 ) as shown in Figure 9.The cumulative sums between the two data sets showed a lesser difference (840 g C m 2 ) at Ben_Kar where 1,622 g C m 2 was recorded for GPP MODIS compared to 2,462 g C m 2 for GPP EC for the same period at Ben_Sav (with two months of missing data between May and June at Ben_Kar).A difference of 483 g C m 2 was observed at Mid_Kar where a total of 1,627 g C m 2 was estimated by GPP MODIS during 2 February 2019 to 6 March 2022 compared to 1,144 g C m 2 by GPP EC .

Discussion
In this study we first set out to test the applicability of 5-day Sentinel-2 VIs retrieved via Google Earth Engine, namely NDVI, EVI and LSWI in computed regression models to validate against EC estimated GPP, and second, to test whether the inclusion of selected meteorological variables, incoming solar radiation, air temperature, SWC and vapor pressure deficit (VPD not available at Mid_Kar) to greenness-based VIs (NDVI, EVI) and a waterbased index (LSWI) can improve the performance of the estimates.Lastly, we compared 8-day sums of a LUE-based GPP model (MOD17AH2) versus EC estimated GPP during the wet and dry seasons in semi-arid ecosystems in South Africa.
We found very strong to moderate correlations between daily GPP EC with NDVI and EVI and LSWI across the sites-indicating potential for modeling GPP in these semi-arid ecosystems using sentinel-2 VIs.Typically, greenness and water based indices are highly correlated with changes in aboveground biomass and canopy structure and have also been shown to be highly correlated to gross primary production (Biudes et al., 2021;Browning et al., 2017;Sjöström et al., 2009).However, the detection of vegetation senescence in dryland systems may be challenging as some drought deciduous species continue photosynthesizing by shedding only some of their leaves during the dry season (Jolly & Running, 2004) while others may shut down photosynthesis during drought while still carrying all their leaves (Smith et al., 2018).

Vegetation Index-Based Models to Predict GPP EC
The validation of VIs with EC estimated GPP showed that both greenness-based and water-based indices were able to generally track the phase and amplitude of GPP during both dry and wet seasons better in the grass-shrub homogenous sites of Nama-Karoo vegetation compared to the mixed Savanna site.Archibald and Scholes (Archibald & Scholes, 2007) found that this is likely due to Savanna trees having a less variable phenological cycle compared to grasses.Therefore an ecosystem footprint would represent a mixed ecosystem scale signal where not all plants would react in unison such as in a mixed savanna site.
In addition, GPP EC is affected by PAR and hydrometeorological conditions which vary highly within a day such as humidity, temperature and soil moisture which affects leaf level gas exchange (Running et al., 2004;Sims et al., 2006).While on the other hand, vegetation greenness depends on factors which vary slower than the abovementioned hydrometeorological conditions, namely, leaf chlorophyll content and canopy structure (Carlson & Ripley, 1997;Gitelson & Merzlyak, 1997), thus having a delayed impact on plant gas exchange regulation.Therefore, when grasses in a semi-arid savanna ecosystem stop photosynthesizing during drought or cold temperatures while the leaves are green, GPP EC may begin to decline while the trees are still relatively active.There might be an apparent lag when satellite sensors detect such changes thus contributing to the apparent lag in comparing fast-response EC estimated GPP with satellite sensor estimated GPP, especially from a savanna site.
Similar instances were observed where MODIS EVI (8-day, 500 m) and flux tower GPP were out of phase during the same phenological stages in grassland, shrubland and mixed savanna sites (Sjöström et al., 2011).However, in our results, for the shrub dominated sites, NDVI and EVI based models were in phase and within amplitude of GPP EC , which is likely due to the higher spatial and temporal resolution of Senintel-2 data sets.LSWI, which is sensitive to leaf water content while retaining a correlation with biomass, was also out of phase with GPP EC most notably at Ben_Sav, and often overestimated GPP EC across all sites during the dry season.This could be due to how LSWI correlates well with vegetation and soil moisture content (Zhao et al., 2009), especially when vegetation is sparse (such as during the dry season).This in contrast to A. Kato et al. (2021) where LSWI better detected phenological stages of dryland ecosystems compared to NDVI and EVI.This can be attribute to the sensitivity of LSWI to also soil moisture changes rather than only vegetation green-up signal thus affecting LSWI and GPP relationships in sparse vegetation situations (Chandrasekar et al., 2022).For example, correlations between LSWI and SWC improved in Nama-Karoo vegetation types compared to a Savanna type vegetation type in the current study.

Combination of Meteorological Variables and Vegetation Indices to Predict GPP EC
The addition of meteorological variables improved the performance of all models except for GPP LSWI at Mid_Kar.The models explained between 62% and 94% of GPP EC variability across the sites when tested compared to between 19% and 85% from the VI models without the inclusion of meteorological variables.While EVI performed better than NDVI overall in terms of R 2 values, NDVI and LSWI showed the most improvement with the inclusion of meteorological data-highlighting the suitability of EVI to perform better as a single VI in estimating GPP.The performance of greenness-based VIs such as NDVI and EVI in estimating GPP has been found to be better in deciduous sites than evergreen sites (Shi et al., 2017;Zhou et al., 2022), with their performance constrained by environmental conditions and features of both climate and vegetation structure playing an important role.These limitations have been highlighted by Smith (Smith et al., 2018) where a phase disagreement in the detection of vegetation senescence was observed and greenness indices estimated a later date of the end of photosynthesis than that observed from ground-based measurements of GPP.
This study also observed that the initial phase disagreements between NDVI, LSWI and GPP EC at Ben_Sav considerably improved with the addition of meteorological variables.The constraining of both NDVI and LSWI by T air , R g and VPD indicates that the addition of meteorological data to fine resolution satellite data provides a convenient approach to estimating site based GPP.However, in situ data may not always be available, thus being a notable limitation in some situations.Nevertheless, these models can be an alternative to LUE models or terrestrial biosphere models which are often reliant on course resolution meteorological inputs and provide a challenge in terms of parameterization (Shi et al., 2017).However, LUE models still present advantages in that while VIs such as NDVI and EVI are highly correlated with the light use efficiency, the structure of multiple linear regressions compared to LUE models lacks the environmental stress reductions on LUE, such as temperature and vegetation water stress and the various ways they are integrated into LUE-based models (Zhang et al., 2015).Second, regarding the nonlinearity of relationships between GPP, R g and VPD and how they are modulated by available water (SWC) within and between seasons, using multilinear regressions for estimating GPP in this case may oversimplify the estimation of GPP even with the inclusion of meteorological variables.

Evaluation of LUE-Based Model Versus GPP EC
GPP MODIS performed reasonably well across the sites in terms of following the temporal patterns of GPP EC , however there were amplitude differences observed between seasons across the sites.In most cases, GPP MODIS underestimated GPP EC especially during the wet season at Ben_Sav resulting in the observed large difference in standard deviation between GPP EC and GPP MODIS during the wet season, while dry season values were generally within range.This is consistent with other sites where GPP MODIS generally underestimated GPP EC by about 34% across 15 Fluxnet sites (L.Wang et al., 2017), with strong underestimation also observed in semi-arid ecosystems in the Sahel (Tagesson et al., 2017) where about 67% variability of GPP was explained by GPP MODIS .About 60% of GPP EC variability was explained by GPP MODIS at scrub site in the Mexican highland (Guevara-Escobar et al., 2021).These are within range with Ben_Sav and Ben_Kar (68%-70%), while at Mid_Kar (83%), GPP MODIS performed considerably better.
While the study was able to relate EC estimated GPP with satellite retrieved estimates and find reasonable agreement between the two data sets over an 8-day scale, there were several differences and limitations to these findings.Notably, the mismatch in spatial extents of the EC footprint and the 500 m pixel representing MODIS GPP product can explain some of the differences in observed amplitude differences in GPP, especially during the dry season at Mid_Kar.Here, the pixel average of 500 m might be too coarse to compare with tower GPP data particularly over a site with grazing paddocks.In addition, this coarse resolution can be challenging in detecting changes in GPP such as sporadic green-up events where surface vegetation is sparse with a matrix of more than one vegetation type within the footprint-characterizing much of semi-arid ecosystems (Chu et al., 2021).In other cases, challenges can be presented in situations where the MODIS footprint is much larger or smaller than the flux tower footprint, the VI signal may be made up of different vegetation sections, introducing the likelihood of different phenology patterns than those captured within the footprint of the EC flux tower.
Limitations regarding the application of the MOD17A2H GPP product have been documented (L.Wang et al., 2017).Amongst these are that existing errors from input data in the MOD17AH2 product which lead to estimation inaccuracies originate from FPAR than the meteorological data.These however can be optimized by using site-specific data where the product improved to explain 91% of GPP variability in arid and semi-arid ecosystems in China (H.Wang et al., 2019).Another issue is that the land cover classification has proved to have frequent misclassification errors, for example, evergreen needleleaf forests were mistook for savannas while grasslands were mistook for open shrublands (L.Wang et al., 2017).The maximum light use efficiency term used in the product was also found to be smaller than the recalculated inferred LUE max for the majority of land cover types.While the evaluation of MOD17AH2 is based on the assumption that the GPP derived from EC flux tower measurements are reliable, there are multiple sources of uncertainty associated with EC derived GPP.These include systematic and random errors associated with instrument precision, calibration and placement as well turbulent transport and statistical errors relating to footprint heterogeneity (Hollinger & Richardson, 2005;Loescher et al., 2006).While aggregating NEE from days to seasons reduces some random errors, the procedures involved in the filtering, filling gaps in the data and the partitioning of NEE into GPP and Re are associated with additional errors.

Conclusion
The greenness-based VIs were able to track the general patterns of GPP EC across all sites, albeit temporal delays between GPP EC and all VIs were observed, particularly at Ben_Sav.The inclusion of meteorological variables to Journal of Geophysical Research: Biogeosciences 10.1029/2023JG007728 support single VI estimates of GPP proved most effective at Ben_Sav in reducing the temporal gaps between GPP EC and each VI (especially NDVI and LSWI); while the predictive strength of NDVI was also substantially improved at Mid_Kar.The least improvement in terms of predictive power with the inclusion of meteorological variables was observed for EVI across all sites-indicating its suitability for estimating GPP in semi-arid ecosystems even without supporting variables.The notable limitation for EVI was a slight temporal delay with GPP EC for the savanna site.While productivity assessments using a single VI may be more favorable, our results suggest that the coupling with meteorological variables could be applied to improve single VIs and accurately detect and characterize productivity transitions especially for NDVI and LSWI.It is this coupling that models including meteorological variables provided better estimates of GPP EC in terms of standard deviation and correlation when tested across the sites.
The use of a single VI such as EVI has however shown to be suitable for both homogenous vegetation types such as Nama-Karoo as well as mixed ecosystems such as the savanna.While multilinear regressions may not be ideal for the integration of VIs with meteorological variables, the study has demonstrated the potential of site-based satellite driven models for scaling up GPP in semi-arid ecosystems and for the likely use of such models as potential site-based GPP gap-filling strategies during lengthy periods of instrument failure.The applicability of the MODIS GPP product in these semi-arid ecosystems was proven to be within phase with EC-derived GPP across all sites but differed in amplitude with issues of overestimation and underestimation for wet and dry season peaks attributed mainly to pixel-footprint mismatches and BPLUT not in agreement with ecosystem properties.Figure 9 illustrates some of the limitations of the MODIS GPP product where for sites that are observed to function differently, are characterized similarly by the product-indicating less variability in cumulative GPP for the period of interest.These biases are highlighted here to illustrate that while this product provides the necessary spatial coverage for vegetation assessments, the BPLUT may still be a notable source of uncertainty when applying this product.This carries significant implications for global vegetation productivity assessments, especially for semi-arid ecosystems where in this study the MODIS GPP product has almost underestimated productivity of these sites by about 50%.In addition, we highlighted the value of in situ meteorological data in improving model performance-a feature that the MODIS GPP product has not yet achieved.
In addressing the in situ data limitations noted in the study on a continental scale, there is a great need to establish and maintain a network of validation sites in semi-arid ecosystems especially due to their growing influence on the interannual variability of the global carbon cycle.This will allow such work to continue reducing the uncertainty of modeling productivity in semi-arid ecosystems.In addition, the relative contribution of trees versus grasses to the measured NEE will need to be partitioned to understand the phenology and functional differences between grasses and trees within the footprint of the Savanna site in the aim to explain some of the lags observed between satellite and EC data.Thus, continuing to study patterns of productivity through incorporating satellite observations and in situ data will improve the ability of planning and management to rely on accurate vegetation assessments to understand the varying controlling effect of ecosystem drivers such as fire, rainfall, incoming radiation and air temperature under varied ecosystem conditions and the reality of climate change.
the notebooks for the analysis and plotting of figures in R software version 4.2.3 are available at https://zenodo.org/records/10650562.

Figure 1 .
Figure 1.Biomes of South Africa and the location of selected study sites (Mucina & Rutherford, 2006).Panel 1 indicates the relative positions of the Benfontein Savanna and Nama-Karoo sites at the Benfontein Nature Reserve.Panel 2 shows the positioning of the Middleburg Karoo site.Panels (a-c) show the footprint extents of each site as well as the positioning of the polygons used to extract satellite data in relation to footprint extent.
Time series of meteorological, eddy-covariance, and satellite data during the periods of interest are shown in Figure2.Seasonal dynamics of GPP, NEE and ER illustrate the wet and dry season periods of productivity across the ecosystems which are associated with seasonal rainfall inputs.The wet season is between October and April -peaking in February and March for Ben_Sav and Ben_Kar with the dry season dominating between May and early September.Ecosystem productivity at Mid_Kar is observed during the wet season periods of November and April-peaking also during February and March.While there are periods of early rains during the dry season across all sites, this does not translate to immediate net ecosystem productivity.While the majority of the rainfall is observed during November and April across the sites, there are periods of rainfall during the dry season, especially during the 2021/2022 dry season (associated with La Niña conditions).

Figure 3 .
Figure 3.Time series from training of vegetation index based models (GPP NDVI , GPP EVI , and GPP LSWI ) and improved regression estimates (GPP NDVI , GPP EVI , and GPP LSWI ) of GPP versus GPP EC estimates for Ben_Sav (left panels), Ben_Kar (middle panels) and Mid_Kar (right panels).In each horizontal panel, the black lines represent GPP EC , the red lines represent vegetation index (VI)-based models and the blue lines represent the improved model estimates.

Figure 4 .
Figure 4. Time series from testing of vegetation index based models (GPP NDVI , GPP EVI , and GPP LSWI ) and improved regression estimates (GPP NDVI , GPP EVI , and GPP LSWI ) of GPP versus GPP EC estimates for Ben_Sav (left panels), Ben_Kar (middle panels) and Mid_Kar (right panels).In each horizontal panel, the black lines represent GPP EC , the red lines represent vegetation index (VI)-based models, and the blue lines represent the improved model estimates.

Figure 5 .
Figure 5.Taylor diagrams to compare all computed vegetation index (VI) model predictions versus GPP estimates from eddy covariance flux towers.Left is Ben_Sav, middle is Ben_Kar and right is Mid_Kar.

Figure 6 .
Figure 6.8-Day sums of GPP MODIS estimates versus 8-day sums of GPP EC across the sites.Red lines indicate GPP EC and blue lines indicate GPP MODIS .Left is Ben_Sav, middle is Ben_Kar and right is Mid_Kar.

Figure 7 .
Figure 7. Linear regression of 8-day sums of GPP MODIS estimates versus 8-day sums of GPP EC across the sites during dry (red) and wet seasons (blue).Left is Ben_Sav, middle is Ben_Kar and right is Mid_Kar.

Figure 8 .
Figure 8.Taylor diagrams to compare GPP MODIS versus observed GPP EC estimates from flux towers during the wet and dry seasons.Left is Ben_Sav, middle is Ben_Kar and right is Mid_Kar.

Table 2
Akaike Information Criterion (AIC) Model Selection Summaries and Coefficients of All Included Variables for Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI), and Land Surface Water Index (LSWI) Models Coupled With Meteorological Variables MALULEKE ET AL.

Table 4
Coefficient of Determination and Residual Standard Error Between GPP EC and GPP Estimates From Vegetation Index (VI)-Based Estimates (GPP NDVI , GPP EVI , and GPP LSWI ) and Model Estimates Improved With Meteorological Variables (GPP NDVI , GPP EVI and GPP LSWI ) for Complete Periods and Between Wet and Dry Seasons During the Testing PeriodThe number of asterisks indicate the level of significance (p < 0.001***, p < 0.001**).