Louise Leroux
email: louise.leroux@teledetection.fr

Christian Baron

Bernardin Zoungrana

Seydou Traore

Danny Lo Seen

Lo Seen

Agnes Begue

Leroux

Crop Monitoring Using Vegetation and Thermal Indices for Yield Estimates: Case Study of a Rainfed Cereal in Semi-Arid West Africa

Keywords: Remote Sensing, Crop yield, NDVI, Land Surface Temperature, Crop model, MODIS, Rainfed cereal, Niger, Harvest Index

come

Introduction

In the Sahelian region of West Africa where traditional rainfed agriculture prevails [START_REF] Marteau | The onset of the rainy season and farmers' sowing strategy for pearl millet cultivation in Southwest Niger[END_REF], over 20 million people suffered from food insecurity in 2014 [START_REF]2014-2016 strategic response plan : Sahel region[END_REF]. Sahelian rainfed farming systems are known to be at high climatic risk due to a high spatio-temporal variability of rainfall and frequent drought events [START_REF] Akponikpè | Spatial fields ' dispersion as a farmer strategy to reduce agro-climatic risk at the household level in pearl milletbased systems in the Sahel : A modeling perspective[END_REF]. Rainfall variability results in large fluctuations in year-to-year crop productivity which leads to episodes of food insecurity. Moreover, the political and socio-economic instability of certain countries in the region also contribute to the variability of agricultural production [START_REF] Traoré | Characterizing and modeling the diversity of cropping situations under climatic constraints in West Africa[END_REF]. These considerations highlight the need for an operational, timely and accurate yield estimation system to assist decision-making [START_REF] Justice | Developing a strategy for global agricultural monitoring in the framework of Group on Earth Observations[END_REF]- [START_REF] White | A continental responses phenology model climatic for monitoring variability vegetation to interannual[END_REF].

Yield estimation systems based on crop modelling allow accurate quantitative assessments (e.g. AGRHYMET in West Africa; the AGRI4CAST action in Europe), but are confronted with input data availability and spatial consistency constraints [START_REF] Huang | Meta-analysis of influential factors on crop yield estimation by remote sensing[END_REF], [START_REF] Son | Prediction of rice crop yield using MODIS EVI-LAI data in the Mekong Delta, Vietnam[END_REF].

For more than two decades, Earth Observation systems have been known to play a significant role in vegetation monitoring by providing synoptic, repetitive, timely, objective and cost-effective information on Earth's surfaces (e.g. [START_REF] Tucker | Satellite remote-sensing of total drymatter production in the Senegalese Sahel[END_REF]- [START_REF] Prince | Monitoring the vegetation of semi-arid tropical rangelands with the NOAA-7 Advanced Very High Resolution Radiometer[END_REF]). They have been acknowledged for their valuable contribution to spatial and temporal monitoring of global vegetation and thus have been used extensively in many parts of the world for crop condition monitoring and yield forecasting [START_REF] Son | Prediction of rice crop yield using MODIS EVI-LAI data in the Mekong Delta, Vietnam[END_REF], [START_REF] Becker-Reshef | A generalized regression-based model for forecasting winter wheat yields in Kansas and Ukraine using MODIS data[END_REF]- [START_REF] Yang | Evaluating high resolution SPOT 5 satellite imagery to estimate crop yield[END_REF]. Combined or not with rainfall data, satellite data is currently being used in early warning systems to assess crop development conditions during the growing season (e.g. FEWS-NET (Famine Early Warning System Network); GIEWS (Global Information and Early Warning System); FOODSEC (Food Security); GEOGLAM (Group on Earth Observation -Global Agricultural Monitoring)), mainly through the production of regularly updated crop growth anomalies' maps based on the Normalized Difference Vegetation Index (NDVI). These systems benefit from timely and synoptic satellite data which compensate the lack of reliable and homogeneous ground data, but which remain mainly qualitative and do not include crop yield monitoring.

The empirical relationship between remote sensing spectral vegetation indices and in situ observations to predict crop yields before harvest has been tested for a long time in many studies (for a review see [START_REF] Rembold | Using Low Resolution Satellite Imagery for Yield Prediction and Yield Anomaly Detection[END_REF]). The simplest approach involves the regression between observed yields and vegetation indices, either on a specific date or through a time integral of vegetation indices between two dates [START_REF] Lobell | The use of satellite data for crop yield gap analysis[END_REF]- [START_REF] Wall | The early explanatory power of NDVI in crop yield modelling[END_REF]. Among the vegetation indices, the NDVI has been widely employed due to its close relationship to several vegetation parameters like the Leaf Area Index (LAI), the fraction of Absorbed Photosynthetically Active Radiation (fAPAR) or the green biomass [START_REF] Tucker | Red and Photographic Infrared linear Combinations for Monitoring Vegetation[END_REF]- [START_REF] Hatfield | Leaf-area estimates from spectral measurements over various planting dates of wheat[END_REF]. Furthermore, several studies have found a good correlation between NDVI and crop yields in many study sites around the world [START_REF] Becker-Reshef | A generalized regression-based model for forecasting winter wheat yields in Kansas and Ukraine using MODIS data[END_REF], [START_REF] Rasmussen | Assessment of millet yields and production in northern Burkina Faso using integrated NDVI from the AVHRR[END_REF], [START_REF] Rasmussen | Operational yield forecast using AVHRR NDVI data: Reduction of environmental and inter-annual variability[END_REF], [START_REF] Maselli | Processing of GAC NDVI data for yield forecasting in the Sahelian region[END_REF]- [START_REF] Groten | NDVI-crop monitoring and early yield assessment of Burkina Faso[END_REF].

Nonetheless, there are intrinsic limitations that prevent an operational use of vegetation indices to estimate crop yield. Apart from technical limitations such as low spatial and temporal resolution leading to mixed pixels and incomplete crop growth descriptions respectively, the main limitation is the indirect link between yield and spectral data. Since the 1980's NDVI is known to be a proxy of the aboveground biomass [START_REF] Tucker | Satellite Remote Sensing of Total Herbaceous Biomass Production in the Senegalese Sahel : 1980-1984[END_REF], [START_REF] Holben | Spectral assessment of soybean leaf area and leaf biomass[END_REF], but the ratio between yield and aboveground biomass (referred hereafter as the harvest index, HI) is also known to be highly variable in space and time. [START_REF] Monteith | Climate and the efficiency of crop production in Britain[END_REF], [START_REF] Monteith | Solar radiation and productivity in tropical ecosystems[END_REF] showed that biomass production is linearly related to fAPAR for crop with no water stress, while [START_REF] Myneni | On the relationship between FAPAR and NDVI[END_REF] showed that a linear relation between NDVI and fAPAR can be assumed since their functional response to leaf orientation, solar zenith angle and atmospheric optical depth is similar.

On the other hand, the harvest index, also known as the reproductive efficiency, is crop-dependent and sensitive to variables that impact the partitioning of the assimilates into grain, such as the genotype, temperature and water/nutrients availability ( [START_REF] Lecoeur | Harvest index increase during seed growth of field pea[END_REF]- [START_REF] Unkovich | Chapter 5 -Variability in Harvest Index of Grain Crops and Potential Significance for Carbon Accounting: Examples from Australian Agriculture[END_REF]). In the Sahelian region, HI is strongly dependent on water conditions during the growing period, in particular during the reproductive stage [START_REF] Misra | Assimilate partitioning in pearl millet (Pennisetum glaucum L.R.Br.)[END_REF]. Crop water conditions can also be derived from remote sensing data, with indices based on the difference between air and surface temperatures, which are useful indicators of water stress for yield estimation. Indeed, since the 1970's various remote sensing-based studies have shown that final yields can be related to thermal indices [START_REF] Idso | Remote-sensing of crop yields[END_REF], [START_REF] Smith | Relationship between wheat yield and foliage temperature: theory and its application to infrared measurements[END_REF]. Based on land surface temperature (LST), the Crop Water Stress Index (CWSI) proposed by [START_REF] Jackson | Canopy temperature as a crop water stress indicator[END_REF] was found useful for yield estimation and crop assessment (e.g. [START_REF] Pinter | Infrared thermometry: A remote sensing technique for predicting yield in water-stressed cotton[END_REF]- [START_REF] Gontia | Development of crop water stress index of wheat crop for scheduling irrigation using infrared thermometry[END_REF]).

