Remote sensing of water use and water stress in the African savanna ecosystem at local scale – Development and validation of a monitoring tool

Savannas are among the most productive biomes of Africa, where they comprise half of its surface. They support wildlife, livestock, rangelands, crops, and livelihoods, playing an important socioeconomic role in rural areas. These water-limited ecosystems with seasonal water availability are highly sensitive to changes in both climate conditions, and in land-use/management practices. Although monitoring programs for African savanna water use have been established in certain areas, most of them are largely restricted to point based measurements or coarse scales, and are not fully capable to provide distributed timely information for planning purposes. In this study we develop a mechanism for monitoring the water used by African savanna from fine scale (meters) to watershed scale, integrating the effects of the water stress. Our hypothesis is that the Ecosystem Stress Index (ESI) is a valuable tool to downscale estimates of actual evapotranspiration at coarse scale, to high resolutions. To monitor savanna water fluxes in a semi-continuous way this study integrates two different ET-estimation approaches: KC-FAO56 model, integrating reflectance-based “crop” coefficients (SPOT 4 & 5 satellites), is used to derive unstressed savanna evapotranspiration (with high spatial resolution), and the two-source surface energy balance model -TSEB, integrating radiometric surface temperature (AATSR satellites) allows the determination of water stress across savannas (ESI, with low spatial resolution). The difference between estimated and observed surface fluxes derived from TSEB (RMSDLE= 53 Wm-2, RMSDH=50 Wm, RMSDRn= 60 Wm, RMSDG=21 Wm) were of the same magnitude as the uncertainties derived from the flux measurement system, being sufficiently accurate to be employed in a distributed way and on a more regular basis. The approach of ESI to downscale ET proved to be useful, and errors between estimated and observed daily ET (RMSD 0.6 mmday−1) were consistent with the results of other studies in savanna ecosystems. The modelling framework proposed provided an accurate representation of the natural landscape heterogeneity and local conditions, with the potential of providing information suitable from local to broader scales.


