Long-Term Grass Biomass Estimation of Pastures from Satellite Data

The general consensus on future climate projections poses new and increased concerns about climate change and its impacts. Droughts are primarily worrying, since they contribute to altering the composition, distribution, and abundance of species. Grasslands, for example, are the primary source for grazing mammals and modifications in climate determine variation in the available yields for cattle. To support the agriculture sector, international organizations such as the Food and Agriculture Organization (FAO) of the United Nations are promoting the development of dedicated monitoring initiatives, with particular attention for undeveloped and disadvantaged countries. The temporal scale is very important in this context, where long time series of data are required to compute consistent analyses. In this research, we discuss the results regarding long-term grass biomass estimation in an extended African region. The results are obtained by means of a procedure that is mostly automatic and replicable in other contexts. Zambia has been identified as a significant test area due to its vulnerability to the adverse impacts of climate change as a result of its geographic location, socioeconomic stresses, and low adaptive capacity. In fact, analysis and estimations were performed over a long time window (21 years) to identify correlations with climate variables, such as precipitation, to clarify sensitivity to climate change and possible effects already in place. From the analysis, decline in both grass quality and quantity was not currently evident in the study area. However, pastures in the considered area were found to be vulnerable to changing climate and, in particular, to the water shortages accompanying drought periods.