In the framework of current early warning systems for food security, crop yield monitoring would certainly benefit from the consistency in space and time of remote-sensing based crop yield estimations.

For this reason, in this study, we investigate the possibility of combining vegetation and thermal indices for crop yield estimation in the Sahelian region, where, to our knowledge, this has not been attempted. Our objective is to build a simple, robust and timely satellite-based model for rainfed cereal yield estimates on the basis that: (i) aboveground biomass can be estimated using NDVI, and that (ii) LST data can provide useful information on the harvest index. Such a model would also provide effective assessments of yearto-year yield variability. The study is conducted in the South-West of Niger (the Niamey Square Degree site) where rainfed pearl millet dominates the agricultural landscape, and soils as well as agricultural pratices are relatively homogeneous. We use the SARRA-H (System for Regional Analysis of Agro-Climatic Risks) crop model [START_REF] Dingkuhn | Decision support tools for smallholder agriculture in Sub-Saharan Africa : A practical guide[END_REF], which has already been validated for pearl millet in the Sudano-Sahelian zone [START_REF] Traoré | Characterizing and modeling the diversity of cropping situations under climatic constraints in West Africa[END_REF], to simulate biomass and the corresponding yield for a period of 11 years (2000-2011), using as input data the rainfall measurements from 28 rain gauges and a meteorological station. We derive the NDVI and the CWSI from MODIS data over the same 11-year period to explore their respective relationship with biomass and the harvest index. The model is then validated using crop statistics data at the scale of the Niger Square Degree site. The proposed approach is finally discussed in the framework of a potential operational yield estimation system that would also include data from the upcoming Sentinel-2 sensors.

Material

Overall approach

We combine satellite-data with agro-meteorological modeling results to analyze the potential of MODIS derived NDVI and LST time series for pearl millet yield assessment in the Niger Square Degree site. The underlying assumptions of our approach are that:

(1) Aboveground biomass can be determined from vegetative indices such as the NDVI [START_REF] Hatfield | Remote Sensing Estimators of Potential and Actual Crop Yield[END_REF],

(2) Harvest indices can be significantly reduced under water-limited conditions [START_REF] Sinclair | Water-Use Efficiency Crop Production[END_REF] due to crop water stress. Land surface temperature (LST) observations can be used as an indicator of crop water stress [START_REF] Jackson | Canopy temperature as a crop water stress indicator[END_REF] and thus be related to the harvest index,

(3) The combination of NDVI and LST provides a better estimate of yields than the NDVI on its own in water-limited regions. Ideally, a remote sensing-based approach has to be calibrated with reliable ground-measurement data. For our study area, the ground truth data currently available are mainly based on farmer's declarative survey and suffer from a lack of consistency both in space and time. Consequently, the choice has been made to use simulated data from SARRA-H crop model to overcome this issue knowing that SARRA-H has been validated for this region [START_REF] Traoré | Characterizing and modeling the diversity of cropping situations under climatic constraints in West Africa[END_REF]. The predictive capacity of the remote sensing-based model is then verified at a regional scale with agricultural statistics.

Study area