Introduction
The aim of this paper is to outline, apply and validate a novel tool to account for water use and water stress in African savannas at fine scale, using Earth Observation data (low to high spatial resolution sensors), to provide accurate and timely information for supporting decisionmaking at different scales. The improvement of rangeland management can increase savanna productivity and resilience, allowing the population depending on these ecosystems to adapt progressively to climate change, and helping to ensure food security by diversifying the current agricultural sector. Savannas are among the most productive biomes of Africa (Scholes and Walker 2004;Beerling and Osborne 2006), where they comprise half of its surface (Sankaran and Ratnam 2013). They support wildlife, livestock, rangelands, crops, and livelihoods (Scholes and Archer 1997), playing an important socioeconomic role in rural areas. Savannas provide valuable ecosystem services at local, regional, and global scales. They greatly influence global land-surface processes and earth water and carbon cycles, maintaining biodiversity, improving soil fertility, and maintaining regional hydrological balance (Kamaljit 2006). In South Africa, where this work was conducted, savannas cover about a third of the land area, providing two vital services: wildliferelated tourism and cattle ranching. It is estimated that around 9.2 million of South Africans benefit directly from these ecosystem services, through the use of savanna resources ( van Wilgen 2010).
There are many definitions for the term "savanna". For the purpose of this study however, they will be defined as grasslands with trees and shrubs, ranging from almost pure pastures to more closed woodlands. The herbaceous layer is almost continuous, with bare soil patches and an uneven cover of trees and bushes, which allows sufficient light penetration to support the understory vegetation (Scholes and Archer 1997). These ecosystems are located between deserts and tropical rainforests (Asner et al., 2009;Asner and Levick, 2012), but this study is focused on semi-arid savannas, which are characterised by a distinct dry season with no rain and high temperatures, where droughts naturally occur. Savanna distribution and vegetation structure (e.g., characteristics and spatial patterns of vegetation, dominant species, water and soil nutrient status, etc.) are mainly determined by anthropogenic effects and; the effect of herbivores; fire; soil characteristics; and water availability, the most decisive factor in semi-arid climates (Sankaran and Ratnam 2013;van Wilgen 2010).
These water-limited ecosystems with seasonal water availability are highly sensitive to changes in both climate conditions and in land-use/ management practices (Schnabel et al., 2013;IPCC 2014;Kueppers et al., 2005;Bond et al., 2008), and water scarcity situations have been aggravated by global warming. Most significant human-induced changes, such as invasive species (changes in vegetation structure), fires (influencing canopy recovery), agriculture pressure and livestock overgrazing (e.g. soil compaction and erosion), indiscriminate hunting, or over-exploitation of resources, have severe implications on the regional water balance, biodiversity and ecosystem productivity (Matongera et al., 2016). Since savannas are highly influenced by human activities, private and institutional practices play a key role in their conservation, while also improving people's living standards and providing employment opportunities. This implies that integrated management strategies for savannas are required, which explicitly consider water resources management in an holistic manner and take the multifunctionality of the ecosystem and needs of all stakeholders into account (Lal, 2015;Schwärzel et al., 2016). There is then a need to develop a mechanism for monitoring water availability and vegetation dynamics in savannas, at both the regional and local scales. Earth Observation (EO) data, particularly data provided by Sentinel 2 and Sentinel 3 satellites, new missions developed by European Space Agency (ESA), allow us to map the water use, and water stress, and the vegetation distribution across the African savannas, as well as to monitor seasonal and long-term temporal variations.
Although monitoring programs for savanna water use have been established in certain areas, most of them are largely restricted to point based/small scale measurements, which are not capable of providing distributed timely and up-to-date information for water resources management and planning purposes. African savanna ecosystem management, with natural vegetation susceptible to water stress, should take advantage of remote sensing advances for estimating water use applying ET as a proxy. The transferability of the findings to African ecosystem remains uncertain and should be tested in-situ. Ramoelo et al. (2014) and Majozi et al. (2016) have evaluated MOD16 global ET product and different ET algorithms over South African savanna sites (equipped with eddy covariance towers). Their results demonstrated the usefulness and potential of EO data for savanna water monitoring. In addition to these results, which used MODIS products with a spatial resolution of 1 km 2 and that did not constrain vegetation water stress, Fig. 1. Proposed modelling framework scheme. On the upper part of the figure K C -FAO56 Model is described. On the down part Two-Source Energy Balance Model is described. In the right, the integration of both models to derived actual ET at high resolution is presented. Symbols: K C is the crop coefficient derived from FAO56, Ks is the stress crop coefficient (equal to ESI), NDVI is the Normalized Difference Vegetation Index, T RAD is the surface radiometric temperature, T S is the soil temperature, T C is the canopy temperature, f C is the vegetation fraction, T AIR is the air temperature, T AC is the air temperature in the canopy -air space, LE is the latent heat flux, H is the sensible heat flux, Rn is the net radiation.
there is a need to downscale ET estimates to farm/plot levels integrating high and coarse resolution images. The strength of the modelling framework presented here is its suitability to cover the different layers of rangeland management, improving the understanding of ET evolution on the complex Southern African context (with mosaics of heterogeneous, partially cover landscapes and land uses).
Our hypothesis is that the Ecosystem Stress Index (ESI), applied with the new generation sensors of Sentinel 2 (S2) and 3 (S3) satellites, is a valuable tool to downscale timely estimates of actual ET derived at coarse resolution to fine scale. Therefore, the main objective of this paper was to model and validate the water used (actual evapotranspiration) by African savanna ecosystem at fine scale (meters), integrating water stress effects. To monitor savanna water fluxes in a semi-continuous way this paper integrates two different ET-estimation approaches, with different conceptual/operational capabilities and limitations. K C -FAO56 (Allen et al., 1998), integrating reflectance-based "crop" coefficients, is used to derive unstressed (potential) savanna evapotranspiration (with high spatial resolution), and the two-source surface energy balance model -TSEB (Norman et al., 1995) integrating radiometric surface temperature allows the determination of stressed (actual) ET (with low spatial resolution). By combining the two approaches the water stress across the savanna can be estimated and used to produce actual ET with high spatial resolution. The choice of the two approaches is based on their proven ability to estimate ET over partially vegetated heterogeneous landscapes (Cammalleri et al., 2010;González-Dugo et al., 2009), as well as over savanna-type ecosystems (Andreu et al., 2018a,b;Campos et al., 2013). Fig. 1 illustrates the modelling framework used in the study. Kc-FAO56 (top of Fig. 1), uses high spatial resolution reflectance data together with rainfall data to derive an empirical crop coefficient K C . The potential ET (ET p ) can then be determined by multiplying K C by reference ET (ET 0 ), which reflects the atmospheric demand and is estimated using meteorological inputs. More details on K C -FAO56 model are presented in section 2.3. For the TSEB scheme (bottom of Fig. 1) the main input is radiometric land surface temperature (LST) at coarse scale. TSEB uses LST to parameterize the energy exchange at the land surface and therefore to estimate the energy used by the actual ET (ET a ). More details on TSEB model are presented in section 2.2. The Ecosystem Stress Index (ESI) is then derived as the ratio of ET a and ET p and used as an additional parameter in the determination of K C , in order to constrain the potential ET originally derived with K C -FAO56 formulation and instead estimate to actual ET at high spatial resolution (right of Fig. 1). This procedure was tested and validated over an eddy covariance experimental site in South Africa during 2011-2012, using SPOT 4 MS & 5 MS satellite data sets (VIS/NIR), with similar spatial characteristics to S2, and AATSR-ENVISAT (TIR) with similar characteristics to S3 thermal sensor, as a preparation for an operational application of this scheme with the new Sentinels sensors, which is currently under development to monitor Kruger National Park savanna evolution. A previous analysis of TSEB was done over Skukuza during December 2015, and the whole 2016, using flux tower data series, integrating radiometric surface temperature derived from field measurements. The assumption made was that possible deviations between Penman-Monteith and Two Source models due to differences in radiation transmission and heat/momentum transfer are lower than the expected errors in satellite retrievals, and that those effects are rather linear, to assume the K Cb -ESI -ET 0 relationship to be true.