Introduction
According to the Intergovernmental Panel on Climate Change Fifth Assessment Report [1], the mean annual temperature is likely to exceed 2 • C over Africa by the end of this century. Heat waves are likely to occur with a higher frequency and longer duration and existing stress on water availability will, moreover, be amplified. Thus, food insecurity brought about by reduced crop productivity and livelihood is identified among the risks for African region. The agricultural economy of developing countries principally relies on crop production but, also, the livestock sector constitutes a very important component giving a contribution that goes beyond food production (e.g., beef cattle, dairy cattle). In received increased attention in the last decade due to its relevancy to climate change and mitigation programs such as REDD+ [17].
Given the usefulness of remote sensing in estimating forest biomass has largely been demonstrated and algorithms have proven to be trustworthy, efficient, and advantageous, few studies have focused on grasslands and also their dynamics in relation to climatic evolution [28][29][30][31][32][33][34]. Even fewer studies have dealt with the lack of information and reference data that characterize remote regions [28] such as Africa. In [35] Schucknecht et al. developed a predictive model based on the spatiotemporal variability of remotely sensed biomass of pastures in Niger for the period 1998-2012 but correlation with climatic variables was not investigated. In [36], Paudel & Andersen aimed at assessing rangeland degradation from precipitation variability in Nepal using Landsat MSS, TM, ETM, and SPOT images covering the years 1976-2008, but without estimating grass yields. Mutanga and Rugege [37] modeled savannah biomass using MODIS data over the Kruger National Park in South Africa, however, no multitemporal analysis on the study area was performed. Research on the topic is still open, and there is a general need for methodologies that could perform a more comprehensive, automatic, and repeatable analysis over longer time windows and variable spatial scales, and that could deal with a lack of reference data.
Biomass estimation of grassland often relies on traditional methods based on field measurements through the harvest method and the pasture meter method [38]. However, they are generally costly, time consuming, imply the destruction of the samples, and can be performed over very limited extents. In this context, satellite remote sensing may constitute an important turning point where ground measurements, if available, can be used for the validation of remotely sensed data. A frequently used method for the estimation of vegetation productivity through remote sensing over large areas is the light use efficiency (LUE)-based model [39], that has been also applied on grasslands [40,41]. Besides this, vegetation indices (VIs) have proven to be effective for grass biomass estimation in terms of both total yield and productivity [29,[42][43][44]. In particular, the most commonly adopted VIs are the normalized difference vegetation index (NDVI) [45,46], enhanced vegetation index (EVI) [47], and the land surface water index (LSWI) [43]. Evidence found in [48] suggests that NDVI is more sensitive to low-density vegetation cover, while EVI is more sensitive to high-density cover, which also means that EVI tends to have lower range of variability on semiarid sites. Moreover, several studies identify the normalized difference vegetation index as the best predictor of grassland ecosystem attributes such as biomass [29,30,49], even if possible poor performances could be observed at high canopy density where the adoption of vegetation indices based on wavelengths located in the red edge is suggested [50]. Observations of biomass by satellite remote sensing VI, such as NDVI, can be repeated in time and in different moments of the phenological season. This approach provides temporal information on growth conditions such as growth rates and maximum/minimum yields, as well as plant responses to dynamic weather conditions [51].
Until recently, Landsat and MODIS have been reference missions for environmental studies focused on vegetation monitoring [5,10]. At present, the European Union (EU) Copernicus program and the Sentinel missions represent a breakthrough due to their ability to provide free products from different sensors with valuable performances in terms of spatial and spectral resolution and revisit time. In the near future, the European Space Agency (ESA) Earth Explorer BIOMASS satellite mission will give a new boost to the topic, providing global estimates of carbon stocked in forests [20].
In this work, we present the results regarding grass biomass estimation provided by an end-to-end automatic EO data processing approach exploiting 20-year time series of historical data in significant regions of Africa where food security is a critical problem. As is well known, the poor and incomplete production of thematic maps and environmental data in some areas limits the development and application of programmatic strategies for the control and management of resources. Furthermore, in the subject matter of this study, scientific data are limited in their temporal and spatial distribution. It has to be noted that given the long time period and the wide area considered, the study tackles and provides solutions to various issues usually neglected when remote sensing is applied to a restricted area and in a limited time span. For example, these involve the combination of various types of images with different spatial resolutions or a high level of automation of the processing steps. In particular, we first developed an automatic land cover classification algorithm for pasture identification by means of neural network techniques. Second, regression analysis was carried out to evaluate pasture biomass from a normalized vegetation index. These two pieces of information can be critical inputs for models estimating future productivity and represent valuable analysis in areas where no previous studies have been carried out.
Zambia has been selected as a significant test area since it is found vulnerable to the adverse impacts of climate change [52]. In Zambia, approximately 84% of the cattle population is maintained by the traditional sector of farming, which is principally diffused in rural areas [53]. According to [54], cattle production is of primary importance, but productivity is low. Moreover, it is reported in [1] that all countries within the Zambezi River Basin could struggle with water shortages even influenced by non-climate factors (e.g., economic growth, population variation). An increase in temperatures and potential variation of rain distribution and water scarcity open up new challenges for livestock and agricultural sectors in the region. Meanwhile, livestock demand will continue to increase significantly in the coming decades [55] along with the increase in global population which is expected to reach 1 billion for the African continent. Finally, climate change is expected to make production more irregular [56][57][58], especially in underdeveloped countries where resources for adaptation are scarce.
Regarding the grass types considered in this study, it is important to remember that several definitions of grasslands exist. For example, to define grasslands, UNESCO [59] counts land areas covered by grass with less than 10 percent of trees and shrubs. Differently, in [60], White et al. include shrubs. Often, in fact, grasslands are mixtures of different plants where shrubs and grasses together can produce several combinations. In this research, grasslands are considered as those areas mainly constituted by grass and characterized by high chlorophyll content during the vegetative stage, suitable for grazing, and able to support cattle for great part of the year. Moreover, in this research study, the terminology for grasslands follows the work of Allen et al. [22], which was published to ensure correct communication during exchange on this topic.

Zambia Strategies in Livestock Sector
Zambia, a republic located in south-central Africa (Figure 1), is found to be particularly vulnerable to climate change, and negative effects of climate change are already challenging its agriculture system [61]. In particular, droughts in 1993-1995, 2001-2002, and 2004-2005 had large impacts on yields and food security [62].
Zambia is currently developing a programmatic approach involving national climate-resilient development strategies. This encompasses the National Climate Change Response Strategy (NCCRS) [63] and the National Adaptation Programme of Action (NAPA) [64]. It also includes the National Adaptation Plan (NAP) and the Intended Nationally Determined Contribution (INDC) which is part of the global climate change agreement signed at the 21th Conference of the Parties (COP-21) in Paris in 2015 [65]. In particular, the INDC focuses its priority actions on building resilience in the agriculture sectors. The Zambian government also facilitated the formulation of the National Policy on Climate Change [52], which was adopted in 2016. Within these agreements, Zambia identifies priority actions as adaptation measures to climate change, such as the diversification and promotion of climate-smart agriculture also in livestock sector, including the improvement of the livestock feed and rangeland management, to which FAO's initiatives relate. The livestock sector is identified by the government as a critical sector for livelihood production, source of revenue, and diversification export strategies. According to Simbaya in [66], more than 35% of the total agriculture output comes from the livestock sector. However, Zambia's agriculture contributes only between 11% and 16% to national GDP (gross domestic product). In particular, the low performances are ascribed to the use of just a small portion of potential arable land and the small-scale farming systems adopted. In fact, the Zambian livestock sector includes both commercial and traditional systems. In particular, in the traditional system, local breeds of cattle and crop farms are generally integrated over 5 hectares of land [67], and 90% of the overall feed resources are constituted by grazing and browsing on natural grasslands mainly constituted by annual and perennial species of grasses [68]. Nevertheless, an increase in small-scale farmers is recognized as the possible way to pursue an increase in the country's economy [66].

Zambia Grazing Areas
Zambia, lying between latitudes 8 • and 18 • S and longitudes 22 • and 34 • E, consists mostly of a highland plateau having an area of about 752,614 km 2 with elevations generally ranging from 915 to 1520 m above sea level. The highest land is located in the eastern part of the country at the Mafinga Hills, where the Mafinga Central is the highest point, reaching 2339 m above sea level. The lowest point of the country is 329 m a.s.l. at the Zambezi River. Two major river basins drain Zambia territory: the Congo basin in the north and the Zambezi Basin, which is spread over the west, center, and south provinces covering about three-quarters of the whole country. Zambia is divided into three main regions, which differ in their agroecological profiles and rainfall characteristics. The first region is the warmest and driest (annual average rainfall of approximately 500 mm), with soils having low nutrients, low organic content, and high pH levels. It comprises the southwestern part of the country. The second region is characterized by a humid subtropical climate with the most fertile soils, and annual precipitation ranging between 800 and 1000 mm. It includes most of the country. Within the second region, the Zambezi floodplain and the Kalahari Sands in the west constitute a subgroup having tropical savannah climate and less fertile soils. The third region is located in the northern part of the country and receives the highest rainfall (>1000 mm/year). Woodlands occur in the high rainfall region in the form of tall and dense forests in the north and northwest, while open forests and scattered trees occur in the drier south. Grasslands range from pure grasslands to grassland with scattered trees [69], and savanna-type-grasslands dominate the low rainfall regions. According to the pilot analysis on global ecosystems by the World Resources Institute [60], Zambia is among the top countries for grassland area with 526,843 km 2 , about 70% of the total country surface (754,676 km 2 ).
Floodplains are generally among the best areas for cattle raising in the context of traditional farming. According to Kulich and Kaluba [70] and among the natural grazing areas, the Kafue and Zambezi floodplains, respectively in the South and West Province, can support large herds of cattle even in the dry season. West Province, in particular, is known as one of the most productive areas of the country for raising cattle due to the Zambezi River and vast floodplains [71,72]. According to Aregheore [69], the potential pastureland in the province is estimated at 10 million hectares and over half a million of head are currently pastured by Lozi population [73], whose settlements are found on the edges of the floodplains. The Barotse Floodplain is known to be one of the most productive cattle areas in the country [72]. There, the Zambezi forms large flat floodplains of grassland at around 1000 m elevation and drained by the river itself and its tributaries. It is located within the Kalahari Sands area (locally known as Barotse Sands), constituted by sands of aeolian origin having very low clay content, 3%-12% silt + clay [74], and extremely dry for great parts of the year [72]. Only in the wet season are the sands waterlogged but are generally nutrient-poor terrains, not very fertile and disturbed by cultivation. In addition to its significance as a grazing area, it is also an important area for its naturalistic value: in the framework of the global Multilateral Environmental Agreements (MEAs), the Barotse Floodplain was designated Ramsar site no. 1662 in 2007, regarded as being a wetland site of significant value. It is also a game management area (GMA), being an area surrounding national parks. It serves as a buffer zone where both animals and humans are supposed to co-exist.
In the Barotseland, floodplains are suitable for supporting cattle productivity and for grazing when floods recede: domestic animals are kept in the floodplains during winter season while, before the floods, they are kept in the adjacent uplands for flood-time grazing. However, especially at the end of the dry winter season, green grasses are substituted by lower quality plants [66].

Study Area
A 10,000 km 2 wide study area was designed within the Barotse Floodplain ( Figure 1). The area, referred to further in this paper as LUE, lies between the Lukulu District and the Mongu District, and covers the Luena flats comprising the Luena River, one of the three tributaries of the Zambezi in the floodplain [74]. The study area falls within the second agroecological zone of Zambia, and in particular, in the subgroup of the Kalahari Sands, characterized by less fertile soils. LUE area was chosen as it falls within the Barotse Floodplain which, as clarified before, is a grassland site of national and supranational importance both as rangeland and a biodiversity hotspot. Moreover, LUE lies on the Luena flats, a major grazing area of the Lukulu and the Mongu Districts as well as the center of life for the Lozi people, who still practice the traditional system of farming.

Climatic Data
The climate in Zambia is characterized by three seasons: a cool and dry season (April to August) with a daily mean of 17 • C or less during the coldest month of July, a hot dry season (August to November) registering mean daily temperatures in the range of 21-27 • C on the plateau and 36 • C or more in the Luangwa and lower Zambezi valleys, and a third warm-wet season (November to April). The latitudinal line at 14 • S roughly divides the country into the high rainfall region to the north and the low rainfall region in the south. Rain is typically confined to the wet season, except for very rare showers in August [69]. The Climatic Research Unit (CRU) dataset [75] was adopted as the reference source for climatic data over the LUE. The principal sources used for the routine updating of the CRU archives are monthly observations at meteorological stations across the world's land areas which come through the World Meteorological Organization (WMO) in collaboration with the US National Oceanographic and Atmospheric Administration (NOAA, via its National Climatic Data Center, NCDC). In particular, CRU Timeseries (TS) v. 4.02 database covers LUE study area in the period 1901-2017 at 0.5 degrees of resolution, and contains monthly estimates of temperature, precipitation, and other climate variables. According to CRU data, LUE experiences mean annual rainfall from 800 to 1100 mm ( Figure 2). To compute yearly averages and seasonal averages, the months of November and December of each year were included in the estimation of the next year as they are part of its wet season. In this way, it is possible to better also understand eventual anomalies in relation to vegetation growth. The rainfall trends relative to each season are shown in Figure 3. The dry seasons have particularly low precipitation levels (close to 1 mm/month for the cool period from April to August, and approximately between 5 and 20 mm/month from August to November), differing from precipitation in the wet season which is abundant and approximately equal to 160 mm/month. As for precipitation, the temperature also changes according to the three seasons ( Figure 4) but with less variation among them: the cool season from April to the mid of August has a mean monthly temperature of 20 • C, while the warm seasons, from late August to April, register approximately 25 • C. By looking at the whole data, the mean maximum temperature is around 27 • C and the mean minimum temperature is between 17 and 20 • C. Minimums and maximums over the time window of interest are shown in Figure 5. Even varying within the time span of the analysis, no significant trends are visible in the data for both precipitation and temperatures. However, anomalous variations are present. In particular, the year 2005 is characterized by particularly low rainfall ( Figure 3) and slightly higher temperatures (Figure 4) in the wet season with respect to the expected means for the period. According to Kaczan et al. [62] in 2004-2005, severe drought occurred in the area, with large impacts on vegetation and productivity. Similar higher temperatures can be also recognized in the wet season of 2016, but with simultaneously high precipitation levels and no drought phenomena.

Land Cover Classification Dataset
Medium-resolution products were considered for land cover classification purposes, to guarantee the highest performances in terms of scalability and processing time. Moreover, the analysis was intended to be extensible to a decadal time span. Due to these considerations, Landsat acquisitions were chosen. In fact, as multispectral missions, Landsat offers continuity of products over 46 years of global land observation, after its first time in 1972. It is also important to mention that Landsat products are free of charge and easy to collect. Landsat scenes are held in the United States Geological Survey USGS archive and can be searched using the Earth Explorer engine. Landsat Collection 1 Level 1 consists of data products generated from Landsat 8 Operational Land Imager (OLI)/Thermal Infrared Sensor (TIRS), Landsat 7 Enhanced Thematic Mapper Plus (ETM+), Landsat 4-5 Thematic Mapper (TM), and Landsat 1-5 Multispectral Scanner (MSS) instruments. In particular, data assigned to "Tier 1" (T1) inventory have the highest available data quality and are considered suitable for time series analysis. This includes Level-1 Precision and Terrain (L1TP) data (i.e., processed through the Level 1 Product Generation System (LPGS) to full precision terrain correction) that have well-characterized radiometric quality and are cross-calibrated among the different Landsat sensors. The extent of the LUE allowed its confinement to one Landsat scene. The temporal range of interest spans the two previous decades . During the period, several Landsat L1TP-T1 products over Zambia were missing, especially during the vegetative period when the vegetation reaches its maximum. Available scenes immediately after were chosen ( Figure 6); in particular, one scene for each year within the time window of analysis was identified and selected for classification purpose.  A total of 19 scenes were downloaded from USGS facility and preprocessed. A comprehensive scheme is shown in Figure 7. It is worth mentioning that starting from 2010, Landsat 5 TM products were no longer available over the LUE. Landsat 7 ETM+ and Landsat 8 OLI imagery were therefore adopted to cover the period until 2017. However, failure of the Landsat 7 occurred on 31 May 2003 and caused no-data stripes which affect product quality. This was corrected here in preprocessing stage. Unfortunately, the mission also did not acquire data on LUE during the vegetative period in 2002 and 2003 and, therefore, products of 2001 and 2004 were respectively used to fill the data gap.

Growth Cycle and Biomass Dataset
To compute grass yield, a large dataset of scenes is required. In fact, products characterized by high revisit time and stability of acquisitions in terms of temporal coverage lead to full description of the quality indicator variations and growth rate predictions over pastures. This is a primary concern especially in regions where the rainy season is the period of maximum grass growth. These regions, like Zambia, experience diffuse cloud cover which limits land monitoring. In this context, MODIS (Moderate Resolution Imaging Spectroradiometer) sensor provide consistent, spatial and temporal time series comparisons of land conditions that can be exploited for further analysis. In particular MODIS VI products estimate vegetation conditions by measuring the vegetation activity at the surface to support an accurate seasonal monitoring. The products are produced and provided at 16-day intervals starting from the daily MODIS surface reflectance, atmosphere-corrected for molecular scattering, ozone absorption, and aerosols. A MODIS compositing method and a constrained view angle approach have been developed to remove low-quality pixels and for pixel selection over the compositing period. Furthermore, quality assessment (QA) flags with thematic information about the quality of the VI product, and provided along with the data, give precise control during analysis, filtering, and application phases [76]. The MODIS sensor Terra was launched on 18 December 1999, and provided the first 16-day VI scene on the 1 January 2000 at its finest spatial resolution of 250 m. Therefore, MODIS was not sufficient to cover the whole period of analysis back to 1996. For this reason, SPOT-Vegetation onboard on SPOT 4 was also considered. SPOT-Vegetation-1 VI (i.e., normalized difference vegetation index) is provided at 1 km resolution. The available acquisitions from 1998 to 2000 have no temporal consistency producing a lack of information from December to April. The latter aspects, together with the coarser resolution, brought us to compute VI also from 30 m Landsat 5 TM. A comprehensive scheme is shown in Figure 8. A total of 482 products were exploited to compute pasture growth cycles and biomass estimations and respectively 411 MOD13 (i.e., NDVI) scenes from 2000 to 2017 while, from 1998 to 2000, 29 SPOT-Vegetation NDVI estimates. Finally, 42 Landsat 5 TM scenes from which NDVI was computed after preprocessing steps using an automatic procedure, implementing its standard formulation, were considered to cover the period from 1996 to 2000. Preprocessing of Landsat regarded radiometric correction, atmosphere correction (BOA), and cloud masking. Original MODIS sinusoidal projection was converted to the UTM reference system by applying reprojection and format conversion. Quality assessment (QA) bands were exploited in each case for low-quality pixel removal in the preprocessing phase of filtering. It is worth mentioning here that MODIS and SPOT-Vegetation preprocessed NDVI data, made available through the sensor's data archives, are particularly useful and convenient for multiple reasons. Firstly, they allow speeding up the processing time, also reducing the number of processing steps to be implemented in the automatic procedure. Secondly, thanks to the extensive validation with multiple different sources of data, they result in particularly stable data including well-documented information on accuracy value and possible uncertainties. In MODIS, for example, as mentioned before, the validation stage produces per-pixel product quality information to mask out poor-quality data before use in research and applications. Product stability is also important when long time analyses are the objective of the research, since this assures homogeneity of outputs during the years.
Regarding the on-field grass yields used to retrieve biomass from satellite vegetation indices, no ground data were available over the study area. To overcome this, biomass field samples collected in the Alfred Nzo district of South Africa ( Figure 9) having similar eco-climatic conditions were used to retrieve an optimal regression formula. The data are part of the HerbMass database by the Grassland Society of Southern Africa (GSSA), a shared database for measurements of wet herbaceous biomass collected within Southern Africa. The used data consist of approximately 200 biomass samples estimated through 10 disc-pasture-meter readings within a~10 m radius (corresponding to a total of approximately 6 ha of grassland), collected on 12 January 2017. GPS coordinates of the specimen locations were recorded at 3 m accuracy.

Methodology
In the following, the proposed methodology for grassland analysis and monitoring in Barotse Floodplain will be presented. In particular, in the first section, the procedure to obtain historical information on grass coverage from vegetation land cover classification of the study area will be described, including details on the design of the multilayer perceptron (MLP) neural network (NN) developed for multispectral satellite sensors. Conversely, in the second section, the method for growth cycle extraction and time series of biomass estimation will be addressed, concerning also the development of the regression model for grass biomass.

Land Cover Classification
The estimation of long-term series of biomass data over pasture areas firstly requires those areas to be identified over the LUE along the time span of interest. The need for longest possible time series of biomass data for statistical analysis is required for continuity of acquisitions of satellite remote sensing data in terms of temporal coverage and product availability, as well as information content: it has to be possible to retrieve the identified information (e.g., index) with continuity even from different satellite products if the condition of satellite interoperability is satisfied or demonstrated. For these reasons, and with the objective to extend the analysis over the two previous decades (1996-2016), Landsat medium-resolution products were considered to guarantee the highest performance in terms of scalability and processing time. Besides standard reprojection and radiometric correction procedures, in order to improve product interoperability, atmospheric correction of the multispectral images was performed using the fast line-of-sight atmospheric analysis of spectral hypercubes (FLAASH) algorithm, which is based on a MODTRAN4 approach for path-scattered radiance, absorption, and adjacency effects [77]. Cloud pixel exclusion was performed by a cloud mask creation on the basis of the quality assessment band of Landsat products. A comprehensive scheme of the developed data processing chain described in this section is depicted in Figure 10.

Multi-Layer Perceptron Classifier
For the detection of grass areas, a land cover classification algorithm was developed. The algorithm had to be capable of producing multiple classifications over the period of the analysis and to be extensible over the whole country. For scalability and generalization capabilities, as well as for the power to accomplish fast automatic processing of massive amounts of data coming from multiple sources, NNs were considered. In particular, an MLP-NN was designed and trained, since MLPs are shown to be the most adopted and competitive approach for image classification [78,79].
MLPs are capable of determining their own input-output discriminating relations given a training dataset representative of each data class to be recognized. This means that it must include samples which are representative of the whole variance of the output classes.
The network was developed and trained using Neumapper free distributed software [80] to perform classification over the whole Zambia area and to identify 4 classes of vegetation, (a) forest/trees, (b) shrubland, (c) grassland, and (d) bare soil as well as a fifth class, (e) water. The samples used for training and testing were collected by visual interpretation of optical very high resolution products from 49 T1-Landsat 5 TM scenes within, or immediately after/before, the vegetation growing period in 2009 ( Figure 11) and then divided for training (80%) and test (20%) datasets. Due to cloud cover and missing scenes, in fact, the year 2009 was considered to ensure the best coverage over Zambia. Despite this, few products from 2009 were still missing and November-December 2008 scenes were used (i.e., beginning of the growing season). All 6 bands of the Landsat MultiSpectral Instrument (MSI) were considered. To improve the network performances further, statistical refinement of the overall dataset was carried out, obtaining the final population as reported in Table 1. Statistical refinements encompass the interquartile range (IQR) test for outlier identification and removal and class prevalence balancing. In our research, for the design of the network, a sigmoid function was used as activation function and applied to the input. The training samples were iteratively inputted to the network at each "epoch". The error backpropagation (EBP) learning algorithm was used as training scheme. A mean square error (MSE) function was iteratively minimized by a gradient descent technique, which changed the weight values according to a chosen learning rate. At the moment, no precise solution exists to finding the optimal topology of the network. Typically, heuristic methods are used to increase the accuracy, such as pruning and growing algorithms [81]. In this research work, a growing method was adopted by progressively increasing the nodes until performance could no longer be improved, and the final network topology was found on the basis of the maximum accuracy (minimum error) achieved during the training and test phases. Final topology counts two hidden layers, resulting in 6-9-7-5, where 6 is the number of inputs or the number of considered spectral bands, 9 is the number of neurons in the first hidden layer, 7 is the number of neurons in the second hidden layer, 5 is the number of output classes. The network reached its optimum learning after about 1500 epochs.

Grass Area Extraction
The trained MLP-NN was applied for land cover classification of Landsat 5 TM, Landsat 7 ETM+, and Landsat 8 OLI products, from 1996 to 2016, for pasture areas retrieval over the LUE, even if, as mentioned before, the classification step has to be intended to be extensible to the whole country area. From the 30 m classifications, pastured areas were then extracted in the masking step to compute vegetation analysis only over grassy parcels. According to cultivation types, Zambian crops during growing season may have had similar spectral response to grassland with no specific geometric patterns to be clearly distinguished and identified. Multitemporal analysis would then be required [82][83][84], but it was not a convenient solution in view of the temporal and spatial coverage of this study. Moreover, crops in the in LUE are not very diffuse (approximately 3% of the total area). Thus, eventual effects related to crop rotation and evolution during the period are not expected to substantially affect the final result. This is also supported by considering the traditional system of farming still characterizing the area, which is generally stable. Then, to recognize and exclude crops from the subsequent biomass analysis on grassland, the ESA Climate Change Initiative (CCI) 20 m land cover [85] was considered for crop mask extraction.

Growth Cycles and Biomass Estimation
Grass biomass has many implications in the management of grazing areas and the adopted farming practices. In fact, it is a measure of grass abundance and can be interpreted as feed yield for livestock. It is worth mentioning that biomass data availability on grassy areas is quite hard to address and, typically, their coverage is limited to precise study areas where the information is manually collected. Multispectral remote sensing allows retrieval of phenological indicators through the computation of several vegetation indices (VIs) which are empirical, and robust measures of the photosynthetic activity of plants is designed to enhance the spectral signal by combining two (or more) wavebands. According to several studies [29,30,49] normalized difference VI (NDVI), which combines the red, 0.6-0.7 µm, NIR, and 0.7-1.1 µm, wavelengths, can be assumed as the best predictor of ecosystem attributes especially in semiarid regions, due to its sensitiveness to low-density cover. From the first optical sensors, NDVI has always been largely adopted for vegetation analysis. Moreover, due to similarities in spectral resolution, NDVI derived from different sources has been used interchangeably in several cases [86][87][88][89] with the aim of computing the analysis over long periods. For example, in their research, Lange et al. [90] reported that MODIS-based vegetation NDVI agreed well with Sentinel-2 NDVI, also compared to ground-based estimations. In their study, Brown et al. [91] compared, among others, Landsat ETM+, MODIS, and SPOT NDVI, finding an encouraging similarity in the curves representing the annual green-up and brown-down of the vegetation. In addition, Boccardo et al. [92] found good homogeneity in comparing MODIS and Landsat NDVI estimates. Moreover, when large and homogeneous surfaces are taken into account, constraint coming from differences in sensor spatial resolution can be generally overcome [88,92,93]. In relation to this, preliminary analysis on NDVI consistency among the different sensors was carried out with a resulting general correlation coefficient of approximately 0.8 between MODIS NDVI and Landsat 5TM NDVI over various vegetation types. Consistency among SPOT-VGT and Landsat 5 TM for the overlapping period was also determined over the patches and the results are shown in Section 5.
It is also important to mention that image co-registration and preprocessing techniques for radiometric and atmospheric corrections were applied in this study to maximize sensor interoperability. A comprehensive scheme of the data processing chain developed for grass growth cycles and biomass is depicted in Figure 12. It is worth remembering that the normalized difference vegetation index from Landsat multispectral scene was computed after the preprocessing phase. Conversely, NDVI products from SPOT-Vegetation and MODIS MOD13 are provided within the sensor's data archives, but still require preprocessing procedures. NDVI from MODIS, SPOT-Vegetation, and Landsat 5 TM were computed against Landsat grass masks to retrieve NDVI-grass related values. To this purpose, grass masks at 30 m resolution were downscaled to MODIS and SPOT resolutions through resampling processing phase using a bilinear interpolation algorithm. For robust and automatic grass growing cycles retrieval, average and standard deviation of NDVI values over pasture areas were computed in Python on a patch-based approach. In our case, the LUE (grassland and no grassland areas) was subdivided into patches 50 × 50 km wide. The implemented algorithm however allows the process scalability since it permits to automatically choose the dimension of the patches to increase or decrease the spatial resolution of the analysis. In this way, NDVI estimates can be obtained only over grassland pixels and, for each patch, mean and standard deviation are estimated in each single time-step, corresponding to the reference-time of the NDVI products.

Regression Model
One of the most common approach for biomass retrieval by multispectral data consists of regression analysis fitted on the available ground samples for the specific geographic region and vegetation type [29]. However, no ground data were available over the study area. To overcome this, biomass field samples collected in the Alfred Nzo district of South Africa were used (Figure 9). Based on the sampling time and geographic coordinates of the field samples, we obtained the NDVI values derived from corresponding multispectral data to establish a database of NDVI vs. biomass values. In particular we used a Landsat 8 OLI product sensed on 01 January 2017 and a Sentinel-2 (S2) multispectral product sensed on 12 January 2017.
Landsat data were preprocessed as previously explained for NDVI calculation. The S2 product was processed in ESA SNAP Toolbox for atmospheric correction, cloud masking (pixel exclusion) and NDVI computation. Biomass dataset was preprocessed for outlier removal and divided into train and validation datasets, on an 80%-20% basis. To increase population dimension of validation dataset and regression consistence, samples from HerbMass database collected on 21 Jan 2016 over the same geographic area were also considered. In particular, after dataset refinement, they consists of 66 sampling points of grass biomass again estimated by disc-pasture-meter technology.
Based on the correlation between NDVI and biomass and on the experience of [29], we compared three regression models: linear (Equation (1)), a power function model (Equation (2)) and an exponential function model (Equation (3)) based on: where B is the biomass, a and c are the regression parameters. We evaluated the precision of the estimation models by using three metrics: the root mean squared error (RMSE), the mean absolute error (MAE) and the mean relative estimation error (REE). The RMSE, MAE and REE were calculated as follows (Equations (4)- (6)): θi is the i th sampled biomass, θ is the mean value of the sampled biomass,θi is the i th estimated biomass, and N is the sample size.

Long-Term Grass Biomass Estimates
As final remark of this section, from remote sensing long-term NDVI estimates, grass biomass values were calculated over the time window of analysis through the regression formula determined from regression analysis as described in the previous sub-paragraph.
Biomass mean and standard deviation are estimated by 1st order Taylor expansion to approximate the non-linear differentiable function by linear function in the neighborhood of a point p (Equation (7)): where X is the variable NDVI, Y is the biomass, g is the non-linear regression formula and g its derivative. If p = E(X) = E(NDVI), it is possible to estimate the mean g and the variance V(Y) of biomass as (Equations (8) and (9)): where V(Y) is the variance of NDVI.
The results are at the spatial resolution of the chosen patch dimension and at the temporal resolution of the NDVI satellite products.
The long temporal series of grassland biomass will be furtherly used to assess the climate change impact on biomass production of Zambian pastures, to provide observations and comments on future estimates of grass yields availability to the livestock sector.
As a conclusion of the section, summary of the considered scenes and sensors among the various stages of the methodology is presented in Table 2. Table 2. Summary of the used scenes and sensors in this study and divided according to their contribution to a specific processing stage of the methodology.

Results
This section is dedicated to the discussion of the obtained results. For a clearer presentation land cover and grass biomass analysis have been reported in the following dedicated subsections.

Land Cover Classification
MLP showed the expected advantages in fast and automatic massive processing, reducing the elaboration time and computational efforts to accomplish the classification task to less than 1 min per product. For the validation of the net a total of 800 reference points were used. In particular, 100 points were determined on a random basis within the dimension of the Landsat scene, forming a non-uniform grid. The grid was then automatically superimposed to each classified MSI image for accuracy evaluation on the whole Zambia region. Given the temporal and coverage extent of the analysis required by this research study, validation was performed using 8 Landsat products randomly selected from the dataset on the whole Zambia region ( Figure 13). Accuracy was calculated on each tile as points correctly classified over the total of the points, and verified through on-screen visual interpretation technique (Table 3).  Confusion matrix of land cover classification is depicted in Table 4, where also user's and producer's accuracies are shown. The first indicates how often a class on the map really exists on the ground, while the second tells us the goodness of the classification considering also the omissions (i.e., how often a feature is correctly classified in the map). We obtained good overall (0.92) and per-class results (≥ 0.87). In the confusion matrix, the kappa coefficient is 0.902. Forest/trees class shows the highest performances while bare soil class gets the lowest value (0.87) since it is sometimes confused with the shrub class. As expected, shrubland obtained 0.84 in user's accuracy meaning that it is frequently assigned to reference points belonging to other classes. Actually, its spectral signature may be influenced by the mutual presence of bare soil and trees. Notwithstanding that, producer's accuracy showed good result. Focusing on grass class, the trained network was able to correctly recognize grassy areas, with producer's and user's accuracy equal to 0.91, ensuring reliability for the subsequent biomass and growth cycles analysis. The neural network classifier was then applied to the yearly (from 1996 to 2016 included) 6-bands Landsat products over the LUE study area. Land cover maps in 1997 and 2016 are shown as examples in Figure 14. Mean extent of grassland class for the period 1997-2016 is approximately 408, 318, 557 and 148 km 2 respectively for Patch 1, 2, 3 and 4. The obtained results have been compared with historical reference studies reporting land forms and floristic presence within the Upper Zambezi River Basin and the Barotse system. The study of Baars & Jeanes [94] delineates the land classification of the West Province, with emphasis on the grasses, to assess the carrying capacity of the rangelands. The authors divide the area in Land Regions (LR) and Land Systems (LS), respectively according to the overall similarity in the pattern of the major woodland types, and to the distinctive vegetation similarity within a LR. LR and LS recognized over the LUE are depicted in Figure 15. LUE study area is mainly described by LR 4, the Luena river lowland, characterized by high presence of linear dambos (i.e., shallow wetlands predominantly grass covered). As confirmed by the latter study, in our classified images (Figures 14 and 15) grasses are mainly found within the annual flooded dambos (LS 4.1) which also support palms and bushes. Moreover, our time series of land cover classifications highlight the variation of the annual extent of the occurring floods. As climatic effect, also floodings restrict the grazing capacity of cattle in the area during vegetation growing season. The yearly not inundated extent of the grassy floodplain for Patches 2 and 3 in the LUE, and calculated from the classified products (Figure 16), is compared with the inundation extents of the Barotse Floodplain estimated on MODIS scenes by Zimba et al. [95]. As stated in the section dedicated to dataset description, the classifications were performed on the available cloud-free images within or in proximity of the growing season (December-May) and which differed from year to year. To increase results comparability, in the graph of Figure 16 we excluded extensions estimated over June and July products (see also Figure 6). The trend reported in Figure 16 is in well agreement with the results of [95].  Around the grassy floodplain, areas with the mutual presence of bare soils and shrubs can be seen in classifications. They are restrained in LS 4.2, which was effectively characterized by waterlogged open savannah and occasional patches of higher open woodlands [94]. In the north and east side of the LUE, forest/trees class dominates the landscape. There, LS 2.1, 2.2 and 3.1 describe respectively the presence of Cryptosepalum-Brachystegia-Julbernardia woodlands, Cryptosepalum forest and Kalahari woodland species. LS 3.1 moreover is characterized by dry plains and/or pans, according to the water filling. Circular and subcircular pans are also in the south-western part of LUE (LS 6.1) while in the northwestern part the very frequently or permanently flooded plain of the Zambezi is characterized by a complex pattern of old and recent channels, swamps and pools. Very good agreement was also found with the landscapes and grassland map of Zambia's Western Province by Kruyff [96].

Regression Analysis
Regression analysis was carried out to establish VI-biomass regression formula. As reported in the previous section, we used normalized difference vegetation index (NDVI) from Landsat 8 OLI and Sentinel-2 sensor, and reference biomass field samples collected by Herbmass dataset at the vegetation growth season (January) to test a power model and an exponential model. Simple linear regression was also considered. It is important to remember that Landsat scenes were not available at the actual sampling date, differently, a Sentinel-2 product was available exactly on the 12 January 2017. Moreover, the NDVI values of Landsat and Sentinel-2 against the reference biomass estimates by Herbmass dataset have a similar meaningful behavior. Certainly the difference in the days of the two acquisitions with respect to the biomass sampling time cannot be ignored. In fact, several changes might have been occurred in the time interval that might have produced deviations among the biomass and NDVI dataset. Considering that Sentinel-2A and Landsat instruments proved to have very good compatibility [97] due to the similarity on both spatial and spectral characteristics, we decided to build regression model considering the Sentinel-2 product only. Hence, correlation coefficients were calculated using all field samples collected on 12 January 2017 and Sentinel-2 data. NDVI was positively correlated to biomass, with a correlation coefficient of 0.71. Both the power and the exponential models well performed over the NDVI-biomass training dataset ( Figure 17). Linear regression and power models had the best performances on calibration dataset (Table 5) with RMSE equal to 572.29 kg ha −1 and 572.54 kg ha −1 respectively. Exponential regression model also performed well, with little difference of performance indicators.  When tested on the validation dataset (Table 6 and Figure 18), the power model showed the best performance with 539.73 in RMSE. Conversely, the linear model is not able to correctly predict biomass values, with larger RMSE, MAE and REE respect to the other models. Power regression model was then adopted for biomass estimation over the LUE area from MOD13, SPOT-VGT and Landsat 5 NDVI.

Growth Cycles and Biomass
The retrieval of grass growth cycles and grass biomass over the LUE along the 21 years of analysis has represented a delicate task as it involved the collection and the integration of heterogeneous data from different multispectral sensors with their own peculiarities. As previously reported, a preprocessing stage was applied to each sensor's dataset in order to maximize the interoperability. Starting from more homogeneous datasets, mean and standard deviation of NDVI were extracted only on grassland pixels within each 50 × 50 km patch ( Figure 19). The NDVI from MOD13 covers the longest time span from 2000 to 2016 and draws complete growth cycles along the years. SPOT-Vegetation NDVI estimates are not able to completely describe the grass evolution during the time interval from 1998 to 2000 with missing data around January 1999. Landsat 5 TM NDVI has missing information from November to April given just partial cycles for 1996, 1997 and 1998 years. NDVI from the different sensors are in well agreement. This is visible in particular for 1998 and 1999 when both Landsat 5 and SPOT-Vegetation estimates are available. Within a cycle (Figure 20) NDVI generally ranges from 0.38 to 0.65. The values appear to be generally low, but they are in accordance with the reported values of biomass for the area. Moreover, it is important to recall that these results are reported as averages within the patches, so they are the result of an average formula applied to both higher NDVI values expected within the channels and lower values in the more arid zones. Lowest NDVI values can be found from July to October, while highest values are from December to April, reflecting the eco-climatic variability within a year for the area. Positive peaks correspond to the wet season and the negative peaks to the dry season. Mean NDVI values (dots) are displayed along to the relative standard deviations (bars). SPOT-Vegetation, MODIS and Landsat NDVI are in the order from the lowest to the highest standard deviation σ. This behavior appears to be correlated to the sensor resolution. Landsat, with its 30 m resolution, is in fact able to discriminate among different types of grasses within the grassland class while, by increasing the pixel size, or rather, by decreasing the product resolution, differences in species and growing stages are less likely to be retrieved producing low deviations from the mean value of NDVI. Standard deviation σ is higher in patches 1 and 2, at the top of the LUE, respect to the patches 3 and 4 at the bottom. Moreover, by looking at the standard deviation of MODIS NDVI within a cycle ( Figure 20 and Table 7), it differs from 0.086 when NDVI is at maximum in the warm season to 0.069 and 0.076 when NDVI is at minimum in the dry winter season. This again can be associated with different grass growth conditions and grass species within a patch, difference obviously less marked when grasses are at their minimum of photosynthetic activity. From NDVI values, biomass is estimated with power regression formula (Equation (2)) and regression parameters as retrieved in the previous section.
Biomass is depicted in Figure 21, where also the precipitation trend from CRU dataset is reported [75]. The measured grass yield generally goes from 1.8 ton ha −1 to 4.1 ton ha −1 . Our estimates agree with the findings by Baars & Jeanes [98]. In their study the authors reported standing crops of grass ranging normally from 1.5 to 4.0 ton ha −1 with occasional higher quantities especially in floodplain channels. Considering the maximums (in January) and minimums (in August) of biomass in each year, all the patches are characterized approximately by the same behavior with patch 2 showing slightly higher values. Period of high biomass presence (≥4 ton ha −1 ) is approximately 2 months long (≈ 50 days), while period of low grass yields (≤2 ton ha −1 ) is approximately 80 days long ( Figure 21). Moreover, the biomass trend is also in agreement with the precipitation variability ( Figure 21, Figure 22 and Table 8). Biomass values rapidly increase in correspondence of the beginning of the wet season in October, and gradually decrease from the end of May when the cold and dry season starts. Peak of precipitation generally occurs in December-January. In the same period peak of biomass is visible, as the two variables are positively correlated (Figure 21). In Figure 22 is also visible that precipitation is more correlated with 1 month later grass yields (correlation coefficient r = 0.88) respect to biomass of the same month (0.75) or of 2 months later (0.76). Moreover, when grass reaches a biomass equal to 4 ton/ha, increase in total precipitation does not produce further increase in yield. From graphs in Figure 21 we can also observe anomalies of biomass trends, as for example longer or shorter period of low grass yields, linked respectively to situations of vegetative stress or of nutrients/water abundance.   In Table 8 and Figure  As we expected, higher precipitation value does not have influence only over maximum and minimum of biomass, but also over the duration of the stages. Shortages of water in 2005 produced a reduction of period of high biomass presence (≥ 4 ton ha −1 ) from 50 days to 15 days approximately, and an increment of 50 days in duration of low available grass yields for grazing (approximately 128 days ≤ 2 ton ha −1 ). Conversely, in 2004 high availability of biomass was possible along 3 months and poor yields lasted only for 1 month and a half. The retrieved growth cycles were also validated on vegetation productivity data from FLUXNET tower in Mongu [99]. Vegetation productivity on flux towers is estimated through the eddy covariance (EC) technique [100][101][102][103][104][105], a micrometeorological method based on turbulent transport theory for the estimation of surface-atmosphere exchanges ecological (e.g., water vapor, carbon dioxide, methane and heat fluxes). The eddy covariance technique is particularly valuable for characterizing natural or managed ecosystems by monitoring their exchanges of gases and energy with the atmosphere. FLUXDATA are considered as one of the primary source for the development and validation of satellite-based productivity models [40,106]. Africa has just few towers belonging to the African regional flux measurements sites CARBOAFRICA, and principally located in forested areas. Fortunately, LUE area is covered by the Mongu flux site located in the south (GPS coordinates -15.4378 N, 23.2528 E) close to the Zambezi River and designated for savannah grass. We used available daytime GPP measurements, between August 2007 and August 2009 ( Figure 24).

Discussion
In the present study optical data from Landsat 5 TM, Landsat 7 ETM+ and Landsat 8 OLI were used to accurately classify a subtropical area in Southern Africa. In particular, the classification algorithm was designed to perform the land cover classification task on the whole Zambia country and over a 21-year time window (from 1996 to 2016). To this purpose, a Multilayer Perceptron Neural Network MLP-NN was used for its generalization capability and trained using the 30 m bands of Landsat products. Datasets for training and testing were populated from multiple products collected over the Zambia country. From the accuracy assessment, MLP-NN demonstrated effectiveness and robustness for multitemporal and wide-spatial analyses as required by studies on climate change effects and impacts. The separation of different vegetation types is a challenging task due to the smooth transition among successional species. The best result was achieved in classifying the forest and water classes (Table 4) principally due to the homogeneity in their spectral response within the resolution cell. Such homogeneous response is also reflected by the distribution of reflectance values in train and test datasets, which is characterized by lower deviation from mean value for forest and water types. Nevertheless, we found a low separability of shrubland and bare soil classes (Table 4). Shrubland is vegetation transition zone, and is principally characterized by shrubs, trees, and grasses at various densities. To this reason, this class is sometimes confused with the others. According to vegetation greenness, shrubland and grassland are eventually attributed to bare soil class if plants are particularly dry or dominated by a woody component. Grassland is confused not only with bare soil and shrubland, but also with water, where riparian vegetation and aquatic plants are present. Despite this, accuracy is high and class reliability is guaranteed.
Extending the classification task to the country Zambia, the retrieved vegetation distribution reflects its actual presence over the region. Forests are found principally in the northern part of the country (Figure 26), in the form of wet Miombo woodlands, while rest of the country is as dry Miombo woods. These woods are partially comprise shrublands where open forests occur and trees are sparse. Bare soils and dry lands mainly occur in the west part of the country, while grasslands suitable for grazing and supporting large herds of cattle are principally located in correspondence of floodplains and water courses and flatlands. In the LUE area, grasses dominate the landscape (Figure 14), where large floodplains are present. The grassland is principally regulated by flooding, which annually covers the grazing areas. Occasionally, flood-time satellite imageries were the only possible exploitable products in the vegetation growing season, making it difficult to properly estimate the extent of interannual variation in grasslands. However, it allowed us to verify the severe reduction of grazing areas during flooding.
To compute grass cycles, we selected NDVI for use in this study, as it is one of the most widely used vegetation indices for biomass estimations, and it has the highest sensitivity for less dense vegetation coverage and is thus more suitable for estimation in semiarid areas. Moreover, its relative simplicity in estimation leads to the possibility to obtain it from multiple sensors and back in time. In this research, we used NDVI from Landsat, SPOT-Vegetation, and MODIS sensors. Simultaneously, the use and integration of such products was made possible by considering the similarities in spectral resolution and also the fact that when large and homogeneous surfaces are taken into account, constraint coming from differences in sensor spatial resolution can generally be overcome. Moreover, a detailed preprocessing stage was applied in each case to maximize sensor interoperability. The integration of NDVI from different sensors allowed the cycles to extend two decades back in time, from 1996 to 2016, with the final result of 21 years of grass biomass data over the study area in Zambia. MODIS allowed covering the longest period, from 2000 to 2016, while Landsat and SPOT-Vegetation were used in the earliest 4 years. The NDVI data in the case of MODIS were produced and collected on a 16-day basis, allowing product stability and availability (no missing temporal windows). This is guaranteed by the compositing method of the MODIS MOD13 products. On the contrary, data from Landsat and SPOT-Vegetation sensors were unavailable, especially at the maximum of the phenological season, corresponding to the maximum grass yield, due to both cloud coverage and data unavailability in the databases. In fact, one of the principal problems of conducting research in the African continent and, in particular, the Sub-Saharan region, is the low amount of collected data; moreover, in tropical-subtropical climate regions like Zambia, a rainy season (semi-permanent cloud coverage) corresponding to summer, matches with the vegetation growth period. This contributes to lowering the suitable optical imagery exactly in the period of the year of greatest interest in vegetation analysis, since clouds mask earth surfaces. Despite the missing information, NDVI from different sources is homogenous. In particular, from 1998 to 2000, both Landsat and SPOT-Vegetation NDVI was obtained to increment biomass estimates and to validate data integration. In Figure 27, the relationship between the NDVI from Landsat 5 TM and the NDVI from SPOT-Vegetation on grasslands is depicted. The two are in good agreement and points are close to the diagonal line (correlation coefficient equal to 0.95). The grass biomass was calculated from NDVI estimates. To retrieve the model for biomass from NDVI data, a regression approach was adopted. We tested a linear and two non-linear regression models: exponential and power formula. All three models performed well over the regression data, but the results of power function were superior to the other models. This is also in agreement with the findings of Jin et al. [29] for steppe in the Xilingol region characterized by semiarid grasslands. In particular, we tested the regression formula derived by the authors for typical steppe region in China on the ground biomass data sampled in Southern Africa, and used in this research work. The outcome is particularly interesting since pretty good agreement was found ( Figure 28): the correlation coefficient between biomass data estimated by Jin et al. using a regression formula [29] and sampled data is equal to 0.72, and comparable to those obtained with our model (0.73). Therefore, even if retrieved for a different geographic area, the regression obtained on Chinese grassland can be extended to African grassland principally due to similarities in vegetation characteristics and climate. Generally speaking, a priori knowledge on vegetation characteristics and climatology allows the use of already developed regression analysis even on different areas of the globe, at least for preliminary studies. This is particularly important when dealing with sites where there is poor availability of data. This notwithstanding, further analysis and dedicated modeling is still required for best-fitting and validation. By considering the biomass estimates over the 21 years of analysis, from 1996 to 2016, no significant trends are visible in the data. Average biomass values at the end of phenological season (when standing crops are at maximum) is constant with time indicating that, in LUE area, grasslands and, thus, the grazing areas are currently stable and do not suffer for adverse climatic trends of variation. If the resulting interannual variability does not seem to reveal the occurrence of critical phenomena, the intra-annual mean biomass reflects grass growth. Grass yields fluctuate according to the seasonality and climatic conditions. Differently, intra-annual standard deviation in biomass reveals the complexity of variation in plant conditions that may be caused by differences in grass species composition, soils, proximity to water, and human-induced changes. No differences in the estimated yields are found among the four patches in the LUE, while the standard deviation is higher in the northern patches. This, again, could be due to differences in grass types from north to south and, by looking at the spatial distribution of grassland within the LUE, variability in spatial configuration can be clearly seen. In fact, it can be noted that the northern patches are characterized by more various environmental conditions respect to the southern patches. In particular, in patch 1 and 2, grassland is found in correspondence with both major rivers and the inner land, while in patch 3, there is a clear prevalence of wetland grasses, and in patch 4, there is a dominance of inner-land species.
Rainfall (soil moisture) is the first limiting factor affecting grass growth. If insufficient precipitation occurs before and during grass phenological season, plants suffer and a negative effect on both duration of stress and total yields are present. In droughty years, such as 2005, both the minimum and the maximum of grass biomass were lower with respect to reference years. However, the strongest impact is on duration of period of minimum and maximum. In fact, both an important reduction in the length of period with high biomass values, and a large increment in the duration of the period with low available grass yields are detected, with further shortcomings in cattle feed resources especially in the winter season. Current water shortages are principally contributing by temporarily stressing and degrading grass conditions by elongating the vegetation dormancy period and retarding regrowth. Plant tolerance to drought varies from species to species. When water shortages occur, photosynthetic activity decreases and total chlorophyll reduces as well as the water potential in cells. This study provides evidence that Barotse grasslands respond quickly to decreases and increases in mean precipitation. In particular, variations in grass photosynthetic activity and yields become more visible after one month from the occurrence of variations in precipitation. If water deficiency stress exceeds the plant resistance to drought, then prolonged exposure and high intensities are likely to result in permanent damage to plant tissues and cells, and in the worst cases, the death of the whole plant may occur [107]. A long-term decline in grass yields according to water availability is expected if the moisture content of soil is not sufficient to support vegetation renewal from one year to another. Longer droughts, involving several years, will have a large impact on grass biomass with important effects on available feed yields. The primary effect of precipitation in limiting and lowering the productivity, and in retarding a new cycle, discloses the crucial contribution of long-term weather forecasts to projecting future grass conditions of grazing areas. This supports the adoption of the MOSAICC tool by FAO to model the climate change impact on biomass. In fact, by deriving statistical relationships between historical biomass data and the observed or simulated water balance variables and weather variables, the model predicts vegetation yields and allows evaluation of possible future climate impacts on biomass. Such results and projections can be considered crucial for government, stakeholders, and farmers, and can then be used in FAO's GLEAM model to enhance ecosystem conservation by testing, for example, the adoption of rotational grazing and deferment strategies and to support livestock management.

Conclusions
Interaction between grazing animals and grassland used as pastures has been addressed by several studies in the past decades [108]. Most of the studies aimed at developing management strategies to avoid long-and short-term effects on vegetation and soils, species composition, and pasture biomass. When grazing occurs on natural perennial grasslands, on the one hand, precise research on the livestock sector and cattle management and practices (such as rotation, transhumance, or deferment strategies) are required [53]; on the other hand, analyses of veld are necessary to identify grass characteristics and sensitivity to different types of pressures.
This study aimed at responding to the need for a methodology that could perform a more comprehensive, automatic, and repeatable analysis of grassland biomass over long time windows and variable spatial scales, and that could deal with the lack of on-field data. The study presented an end-to-end automatic EO data processing approach exploiting 20-year time series of historical data in significant regions of Africa where food security is a critical problem. In fact, this research work focused on the analysis of grassy areas in Barotseland, Zambia, to support the livestock sector and increase its resilience to climate change. In particular, the materials and methods adopted for the characterization of grass distribution, grass growth cycle, and biomass were described. In general, the main challenge of climate change-related studies is the necessity to rely on long historical datasets to enable the identification of long-term trends, changes, and other statistical indicators, such as the correlation between climatic and non-climatic factors. This study addresses the challenge by performing the analysis over a very long time window (21 years), which enabled investigation of the correlations of productivity with climate variables to clarify sensitivity to climate change and possible effects already in place.
To process and integrate the huge amount of data, this research work exploited the automatism and generalization capability of the MLP-NN algorithm for land cover classification, which has been integrated in a processing chain capable of multitemporal analysis, producing biomass trends as a final output. The results obtained on the 21-year long time series highlight the effectiveness of MLP-NN in accomplish massive processing of satellite multispectral data coming from different sources. Moreover, the research provides evidence that medium-quality multispectral data are suitable for subtropical vegetation identification and extraction of grazing areas in remote regions of Africa. The established workflow can provide a robust methodology which is capable of retrieving biomass at various resolutions and extensions as it is characterized by a high degree of automatism and scalability. According to the results of this research, Zambia pasturelands were found to be vulnerable to changing climate and, in particular, to water shortages that come with drought periods. Despite this vulnerability, the estimations do not suggest that grass quality decline is occurring. The yearly reduction in the extent of grassland principally corresponds to variation in the extent of floods in Zambezi, while the reduction in grass quality and availability comes from a temporary increase in the length of a water stress period and a retard in the beginning of regrowth. At the moment, grazing on the floodplain is impractical from around January to June [98], due to floods, and in that period, cattle are moved inland. Future reduction in grass quality and extent is prone to be more limiting for local livestock productivity. Although exploitation of long time series of grass biomass values increases the possibility of regression capability on climatic variables and the certainty of the determined increments and decrements, field sampling remains necessary to validate the estimations.
Current previsions on Zambia climate suggest that what actually is still an eventuality is likely to become more frequent (e.g., droughts and severe floods), creating urgency for shared and sustainable approaches. Such estimations are clearly important for the definition and establishment of grassland conservation strategies and grazing management to overcome and reduce negative effects derived from a changing climate. Suggested initiatives include rotational grazing to reduce the pressure on preferred species, community participation, shared management programs to encourage farmers in matching cattle numbers to feed resources, and the introduction of leguminous crop.
For more detailed analysis, the possible effect of grassland degradation deriving from changes in grass composition should be investigated and included in future projections of feed yields. In fact, the disappearance of leguminous crop caused by overgrazing reduces the potential for nitrogen cycling and leads to lower nitrogen-use efficiency.
This research work emphasizes the important contribution coming from Earth observation in conducting large-scale analyses of vegetation with the added value of multitemporal and historical evaluation of changes and trends. The latter aspect is even more valuable in remote areas where ground data are absent, scarce, or have poor homogeneity and continuity with regard to acquisition in addition to low availability of technical resources and personnel. The use of SAR images should also be considered in future. Indeed, these data have already been shown to contain useful information, complementary to that provided by optical systems, for biomass estimation of vegetation [109]. Hence, an update of the whole methodology including these data may be appropriate.
The ongoing increase in provided satellite remote sensing products and the involvement of the scientific community in designing and developing Earth observation missions more focused and functional to specific objectives and themes give new and great opportunities for analysis to interpret climate change at different scales. This, moreover, allows the establishment and definition of mitigation strategies which start from broader and more accurate awareness and knowledge of the problems. Of note, the approach is extensible and applicable to wide areas at different scales to support different entities, from small farmers to whole countries, continents, and world regions.  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.

Conflicts of Interest:
The authors declare no conflict of interest.