The study area (12.9° -13.9°N; 1.6° -3.1°E, hereafter referred to as the NSD site) which includes the Niamey Square Degree site covers about 18,000 km² and is located in South-West Niger (Fig. 2a). The site is part of the AMMA-CATCH observatory (African Monsoon Multidisciplinary Analysis-Coupling the Tropical Atmosphere and the Hydrological Cycle; http://www.amma-catch.org/) and has been chosen for two reasons: (1) rainfall is considered as the main driver of crop yield [START_REF] Sultan | Agricultural impacts of large-scale variability of the West African monsoon[END_REF], and (2) the site is instrumented since the early 90's including a dense network of rain-gauges which are continuously recording rainfall. The climate is typically Sahelian. Annual ambient temperatures are high and rainfall distribution is monomodal during June-September. Rainfall is highly variable spatially [START_REF] Barbé | Rainfall climatology of the HAPEX-Sahel region during the years 1950-1990[END_REF] and temporally [START_REF] Hiernaux | Trends in productivity of crops , fallow and rangelands in Southwest Niger : Impact of land use , management and variable rainfall[END_REF] with 10-22% inter-annual variations between 2000-2010 (Fig. 2b). In addition, despite the small size of the study area (about 160 km x 110 km; Fig. 2c), the regional rainfall pattern shows a high latitudinal gradient from 480 mm/year (north of the study site) to 630 mm/year (south).

The production system is rainfed, dominated by pearl millet [START_REF] Traoré | Characterizing and modeling the diversity of cropping situations under climatic constraints in West Africa[END_REF] which is drought-resistant and well adapted to the sandy soils predominant in the study area [START_REF] Rockström | Water , nutrients and slope position in on-farm pearl millet cultivation in the Sahel[END_REF]. It is characterized by low inputs [START_REF] Baron | From GCM grid cell to agricultural plot: scale issues affecting modelling of climate impact[END_REF] and low yields (generally lower than 700 kg/ha; [START_REF] Rockström | Water , nutrients and slope position in on-farm pearl millet cultivation in the Sahel[END_REF]).

Satellite data 2.3.1. The MODIS Vegetation Indices product (MOD13Q1)

The MODIS Vegetation Indices (VI) product (MOD13Q1 collection 5) was used in this study because of its data consistency, providing spatial and temporal information on vegetation conditions every 16 days at 250-m spatial resolution since 2000 [START_REF] Justice | The Moderate Resolution Imaging Spectroradiometer (MODIS): land remote sensing for global change research[END_REF]. Even if the MODIS data are pre-processed with the CV-MVC (Constrained view angle-maximum value composites) algorithm, noise still exists in the time series due to cloudiness, sensor problems or Bidirectional Reflectance Distribution Function (BRDF) effects [START_REF] Chen | A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky-Golay filter[END_REF]. In consequence, we applied a Savitzky-Golay filter to reduce noise and improve the quality of the NDVI time series towards a more efficient crop yield monitoring [START_REF] Wan | Modis Land Surface Temperature Products Users' Guide[END_REF]. After testing different smoothing parameters, a filter width of 4 and a degree of smoothing polynomial of 6 were retained, which allowed to match the upper envelope of the NDVI time series

The MODIS Land Surface Temperature product (MOD11A2)

The MODIS LST product (MOD11A2, collection 5) is composed of the average value of daily 1kilometer LSTs under clear sky conditions for an 8-day period [START_REF] Wan | Quality assessment and validation of the MODIS global land surface temperature[END_REF]. The MODIS LST product was validated with in situ temperature measurements recorded at various places and under various surface and atmospheric conditions [START_REF] Friedl | MODIS Collection 5 global land cover : Algorithm refinements and characterization of new datasets[END_REF]. According to [START_REF] Friedl | MODIS Collection 5 global land cover : Algorithm refinements and characterization of new datasets[END_REF] the MODIS LST accuracy is better than 1 Kelvin. The LST data has been converted to degrees Celsius. As for the MODIS NDVI data, noisy pixels affected by clouds or other atmospheric disturbances were removed when temperatures were below 0°C and the neighboring values in the time series have been linearly interpolated.

The MODIS Land Cover Type product (MCD12Q1)

The MODIS LCP (MCD12Q1, version 51) contains the International Geosphere Biosphere Program's (IGPB) classification, describing 17 land cover classes on a yearly basis at a spatial resolution of 500-m [START_REF] Friedl | Global land cover mapping from MODIS: Algorithms and early results[END_REF], [START_REF] Sultan | Assessing climate change impacts on sorghum and millet yields in the Sudanian and Sahelian savannas of West Africa[END_REF].Two classes are related to agriculture: cropland (class number 12) and cropland/natural vegetation mosaic (class number 14). Assuming that cultivated land cover area did not vary considerably during the 10-year period of study, only "consistent" pixels (i.e. pixels classified as cropland for more than six years between 2001 and 2010) were kept as cropland and the rest masked out. This crop mask was tested against a land cover map based on Landsat images in 2013 and displayed a user accuracy of 73% and a producer accuracy of 50% for the crop classes (not shown here). Because of its availability at a regional scale, we chose to conduct the analysis with the MODIS LCP to ensure the reproducibility of the methodology elsewhere. In this study, we considered that the resulting cropland was approximately equivalent to the pearl millet cultivated area (since pearl millet represents over 70% of the total agricultural production in the study area; [START_REF] Traoré | Characterizing and modeling the diversity of cropping situations under climatic constraints in West Africa[END_REF], [START_REF] Maselli | Processing of GAC NDVI data for yield forecasting in the Sahelian region[END_REF]).

Climate data

A set of daily rainfall data recorded throughout the period 2000-2010 at 28 rain-gauges (corresponding to 28 villages) distributed across the study area (Fig. 2c) was used. This dataset was provided by the AMMA-CATCH observing system. Other weather data including daily minimum and maximum air temperature, wind speed, solar radiation and minimum and maximum air relative humidity measurements were obtained from a weather station located south of Niamey (Fig. 2c). According to [START_REF] Baron | From GCM grid cell to agricultural plot: scale issues affecting modelling of climate impact[END_REF],

the variability of other meteorological data is very low compared to rainfall in this area, such that only one weather station was considered necessary.

Agricultural statistics

Agricultural statistics were used in the validation process of the remotely sensed-based yield model.

Pearl millet yield data, collected from ground surveys of major staple crops in Niger, was used. These ground surveys are conducted every year by the Niger Agricultural Statistics service at department level and were therefore available for the 2000 and 2010 period. In this study, yield data for the Kollo department and the four surrounding departments were considered (Fig. 2a).

Methods

Crop model simulations

The SARRA-H crop model

SARRA-H V3.3 [START_REF] Dingkuhn | Decision support tools for smallholder agriculture in Sub-Saharan Africa : A practical guide[END_REF], [START_REF] Baron | From GCM grid cell to agricultural plot: scale issues affecting modelling of climate impact[END_REF] was used in this study to simulate attainable pearl millet yields under climatic constraint in the NSD site at village level. This model is particularly suited for the analysis of climate impacts on cereal growth and yield in dry environments. It is currently used by AGRHYMET for operational agro-meteorological forecasting across West Africa. It simulates attainable yield under waterlimited conditions taking into account potential and actual evapotranspiration, phenology, potential and water-limited assimilation and biomass partitioning (for more details about the SARRA-H crop model see, http://sarra-h.teledetection.fr). The crop model SARRA-H has been calibrated and validated for local photoperiod sensitive pearl millet cultivars using ground surveys conducted in various location across West Africa such as in Senegal, Burkina Faso, Mali or Niger [START_REF] Traoré | Characterizing and modeling the diversity of cropping situations under climatic constraints in West Africa[END_REF]. The model was found to perform well over West Africa through comparison with FAO statistics [START_REF] Sultan | Assessing climate change impacts on sorghum and millet yields in the Sudanian and Sahelian savannas of West Africa[END_REF]- [START_REF] Ramarohetra | How satellite rainfall estimate errors may impact rainfed cereal yield simulation in West Africa[END_REF] or in comparison with other crop models in the framework of the Agricultural Model Intercomparison and Improvement Project (AgMIP, [START_REF] Bassu | How do various maize crop models vary in their responses to climate change factors?[END_REF]).

Aboveground biomass, harvest index and yield simulations

Attainable pearl millet aboveground biomass, harvest index and yield were simulated with the SARRA-H crop model for each of the 28 rainfall stations of the NSD site between 2000 and 2010, according to soil type, rainfall regime and agricultural practices (crop varieties and sowing dates). A total of 1276 simulations were conducted. The range of parameters used for the simulation was derived from previous studies and expert knowledge:

 Crop varieties: Two local pearl millet photoperiodic cultivars are found at the NSD site; Hainy Kirey (90-120 days cycle duration) and Somno (120-150 days cycle duration). These two photoperiodic varieties are particularly adapted to spatial and temporal variability of the length as well as the onset of the rainy season of the Sahelian zone [START_REF] Marteau | The onset of the rainy season and farmers' sowing strategy for pearl millet cultivation in Southwest Niger[END_REF], [START_REF] Sultan | Robust features of future climate change impacts on sorghum yields in West Africa[END_REF]. In the NSD site pearl millet HK represents among 80% of the crop [START_REF] Ramarohetra | How satellite rainfall estimate errors may impact rainfed cereal yield simulation in West Africa[END_REF], [START_REF] Bassu | How do various maize crop models vary in their responses to climate change factors?[END_REF]. Pearl millet aboveground biomass, harvest index and yields were simulated considering neither fertilization nor irrigation.

 Sowing dates: In Sahelian regions, farmer's agricultural practices choice is highly determined by the climatic constraints. Farmers generally start sowing photoperiodic millet varieties as soon as possible after the first significant rain, to benefit from the flush of available nitrogen associated with early rains, in spite of a high risk of failure and subsequent need of re-sowing [START_REF] Sultan | Robust features of future climate change impacts on sorghum yields in West Africa[END_REF], [START_REF] Kouressy | Adaptation to diverse semiarid environments of sorghum genotypes having different plant type and sensitivity to photoperiod[END_REF]. In the model, the beginning of the time window considered for the search of the satisfying conditions for sowing was set on the 1 st May, and the sowing date was automatically generated by the model as the day when simulated soil water available for the plant is greater than 10 mm at the end of the day.

 Soil type: According to the Harmonized World Soil Database [START_REF] Bezançon | Changes in the diversity and geographic distribution of cultivated millet (Pennisetum glaucum (L.) R. Br.) and sorghum (Sorghum bicolor (L.) Moench) varieties in Niger between 1976 and 2003[END_REF], 75% of soils in the NSD site are sandy and 25% are sandy clay loam. Since there is no existing data presenting the proportion of each soil type in each of the NSD site's villages respectively, we assumed the proportions proposed for the whole NSD site as being equivalent to the proportion in each village. Yields, aboveground biomass and harvest index were simulated for these two types of soils, weighted according to these proportions and considering two rooting depths (600 mm and 1800 mm) per type of soil.

An example of the aboveground biomass output obtained for Torodi village in 2008 is presented in Fig. 3.

Relationships between crop model simulations and remote sensing indices for pearl millet

Processing of remote sensing indices

MODIS NDVI time series

In Niger cultivated areas are principally gathered around villages within a distance of less than 10-km [START_REF] Maselli | Processing of GAC NDVI data for yield forecasting in the Sahelian region[END_REF]. To compare the NDVI with simulated aboveground biomass, NDVI median values within a square of 10-km by 10-km (corresponding to 1600 MODIS pixels) around each village were extracted in order to limit the analysis to areas with the higher density of crop surfaces. The median value was used to represent the average situation while minimizing the effect of pixels with a significant proportion of natural vegetation as can be expected when working with a broad-scale crop mask. Mean values were also not appropriate because the NDVI values were found not to be normally distributed. In this study three NDVI time integrals (cumulative values) were defined (Fig. 4):

 The rainy season (NDVI_RS) extends from its onset to its retreat. In order to take into consideration the spatial and temporal variability of the length of the rainy season, the onset and retreat of the rainy season was computed for each year and for each village of the NSD site following Sivakumar's definition [START_REF] Roudier | An exante evaluation of the use of seasonal climate forecasts for millet growers in SW Niger[END_REF].

 The growing period (NDVI_GP) extends from the onset of the rainy season to the end of September (Fig. 4). The end of the crop growing period corresponds approximatively to the harvesting period which was fixed here to the end of September (270 th day of the year) since it generally occurs during September [START_REF] Vaksmann | Adéquation entre risque climatique et choix variétal du mil : Cas de la zone de Bankass au Mali[END_REF].

 The productive period (NDVI_PP) of the crop growing period corresponds to phenological stages, such as the reproductive or the first maturing stages, which are especially sensitive to water stress.

Consequently, yield loss becomes significant under water stress conditions during these droughtsensitive stages. The NDVI between the beginning of August (213 th day of the year) and the end of September (including the reproductive and the maturation phases as well as the harvesting period) were used to calculate the NDVI integral during the crop productive period.

The cropland extent obtained from the MODIS Land Cover Product was used to keep only cropland classified pixels in the NDVI integral calculation which allows minimizing the influence of natural vegetation signals.

The Crop Water Stress Index (CWSI)

The CWSI, commonly used as a plant stress detection index, is originally based on canopy-air temperature difference and their relation to air vapor pressure deficit. It ranges from 0 (ample water) to 1 (maximum stress) [START_REF] Jackson | Canopy temperature as a crop water stress indicator[END_REF]. [START_REF]Harmonized World Soil Database (version 1.2)[END_REF] suggest an equivalent approach based only on canopy-air temperature differences. The CWSI used in this study can be expressed as:

(

where is the canopy-temperature from MODIS LST data and is the air temperature measurement from the meteorological AGRHYMET station. Subscripts , and refer to the minimum (nonstressed crop), maximum (cover no longer transpiring), and observed canopy-air temperature differences respectively, computed for each date within the crop mask over the study area. Since the HI of pearl millet is more sensitive to water stress during the crop productive period of the growing season [START_REF] Misra | Assimilate partitioning in pearl millet (Pennisetum glaucum L.R.Br.)[END_REF], an integral of CWSI was calculated over the productive period as defined previously (CWSI_PP).

Model development for aboveground biomass, HI and yield estimation

In Niger, pearl millet is characterized by a LAI generally lower than 2, which suggests that the relationships between NDVI and LAI are below the saturation level explained in [START_REF] Sivakumar | Predicting rainy season potential from the onset of rains in Southern Sahelian and Sudanian climatic zones of West Africa[END_REF]. The relationship between simulated aboveground biomass and each of the three NDVI time integrals was modeled with an Ordinary Least Square regression (OLS) through the following expression:

(

where represents the simulated aboveground biomass in year t and village n with the crop model SARRA-H, is the NDVI variable for the same year and village, are the parameters to be estimated and is the error term. An OLS was run at village level for the three NDVI time integrals.

As for the aboveground biomass estimation, an OLS regression was applied to derive HI from the CWSI, while the crop model output was used to calibrate the remote sensed based model:

where represents the simulated HI in year t and village n with the crop model SARRA-H, is the CWSI variable for the same year and village, are the parameters to be estimated and is the error term.

The basic equation to estimate yield, is:

Thus, by replacing each term of Eq. ( 4) by Eq. ( 2) and Eq. ( 3), the following model for yield estimation can be derived (Eq.5):

(5)

Results

Crop model simulation results

The crop model SARRA-H was run for the 28 villages of the NSD site for a period of 11 years (from 2000 to 2010). In these simulations, the mean annual simulated yields at village scale vary from 100 kg ha -1 to 1400 kg ha -1 (not shown). The yields are in the same order of magnitude that the ones measured by CIRAD (French agricultural research center for development) and AGRHYMET in the NSD site between 2004 and 2008 (400 to 1100 kg ha -1 ; [START_REF]Farmers yield variability assessment and validation of crop model to predict 'average regional' farmers yield for the main cropped varieties of millet, sorghum and maize[END_REF]). The temporal and spatial variability of the outputs of the simulation protocol are presented in Table 1 and Table 2 respectively. Table 1 shows a general high temporal variability of simulated pearl millet aboveground biomass for the 28 villages with a coefficient of variation (CV) ranging from 31% for Gorou Goussa to 63% for Kollo. Compared to the high year-toyear variability of the aboveground biomass, the temporal variability of the simulated yields (CV ranged from 19% to 46% between 2000 and 2010) and harvest indices (CV below 40% and mean HI = 0.29) are moderate. Given the size of the study area, the aboveground biomass, the HI and the yield's spatial variability could be considered relatively high (CV between 9% and 59% Table 2). The years 2000, 2002,

2007 and 2010 are those showing the highest spatial variability between the villages (e.g. 30%, 36%, 30% and 52% respectively, for simulated yields). The analysis of the crop model output during these years (not shown) reveals high water stress conditions at the beginning of the growing period (during the vegetative stage), affecting locally some of the villages and resulting in very low simulated aboveground biomass and yields for those years. We also validate the SARRA-H crop model against agricultural statistics by averaging simulated yield at the NSD site level (Fig. 5). The yields are overestimated which is one of the main drawbacks of many crop models since they simulate potential yields limited by water supply which 324 could be different from the actual yields attained in the field [START_REF] Challinor | Design and optimisation of a large-area process-based model for annual crops[END_REF]. 325

Biomass estimation based on NDVI data 4.2.1. Results at village scale

In order to support the choice of using median values to extract NDVI around each villages, different descriptive statistics have been extracted for each of the NDVI-integrated variables in order to determine the best NDVI-integrated x descriptive statistics combination for aboveground biomass estimation: the median value, the maximum value, the range (the difference between the maximum and the minimum)

and the standard deviation. The results are illustrated in Table 3. The four descriptive statistics x the three NDVI-integrated variables were compared to the simulated aboveground biomass using an OLS regression (Table 3). For all the combinations tested the correlation coefficients are low (below 0.40 but all highly significant). The Root Mean Square Errors (RMSE) is high with an RMSE equal to 989 kg ha -1 (RRMSE=42%) for the best combination (NDVI median x NDVI_PP), and an RMSE equal to 1060 kg ha - 1 (RRMSE=46%) for the less performing combination (NDVI range x NDVI_RS). Fig. 6 shows the resulting scatterplot of NDVI_PP versus simulated aboveground biomass. The dispersion of the points along the regression lines suggests the low ability of MODIS NDVI to reveal spatial and temporal aboveground biomass variability at a village scale. According to the Table 3, the best results are observed for the median NDVI values extracted around villages, thus only the combination NDVI median x NDVIintegrated variables are considered in the remainder of the study.

Results at the NSD site scale

Since neither NDVI observations nor simulated aboveground biomass follow normal distributions, median values were preferred to mean values to compare NDVI and simulated aboveground biomass over the 28 villages (NSD site scale). The aggregated NDVI value at the NSD site scale was computed considering the NDVI median value for all cropped pixels of the 28 villages. Fig. 7 shows that overall NDVI observations represent well the magnitude of the simulated aboveground biomass variability (Fig. 7a, 7b and 7c) as well as the global trends and extreme events (Fig. 7d). Among the three NDVI variables, the NDVI_PP presents the best indicator of pear millet aboveground biomass with a correlation coefficient 0.60 (significant at 10%) and a RMSE of 654 kg ha -1 which is equivalent to a RRMSE of 28%, whereas NDVI_RS appears to be the less reliable indicator (Fig 7c and Fig. 7a, respectively). The year-to-year variability is correctly displayed, with a positive trend between 2000 and 2005, a negative trend between 2005 and 2010, and NDVI observations differing from simulated aboveground biomass by less than one standard deviation (Fig. 7d). These NDVI trends coincide with the observed rainfall anomalies at the NSD site scale (Fig. 2b). At the site scale, the remote sensing based model for aboveground biomass estimation is expressed as follows:

20 (6)
where is the production of pearl millet aboveground biomass estimated at the harvest period in kg ha -1 , and _PP is the NDVI integral during the productive period at the NSD site scale.

Harvest index estimation based of LST data

Since aboveground biomass is estimated at the NSD site scale, the model for HI estimation was developed at this same scale by taking the median value of the CWSI_PP derived from the LST data, and integrated over the crop productive period. The resulting model is presented in Fig. 8, which shows that the HI and the CWSI_PP are linearly and negatively correlated, with a correlation coefficient of -0.68 (significant at 5%) and a RMSE of 0.07 (Fig. 8a). This relationship may be explain by a new biomass production allocated to grain decreasing as crop water stress increases, leading to a consequent decrease in yield. In order to better visualize the year-to-year variability of both simulated HI and CWSI_PP, we have plotted the (1-CWSI_PP) value (Fig. 8b). The year-to-year variability is generally well represented by the CWSI_PP except for 2005. The model derived for the HI estimation is expressed as follows: [START_REF] White | A continental responses phenology model climatic for monitoring variability vegetation to interannual[END_REF] where is the estimated harvest index and is the Crop Water Stress Index's integrated over the productive period at the NSD site scale.

Yield estimation based on NDVI and LST data and evaluation

Pearl millet yields at the NSD site scale were obtained by multiplying the estimated aboveground biomass (Eq.6; Fig. 9a) by the estimated HI (Eq.7; Fig. 9b). The estimated yields vary from 390 kg ha -1 to 1294 kg ha -1 (Fig. 9c). The estimated yields show an overall stable trend between 2000 and 2010 and a decline between 2005 and 2007 (Fig. 9c). The predictive capacity of the remote-sensing-based model for pearl millet yield estimation is shown in Fig. 10. The combined model based on NDVI and LST data is first evaluated by comparing simulated crop yield (from SARRA-H) to estimates based on the remote sensing-based model (Fig. 10a).

The combined model is in moderate agreement with simulated yields (r=0.50, RMSE of 219 kg ha -1 and a Mean Signed Difference [MSD] of 74 kg ha -1 ; Fig. 10a).

In order to show the contribution of thermal indices in crop yield estimation, we compared the results with the estimated yields based only on NDVI data. The model (based on both NDVI and LST data) results are in good agreement with the official yield statistics (r=0.82 significant at 5%, Fig. 10b).

Furthermore, the combination of NDVI and LST data clearly contributes to improve yield estimation compared to NDVI data alone (r=0.59, Fig. 10c). However, like the crop model used for the calibration, the remote sensing-based models clearly overestimate yields (Fig. 10b and Fig. 10c) which leads us to consider the ability of these models to render the yield's year-to-year variability observed by the agricultural statistics. To do so, both estimated and observed yields were normalized. For each year, the absolute differences between agricultural statistics z-score values and those of the models were computed (Fig. 11a). In order to provide an overall indication on the performance of each of the models, the sum of the absolute differences is also assessed. Yield's year-to-year variability from 2000 to 2010 is quite well rendered in both models in Fig. 11a, particularly for the second half of the period (between 2005 and 2010). The combined model based on NDVI and LST data is the closest of the agricultural statistics temporal profile (absolute difference sum = 5.61), particularly in extreme dry years such as in 2000 (Fig. 11b). Nevertheless, the overall trend is also well transcribed, split in a stable period between 2000 and 2005, followed by a decrease trend in yields between 2005 and 2010 (Fig. 11b).

To test the robustness of the remote sensing-based model, yields for the four surrounding departments were computed and compared with the corresponding official yield statistics (Table 4).

Overall, computed yields coincide with the yield statistics, with correlation coefficients above 0.50 (significant at 10%) for 3 departments (Table 4). As for the NSD site the remote-sensing based model systematically overestimates yields (RMSE ranging from 237 kg ha -1 to 742 kg ha -1 ).

Discussion

Aboveground biomass estimation based on NDVI time series

The first stage of the remote sensing-based model consisted in developing an empirical relationship between NDVI time series and pearl millet aboveground biomass simulated by the crop model SARRA-H.

The study first highlighted that the ability of the MODIS NDVI time series to estimate aboveground biomass depends on the scale considered. At the village scale (considering the whole dataset:

28 villages, 11 years) , the study found out that the MODIS NDVI time series are not able to reveal both the spatial and temporal variability of the simulated aboveground biomass (RRMSE > 40%; Table 3 and Fig. 6). As previously shown by [START_REF] Sultan | Agricultural impacts of large-scale variability of the West African monsoon[END_REF], in the semi-arid zone of Niamey, aboveground biomass and final yields are mainly influenced by the spatio-temporal distribution of rainfall, and so a high variability of aboveground biomass can be observed between villages which are only a few kilometers apart. Thus, the low correlation between NDVI and aboveground biomass at the village scale implies that the spatial variability of NDVI is not as strongly associated with the spatial variability of rainfall. Further analyses are required on other potential factors that could influence NDVI at this scale. We could assume for instance that, in semi-arid regions where vegetation cover is relatively sparse, soil may cause high variations in the NDVI values at such a small scale, causing NDVI values artifacts [START_REF] Huete | Investigation of soil influences in AVHRR red and near-infrared vegetation index imagery[END_REF] and therefore reducing the correlation between NDVI and aboveground biomass. [START_REF] Rasmussen | Operational yield forecast using AVHRR NDVI data: Reduction of environmental and inter-annual variability[END_REF], considering a direct relation between NDVI and yield, found that including soil information improved yield prediction in the Peanut Basin in Senegal. On the other hand, at the NSD site scale (temporal analysis), a good correlation was found between simulated aboveground biomass and NDVI_PP (r=0.60). This improvement could be explained by [START_REF] Marteau | The onset of the rainy season and farmers' sowing strategy for pearl millet cultivation in Southwest Niger[END_REF] the reduction of the noise in the NDVI time series when aggregating at a coarser level and (2) a better representativeness of the overall crop growth conditions over the NSD site that is mainly driven by rainfall variability.

The capacity of the MODIS NDVI time series to estimate aboveground biomass depends also on the time period used for the integration. On that point our results are different from [START_REF] Tucker | Satellite Remote Sensing of Total Herbaceous Biomass Production in the Senegalese Sahel : 1980-1984[END_REF], [START_REF] Diouf | Monitoring land-cover changes in semi-arid regions: remote sensing data and field observations in the Ferlo, Senegal[END_REF] who found a good correlation between NDVI integrated over the whole growing season and aboveground biomass in Senegal. In these studies, only natural herbaceous vegetation was considered, for which final aboveground biomass is not much different from vegetative biomass, thus justifying NDVI integration over the entire length of the growing season. Our study focuses on a final aboveground biomass that depends on both vegetative biomass and grains. NDVI values were integrated over the crop productive period to account for grains, since it corresponds to the reproductive period and maturation phases, which include grain filling when plants reach their maximum development [START_REF] Maselli | Processing of GAC NDVI data for yield forecasting in the Sahelian region[END_REF]. Our results corroborate other studies that directly relate NDVI to yields such as [START_REF] Tucker | Relationship of spectral data to grain yield variation[END_REF] who found that the strongest correlation of NDVI with wheat yields is achieved when taking into consideration NDVI values around their maximum which includes the sensitive stages of grain production. [START_REF] Rasmussen | Assessment of millet yields and production in northern Burkina Faso using integrated NDVI from the AVHRR[END_REF] then tested the influence of different NDVI integration periods and found a coefficient of determination R²=0.50 (i.e. r=0.70) for the productive period. In another analysis, using NOAA AVHRR data between 1982 and 1990 for Niger, [START_REF] Maselli | Processing of GAC NDVI data for yield forecasting in the Sahelian region[END_REF] concluded that the best time integration period for millet and sorghum yield assessment is from August to September. Finally, more recently in a study conducted in China [START_REF] Huang | Analysis of NDVI Data for Crop Identification and Yield Estimation[END_REF] it was also found that the productive and maturing stages including the heading, flowering and filling of the crops are the best suitable periods for yield estimation of paddy rice, corn and winter wheat due to the stress sensitivity of these periods that would lead to biomass reduction and thus potentially yield losses. In our study, NDVI_RS and NDVI_GP (both determined by the onset of the rainy season) appear to be less correlated to aboveground biomass. A potential explanation for this could be the delay between the NDVI onset of the growing season and the calculated start-of-season, which occur one month apart, as previously shown by [START_REF] Vintrou | A Comparative Study on Satellite and Model-Based Crop Phenology in West Africa[END_REF]. At the beginning of the growing season in a MODIS pixel the proportion of the millet cover is probably lower than the proportion of the surrounding natural vegetation. The latter reacts immediately to the first significant rainfall, whereas crops are sown later, when sufficient water (>10 mm) is available in the soil [START_REF] Marteau | The onset of the rainy season and farmers' sowing strategy for pearl millet cultivation in Southwest Niger[END_REF] and have a growth rate lower than natural vegetation.

On the year-to-year variability analysis, a decrease in both the simulated aboveground biomass and the NDVI was observed from 2005 to 2010, with an important decline between 2005 and 2006 (Fig. 7d.). When comparing this result with annual rainfall anomalies (Fig. 2b) and it can be concluded that both aboveground biomass and NDVI follow the major trends of rainfall anomalies (as seems particularly evident between 2005 and 2006). This comes in support of the previous assumption that rainfall remains the main determinant of NDVI variability at the NSD site scale.

Harvest index estimation based on an indicator of crop water stress: the CWSI

For most crop models, including SARRA-H, DSSAT and CROPWAT [START_REF] Dingkuhn | Decision support tools for smallholder agriculture in Sub-Saharan Africa : A practical guide[END_REF], [START_REF] Jones | The DSSAT cropping system model[END_REF], [START_REF] Fao | A computer program for irrigation planning and management[END_REF], water stress during the reproductive and maturation phases is considered a crop yield limiting factor. In the remote-sensing model, we take into account the crop water stress effect on yield through the use of the CWSI, an indicator based on LST. To our knowledge, it is the first time that a link is sought between an indicator of crop stress and HI. An overall good correlation (r=-0.68) was found between HI and CWSI_PP at the NSD site scale, meaning that the HI decreased linearly as the water supply became more limited for plants. However, as for the use of vegetation indices in semi-arid zones, the main issue with thermal indices based on canopy temperature is the spatial heterogeneity due to the soil influence when the canopy does not completely cover the ground. Because bare soil is often much warmer than the air, the soil background temperature included in the LST can lead to false detections of crop water stress [START_REF] Hatfield | Agriculture and Remote Sensing[END_REF]. To overcome this limitation, a possibility may be to use the Water Deficit Index developed by [START_REF] Moran | Estimating Crop Water Deficit Using the Relation between Surface-Air Temperature and Spectral Vegetation Index[END_REF], which considers both the difference between air and surface temperatures and the fraction of crop cover derived from vegetation indices, to estimate the water status. This method was not tested in this study, as some adaptations are ongoing to test the construction of the vegetation indextemperature trapezoid from satellite time series.

Estimation of pearl millet yields

The two previous approaches for aboveground biomass and HI estimation were combined into a simple, robust and timely satellite-based model of rainfed cereal yield, applicable at the department level.

If in absolute values, yields are overestimated compared to official agricultural statistics of the Kollo department, the analysis of the standardized values has shown a good agreement in terms of year-to-year variability reproduction, translating into a high correlation with statistics. In their recent meta-analysis [START_REF] Huang | Meta-analysis of influential factors on crop yield estimation by remote sensing[END_REF] found that for four studies conducted in Senegal, Burkina Faso and Niger using NOAA AVHRR data, the correlation coefficients between NDVI alone and millet yield were comprised between 0.75 and 0.94 which is comparable to the present work (r=0.82). However, caution in the interpretations has to be taken particularly because [START_REF] Marteau | The onset of the rainy season and farmers' sowing strategy for pearl millet cultivation in Southwest Niger[END_REF] although the size of the study area considered in these studies is similar to that of the present study (i.e. results aggregated at a department level), the time period considered was much shorter (2 years in [START_REF] Rasmussen | Assessment of millet yields and production in northern Burkina Faso using integrated NDVI from the AVHRR[END_REF]) and (2) when the time period considered is comparable to ours, results were aggregated at higher administrative levels than for us (several departments or country level; e.g. [START_REF] Rasmussen | Operational yield forecast using AVHRR NDVI data: Reduction of environmental and inter-annual variability[END_REF], [START_REF] Groten | NDVI-crop monitoring and early yield assessment of Burkina Faso[END_REF]).

The comparison with a model based only on NDVI has highlighted the usefulness of combining vegetation and thermal indices (NDVI and CWSI) for yield estimation. The ability to render the year-toyear variability of pearl millet yield was clearly improved through this combination, with a correlation coefficient increasing from 0.59 to 0.82 and the z-score absolute difference sum decreasing from 7.28 to 6.21. Indeed, because of the spatial variability of management practices, soil water capacity or nitrogen availability, different yields could be observed for the same amount of biomass. In addition, events such as droughts during the reproductive stage, with potentially drastic yield reduction but negligible effects on vegetative biomass, are certainly poorly detected by a model based only on vegetation indices. Thus the direct relation NDVI/yield mostly allows assessing potential harvestable yields when assuming nonlimiting conditions (i.e. when yield is proportional to aboveground biomass). These potential yields could the localization of cropland in the MODIS LCP over West Africa and showed a high spatial variability with user accuracy varying between 17% to 70% according to farming systems. Thus to extrapolate our methods in other locations, further efforts are needed to develop at least a map locating cultivated zones and if possible the main crop type at a regional scale.

The use of a crop model instead of ground measurements to calibrate the remote-sensing model can also be questioned. SARRA-H as most crop models tends to overestimate yields (Fig. 5) since it simulates attainable yields according to agro-meteorological constraints but does not integrate all biotic (e.g. birds, pests, and diseases due to excess moisture) or other non-environmental factors that influence crop management which can lead to yield variations [START_REF] Prasad | Crop yield estimation model for Iowa using remote sensing and surface parameters[END_REF], [START_REF] Sultan | Agricultural impacts of large-scale variability of the West African monsoon[END_REF], [START_REF] Lobell | Global scale climatecrop yield relationships and the impacts of recent warming[END_REF] . Remote sensing indices do integrate biotic and non-environmental factors, and because they are calibrated using crop model outputs, an overestimation of yields by the remote sensing-based model could be expected. In addition, since the simulated yields from SARRA-H are overestimated, does that mean that the aboveground biomass and the harvest index are also overestimated? For the latter, the simulated HI as well as the estimated HI are within the range of those measured by [START_REF] Bégué | Radiation use efficiency of pearl millet in the Sahelian zone[END_REF] over 168 pearl millet plots in the Niamey area. Authors found a mean HI of 0.22, we found a mean simulated and estimated HI of 0.29. For aboveground biomass, reliable measurements in on-farm situations are not available. However, under controlled conditions it has been shown in [START_REF] Sultan | Agricultural impacts of large-scale variability of the West African monsoon[END_REF] in Senegal for pearl millet and in [START_REF] Kouressy | Adaptation to diverse semiarid environments of sorghum genotypes having different plant type and sensitivity to photoperiod[END_REF] for Sorghum in Mali that the aboveground biomass (both yields and growth dynamics) were well simulated by SARRA-H. The same conclusion can be drawn from the study of [START_REF] Traoré | Characterizing and modeling the diversity of cropping situations under climatic constraints in West Africa[END_REF] based on on-farms survey near Niamey (Niger) for pearl millet.

Nonetheless, beyond the yields overestimation, our study show that the year-to-year variability is quiet well simulated by the remote sensing based model.

Remote sensing indices also present intrinsic limitations. Despite the fact that a filter was applied to reduce noise in the NDVI and LST time series, the presence of clouds, aerosols or dust residues may lead to noise and the downgrading of data quality [START_REF] Hird | Noise reduction of NDVI time series : An empirical comparison of selected techniques[END_REF]. Thus, the poor performance of the NDVI/aboveground biomass relation at local scale may also be explained by the 250m x 250m pixel size of MODIS images that integrates a mixture of elements (crops, natural vegetation, bare soils) particularly in the semi-arid region with low and sparse vegetation and where crop fields are often smaller than the pixel size.

Finally, our study is limited to a period of eleven years and to 28 sites due to the unavailability of more climatic data from ground observations to run the crop model. Agro-meteorological variables derived from satellite could also be considered as an alternative. However, the correct estimation of these variables from satellite, especially rainfall, remains an open issue. For instance, [START_REF] Leroux | Testing satellite rainfall estimates time series for crop yield simulation of a rainfed cereal in West Africa[END_REF] found in the same area that the TRMM 3B42 product, which delivers rainfall estimates at a daily time step, was not able to accurately detect rainfall temporal pattern at the station level, and particularly the intra-seasonal rainfall distribution. We hope that in a few years, the statistical relationships between aboveground biomass and NDVI, and between HI and CWSI, can be updated and made more robust when more climatic data are available.

Conclusion and perspectives

The difficulty to access ground measurements in West Africa and to estimate yields over large areas using other monitoring methods such as agrometeorological modelling makes remote sensing observation a good alternative or addition to consider for early warning systems. In this study, we investigated a new approach based on the combination of vegetation and thermal indices for rainfed cereal yield assessment in the Sahelian region. Empirical statistical models were developed between remotesensing indices (MODIS NDVI and LST), and SARRA-H simulated aboveground biomass and harvest index respectively, and combined for the assessment of crop yield. We demonstrated that the combined model performed better than the one using vegetation index alone. The inclusion of LST improves yield estimations by accounting for the harvest index which is an indicator of the proportion of total aboveground biomass really transformed into grains. In addition, it allows using NDVI as an estimator of aboveground biomass, which is its primary function, rather than an indirect estimator of yield.

Fig. 1

 1 Fig. 1 summarizes the overall methodology. Empirical statistical relationships are sought (1) between

Figure 1 :

 1 Figure 1: Flowchart of the approach adopted corresponding to the three stage of the remote sensed based model development.

Figure 2 :

 2 Figure 2: The NSD site (Niger) : a) location of the NSD site (red square) and the five departments considered in the study; b) anomalies of annual rainfall (deviation from the mean between 2000 and 2010); c) mean NDVI annual amplitude between 2000 and 2010, and location of the 28 rain gauge (red circles) and meteorological (star) stations.

Figure 3 :

 3 Figure 3: Example of Somno millet simulated aboveground biomass with the crop model SARRA-H for the village of Torodi in 2008. The dark curves represent aboveground biomass for a sandy clay loam soil (thin soil in dashed line and deep soil in solid line), the gray curves represent a sandy soil (thin soil in dashed line and deep soil in solid line) and red curve represent the resulting weighted average.

Figure 4 :

 4 Figure 4: The three NDVI time integrated periods, example of the Kollo village in 2005. NDVI profile (black line) during the season is compared to the simulated aboveground biomass (red line).

Figure 5 :

 5 Figure 5: Observed pearl millet yields from agricultural statistics for the department of Kollo vs simulated yields obtained with SARRA-H aggregated at the NSD site level.

Figure 6 : 1 )

 61 Figure 6: Scatterplot of the simulated aboveground biomass (kg ha -1 ) and the NDVI integrated over the productive period for the 28 villages of the NSD site and over the 11 years of data. The RMSE of the aboveground biomass is 989 kg ha -1 which is equivalent to a RRMSE of 42%, and the correlation coefficient is 0.38. The solid line is the linear regression line and the blue area is the confidence interval for pvalue<0.1.