Two Source energy balance model (TSEB) -basic formulations
For regional estimations, considering the effects of a partial vegetation cover, the two-source EB model used in this study (Norman et al., 1995;Kustas and Norman 1999) is of great interest; because it separately formulates the flux energy exchange between the atmosphere and the soil, and between the atmosphere and the vegetation. The model used in this study was the updated version of the Two-Source Energy Balance (TSEB) model as described by Kustas and Norman (1999) and Li et al. (2005). The Python TSEB code from H. Nieto and R. Guzinski (https://github.com/hectornieto/pyTSEB, last accessed 22.01.2018) was implemented in this study. The model assumes that the surface radiometric temperature (T RAD ) is a combination of soil (T S ) and canopy (T C ) temperatures [K], weighted by the vegetation fraction (f C ) derived from the sensor. (1) The model also assumes surface energy-balance expressed by the equation below, where Rn [Wm −2 ] is the net radiation reaching the surface, H [Wm −2 ] the sensible heat flux used to heat the atmosphere, LE [Wm −2 ] the latent heat flux to evapotranspire the water, and G [Wm −2 ] the soil heat flux heating the first layers of the soil: It can be formulated for the entire soil-canopy-atmosphere system, or for the soil (subfix s) and canopy (subfix c) components separately. Since the radiation formulation follows the "layer-approach" (Lhomme and Chehbouni 1999), a simple summation of the soil and canopy components yields the total flux.
The soil heat flux (G) was estimated as a ratio (Є = 0.35) of the net radiation reaching the soil, as follows: Within a series resistance scheme, the sensible heat fluxes H C , H S , and H were expressed as: where T AC [K] is the air temperature in the canopy -air space, R X [s m −1 ] is the resistance to heat flow of the canopy boundary layer, R S [s m −1 ] is the resistance to heat flow in the boundary layer above the soil, and R A [s m-1] is the aerodynamic resistance calculated from the stability-corrected temperature profile equations (Brutsaert 2010), using Monin-Obhukov Similarity Theory (MOST).
Since Eq. (1) contains two unknowns (T C and T S ), the system of equations (1)-(9) is initialized assuming a potentially transpiring canopy (LE C ) following the Priestley-Taylor equation (Priestley and Taylor 1972): where α PT is the Priestley-Taylor coefficient, usually taken as 1.26 [-], fg [-] is the green vegetation fraction (1), Δ [kPa K −1 ] is the slope of the saturation vapour pressure versus temperature, and γ [kPa K −1 ] is the psychrometric constant. If the vegetation is stressed, the Priestley-Taylor approximation, i.e. Equation (10), overestimates the transpiration of the canopy and negative values of LE S are computed by the model (through equations (1)- (9)). This unlikely condensation over the soil during daytime indicates the existence of vegetation water stress, which is then solved by an iterative process that reduces α PT in Eq. (10) until Eqs. (1)-(9) yield a non-negative value for LE S . If α PT reaches 0 A. Andreu, et al. Physics and Chemistry of the Earth 112 (2019) [154][155][156][157][158][159][160][161][162][163][164] and LE S is still negative, then it is assumed that there is no ET and backup model formulations are used.

Crop coefficient (K C )-FAO56 model -basic formulation
In the FAO56 method the ratio of non-stress canopy ET (ETp) to reference ET (ET 0 ), the crop coefficient (K C ), is determined using VI-NDVI (Normalized Difference Vegetation Index) obtained from multispectral imagery (Annex 1). This provides robust and continuous estimation of the ET under unstressed conditions. Due to the high spatial resolution of the VIS/NIR information (10 m), this approach provides the resolution required for farm/local scale applications. In the dual crop coefficient method (Allen et al., 1998), K C is divided into two coefficients, one that represents the transpiration from the crop (K Cb ) and the (K S ) water stress, and another one as indicator of the soil evaporation (K e ).
K Cb was derived using Campos et al. (2013) proposal for savanna systems. K e is direct influenced by the climate conditions and the precipitation range, being higher when the evaporative demand, determined by the available energy rates, could be satisfied with the available water in the soil. It was computed following González-Dugo et al. (2013) by means of a simple water balance, where K e is derived using precipitation data and a reduction coefficient. Finally, in order to take into account the water stress, K Cb will be multiplied by a transpiration reduction factor, K S (in this case approximate by the ESI) which ranges from 1 to 0 depending on the available water in the ecosystem.

Study area
Data from two different savanna sites managed by the Centre for Scientific and Industrial Research (CSIR), Skukuza and Malopeni (CarboAfrica), were used to evaluate the vegetation and ecosystem parameters. Skukuza experimental site, initiated in 2000 as part of the Southern Africa Regional Science Initiative (SAFARI 2000) and within the FluxNet network (Majozi et al., 2016), equipped with energy flux measurement systems and meteorological stations, was used to validate the modelling framework and for input meteorological forcing for the analysis. These sites are part of the Kruger National Park (Fig. 2), in the north-eastern part of South Africa. The Kruger National Park constitutes one of the largest game reserves in Africa, covering approximately 19485 km 2 . The park dissects the Limpopo and Mpumalanga provinces.
The Skukuza flux tower (25.0197°S, 31.4969°E) lies at 365 m above sea level at the boundary of broad-leaved Combretum savanna and fineleaved Acacia savanna (Ramoelo et al., 2014). The climate of the area is semi-arid, receiving annual rainfall of about 547 mm, which falls between November and April (Scholes et al., 2001). Annual minimum average is ∼14°C and maximum average ∼29°C. The Malopeni flux tower (23.8325°S, 31.2145°E, 384 masl) is located in the northern part of the park, approximately 130 km north-west of Skukuza, and it is dominated by broad-leaf Colophospermum mopane. Rainfall of the area ranges between 0 mm during the dry season to 472 mm in the rain season. The average temperatures are similar to Skukuza, with minimums of ∼12 and maximums of ∼30°C. The areas are classified as hot and dry savanna (Ramoelo et al., 2014), with annual variability and irregularity of rainfall. These contrasting savanna occur on soils with different texture, water holding capacity and nutrient levels, as well as different physical structure, physiology and phenology (Scholes et al., 2001). They are homogeneous landscapes with smooth topography and low slopes. Both areas are in a natural state and act as wildlife reserves, with the main land use being tourism (natural safaris). Elephant, buffalo, impala, kudu, wildebeest, rhinoceros, among other wild-animals, constitute the grazers and browsers of the park. Kruger National Park is surrounded by both communal and urban areas, as well as private conservancies.
All the energy balance components, Rn, G, H and LE, were measured directly over the field in Skukuza, at 16 m above ground level with an eddy covariance tower (ECT). Wind speed was measured with a 3D sonic anemometer (CSAT3 3-D sonic anemometer, Campbell Scientific Inc, USA), and water vapour and carbon dioxide concentrations were measured with an open-path infrared gas analyser (IRGA; model LI-7500A, Li-Cor, Lincoln, NE, USA). The anemometer was oriented in the prevailing wind direction (W). Air temperature and humidity were measured independently at 15 m agl using a HMP50 probe (Campbell Scientific Inc., Logan, UT, USA). A NRLite Kipp and Zonen net radiometer was installed (z = 15 m agl) from 2009, and incoming and outgoing shortwave and longwave radiation was measured with a Kipp and Zonen CNR4 Net Radiometer (z = 15 m agl) since 2015. Post-processing of eddy covariance datasets from 2010 to 2012 and 2015 to 2017 were done using two different eddy-covariance software packages, which include the necessary corrections. Dataset from 2010 to 2012 was processed using EddySoft developed by Max Planck Institute for Biogeochemistry, Germany (Kolle and Rebmann 2007) and dataset 2015 to 2017 was processed using EddyPro developed by LI-COR Bioscience, USA (Burba and Anderson 2010). After processing, halfhourly average values of the turbulent energy balance components were obtained. Soil heat flux was determined by averaging the measurements of three ground heat flux plates (HFT3, Campbell Scientific) installed over the broad leaved combretum savanna area. Meteorological stations were also located inside Skukuza, which measured variables of interest such as precipitation (Texas TR525M tipping bucket rain gauge).
Two field campaigns were made over Skukuza and Malopeni, together with the CSIR, to collect ecosystem vegetation parameters. The first one was conducted over the dry period, from the 22nd to the 28th of June 2016, and the second over the wet season, from the 22nd of January 2016 to the 1st of February 2017. Ecosystem and local leaf area index (LAI) measurements were made in the vicinity of the tower using the LICOR 2200 Plant Canopy Analyser, following the distribution of the ecosystem, and taking isolated measurements of different vegetation. For the ecosystem LAI, a total of 120 sample points were measured during each season. Measurements of trees and grass spectral responses were made in the field during each season using a portable system, to study the seasonal variability of the vegetation. Reflectance measurements were made using the ASD FieldSpec 3 spectroradiometer (ASDInc.), which registers radiance data in the range 350-2500 nm, a reference panel calibrated in a laboratory was used (Spectralon, Labsphere, North Sutton, NH) to measure incident shortwave radiation for calibration purposes. Around 100 samples with a bare fiber (FOV = 25°) were taken over different mature and young trees and scrubs, at 0.5 m height from the leaves. The understory was more variable, with different species, canopy heights, fractional cover and green fractions. Fifty measurements were taken at 1 m height above the soil (grass), with a visible area diameter equal to 0.93 m. Significant parameters for the description of canopy structure were determined in the field or collected from the literature, such as tree leaf size (s = 0.02 m), and canopy height -h C as a constant average tree height A. Andreu, et al. Physics and Chemistry of the Earth 112 (2019) 154-164 of 6.5 m (Majozi et al., 2016). Local LAI from grass and oak and ecosystem LAI were derived from field reflectance measurements in order to determine whether the derivation of LAI using the broad bands from satellites could be used as a proxy, following Choudhury (1987). Data from December 2015, and 2016 of Skukuza flux tower was used to study the behaviour of TSEB, due to the possibility to derive thermal information from the 4-way radiometer located over the site. Data from 2010 to 2012 of Skukuza ECT was used to validate the modelling framework using remote sensed image series (SPOT 4/5 and AATSR). The energy balance closure (EBC) error (Twine et al., 2000;Wilson et al., 2002;Foken, 2008;Franssen et al., 2010) was used to describe the accuracy of the data-flux series collected. To determine the area that contributed most to the measured fluxes at the tower, and assure sufficient fetch for remote sensing integration, an approximate solution for the contribution to the vertical flux (Schuepp et al., 1990) was computed. Field data collected over Skukuza and Malopeni on 2016/2017 were used to study the ecosystem vegetation parameters, the evolution and structure of the canopy, and the annual cycle of vegetation phenology (leaf area index and fractional cover).

Earth Observation data and thermal field data
The dates used in the analysis have been selected from the data series from Skukuza ECT system and the remote sensed data, discarding days according to the following criteria: (a) lack of satellite thermal/ optical imagery over the ECT site due to clouds or temporal sensor resolution and (b) periods with gaps due to instrument failure. 10 images of SPOT 4/5 from Skukuza were analyzed (Table 1) to compute K Cb , but 2 of them had to be excluded due to cloud coverage. 200 days were used to evaluate TSEB with AATSR-ENVISAT images from 2010 to 2012, using radiometric surface temperature from the sensor over Skukuza and NDVI derived from the satellite.

SPOT 4 MS & 5 MS as a VIS/NIR source for K C -FAO56
The sun-synchronous SPOT 4 (Satellite pour l' Observation de la Terre) launched in March 1998 and SPOT 5 HRG multispectral sensor data was used in this study. SPOT 4/5 images were acquired between April 2010 and May 2012, from Zone 36 (Skukuza coordinates). SPOT mission was specifically designed to monitor land surface parameters such as vegetation productivity at high spatial resolution (m). SPOT 4 and SPOT 5 have multiple spectral bands including red and near infrared, and were collecting data from an orbital altitude of 822 km. Based on these spectral bands, f C and LAI were derived after the method by Choudhury et al. (1994), described in Annex 1. Both satellites sensors have a pixel size of 10 m, making them suitable for local scale monitoring of critical environmental processes. Atmospheric and surface emissivity effects were corrected using an atmospheric radiative transfer model MODTRAN4 (Berk et al., 1999), from the FLASSH ENVI module.

AATSR-ENVISAT as a T RAD source for TSEB
To take into account the savanna ecosystem water stress, we applied TSEB with thermal information from the Advanced Along-Track Scanning Radiometer (AATSR), acquiring images between 2010 and 2012 over Skukuza. AATSR was one of the instruments on board the European Space Agency (ESA) satellite ENVISAT. The instrument included an along track scanning technique, continuous on-board calibration of the thermal bands, low-noise, infrared detectors and cooled to near-optimum operating temperatures. All these technical advancements makes the sensor suitable for the provision of accurate thermal measurements of the earth's surface. Global land surface temperature has been produced operationally at 1 km resolution since 2004, with an algorithm accuracy within the target specification (Prata, 2002;Coll et al., 2009Coll et al., , 2012. The product used in the study was the UOL LST product, which contains full resolution nadir-view land surface temperature generated by the University of Leicester (UOL), UK. The same LST retrieval method as for the Gridded Surface Temperature (Prata, 2002) products is used, although improved auxiliary datasets for land cover, green vegetation fraction and total column water vapour are applied.   Choudhury et al. (1994) model to produce LAI for the application of TSEB during 2015-2016. MOD13Q1 MODIS product, with 250 m of resolution and provided every sixteen days, was regarded as accurate enough for our requirements and the vegetation evolution (Andreu et al., 2018a,b). This product selects an NDVI value representative of the 16-day period, as the average of the two days with maximum NDVI and highest-quality information.

T RAD derived from field data -TSEB analysis during 2015/2016
In situ land surface temperature was derived from converting the upwelling (LW OUT ) and downwelling (LW DN ) longwave radiation measurements on the field using the following equation: where σ is the Stefan-Boltzmann constant and ε is the surface emissivity (derived from standard previously published values on Campbell and Norman, 2009).

Study area analysis
The energy balance closure was evaluated over Skukuza experimental site, yielding average closures from 2011 to 2016 of 17% (Fig. 3), within the error range found by other authors (Foken, 2008;Franssen et al., 2010;Majozi et al., 2016). Generally, the 4-ways radiometer measurements provide better quality and precision, regarding the energy balance closure, but the sensor was installed in 2015.
LE presents a peak during the wet season (from November to April), and a minimum during the dry season (from May to October), related to the typical rainfall temporal distribution. However, the magnitude and timing of precipitations can vary significantly from year to year. Incoming shortwave radiation does not have high variations (∼250 Wm -2 ), with maximums from November to February and minimums over the dry season (April to September), also depending on the daily cloud coverage. Due to the relative smaller magnitude of G compared to the other components, this flux is usually simplified and discarded on daily scales. However, in savanna type ecosystems, G may represent more than 20% of the net radiation during the dry period, similar to LE (Andreu et al., 2018a,b), due to low fractional covers. These fluxes are consistent with Ramoelo et al. (2014) and Majozi et al. (2016) and the climatic conditions of the area, with high temperatures and no precipitation during the dry period.
A cumulative normalized contribution to flux measurement curve (footprint) was estimated as a function of the distance from the measurement point, as well as the relative contribution to the measured values. The area that most contributed to the energy fluxes measured at Skukuza ECT using the approximate solution of Schuepp et al. (1990) and data from 2015 to 2016 was located between 0 and 500 m in the predominant wind direction, where the 60% of the fluxes captured were originated, reaching the 80% within 1000 m, with the maximum located in the first 200 m (Fig. 3).
During 2016, 2017, ecosystem LAI measured over Skukuza during dry season (May) was equal to 0.3, with a fractional cover of 0.4. Similar to the tree cover computed by Scholes et al., 2001) (0.3), due to the nonexistence of grass during that period. With the ASD LAI from the ecosystem reached 0.5. During wet season (February) ecosystem LAI, integrating tree and grass cover, reached 2.65 (ecosystem fractional cover of 0.9). Malopeni LAI during dry season was higher, 0.63, due to the different tree vegetation and fractional cover. During the wet season was equal to 1.21, probably due to the lower tree and grass cover and leaf structure. Fractional cover was also estimated using SPOT 4 and 5 images during 2010-2012. Usually values found were higher than the ones computed on the field, with maximums of 0.6-0.86 ranging from December to May and minimum during the period from May to October (0.1-0.5), following a similar trend. Fig. 4 compares estimated and observed half-hourly fluxes applying TSEB for all days with radiometric surface temperature derived from the tower longwave outgoing radiation measurements. As model input (air temperature and humidity, solar incoming radiation, wind speed, and T RAD ) and validation data (four surface energy fluxes) the data-set collected over 2015-2016 by the ECT over Skukuza was used. Enegy balance closure was solved using the Bowen Ratio. Root Mean Square Differences (RMSD) between estimated and observed values were 24 Wm -2 for Rn measured with the 4-ways radiometer and 114 Wm -2 when compared with the net radiometer, 28 Wm -2 for G, 82 Wm -2 for H, and 43 Wm -2 for LE. Sensible heat flux error is high compared with other authors results, although LE error is within the limits found for more uniform and homogeneous canopies (Norman et al., 1995;Kustas and Norman, 1999;Timmermans et al., 2007;Sánchez et al., 2008;  Modeled LE tend to overestimate measured LE for the whole season. RMSD for the dry season are lower for all fluxes, with less error dispersion. Since H is the main contributor to fluxes, simpler -non thermal based -models may greatly overestimate LE. A better characterization of the ecosystem aerodynamic components and vegetation structure (e.g. green fraction) can improve the final estimations for this systems (Andreu et al., 2018a,b;Kustas et al., 2016).

Evaluation of the performance of TSEB model using T RAD and NDVI from AATSR 2010-2012
The TSEB model was applied and evaluated over Skukuza with data from AATSR satellite. As model input (air temperature and humidity, solar incoming radiation and wind speed) and validation data (four surface energy fluxes) we used the data-set collected over the same period by the eddy covariance tower. RMSD for the net radiation was found equal to 80 Wm -2 , estimated Rn higher than measured values for almost all points (Fig. 5a). The great error of Rn was probably due to the field measurements, since during the period just the net radiometer was available, and, as it is possible to see on Fig. 5b, available energy over the area was always lower than the turbulent fluxes measurements. Measurements from the net radiometer and the 4-ways radiometer were compared during 2015/2016, and a calibration curve was develop and applied to the Rn measurements from 2010 to 2012, and the resulting Rn values were used for the model validation (Fig. 6). Fig. 6 compares estimated and observed fluxes applying TSEB with data from AATSR satellite. The following RMSD values for the energy fluxes were found, closing the balance by the Bowen ratio, RMSD LE = 53 Wm -2 , RMSD H = 50 Wm -2 , RMSD Rn = 60 Wm -2 , RMSD G = 21 Wm -2 . RMSD for the turbulent fluxes are within the range found by other authors (Norman et al., 1995;Kustas and Norman, 1999;Timmermans et al., 2007;Sánchez et al., 2008;Gonzalez-Dugo et al., 2009), consistent with typical uncertainties derived for the flux measurement system, and average balance closures error of the area for that days, representing 70 Wm -2 . As in the previous section, the model tend to overestimate LE, as well as the net radiation. RMSD for the dry season is lower for all fluxes. Cammalleri et al. (2010) found in a similar semi-arid ecosystem, RMSD values for Rn of 28 Wm -2 , for G of 17 Wm -2 , and 40 and 43 Wm -2 for H and LE, respectively. Morillas et al. (2013) found in a more arid steppe environment RMSD values for Rn of 58 Wm -2 , 64 Wm -2 for H and 105 for LE Wm −2 , with canopy and soil radiometric temperatures ground data. In Kustas et al. (2016), a revision of Morillas et al. (2013) data analysis was done, confirming the results but significantly reducing the bias using both key vegetation inputs to TSEB and providing a more realistic parametrization to Rs for rough soils.

ET O versus actual from TSEB and derivation of ecosystem stress index (ESI)
Using AATSR ET estimated values and potential ET derived from meteorological information and reflectance values, monthly mean ETa and ET p values were derived and plotted in Fig. 7, together with monthly mean accumulate precipitation. Although days with data points are linearly connected in Fig. 7, to illustrate the evolution of monthly mean values, the interpolation between dates is not linear due to sharp changes related to weather conditions (cloud cover, vapour pressure, air temperature and wind) and the averaging process. The transpiration rates are reduced during the dry season due to lower atmospheric demand and water deficiency. Higher atmospheric demand is observed during the wet season, where solar incoming radiation is also higher. During these periods rainfall rates reach their peak, so the water demand will be mostly/partially covered by the soil water content. In this situation, the natural vegetation will be slightly stressed. From May to October, when the atmospheric demand is still high but the water availability does not allow the canopy to transpire potentially, the vegetation would be water stressed, with a dormant period where the demand is reduced. Monthly mean ESI index is derived from the mean potential and actual ET values, and this index was later used as a K S proxy to integrate the water stress of the vegetation into the K Cb model.

ET at high resolution (K C -FAO56 + ESI)
For this first attempt to evaluate evapotranspiration over this savanna ecosystem, 8 days from 2010 to 2011 were selected and analyzed over Skukuza, depending on the cloud cover. K Cb was derived using NDVI from SPOT images, K e was derived doing a simple water balance over the soil root-zone, and K S was approximated as monthly mean ESI computed previously. Differences between measured and estimated values (Table 2), RMSD 0.6 mm day −1 , were consistent with other studies using K C -FAO56 and TSEB for similar ecosystems (Andreu et al., 2018a,b, Campos et al., 2013. It can be observed that the estimated ET results using K C -FAO56 + ESI were slightly lower than the measured results for the entire period under study, while using the usual K C -FAO56 ET is usually overestimated, except when the values are close to A. Andreu, et al. Physics and Chemistry of the Earth 112 (2019) 154-164 zero. For the months of September and October both models produced ET estimates of 0 mm day −1 whereas measured ET ranged between 0.3 and 0.5 mm day −1 . However, the relative error of this estimations is lower due to the low relative proportions of this months ET over the period. As it is possible to see, using a stress crop coefficient restrain the ET estimations, lowering the range, modelling the natural partially covered ecosystem in a more realistic way. Distributed real ET over Skukuza was mapped as a first approach to monitor the ecosystem status on a regular basis, with the objective of assessing the applicability of S2 and S3 satellites. It was possible to derive reliable ET values that reflect the local conditions and climate, and the evolution of the vegetation cover (Fig. 8). During dry season ET rates ranged in between 0 and 1.5 mm day −1 , and during wet season can reach 6 mm day −1 . Distributed ET variations can be seen over the studied area due to the high spatial resolution of the satellite images. As an example, near the river, were the water availability to the root zone is maintained throughout the year, vegetation is always green and transpiring. Some areas are perceptible dryer during all the studied period (upper left zone from the image), and this information can be used in rangeland management, keeping the grazing animals outside those vulnerable zones.

Conclusions
Mapping actual ET (integrating vegetation water stress) over savanna on a plot scale was possible using high spatial resolution reflectance data, such as SPOT, and radiometric surface temperature, derived from AATSR. Combined TSEB and K C -FAO56 models provided an accurate representation of the heterogeneity over natural landcover, and the local meteorological conditions. ESI values and the associated actual ET values were derived at medium spatial resolution (km) using TSEB, and actual ET at high spatial resolution was estimated using K C -FAO56 model. ET estimates reflected the vegetation water stress conditions caused by the local climate, whereas the evolution trend of the vegetation was perceptible.
The differences between estimated and observed surface fluxes derived from TSEB were of similar magnitude as the uncertainties derived from the Skukuza ECT measurement system. Therefore, the estimated fluxes were assessed as being sufficiently accurate to be employed in a distributed way and on a more regular basis. The approach of ESI downscaling proved to be useful for natural semi-arid vegetation. Thermal-based models are well suited to detect short-term water stress, frequent in these semi-arid areas. Errors between estimated and observed daily evapotranspiration were consistent with the errors found by other studies over savanna ecosystems (Campos et al., 2013).
An important source of error over the area for the period studied (before 2015), responsible for the discrepancies between net radiation observed and estimated fluxes, are the field ground data from the net radiometer, although the calibration or Rn net radiation measurements improve the final results. The mismatch between measured and estimated values footprint can also be an important source of error, as can be the soil heat flux, which is difficult to measure in the field due to land heterogeneity. G influences the system available energy, limiting the amount available for the latent and sensible heat fluxes. A better characterization of the ecosystem aerodynamic components and vegetation structure (leaf are index, clumping factor and green fraction) can improve the final estimations for these heterogeneous, partial cover systems.
Sentinel 2 & Sentinel 3 (also MODIS with 1 km spatial resolution,  and Landsat with higher thermal spatial resolution) images will be used for the operational application of the modelling framework, to monitor Kruger National Park savanna evolution. Distributed maps of meteorological variables will be used as input, to integrate the heterogeneity of the region where the differences in orographic, meteorological, abiotic and biotic conditions create different subtypes (bioclimatic levels). Further publications exploring this operational application are planned (e.g. , to determine the possible relevant times and spatial scale resolutions needed for each management objective, and the implications of our findings for savanna water resources management. The precision/resolution/accuracy of the information required for the Park management will differ at each scale: farm-local (e.g. evaluating the effect of management practices, livestock and wild animals densities, grazing intensity), to watershed (e.g. detection of vulnerable areas) and regional (e.g. early prediction of drought). The strength of the modelling framework proposed is its suitability to cover rangeland management from local to broader scales. To characterize local effects and practices, the 10 m spatial resolution determined by Sentinel 2 will probably provide the detail required by the authorities (e.g. SANParks), and water managers, as some surveys suggested (Nieto et al., 2018). For larger scales, lower temporal and spatial resolution will suffice in order to analyze the evolution of vulnerable stressed savanna areas. However, this methodology is highly hindered by clouds, and gaps will cause losses of key information especially during the wet season. The integration of other data sources (through modelling) not sensitive to clouds (e.g. radar) can be a possibility to overcome this limitation.

Declarations of interest
None.

Funding
This work was supported by the European Space Agency, under the TIGER Initiative (Tiger Bridge project #410). entire archive of AATSR data products is available free of charge from ESA. The authors do not have permission to share SPOT 4/5 data. This work partly used eddy covariance data acquired and shared by the FLUXNET community, including these networks: CarboAfrica. The ERA-Interim reanalysis data are provided by ECMWF and processed by LSCE. The FLUXNET eddy covariance data processing and harmonization was carried out by the European Fluxes Database Cluster, AmeriFlux Management Project, and Fluxdata project of FLUXNET, with the support of CDIAC and ICOS Ecosystem Thematic Center, and the OzFlux, ChinaFlux and AsiaFlux offices.

Annex 1. NDVI -Normalized Difference Vegetation Index
Near infra-red (NIR) and Red bands are used to retrieve Normalized Difference Vegetation Index, following the basic equation: The scaled NDVI approach (Choudhury et al., 1994) method is used to retrieve f C (0) with NDVI remote data: where NDVI MAX and NDVI MIN , represent a surface fully covered by vegetation (∼0.9) and completely bare (∼0.08), respectively. The parameter p A. Andreu, et al. Physics and Chemistry of the Earth 112 (2019) 154-164 represents the ratio of a leaf angle distribution term (k) to canopy extinction (k'), where p = k/k'. k is the leaf angle distribution function, which appears to range between 0.5 and 0.7 (Ross, 1975) depending on the leaves having a spherical (k = 0.5), vertical (k < 0.5) or horizontal (k > 0.5) distribution. We assume a random distribution as our savanna ecosystem contains mixture of different plant types. k' is the damping coefficient, ranging between 0.8 and 1.3 for the NDVI (Baret and Guyot 1991;Delegido et al., 2011). We used a weighted p parameter of ∼0.9, with k ∼0.5 and k' ∼0.55. Leaf area index is derived from fractional cover by an exponential function as (Choudhury 1987) suggested: