Integrating Sentinel-2 Imagery with AquaCrop for Dynamic Assessment of Tomato Water Requirements in Southern Italy

: A research study was conducted in an open ﬁeld tomato crop in order to: (i) Evaluate the capability of Sentinel-2 imagery to assess tomato canopy growth and its crop water requirements; and (ii) explore the possibility to predict crop water requirements by assimilating the canopy cover estimated by Sentinel-2 imagery into AquaCrop model. The pilot area was in Campania, a region in the south west of Italy, characterized by a typical Mediterranean climate, where ﬁeld campaigns were conducted in seasons 2017 and 2018 on processing tomato. Crop water use and irrigation requirement were estimated by means of three di ﬀ erent methods: (i) The AquaCrop model; (ii) an irrigation advisory service based on Sentinel-2 imagery known as IRRISAT and (iii) assimilating the canopy cover estimated by Sentinel-2 imagery into AquaCrop model Sentinel-2 imagery proved to be e ﬀ ective for monitoring canopy growth and for predicting irrigation water requirements during mid-season stage of the crop, when the canopy is fully developed. Conversely, the integration of the Sentinel-2 imagery with a crop growth model can contribute to improve the irrigation water requirement predictions in the early and development stage of the crop, when the soil evaporation is not negligible with respect to the total evapotranspiration. tool sen2cor_v2-5-5 ). Sen2Cor tool performs the atmospheric, terrain and cirrus correction of Top-Of-Atmosphere Level 1C input data, and creates Bottom-Of-Atmosphere, optionally terrain and cirrus corrected reﬂectance images; additional, aerosol optical thickness, water vapor, scene classiﬁcation maps and quality indicators for cloud and snow probabilities. (CDC) coe ﬃ cients; full canopy crop transpiration coe ﬃ cient (Kc); biomass WP and soil water depletion thresholds. Non-conservative parameters vary depending on crop and ﬁeld management, soil type, and climate (sowing date and density, length of crop cycle and phenological stages, maximum canopy cover, etc.). Such parameters can be either retrieved from AquaCrop literature or calibrated by the user (e.g., by means of ﬁeld experiments). Main AquaCrop outputs are crop production (biomass and yield) and crop water use.


Introduction
Worldwide significant progress has been made to utilize precision agriculture for irrigation as a mean to increase water use efficiency or decrease the water footprint in irrigated agriculture [1,2]. The progress is mainly restricted to advances at the plot scale and individual systems such as installations for drip irrigation or central pivots. As well known, actual crop evapotranspiration (ET) is a major term of water budget in agriculture and it is the main variable used to determine crops water requirement. Beside their massive progress during the recent years, accurate field measurements (soil moisture, plant-based sensors, etc.) are very scarce because sensors and measuring devices are expensive, their use requires specific expertise and complex maintenance and are typically limited to experimental stations. For this reason, many attempts have been made for developing indirect ET estimation methods based on crop data at a field scale easily available across large regions. Concordantly, the use of remotely sensed data has become more common to monitoring and controlling activities at different