Figure 7 :

 7 Figure 7: SARRA-H simulated aboveground biomass (kg ha -1 ) vs a) MODIS NDVI integrated during the rainy season, b) MODIS NDVI integrated during the growing season, and c) MODIS NDVI integrated during the productive period. The regression line is in black solid line. d) Comparison of the interannual variability of simulated aboveground biomass and NDVI observations, expressed in z-score values. The grey area is the ± standard deviation computed from simulated aboveground biomass.

Figure 8 :

 8 Figure 8: a) SARRA-H simulated harvest index vs CWSI_PP estimated from MODIS LST data, over the 2000-2010 period (the regression line is in black solid line); b) comparison of the interannual variability of SARRA-H simulated harvest index and (1-CWSI_PP) values, expressed in z-score values. The grey area is the ± standard deviation computed from simulated harvest index.

Figure 9 :

 9 Figure 9: Evolution of a) the aboveground biomass estimated from the MODIS NDVI model (Eq. 6), b) the harvest index estimated from the MODIS-derived CWSI model (Eq. 7), and c) the resulting pearl millet yield derived from the combination of Eq. 6 and Eq. 7, over the study site.

Figure 10 :

 10 Figure 10: a) Simulated yields from SARRA-H vs estimated yields from the combination of NDVI and LST data, and vs estimated yield from remote sensing with (b) or without (c) LST data. The 1:1 line is given in grey dashed line.

Figure 11 :

 11 Figure 11: a) Year-to-year yield variability (SARRA-H, NDVI data, NDVI x LST data) comparison with agricultural statistics. The y-axis indicates the absolute difference between yields anomalies (expressed in z-score) estimated and yield anomalies from agricultural statistics. In brackets are specified the sum of absolute differences. b) Agricultural statistics and simulated yields (NDVI x LST data) standardized anomalies (in z-score).

Table 1 : Temporal variability of simulated aboveground biomass, harvest index (HI) and yield. The mean values and the 326 coefficients of variation (CV) are calculated on the 2000-2010 period, and are given for each village. In bold, the values 327 averaged of means and CV over the dataset are given.

 1

	328
		Aboveground biomass	Harvest Index	Yield
		Mean (kg ha -1 ) CV (%) Mean CV (%) Mean (kg ha -1 ) CV (%)
	Alkama	2063	46	0.29	27	813	33
	Banizoumbou	2285	51	0.27	26	768	34
	Beri Koira	2290	48	0.31	33	911	19
	Berkiawel	2365	52	0.30	30	920	32
	Bololadie	2138	58	0.29	28	789	46
	Boubon Golf	2387	44	0.31	31	972	17
	Darey	2012	40	0.32	21	914	22
	Debere Gati	2381	55	0.29	24	888	41
	Fandou Beri	2001	43	0.31	22	903	33
	Gardamakora	2066	51	0.30	33	808	36
	Gorou Goussa	2653	31	0.26	27	956	18
	Guilahel	2416	52	0.28	35	855	30
	Harikanassou	2732	33	0.27	22	1033	9
	IH Jachere	2254	49	0.30	29	902	23
	Kaligorou	2349	35	0.28	27	896	25
	Kare	2318	51	0.30	25	922	30
	Kokorbe Fandou	1936	62	0.33	31	829	37
	Kollo	2074	63	0.30	36	754	41
	Koure Sud	2321	42	0.29	21	940	25
	Koyria	2350	38	0.29	26	924	18
	Massi Koubou	2155	49	0.30	35	857	34
	Niamey Aeroport	2386	52	0.30	31	892	25
	Niamey Orstom	2103	51	0.32	26	907	24
	Sandideye	2573	42	0.28	26	963	25
	Tanaberi	2302	38	0.29	22	949	25
	Torodi	3271	43	0.24	39	934	37
	Wankama	1915	49	0.32	23	844	37
	Yillade	2674	40	0.27	24	994	19
	Mean	2313	47	0.29	28	894	28
	329