Satellite-Based Assessment of Crop Water Requirement: IRRISAT
The remote sensing technique employed is the one employed by IRRISAT, a satellite-based irrigation advisory service developed in Italy and operational since 2007 in Southern Italy [7]. The service aims at providing farmers and water managers with real time information on crop water needs, which are estimated by combining high resolution data from Earth Observation satellites with weather data for the calculation of crop water requirements. More recently, IRRISAT has been combined with numerical weather predictions for forecasting crop water needs up to five days in advance [19,20]. Information is delivered in near-real time (24 h) to users (farmers, Water User Association and water agencies) by means of a dedicated WebGIS accessible from PC, tablets and smartphones (https://www.irrisat.com/en/). In 2016, the operational irrigation advisory service of the Campania region has reached about 2000 farmers with a total irrigation area larger than 80,000 hectares: The achieved water saving has been estimated to be larger than 30% [7]. In 2016 the Italian Ministry of Agriculture has listed IRRISAT among the applicable methodologies for the estimation of the irrigation volume, complying with the EU Water Framework Directive.
The IRRISAT methodology is summarized in Figure 2. Crop potential evapotranspiration (ETp) is computed with the Penman-Monteith equation, with crop parameters albedo (α) and leaf area index (LAI) derived from processing Sentinel-2A/B images in the visible and infrared regions [7] while assuming fixed values for the stomatal resistance (sr ≈ 100 sm −1 ) and crop height (hc = 0.4 m) for herbaceous crop [8]. Following this approach, the calculation of ETp requires standard meteorological data (daily air temperature, relative humidity, solar radiation, wind speed and precipitation), LAI and surface α. Processing tomato (Solanum lycopersicum L.) was cultivated in years 2017 (41 • 00 26.33" N 14 • 10 13.93" E) and 2018 (41 • 01 39.60" N 14 • 10 32.87" E) in two parcels of 4 ha each. Soil samples (three samples were averaged for analysis) were examined for assessing the main physical-chemical properties, including gravel (%), soil texture, bulk density (t/m 3 ), pH and organic matter (%), at a depth of 15 cm. In 2018 soil characteristics were obtained by the soil map of Region Campania (pedological maps 1:50,000; http://agricoltura.regione.campania.it/pedologia/suoli.html). SPAW software (Soil-Plant-Air-Water; v6.02.75, United States Department of Agriculture-USDA, Washington, DC, USA) was used for assessing soil hydraulic properties, such as soil water retention, soil hydraulic conductivity, field capacity and plant available water.
Soil was ploughed at 40 cm depth, and tomato seedlings were transplanted on April 9 in continuous double rows with 33.5-40.0 cm space between plants, 50 cm between rows, 1.10-1.20 cm between double rows, with a final plant density of 32,000 and 33,500 plants/ha in year 2017 and 2018, respectively.
Meteorological daily data of maximum and minimum temperature ( • C), relative humidity (%), wind speed (m s −1 ) and precipitation (mm) were collected at a complete weather station located in the study area for the two cropping seasons.
Plants were watered by light driplines, with 30 cm dripper spacing and 2 L/h flow rate at 1 bar.
Based on the farmer best knowledge, a total volume of 130 L/plant of water, corresponding to 4160 m 3 /ha, and 120 L/plant corresponding to 4020 m 3 /ha of irrigation was given throughout the growing period in 2017 and 2018, respectively. Such volumes were then considered as reference values against which the different estimation methods (AquaCrop, IRRISAT, integration of Sentinel 2A imagery into AquaCrop model) were tested.

Satellite-Based Assessment of Crop Water Requirement: IRRISAT
The remote sensing technique employed is the one employed by IRRISAT, a satellite-based irrigation advisory service developed in Italy and operational since 2007 in Southern Italy [7]. The service aims at providing farmers and water managers with real time information on crop water needs, which are estimated by combining high resolution data from Earth Observation satellites with weather data for the calculation of crop water requirements. More recently, IRRISAT has been combined with numerical weather predictions for forecasting crop water needs up to five days in advance [19,20]. Information is delivered in near-real time (24 h) to users (farmers, Water User Association and water agencies) by means of a dedicated WebGIS accessible from PC, tablets and smartphones (https://www.irrisat.com/en/). In 2016, the operational irrigation advisory service of the Campania region has reached about 2000 farmers with a total irrigation area larger than 80,000 hectares: The achieved water saving has been estimated to be larger than 30% [7]. In 2016 the Italian Ministry of Agriculture has listed IRRISAT among the applicable methodologies for the estimation of the irrigation volume, complying with the EU Water Framework Directive.
The IRRISAT methodology is summarized in Figure 2. Crop potential evapotranspiration (ETp) is computed with the Penman-Monteith equation, with crop parameters albedo (α) and leaf area index (LAI) derived from processing Sentinel-2A/B images in the visible and infrared regions [7] while assuming fixed values for the stomatal resistance (sr ≈ 100 sm −1 ) and crop height (hc = 0.4 m) for herbaceous crop [8]. Following this approach, the calculation of ETp requires standard meteorological data (daily air temperature, relative humidity, solar radiation, wind speed and precipitation), LAI and surface α. The irrigation water requirement (IWR) is then calculated as the difference between ETp and the crop effective rainfall (Pn), according to the following equation: (1) The irrigation water requirement (IWR) is then calculated as the difference between ETp and the crop effective rainfall (Pn), according to the following equation: It is assumed that capillary rise does not contribute to root zone soil moisture in the summer season, as usually occurs in the southern European regions. Runoff and deep drainage are assumed to be negligible considering the low amount of rainfall during the two growing seasons. Pn is obtained by reducing the precipitation above canopy (P) by a quantity that depends on canopy development, according to an empirical function of the LAI and the fractional vegetation cover fc, as reported in Vuolo et al. [7].
In this study, 10 and 21 multispectral high-resolution images from Sentinel-2A and 2B have been acquired during 2017 and 2018, respectively. The multi spectral instrument (MSI) on board of Sentinel-2A/2B captures data at 10, 20 and 60 m of spatial resolution over 13 spectral bands with a very high temporal resolution of five days at the equator. Individual Sentinel-2 granules Level-1C (processed at the top-of-atmosphere reflectance) were acquired from Copernicus Open Access Hub (https://scihub.copernicus.eu/), already ortho-rectified in UTM/WGS84 (image tiles of 100 × 100 km 2 ). The information gathered by Sentinel-2 system (orbit, attitude, date accuracy and viewing directions of all detectors) are exploited for geolocating all Sentinel-2 pixels with an accuracy of about 11 m for about 97% of the cases, which is about the size of one Sentinel-2 pixel. The standard need for multi-temporal registration errors is 0.3 pixels, and the current performances show that for more than 50% of the cases, the performance does not meet that requirement. The resolution is estimated to be three times the registration error, thus the resolution Sentinel-2 time series is around 30 m.
In order to obtain homogeneous and comparable products as time series, all value-added products (LAI, α and fc) are calculated based on atmospherically corrected Level-2A data. LAI and fc are calculated by S2ToolBox [21], an artificial neural network (ANN) algorithm, trained by using radiative transfer simulations from PROSPECT [22] and SAIL [23] models, and tailored for Sentinel-2 data. The algorithm requires eight Sentinel-2 spectral bands (B3-B7, B8a, B11 and B12) at 10 and 20 m (pixel size), which are all resampled to 10 m to derive LAI and fc. Experimental studies have shown the accuracy of this approach for LAI estimation in different environments and crops [24,25]. In this study, average and variance of LAI and fc at parcel scale were assessed by taking a minimum of 50 pixels falling within each parcel, after excluding pixels affected by boundary effects or cloudiness, according to the quality indicator provided by S2ToolBox.
The broadband surface albedo has been calculated, when the observed surface is considered as Lambertian, as the integration of at-surface reflectance across the shortwave spectrum [26], as shown in equation: where α is albedo, ρ bi is surface reflectance for a given band bi at Level-2A Sentinel-2 surface reflectance, ω bi is the weighting coefficient representing the solar radiation fraction derived from the solar irradiance spectrum [26] within the spectral range (spectral response curves) for bands bi and is calculated as equation: where R sλ is extra-terrestrial irradiance for wavelength λ (µm); and UP bi and LO bi are upper and lower wavelength bounds for Sentinel-2A/B band bi, respectively.

Model-Based Assessment of Crop Water Requirement: AquaCrop
AquaCrop simulates crop yield in four steps: Crop development, crop transpiration, biomass production and yield formation. It calculates the daily soil water balance and divides evapotranspiration into soil evaporation and crop transpiration. AquaCrop describes the foliage development of the crop by the canopy cover (CC), that is formally equivalent to the fractional cover (fc) estimated by Sentinel-2 imagery, i.e., it is the fraction of soil surface covered by the green canopy. Hereinafter, we use the two terms canopy cover (CC) and fractional cover (fc) just to distinguish the two variables, respectively derived with AquaCrop and Sentinel-2 imagery.
Transpiration is a function of CC, while evaporation is proportional to the area of soil not covered by vegetation. The CC is multiplied by reference evapotranspiration (ETo), determined by the FAO Penman-Monteith equation, and the crop coefficient (Kc) to calculate potential crop transpiration. Actual transpiration (Ta) is calculated starting from potential one by accounting for water stress. Then, Ta is used for the calculation of crop biomass though its multiplication with water productivity normalized for climate. By using a harvest index (HI), crop yield is obtained by the biomass. To describe the effect of water stress, the model considers different thresholds of water available to the root zone. The first affects leaf canopy expansion (slowing down); the second threshold affects canopy senescence (quickening); the third is referred to as stomata closure (increase) and so to transpiration. Stress coefficients (Ks) range between 1 (no stress) and 0 (complete stress) and are multiplicative factors of the target process.
Model parameters are grouped into two classes: Conservative and non-conservative. Conservative parameters are not dependent on local and management conditions: Canopy growth (CGC) and canopy decline (CDC) coefficients; full canopy crop transpiration coefficient (Kc); biomass WP and soil water depletion thresholds. Non-conservative parameters vary depending on crop and field management, soil type, and climate (sowing date and density, length of crop cycle and phenological stages, maximum canopy cover, etc.). Such parameters can be either retrieved from AquaCrop literature or calibrated by the user (e.g., by means of field experiments). Main AquaCrop outputs are crop production (biomass and yield) and crop water use.

AquaCrop Implementation
In this study, a limited number of AquaCrop parameters were partly calibrated with field observations, including management information: Transplant dates and densities, flowering date and duration, starting of senescence, maturity, and final yield were used for local calibration of the model. For simulating irrigation, the model was set in net irrigation requirement mode, which estimates the crop water requirement based on a selected threshold of allowed root zone (water) depletion (RZD). In order to reproduce the irrigation method adopted by the farmer, drip irrigation was simulated to ensure that RZD was always above 50% of the readily available water (RAW).

Assimilation-Based Assessment of Crop Water Requirement
The third method for assessing crop water requirements was based on the integration of Sentinel-2 crop derived data with AquaCrop. The fractional cover (fc) estimated by Sentinel-2 has been sequentially assimilated into AquaCrop, by direct insertion, in place of the canopy cover (CC) simulated by the model. The sequential direct insertion is applied under the assumption that a continuous update of one crop model state based on remote observations can reduce the biases induced by the model simplifications of the processes and environmental conditions influencing the crop growth dynamics.
Crop CC simulated by AquaCrop along the growing season and the fc values measured by satellite were compared and the differences were statistically analyzed by means of the Pearson correlation

Test Site Characteristics
The parcel cultivated in 2017 had a Vitric Phaeozems (Eutric) soil (WRB classification), with a loamy texture (USDA classification), while soil was Vitric Cambisols with a loamy-sand texture in 2018. Soils were both deep, well drained and with no fertility constraints (Table 1). Considering the climatology of the study area, average temperature of 15.

AquaCrop Calibration and Implementation
AquaCrop parameters were calibrated to obtain the best fit between field observations and simulations, both in 2017 and 2018 (Table 2). AquaCrop simulated yields were 7.23 and 7.60 t/ha (dry weight) in 2017 and 2018, with an error of 0.42% and 3.40%, respectively (Table 3). The growth of the green canopy simulated by the model was compared with the fc values observed by IRRISAT. In both seasons the simulated growing curve fitted well with the satellite observations, although an underestimation for the initial canopy cover (late April-early May) and an overestimation during the last part of the growing season (July) was observed ( Figure 4). Nevertheless, the statistical analysis (RMSE, EF, d) showed a very good agreement between simulated CC and fc, which correlation was highly significant in both years (α = 0.001; Table 4; Figure 5). In 2017 the spatial variability of the crop canopy within the study parcel was larger than in 2018, as testified by the larger variance of the fc values retrieved in 2017.

AquaCrop Calibration and Implementation
AquaCrop parameters were calibrated to obtain the best fit between field observations and simulations, both in 2017 and 2018 (Table 2). AquaCrop simulated yields were 7.23 and 7.60 t/ha (dry weight) in 2017 and 2018, with an error of 0.42% and 3.40%, respectively (Table 3). The growth of the green canopy simulated by the model was compared with the fc values observed by IRRISAT. In both seasons the simulated growing curve fitted well with the satellite observations, although an underestimation for the initial canopy cover (late April-early May) and an overestimation during the last part of the growing season (July) was observed ( Figure 4). Nevertheless, the statistical analysis (RMSE, EF, d) showed a very good agreement between simulated CC and fc, which correlation was highly significant in both years (α = 0.001; Table 4; Figure 5). In 2017 the spatial variability of the crop

Estimation of Crop Water Irrigation Requirements
During the 2017 crop growing season, spanning from June 9th to July 20th, the accumulated rainfall was 40.4 mm, the estimated reference evapotranspiration (ETo) was 553 mm, while crop transpiration (Tr) and soil evaporation (E) estimated by AquaCrop were 345 and 192 mm, respectively. In such conditions, AquaCrop modeled an irrigation water requirement (IWR) of 461 mm corresponding to 144 L/plants.
In year 2018, wetter conditions occurred, with an accumulated rainfall of 125 mm and an estimated ETo of 439 mm. Under such conditions, both crop Tr (291 mm) and E (137 mm) were lower. Accordingly, a lower IWR, 332 mm, was estimated by the model. Being the final yields very similar, the consequent evapotranspiration (WP ET ) was higher in 2018 than 2017 (Tab 3).
In Figure 6, three different daily ETp series were compared, respectively produced by calibrated AquaCrop, IRRISAT and AquaCrop after replacing CC with fc estimated by IRRISAT. During the 2017 crop growing season, spanning from June 9th to July 20th, the accumulated rainfall was 40.4 mm, the estimated reference evapotranspiration (ETo) was 553 mm, while crop transpiration (Tr) and soil evaporation (E) estimated by AquaCrop were 345 and 192 mm, respectively. In such conditions, AquaCrop modeled an irrigation water requirement (IWR) of 461 mm corresponding to 144 L/plants.
In year 2018, wetter conditions occurred, with an accumulated rainfall of 125 mm and an estimated ETo of 439 mm. Under such conditions, both crop Tr (291 mm) and E (137 mm) were lower. Accordingly, a lower IWR, 332 mm, was estimated by the model. Being the final yields very similar, the consequent evapotranspiration (WPET) was higher in 2018 than 2017 (Tab 3).
In Figure 6, three different daily ETp series were compared, respectively produced by calibrated AquaCrop, IRRISAT and AquaCrop after replacing CC with fc estimated by IRRISAT. The agreement between AquaCrop and IRRISAT ETp estimates is excellent in the mid-season stage, between about the 50th and 85th day after transplant, when the crop canopy is fully developed. This result is remarkable considering that IRRISAT is fully based on remote observations weather data, and not relying on field or crop data. Moreover, IRRISAT predictions are not based on a continuous assessment of the crop dynamics, rather it provides ETp estimates based on crop parameters retrieved from the last available image, thus with a delay of five or more days, depending The agreement between AquaCrop and IRRISAT ETp estimates is excellent in the mid-season stage, between about the 50th and 85th day after transplant, when the crop canopy is fully developed. This result is remarkable considering that IRRISAT is fully based on remote observations weather data, and not relying on field or crop data. Moreover, IRRISAT predictions are not based on a continuous assessment of the crop dynamics, rather it provides ETp estimates based on crop parameters retrieved from the last available image, thus with a delay of five or more days, depending on the sky conditions (i.e., cloudiness).
In the initial stage and in the crop development stage (first 50 days after transplant), before canopy is fully developed, the estimation of ETp in AquaCrop is higher than the one determined by IRRISAT. This is due to the conceptually different calculation schemes: Whilst IRRISAT has a resistance-based approach in the Penman-Monteith linked to the LAI at the pixel scale (according to the "big leaf" schematization), AquaCrop has a crop coefficient approach for soil evaporation depending on soil properties. Diversely, with growing LAI and soil cover, the two estimates converge toward very similar values of ETp, hence irrigation requirements.
In the senescence stage, after about the 85th day from transplant, the values of ETp derived by means of the two methods slightly diverge again, with higher ETp for IRRISAT respect to AquaCrop, again due to the difference between the LAI-based approach of the first one and the crop coefficient of the second one. Furthermore, AquaCrop explicitly consider the transpiration reduction associated with the senescence, which cannot be fully predicted by the reduction of the observed LAI applied into the Penman-Monteith equation.
The canopy cover fc maps retrieved from Sentinel-2 imagery can be more effectively employed for constraining AquaCrop predictions to the actual observed crop growth and for assessing its spatial variability. In this study we applied a simple assimilation technique, known in the literature as "direct" insertion [27], consisting in replacing the model state variable CC with the remotely retrieved variable fc. As displayed in Figure 6, AquaCrop ETp estimates after direct insertion are essentially equal to the calibrated AquaCrop, except for the development and senescence stages, when the average fc sensibly deviates from the model calibrated growth curve. As illustrated in Table 3, fc direct insertion does not affect the cumulative ETp, rather it affects its partitioning in transpiration (Tr) an evaporation (E), with a slight increase in Tr and a decrease in E. In 2018, instead, fc direct insertion implies a slight reduction of ETp, mainly associated with a reduction of Tr. Correspondingly, tomato yield was overestimated compared with the observed and the calibrated values in 2017 (+14.3%), while slightly underestimated in 2018 (−1.03%). This result confirmed that sequential assimilation of one state variable does not necessarily improve the prediction performance of all model state variables. In this sense, more complex sequential assimilation techniques are needed, in order to account for the structure of the observation and prediction errors [28], as well of the atmospheric forcing [29].
The impact of the direct insertion in terms of irrigation water requirement (IWR) is null in 2017, while it determines a 5% decrease in 2018. As illustrated in Figure 7, also the temporal pattern of IWR is almost unchanged in the two simulated seasons. A potential water saving of 70 mm could be achieved in 2018 according to AquaCrop, by accounting for the effective contribution of the summer rainfall (125 mm) to the soil water deficit.
Looking again at Table 3, it is interesting to note that IRRISAT, despite the significant reduction of the cumulative ETp both in 2017 (16%) and 2018 (19%) compared with the calibrated AquaCrop, predicts cumulative IWR just 3% smaller in 2017 and 11% smaller in 2018 than the calibrated AquaCrop. As illustrated in Figure 7, this reduction is essentially due to the initial and development stage, when AquaCrop, differently form IRRISAT, accounts for the soil evaporation in the IWR assessment. After the crop is fully developed, the cumulative IWR curves are almost parallel, testifying the good agreement of the remotely assessed daily IWR with AquaCrop.
Looking again at Table 3, it is interesting to note that IRRISAT, despite the significant reduction of the cumulative ETp both in 2017 (16%) and 2018 (19%) compared with the calibrated AquaCrop, predicts cumulative IWR just 3% smaller in 2017 and 11% smaller in 2018 than the calibrated AquaCrop. As illustrated in Figure 7, this reduction is essentially due to the initial and development stage, when AquaCrop, differently form IRRISAT, accounts for the soil evaporation in the IWR assessment. After the crop is fully developed, the cumulative IWR curves are almost parallel, testifying the good agreement of the remotely assessed daily IWR with AquaCrop.

Conclusions
Sentinel-2 imagery can be effectively exploited for monitory canopy growth of tomato crops in open field. Irrigation advisory services, such as IRRISAT, which are only based on crop data retrieved by Sentinel-2 and weather data, can provide a reliable assessment of crop water requirements of tomato field crops especially when crop canopy is fully developed. It should be considered that IRRISAT does not require input data concerning the soil or crop phenology, since it is entirely based on crop growth monitoring from space. Hence, integrating Sentinel-2 imagery with a crop growth model such as AquaCrop can be an effective strategy for assessing crop water requirement in the initial and development stage of the crop, as well as for identifying the senescence stage. Further, being the satellite imagery a spatial information, the integration into a crop model can help in assessing crop water requirement at field or higher scales, i.e., at territorial level. Thus, a sequential assimilation can be used to support irrigation planning by Irrigation and Land Reclamation consortia. In this study a simple direct insertion method has been applied for assimilating canopy cover retrieved by Sentinel-2 imagery into AquaCrop, which does not guarantee an optimal model-data