Table 2 : Spatial variability of simulated aboveground biomass, harvest index (HI) and yield. The mean coefficients of variation

 2

		Mean (kg ha -1 ) CV (%) Mean CV (%) Mean (kg ha-1) CV (%)
	2000	2332	24	0.22	21	719	30
	2001	2536	23	0.25	23	943	18
	2002	1501	52	0.35	15	768	36
	2003	3054	24	0.24	17	1050	14
	2004	2386	34	0.28	16	949	22
	2005	3967	23	0.21	17	1082	12
	2006	1706	26	0.36	16	911	15
	2007	1989	41	0.31	17	879	30
	2008	2781	28	0.27	21	1029	11
	2009	2365	41	0.31	24	990	19
	2010	828	59	0.43	9	518	52
	Mean	2313	34	0.29	18	894	24

Table 3 : Elements of the regression analysis obtained between the simulated aboveground biomass and the descriptive statistics x NDVI variables (NDVI integrated during the rainy season, the growing period and the productive period) obtained at the village scale for years 2000-2010.

 3

	Descriptive statistics	NDVI Variables	Intercept Slope	r	p-value	RMSE (kg ha-1)	RRMSE (%)
		NDVI_RS	336	0.07 0.32 6.08E-09	1012	43
	Median	NDVI_GP	255	0.10 0.34 1.20E-09	1006	43
		NDVI_PP	-704	0.24 0.38 5.80E-12	989	42
		NDVI_RS	893	0.03 0.26 5.07E-06	1033	45
	Max	NDVI_GP	437	0.06 0.33 5.11E-09	1011	43
		NDVI_PP	-99	0.13 0.31 1.86E-08	1017	44
		NDVI_RS	1886	0.01 0.13	0.02	1060	46
	Range	NDVI_GP	1585	0.04 0.21	0.0002	1046	45
		NDVI_PP	1634	0.06 0.18	0.002	1052	46
		NDVI_RS	1810	0.14 0.16	0.005	1056	45
	Standard Deviation	NDVI_GP	1389	0.34 0.24 3.07E-05	1039	45
		NDVI_PP	1280	0.59 0.21	0.0001	1045	45

Table 4 : Estimated yields from the remote-sensed based model vs the agricultural statistics yields

 4

		r	p-value RMSE (kg ha -1 )
	Fillingue 0.45	0.15	238
	Kollo	0.82	0.01	423
	Ouallam 0.23	0.48	237
	Tera	0.58	0.06	505
	Tillaberi 0.64	0.03	742

(CV) are calculated on the 28-village data set, and are for each year. In bold, the values averaged of means and CV over the

dataset are given.

Acknowledgments

Leroux L. is supported by the CIRAD and by the Centre National d'Etudes Spatiales (Project CNES-TOSCA "Dynafrique"). The authors thank the MODIS team for sharing the MODIS Vegetation Indices, Land Surface Temperature and Land Cover products. Rainfall data were available thanks to the AMMA program and especially Guillaume Quantin and Théo Vischel (Laboratoire d'Etude de Transferts en Hydrologie et Environnement, Grenoble, France). Agricultural statistics data were kindly provided by the Niger Agricultural Services. We thank Beatriz Belón for a helpful review of the manuscript. Finally, we would like to thank anonymous reviewers for their valuable comments.

Aboveground biomass

Harvest Index Yield however be reduced by crop water stress during the reproductive stages as shown in this study.

Consequently the direct relation NDVI/yield should be considered valid only for specific areas or years without major limiting factors affecting yield.

Limitations of the method

The remote sensing-based model was applied directly to four surrounding departments and the correlation coefficients were globally good despite an overall tendency to yield overestimation by the model. The four departments are situated at the North of Kollo. They are mainly dominated by agropastoral activities, with a mixture of livestock and crop cultivation [START_REF] Fews Net | LIVELIHOODS ZONING ' PLUS ' ACTIVITY IN NIGER[END_REF]. Therefore, the probability to have a mixture of crop vegetation and grasslands within a MODIS cropped pixel is high, which may explain a lower performance of the model. Moreover, in these mixed zones of pasture, the seeding rates are also very low leading to a sparse vegetation cover that causes high NDVI variations due to soil effects.

This highlights the main limitation of such models, based on empirical relationships between remote sensing indices and yields: they depend on the environmental characteristics of the study area, which restricts their application elsewhere without recalibration. In addition, such models also depend on the farming system considered. For this reason, the model we developed in this study is only valid for a system based on a single crop and should be tested or adapted for other farming systems such as in the cereal-root crop mixed system where a wide range of different cereals is grown (maize, millet, sorghum or cassava among other) including cases of intercropping.

Another consideration to take into account concerning our methodology is the need of a crop mask to isolate cropped pixels. Since a pearl millet crop type map is not available for the NSD site, a crop mask from the MODIS LCP was used here. The same approach was also applied to NDVI and CWSI values extracted from the Landsat Crop mask. A coefficient of correlation of 0.80 is obtained when the resulting estimated yields are compared to official statistics (not shown) which is close to the one obtained with the MODIS LCP product. This confirms the relevance of the approach for the NSD site. However, while the MODIS LCP has been validated for our study area, [START_REF] Leroux | How Reliable is the MODIS Land Cover Product for Crop Mapping Sub-Saharan Agricultural Landscapes?[END_REF] recently spatialized the uncertainties in Furthermore, by using a crop model validated over the study area, this study showed that the combination of satellite data with crop modelling is a good option for yield estimation and its year-to-year variability based on remote sensing, especially for areas where ground measurements, required for the calibration of the remote sensing-based model, are not available.

Our study confirms that even in small-holder agriculture such as those of the Sahelian region, the use of coarse resolution satellite information for yield monitoring is possible. As the model proposed is simple, robust and based on empirical relations with vegetation and thermal MODIS indices, there is scope for operational implementation of yield estimation at regional scale in a food security early warning system, in particular for the assessment of the year-to-year yield variability in regions with agronomic and climatologic characteristics close to those of the NSD site. In addition, such a system could provide an early estimation of yield shortly after harvest for an area equivalent to an administrative unit unlike agricultural statistics that are currently available from three to six months after harvest. But that would require addressing the issue of multi-crop type systems on which, to the best of our knowledge, no studies