Giovanni Paolini
email: giovanni.paolini@isardsat.cat

Maria Jose Escorihuela
email: mj.escorihuela@isardsat.cat

Magí Olivier Merlin

Pamies Sans

Joaquim Bellvert

M Pamies

Classification of Different Irrigation Systems at Field Scale Using Time-Series of Remote Sensing Data

Keywords: irrigation systems, actual evapotranspiration, soil moisture, time-series classification, machine learning, remote sensing, field scale. I

Maps of irrigation systems are of critical value for a better understanding of the human impact on the water cycle, while they also present a very useful tool at the administrative level to monitor changes and optimize irrigation practices. This study proposes a novel approach for classifying different irrigation systems at field level by using remotely sensed data at sub-field scale as inputs of different supervised Machine Learning (ML) models for time-series classification. The ML models were trained using ground-truth data from more than 300 fields collected during a field campaign in 2020 across an intensely cultivated region in Catalunya, Spain. Two hydrological variables retrieved from satellite data, actual evapotranspiration (ETa) and soil moisture (SM ), showed the best results when used for classification, especially when combined together, retrieving a final accuracy of 90.1±2.7%. All the three ML models employed for the classification showed that they were able to distinguish different irrigation systems, regardless of the different crops present in each field. For all the different tests, the best performances were reached by ResNET, the only Deep Neural Network model among the three tested. The resulting method enables the creation of maps of irrigation systems at field level and for large areas, delivering detailed information on the status and evolution of irrigation practices.

semi-arid areas where freshwater resources are limited such as the Mediterranean region, where irrigation currently accounts on average for 69% [START_REF] Fao | Aquastat -fao's global information system on water and agriculture[END_REF], [START_REF] Malek | Adaptation of land management in the mediterranean under scenarios of irrigation water use and availability[END_REF], and it goes up to over 84% for south and east-Mediterranean countries [START_REF] Unep | United nations environment programme, mediterranean action plan and plan bleu, state of the environment and development in the mediterranean[END_REF].

These percentages are projected to increase between 4 and 18% by the end of the century [START_REF] Guiot | Climate and Environmental Change in the Mediterranean Basin-Current Situation and Risks for the Future[END_REF], [START_REF] Cramer | Climate change and interconnected risks to sustainable development in the mediterranean[END_REF], but important uncertainty factors such as population growth and food demand may raise this estimate to between 22% and 74% [START_REF] Fader | Mediterranean irrigation under climate change: More efficient irrigation needed to compensate for increases in irrigation water requirements[END_REF].

Monitoring, understanding, and improving the efficiency of irrigation practices is a fundamental step toward controlling and mitigating water demands. Precise knowledge of irrigation practices is needed to better constrain and reduce uncertainty in hydrological models that predict future trends for water use and account for the anthropological impact on the water cycle. At the local scale, it is important to have an overview of irrigation practices at fine spatial resolution for administrative and irrigation management purposes, such as monitoring irrigation water usage and optimizing the efficiency of irrigation. Field irrigation efficiency mainly depends on the irrigation system and the level of field modernization. Generally, irrigation systems used in semi-arid areas include surface systems such as flood and furrow and pressurized systems such as sprinkler and drip [START_REF] Howell | Irrigation efficiency[END_REF]. Efficiencies vary between 90% for drip, 75% for sprinkler, and 60% for flood irrigation [START_REF] Brouwer | Irrigation water management: irrigation scheduling, training manual no. 4. annex i: Irrigation efficiencies. fao land and water development division[END_REF]: clearly, pressurized systems are more efficient. The low application efficiencies obtained by surface systems are mainly due to water losses associated with deep infiltration, soil evaporation, and flooding in some parts of the soil because of a poor flow design at the entrance of the furrow [START_REF] Chavez | Design and evaluation of surface irrigation systems applying an analytical formula in the irrigation district 085, la begoña, mexico[END_REF], [START_REF] Chávez | Water-use efficiency and productivity improvements in surface irrigation systems[END_REF].

Remote sensing is a unique and valuable tool, capable of addressing the lack of large-scale precise information over irrigation practices, and overcoming the limitations of analyses based on in-situ observations, which are often prone to inconsistencies and gaps in the information collected. Current results in the field of remote sensing for irrigation practices featured the creation of global or regional scale maps of irrigated areas [START_REF] Gao | Irrigation mapping using sentinel-1 time series at field scale[END_REF], [START_REF] Bazzi | Mapping irrigated areas using sentinel-1 time series in catalonia, spain[END_REF], [START_REF] Dari | Detecting and mapping irrigated areas in a mediterranean environment by using remote sensing soil moisture and a land surface model[END_REF], [START_REF] Peña-Arancibia | Dynamic identification of summer cropping irrigated areas in a large basin experiencing extreme climatic variability[END_REF], [START_REF] Deines | Annual irrigation dynamics in the u.s. northern high plains derived from landsat satellite data[END_REF], [START_REF] Pun | Land use classification: A surface energy balance and vegetation index application to map and monitor irrigated lands[END_REF], irrigation timings [START_REF] Chen | Detecting irrigation extent, frequency, and timing in a heterogeneous arid agricultural region using modis time series, landsat imagery, and ancillary data[END_REF], [START_REF] Page | Potential for the detection of irrigation events on maize plots using sentinel-1 soil moisture products[END_REF], [START_REF] Bazzi | Near real-time irrigation detection at plot scale using sentinel-1 data[END_REF] and quantification of irrigation amounts at variable resolutions [START_REF] Dari | Exploiting high-resolution remote sensing soil moisture to estimate irrigation water amounts over a mediterranean region[END_REF], [START_REF] Ouaadi | Irrigation amounts and timing retrieval through data assimilation of surface soil moisture into the fao-56 approach in the south mediterranean region[END_REF], [START_REF] Jalilvand | Quantification of irrigation water using remote sensing of soil moisture in a semi-arid region[END_REF], [START_REF] Zaussinger | Estimating irrigation water use over the contiguous united states by combining satellite and reanalysis soil moisture data[END_REF], [START_REF] Zappa | Detection and quantification of irrigation water amounts at 500 m using sentinel-1 surface soil moisture[END_REF], [START_REF] Modanesi | Optimizing a backscatter forward operator using sentinel-1 data over irrigated land[END_REF], [START_REF] Olivera-Guerra | Irrigation retrieval from landsat optical/thermal data integrated into a crop water balance model: A case study over winter wheat fields in a semi-arid region[END_REF], [START_REF] Abolafia-Rosenzweig | Soil moisture data assimilation to estimate irrigation water use[END_REF]. In particular, for studies oriented on the mapping of irrigated areas, remote sensing data are often coupled with Machine Learning (ML) models, proving to be successful with both supervised [START_REF] Gao | Irrigation mapping using sentinel-1 time series at field scale[END_REF], [START_REF] Bazzi | Mapping irrigated areas using sentinel-1 time series in catalonia, spain[END_REF] and unsupervised approaches [START_REF] Dari | Detecting and mapping irrigated areas in a mediterranean environment by using remote sensing soil moisture and a land surface model[END_REF].

Among ML models for classification, one of the most popular algorithms for classification is Random Forest [START_REF] Breiman | Random forests[END_REF],

widely used in the field of earth observation [START_REF] Belgiu | Random forest in remote sensing: A review of applications and future directions[END_REF] both for land cover mapping [START_REF]Random forest classifier for remote sensing classification[END_REF], [START_REF] Rodriguez-Galiano | Random forest classification of mediterranean land cover using multi-seasonal imagery and multi-seasonal texture[END_REF], [START_REF] Noi | Comparison of random forest, k-nearest neighbor, and support vector machine classifiers for land cover classification using sentinel-2 imagery[END_REF] and irrigation mapping [START_REF] Deines | Annual irrigation dynamics in the u.s. northern high plains derived from landsat satellite data[END_REF], [START_REF] Gao | Irrigation mapping using sentinel-1 time series at field scale[END_REF], [START_REF] Bazzi | Mapping irrigated areas using sentinel-1 time series in catalonia, spain[END_REF]. Recently, an adaptation of the Random Forest model for time-series classification was proposed by Deng et al., 2013 [36] and it is particularly useful for classification problems involving temporal changes. This model is often the baseline to test for new more performing ML models [START_REF] Bagnall | The great time series classification bake off: a review and experimental evaluation of recent algorithmic advances[END_REF] such as Rocket [START_REF] Dempster | Rocket: exceptionally fast and accurate time series classification using random convolutional kernels[END_REF], which has proven to be competitive with the state-of-the-art algorithms for classification of multivariate time-series [START_REF] Ruiz | Benchmarking multivariate time series classification algorithms[END_REF], but with a much faster computational speed. Deep learning models have also been recently adapted to time-series classification [START_REF] Fawaz | Deep learning for time series classification: a review[END_REF]. ResNET [START_REF] Dempster | Rocket: exceptionally fast and accurate time series classification using random convolutional kernels[END_REF] is a Deep Neural Network model that presents a more complex architecture than traditional models, which allows it to adapt to complex tasks at the expense of computational power. [START_REF] Fawaz | Deep learning for time series classification: a review[END_REF] [START_REF] Fawaz | Deep learning for time series classification: a review[END_REF] proved how ResNET outperforms traditional models for classification tasks in a comparison study performed using This article has been accepted for publication in IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing. This is the author's version which has not been fully e content may change prior to final publication. Citation information: DOI 10.1109/JSTARS.2022.3222884 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ a large variety of multi-disciplinary datasets. Moreover, in the field of irrigation mapping, Bazzi et al., 2020 [START_REF] Bazzi | Near real-time irrigation detection at plot scale using sentinel-1 data[END_REF] confirmed how deep learning models outperform traditional ML techniques.

The study of irrigation practices employs one or a combination of remotely sensed datasets from different sensors, from optical/thermal to active/passive microwave, at different spatial resolutions, but generally, two main hydrological variables are recognized to be key for the monitoring of irrigation practices: Actual Evapotranspiration (ET a ) and Soil Moisture (SM ) [START_REF] Massari | A review of irrigation information retrievals from space and their utility for users[END_REF]. ET a has been widely assessed during the last decades through Surface Energy Balance (SEB) models [START_REF] Norman | Source approach for estimating soil and vegetation energy fluxes in observations of directional radiometric surface temperature[END_REF], [START_REF] Bastiaanssen | A remote sensing surface energy balance algorithm for land (sebal).: Part 2: Validation[END_REF], [START_REF] Allen | Satellite-based energy balance for mapping evapotranspiration with internalized calibration (metric)-model[END_REF]. To be applied at field scale, SEB models require accurate Land Surface Temperature (LST) data with sufficient spatial resolution. Novel methods have been recently developed to derive ET a and potential evapotranspiration (ET p ) at 20 m resolution using data from the European Commission's Copernicus program [START_REF] Guzinski | Evaluating the feasibility of using sentinel-2 and sentinel-3 satellites for high-resolution evapotranspiration estimations[END_REF]. In particular, the method employed for this study consists in applying the Two-Source Energy Balance (TSEB) model [START_REF] Kustas | A two-source energy balance approach using directional radiometric temperature observations for sparse canopy covered surfaces[END_REF] with data from Sentinel-2 and Sentinel-3 in combination with meteorological data forcing from the Copernicus Climate Data Store (CDS) [START_REF] Guzinski | Modelling high-resolution actual evapotranspiration through sentinel-2 and sentinel-3 data fusion[END_REF], [START_REF] Bellvert | Feasibility of using the two-source energy balance model (tseb) with sentinel-2 and sentinel-3 images to analyze the spatio-temporal variability of vine water status in a vineyard[END_REF]. This approach relies on downscaling Sentinel-3 thermal bands to Sentinel-2 spatial resolution using a data mining sharpener (DMS) approach [START_REF] Gao | A data mining approach for sharpening thermal satellite imagery over land[END_REF].

Similarly, SM data from remote sensing has also been obtained at 20m resolution through disaggregation techniques. Passive L-band microwave sensors are usually preferred as low resolution SM input since they are recognized to have the highest sensitivity to SM and lower signal-to-noise ratio with respect to active radar or optical sensors [START_REF] Sabaghy | Spatially enhanced passive microwave derived soil moisture: Capabilities and opportunities[END_REF], [START_REF] Ulaby | Microwave Remote Sensing, Active and Passive[END_REF], at the expense of spatial resolution, which is in the order of tens of kilometers. A common solution to overcome this limitation in terms of spatial resolution is the downscaling of the SM products with optical/thermal data, which provides land surface parameters at higher spatial resolution than their radar counterparts [START_REF] Peng | A review of spatial downscaling of satellite remotely sensed soil moisture[END_REF]. DisPATCh (Disaggregation based on Physical And Theoretical scale Change) [START_REF] Merlin | A simple method to disaggregate passive microwave-based soil moisture[END_REF], [START_REF] Merlin | Disaggregation of smos soil moisture in southeastern australia[END_REF] has been applied numerous times to disaggregate SM from passive microwave sensors to higher spatial resolution, through the use of optical/thermal products from MODIS and Sentinel-3 at 1 km [START_REF] Merlin | Disaggregation of smos soil moisture in southeastern australia[END_REF] or Landsat at 100 m [START_REF] Merlin | Self-calibrated evaporation-based disaggregation of smos soil moisture: An evaluation study at 3 km and 100 m resolution in catalunya, spain[END_REF], [START_REF] Oiha | Sequential downscaling of the smos soil moisture at 100 m resolution via a variable intermediate spatial resolution[END_REF]. A disaggregation at 20 m has been recently proposed [START_REF] Paolini | Disaggregation of smap soil moisture at 20 m resolution: Validation and sub-field scale analysis[END_REF] with the use of SMAP daily low resolution SM gridded at 9 km, Sentinel-2 Normalized Difference Vegetation Index (N DV I) at 20 m, and the enhanced Sentinel-3 LST products disaggregated at 20 m.

As noted by Massari et al., 2021 [41] there is still an open question on obtaining maps of irrigation systems from satellite data. So far, studies have been limited to mapping irrigated from non-irrigated areas or singular irrigation systems. Numerous studies have been performed to distinguish irrigated from rain-fed areas. At field level, Gao et al., 2018 [14] proposed to directly use a Sentinel-1 backscatter product to train two ML models (random forest and support vector machine) and detect differences in the satellite signal between irrigated and not-irrigated fields, reaching an overall classification accuracy of 81%. Similarly, Bazzi et al., 2019 [15] trained different ML models (random forest and a convolutional neural network) using Sentinel-1 backscatter signal and Sentinel-2 NDVI timeseries, and reaching an overall accuracy of 94%. Passive microwave sensors also showed promising results in the detection of irrigation signals [START_REF] Lawston | Assessment of irrigation physics in a land surface modeling framework using non-traditional and human-practice datasets[END_REF], [START_REF] Escorihuela | Comparison of remote sensing and simulated soil moisture datasets in mediterranean landscapes[END_REF] but their coarse resolution does not allow for detection at field-level. [START_REF] Dari | Detecting and mapping irrigated areas in a mediterranean environment by using remote sensing soil moisture and a land surface model[END_REF] produced irrigation maps at 1 km resolution from SMOS and SMAP products disaggregated at 1 km, using an unsupervised clustering ML model and suggesting the need for high spatial resolution product to resolve the high spatial variability of irrigated areas.

Regarding the detection of unique irrigation systems, different studies have used deep learning to recognize the rounded shape of center pivot systems [START_REF] Rundquist | The nebraska center-pivot inventory: An example of operational satellite remote sensing on a long-term basis[END_REF], [START_REF] Yan | Automated crop field extraction from multi-temporal web enabled landsat data[END_REF], [START_REF] Saraiva | Automatic mapping of center pivot irrigation systems from satellite images using deep learning[END_REF]. Moreover, a deep learning approach was recently employed by Liang et al., 2021 [63] to map contour-levee (flood) irrigation from aerial pictures. Despite these techniques being very effective in identifying a unique irrigation system, a general approach for creating maps of irrigation systems is still missing. To this date, there is no study (to the best of our knowledge) on creating maps of irrigation systems (i.e. sprinkler, drip, flood) at field level. These maps can have a wide application in the scientific community since they could replace the simplistic assumption of irrigation scenarios used in many Land Surface Models (LSM), e.g. Noah [START_REF] Ozdogan | Simulating the effects of irrigation over the united states in a land surface model based on satellite-derived agricultural data[END_REF], [START_REF] Evans | Modeling the large-scale water balance impact of different irrigation systems[END_REF], [START_REF] Lawston | Assessment of irrigation physics in a land surface modeling framework using non-traditional and human-practice datasets[END_REF], [START_REF] Vahmani | Incorporating an urban irrigation module into the noah land surface model coupled with an urban canopy model[END_REF], or ORCHIDEE [START_REF] Yin | Improvement of the irrigation scheme in the orchidee land surface model and impacts of irrigation on regional water budgets over china[END_REF] models. Additionally, maps of irrigation systems could provide a useful tool for local policies, given that a complete and continuous overview of irrigation systems is lacking in many areas. These maps could give unprecedented insights to monitor the evolution of irrigation practices and promote and supervise the shift towards more sustainable and efficient irrigation methods.

In this context, we propose a novel methodology to produce maps at the field scale distinguishing between the three main irrigation systems, drip, sprinkler and flood, and also not irrigated fields. The hypothesis of this study is that differences between irrigation systems should be detectable by analyzing temporal patterns of Actual Evapotranspiration (ET a ) and Soil Moisture (SM ) at the field or sub-field scale, through the use of a supervised ML model. It is expected that time series of remote sensing data reveal distinctive temporal patterns among different irrigation systems, given the large variation in the amount and timing of water applied for the different systems.

The proposed methodology will be applied to a semi-arid area of the Ebro basin (Spain) that is characterized by high variability of crop types and irrigation systems. The methodology will be then validated against both in-situ data and independent administrative data, retrieved from statistical estimations at the district level and large surveys among farmers.

The study is organized as follows: Section II presents the study area, the field campaign, the administrative dataset, and remote sensing data used. The first part of Section III presents how time-series of remote sensing data are prepared, the additional crop classification of the ground truth dataset, the ML methods used, the final post-processing applied after the ML classification and the metrics used for evaluation of the performances. A second part of section III then introduces the different experiments performed in this study. Section IV presents the results in classifying irrigation systems using different remote sensing variables and ML methods. Finally, section V presents a final summary of the results and discussion.

II. MATERIALS

A. Study Area

The selected study area (41.28-42.02 N, 0.27-1.3 E) is located in the northeast of the Iberian Peninsula, in the province of Lleida (Spain). The climate in the region is typically Mediterranean with an average annual precipitation and reference evapotranspiration (ET 0 ) of 350 mm and 1100 mm, respectively. Irrigation in the area usually starts in mid-March and lasts until the end of November. The area is densely irrigated, with a variety of different techniques depending on the degree of modernization and the water allocation for the respective irrigation district: from traditional irrigation systems based on flooding techniques to more recent and efficient techniques that use sprinkler or drip irrigation.

The study area is divided into eight irrigation districts, covering a total surface of around 3000 km 2 . An overview of the study area is provided in Fig. 1. Irrigation practices vary depending on the seasonal water allocation, the different crop types, and the modernization level of each irrigation district. As an example, the "Canals d'Urgell" district is one of the oldest districts in the area, and irrigation is mainly performed through flooding. Farmers have full water availability throughout the growing season, but irrigation is performed in turns every 15-20 days. On the other hand, Algerri-Balaguer is a modernized pressurized irrigation district with a water allocation of 6000 m 3 /ha. Crops are mainly irrigated by sprinkler or drip irrigation. The Garrigues Sud district, located in the southern part of the study area, has a seasonal water allocation of around 1300 m 3 /ha, affecting, therefore, the type of crops grown (mostly olives and almonds) and irrigation practices (sustained deficit irrigation). The region has various types of cultivated crops, that can be grouped into: winter cereals (accounting for around 34 % of the total area according to administrative databases), maize (accounting for around 7 %, but not considering the percentage of maize growing as a second crop after winter cereals), fruit and nut trees (14 %), vineyards (1 %) and olives (9 %).

This article has been accepted for publication in IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing. This is the author's version which has not been fully e content may change prior to final publication. Citation information: DOI 10.1109/JSTARS.2022.3222884

B. Field Campaign

A field campaign was performed inside the study area during 2020, in order to collect a dataset of ground truth samples of fields with different irrigation systems and crop types. More than 300 fields were classified using 4 different labels: sprinkler, flood, drip, and non-irrigated. Fig. 1 shows how the samples were randomly distributed across all the irrigation districts considered, in order to have a representative dataset for the area. Table I summarizes the different irrigation systems of all the fields collected. It is possible to notice how approximately the same number of samples was collected for each irrigation system, ensuring a balanced dataset. In addition, each field was initially classified also by crop type by visual inspection and grouped in different classes: winter cereals, maize, alfalfa, olive, vineyards, or fruit and nut trees. Collecting these two variables provided an overview of the relationship between irrigation systems and crop types, which was needed to ensure the collection of a complete dataset representative of the different typologies of fields in the area.

C. SIGPAC-DUN administrative dataset

An administrative database, SIGPAC-DUN (Sistema d'informació geogràfica de parcelles agricoles), was used in order to verify and expand the information about crop types collected during the field campaign, and to extract the exact shape contour of each field. SIGPAC-DUN is provided by the Catalan Ministry of Climate Action, Food and Rural Agenda and it contains a large variety of spatial and alphanumeric information about agricultural practices at parcel level, with yearly updates. Most of the information contained in SIGPAC-DUN is directly submitted by the farmers through an annual agrarian declaration of cultivated crops and it gathers multiple details over the usage of the fields such as crop type. Parcel's shapes are also contained in SIGPAC-DUN, created from cadastral maps and image interpretation [START_REF] Sigpac | sistema de información geográfica de parcelas agrícolas), ministry of climate action, food and rural agenda[END_REF]. The dataset also contains information on presence of irrigation, indicated as a percentage of the field subject to irrigation (from 0 to 100), exclusively based on administrative data that indicates if irrigation is installed in the field, not its actual usage. the year and what irrigation system (if present) was installed in the field. For this reason, when using SIGPAC-DUN it is not possible to catalog the fields with double crops when a second crop type is present during the year. When comparing the systems of irrigation declared in the latest SIGPAC-DUN database against the ground-truth dataset collected during the field campaign, a discrepancy of around 10% was found between the two datasets (33 out of 332 collected fields for 2020, as shown in Fig. 2). This discrepancy seems to suggest that SIGPAC-DUN reflects an outdated catalog of irrigation systems, since most of the misclassification between SIGPAC-DUN and the groundtruth dataset are between traditional flooding systems or not irrigated fields for SIGPAC-DUN and modern irrigation systems for the ground-truth dataset. This suggests that a process of modernization of the irrigation systems is taking place in the area, but it is not registered. As a matter of fact, Fig. 2 clearly shows that the highest discrepancy is found for the 10 fields mis-classified as flood by SIGPAC-DUN, which are in reality sprinkler systems.

D. Remote Sensing Data

Various remote sensing products were evaluated as potential inputs for the classification task, derived from different satellite data. ET a and SM at the sub-field level (20 m) were the main hydrological variables considered, but additional variables were also considered in order to evaluate their feasibility in the classification task. These inputs can be broadly grouped into two categories, called Level 4 (L4) and Level 2 (L2) variables. A general overview of these variables is presented in Table II. L4 variables were estimated by the combination of multiple satellite data into different models, in order to obtain a set of more detailed hydrological information with a unified spatial resolution of 20 m, while L2 variables represent data directly retrieved from the satellites at their original processing level. L4 variables are: Actual Evapotranspiration (ET a ), DisPATCh Soil Moisture (SM ), crop water stress coefficient (K s ) and Sentinel-2 Leaf Area Index (LAI).

K s was considered given the proven strong link between this stress index and root-zone water depletion [START_REF] Colaizzi | Estimating soil moisture under low frequency surface irrigation using crop water stress index[END_REF], which could provide new additional information about the field water content at a different depth level than surface SM .

LAI was instead selected to test if a variable only related to vegetative growth could perform well in the task of classifying irrigation systems.

ET a estimates were obtained with the Priestley-Taylor Two-Source Energy Balance (TSEB-PT) model using

Copernicus-based inputs [START_REF] Norman | Source approach for estimating soil and vegetation energy fluxes in observations of directional radiometric surface temperature[END_REF], [START_REF] Kustas | A two-source energy balance approach using directional radiometric temperature observations for sparse canopy covered surfaces[END_REF], [START_REF] Morillas | Using radiometric surface temperature for surface energy flux estimation in mediterranean drylands from a two-source perspective[END_REF] and following the methodology described by Guzinski et al. 2020 [START_REF] Guzinski | Modelling high-resolution actual evapotranspiration through sentinel-2 and sentinel-3 data fusion[END_REF] which produced and validated a 20 m ET a product derived by applying TSEB-PT to remotely sensed data. The main input data required to run the TSEB-PT were retrieved from Sentinel-2 shortwave observations, Sentinel-3 LST , and ERA5 meteorological reanalysis data. Sentinel-2 images were used to retrieve the biophysical parameters of the vegetation through the Biophysical Processor [START_REF] Weiss | Sentinel2 toolbox level2 products s2toolbox level 2 products: Lai, fapar, fcover version 1.1[END_REF]. These biophysical parameters were used to derive inputs needed in the TSEB-PT such as leaf area index (LAI), leaf optical properties, and transmittance. High-resolution thermal data at 20 m was retrieved by applying a Data Mining Sharpening algorithm [START_REF] Gao | A data mining approach for sharpening thermal satellite imagery over land[END_REF] to Sentinel-3 SLSTR LST images at 1 km, using shortwave multi-spectral data from Sentinel-2 [START_REF] Guzinski | Modelling high-resolution actual evapotranspiration through sentinel-2 and sentinel-3 data fusion[END_REF], [START_REF] Bellvert | Feasibility of using the two-source energy balance model (tseb) with sentinel-2 and sentinel-3 images to analyze the spatio-temporal variability of vine water status in a vineyard[END_REF] as a higher resolution proxy. Meteorological parameters were retrieved from the ERA5 meteorological reanalysis, which delivers an hourly product gridded at 0.1º. The required meteorological parameters are: air and dew point temperature at 2m, wind speed at 100 m, surface pressure, and total column water vapor. Finally, ancillary data such as vegetation height, and leaf inclination angle were set based on a land cover map obtained from the Copernicus Global Land Service and a look-up table.

Instantaneous energy fluxes at the satellite overpass were upscaled to daily water fluxes, expressed in units of mm/day, by multiplying the instantaneous ratio of latent heat flux over solar irradiance by the average daily solar irradiance [START_REF] Cammalleri | Upscaling of evapotranspiration fluxes from instantaneous to daytime scales for thermal remote sensing applications[END_REF].

Crop stress coefficient (K s ) was retrieved from the ratio of ET a and potential evapotranspiration (ET pot ). In order to be consistent with the TSEB-PT, the two-layer Shuttleworth-Wallace (SW ) model [START_REF] Shuttleworth | Evaporation from sparse crops-an energy combination theory[END_REF] The considered L2 variables were: SMAP SM , Sentinel-3 LST and Sentinel-2 N DV I. SMAP SM is the enhanced L2 passive SSM product (L2_SM _P _E) from the SMAP mission, gridded at 9 km [START_REF] Entekhabi | The soil moisture active passive (smap) mission[END_REF]. LST is the L2 product from the SLSTR instrument on-board the Sentinel-3 satellite, which delivers daily 1 Km data [START_REF] Team | Sentinel-1 user handbook[END_REF].

N DV I was produced by combining band 4 (visible) and band 8A (near-infrared) from the Level 2 product of the Sentinel-2 satellite, with a 20 m resolution and a temporal resolution of around five days.

III. METHODS

A. Time Series Data Preparation

Annual time-series were extracted for each pixel of each field of the ground truth database. Three different years, 2018, 2019, and 2020 were used to create three annual time-series per pixel. While data on irrigation systems were only collected during the field campaign of 2020, no changes in irrigation systems were assumed for the two previous years: an assumption that was confirmed for the majority of the fields by inquiry with farmers and/or professionals working in the area. Considering multiple years is beneficial to 1) substantially increase the ground truth dataset and 2) allow the models to learn and generalize from a larger dataset, more diverse in terms of meteorological and crop growing conditions. This increased variability in the dataset allows the ML models to be more robust to changes. After the extraction of each time-series, a strategy was selected in order to fill the gaps whenever the data was unavailable for a particular day or pixel. For the case of ET a time-series, gaps were filled following the methodology proposed by Jofre-Čekalović et al., 2022 [START_REF] Jofre-Čekalović | Accounting for almond crop water use under different irrigation regimes with a two-source energy balance model and copernicus-based inputs[END_REF]: when not available, ET a was retrieved by multiplying reference evapotranspiration (ET 0 ) with the crop coefficient (K cs ). K cs was obtained as the ratio between ET a and ET 0 for those days with available data, while temporally interpolated for the missing dates. For the case of the DisPATCh SM time-series, the filling was performed using the original SM values from SMAP SM . For the rest of the variables used in this study, a simple linear interpolation was implemented as a gap-filling methodology.

As an additional pre-processing step, every time-series from each pixel and each variable are scaled through z-normalization, a standard technique that can speed up ML model convergence and improve performances [START_REF] Bagnall | The great time series classification bake off: a review and experimental evaluation of recent algorithmic advances[END_REF], [START_REF] Zhang | Understanding deep learning (still) requires rethinking generalization[END_REF]. From this dataset at the pixel level, a field level dataset was created by calculating the median of all pixels contained in each field. The dataset at field level was used for the experiment with simpler ML models, while for deep neural network models the large dataset at pixel level was needed in order to tune all the parameters and avoid overfitting. Finally, these datasets were split into two equal parts, 50% for training and 50% for testing of the classical ML models. keeping an equal distribution of irrigation systems and crop types in the two groups in order to avoid imbalances towards a particular irrigation or crop type in the training or testing of the classification. Moreover, time-series from the same field were used consistently for only one task, training or testing, in an attempt to avoid undesired correlation between the two datasets. Finally, for each ML model, 10 different runs were performed in order to extract more reliable performance metrics. During the different runs, train and test datasets were shuffled each time in a random fashion, but keeping the same constraints on the distribution of crop types, irrigation systems and using same-field time-series for training or testing only.

B. Classification of Crop Types

Regarding specific information about crop types, using the SIGPAC-DUN dataset only partially completed the missing information about the years previous to 2020. SIGPAC-DUN does not contain information on the presence of secondary crops. For this reason, an additional analysis was performed on the fields with annual crops. This analysis consisted of a simple crop classification algorithm applied to the LAI time-series to detect the number of peaks occurring during the growing season and check for the presence of multiple crops along the same year. More specifically, LAI time series from Sentinel-2 at 20 m resolution were collected for each ground truth field and for each considered year. The median value was extracted among the pixels inside every field, for every available date. The resulting time series were first processed in order to remove outliers with the Hampel filter algorithm (using a threshold value of 3 and a window length equals to 3) and then linearly interpolated to daily intervals.

Each time-series was classified based on the number of peaks present during the year using a simple peak detection algorithm which distinguished among winter crops (a single peak in the winter/spring period), summer crops (a single peak in the summer period), double crops (double peaks), or alfalfa, grown and harvested multiple times during the spring-summer period (multiple peaks in the summer period). Fig. 3 shows the LAI time-series for the different classes of annual crops that are grouped based on the number and position of peaks. Table III summarizes the number of fields collected in the campaign and used in this study as a dataset to train and test the ML models. These fields are divided per crop type and irrigation system, for each one of the three years considered. It is possible to notice how fields with annual crops vary in number throughout the three years since they showed changes in crop type from one year to another. Additionally, the total number of fields varies across the years varies since for a small number of fields the crop type could not be identified clearly, hence they were not considered for the specific year. to distinguish different systems of irrigation [START_REF] Massari | A review of irrigation information retrievals from space and their utility for users[END_REF]. Three different ML models were used, slightly adapted from the 283 model presented in [START_REF] Löning | sktime: A unified interface for machine learning with time series[END_REF] and [START_REF] Fawaz | Deep learning for time series classification: a review[END_REF].

MODELS Reason for selection Time Series Forest

Random forest is widely used in similar research [START_REF] Gao | Irrigation mapping using sentinel-1 time series at field scale[END_REF]. Used as a benchmark.

ROCKET

State-of-the-art accuracy [START_REF] Ruiz | Benchmarking multivariate time series classification algorithms[END_REF] with low computational requirements.

ResNet

Deep neural network with flexible structure, it outperforms other models in general reviews [START_REF] Fawaz | Deep learning for time series classification: a review[END_REF].

Time series forest [START_REF] Deng | A time series forest for classification and feature extraction[END_REF] comes from the family of (decision) tree ensemble classifiers: it is a random forest [START_REF] Breiman | Random forests[END_REF] adapted to detect temporal features. Random forest is a set of classification trees, where each tree is trained on a random but independent portion of training data, using bagging or bootstrapping to select these training subsets [START_REF] Breiman | Bagging predictors[END_REF]. Random forest algorithm are widely used [START_REF] Gao | Irrigation mapping using sentinel-1 time series at field scale[END_REF], [START_REF] Bazzi | Mapping irrigated areas using sentinel-1 time series in catalonia, spain[END_REF], [START_REF] Deines | Annual irrigation dynamics in the u.s. northern high plains derived from landsat satellite data[END_REF], [START_REF] Belgiu | Random forest in remote sensing: A review of applications and future directions[END_REF], [START_REF]Random forest classifier for remote sensing classification[END_REF], [START_REF] Rodriguez-Galiano | Random forest classification of mediterranean land cover using multi-seasonal imagery and multi-seasonal texture[END_REF], [START_REF] Noi | Comparison of random forest, k-nearest neighbor, and support vector machine classifiers for land cover classification using sentinel-2 imagery[END_REF]: the reason for their success lies in the low computational power required when compared to similar ML techniques, its stability (by design) against over-fitting, and its robustness against mislabelled training data. Additionally, this algorithm has a notable advantage in terms of interpretability of its prediction: each prediction has a confidence level that is retrieved from the percentage of trees that voted the same class. Time-series Forest is a variation of Random Forest where each tree is split using a combination of distance and entropy gain: an approach that captures well temporal characteristics of the inputs. The model also allows for an insightful inspection of the classification process thanks to the possibility of producing temporal importance curves. These curves underline the parts of the time-series that contain the most useful information and reveal which is the most important statistical feature among the ones extracted. Given the wide diffusion of this ML method and its interpretability, Time Series Forest was chosen as the benchmark model to run the classification experiment presented in this study.

ROCKET (RandOm Convolutional KErnel Transform) is another algorithm designed for time-series classification.

It was proposed by Dempster et al., 2020 [START_REF] Dempster | Rocket: exceptionally fast and accurate time series classification using random convolutional kernels[END_REF] and proven to be competitive with the state-of-the-art algorithms for classification of multivariate time-series [START_REF] Ruiz | Benchmarking multivariate time series classification algorithms[END_REF], but with a much faster computational speed. It is a kernel-approach classification inspired by convolutional neural networks. After producing a large number of kernels, two main features are extracted (maximum value and portion of positive values) that are then used to train a linear classifier (a ridge regression algorithm, as proposed in the original paper). An innovative aspect of ROCKET is the existence of a single hyper-parameter, corresponding to the number of kernels (set to 10, 000 by default), which avoids computationally intensive hyperparameters's tuning, a task required by many other classifiers. [START_REF] Fawaz | Deep learning for time series classification: a review[END_REF], [START_REF] Wang | Time series classification from scratch with deep neural networks: A strong baseline[END_REF], given their success in terms of performances and wide diffusion for similar classification tasks. ResNET was shown to consistently outperform traditional ML models when applied to a large variety of different datasets [START_REF] Fawaz | Deep learning for time series classification: a review[END_REF]. ResNET main architecture is used as in Wang et al., 2017 [83]. This model was adapted in this study to be run with multiple variables (multivariate model), through the use of a late fusion of parallel networks [START_REF] Bazzi | Near real-time irrigation detection at plot scale using sentinel-1 data[END_REF], [START_REF] Audebert | Beyond rgb: Very high resolution urban remote sensing with multimodal deep networks[END_REF], [START_REF] Gaetano | A two-branch cnn architecture for land cover classification of pan and ms imagery[END_REF], [START_REF] Feng | Integrating multitemporal sentinel-1/2 data for coastal land cover classification using a multibranch convolutional neural network: A case of the yellow river delta[END_REF], [START_REF] Interdonato | Duplo: A dual view point deep learning architecture for time series classification[END_REF].

The more complex architecture of this Deep Learning model requires the use of a large training dataset in order to prevent the rapid over-fitting of the model. For this reason, pixel-level time-series were used. Time-series contained in this dataset could be redundant, since adjacent pixels are expected to show similar values, given that spatial variability at 20 m resolution is limited. Nevertheless, the small variations that they present could still potentially improve the model accuracy, similarly to the improvement produced by most of the data-augmentation techniques, often used in deep learning, where slight changes are introduced in the dataset to create new training data [START_REF] Fawaz | Deep learning for time series classification: a review[END_REF], [START_REF] Forestier | Generating synthetic time series to augment sparse datasets[END_REF].

Training and testing datasets are still selected following the criteria of the other ML models, keeping all time-series from the same fields either for training or testing, thus avoiding data-leaking effects. Final performance metrics are still evaluated at field-level (after a spatial aggregation of the irrigation systems).

D. Post-Processing

After the model's training, annual maps of irrigation systems were produced and two main post-processing steps are implemented in order to correct for possible misclassification. Correcting the classification output with a statistical or knowledge-based approach is commonly used for multi-temporal geo-spatial classifications [START_REF] Deines | Mapping three decades of annual irrigation across the us high plains aquifer using landsat and google earth engine[END_REF], [START_REF] Elhag | Land use changes and its impacts on water resources in nile delta region using remote sensing techniques[END_REF], [START_REF] Herrault | Automatic extraction of forests from historical maps based on unsupervised classification in the cielab color space[END_REF]. The first post-processing steps involved spatial aggregation at field-level for the irrigation systems maps produced at pixel-level. The SIGPAC-DUN field shapes were used as a mask and spatial aggregation was performed in order to select only one irrigation system for each field, based on the most recurring irrigation system predicted among all the pixels contained in the field. This approach is also realistic since only one system of irrigation is expected to be found for each field. The second post-processing step involved a temporal analysis to detect and filter unlikely single-year changes of irrigation for a specific field. As a general rule, the presence of a single-year anomaly in the irrigation systems was always corrected, except for cases where the anomaly was found in the first or last year analyzed and the change could be explained by a modernization of the irrigation system (from not irrigated to irrigated or from traditional irrigation as flood to more modern irrigation systems such as drip and sprinkler). This knowledge-based temporal correction assumes that irrigation practices are not interrupted from one year to another and that modern irrigation systems are never replaced by traditional irrigation, given the significant infrastructure cost and no real production benefits of such change.

E. Evaluation Metrics

In order to compare the results from different ML models and variables, standard evaluation metrics are calculated by comparing the predicted versus actual irrigation system of the fields present in the test dataset. The Accuracy is computed as the number of correct predictions over the total number of samples. The Precision, also called User Accuracy, is instead calculated for each irrigation system category as the number of correct predictions over the number of samples predicted with the same category. The Average Precision is expressed as the average of all the precisions calculated for all the categories. The Recall, also called Producer Accuracy, is instead calculated for each irrigation system category as the number of correct predictions over the number of samples with the same category in the ground truth. The Average Recall is expressed as the average of the recalls calculated for all the categories.

F. Experiment Design

Tree different experiments were designed in order to explore different approaches in classifying irrigation systems.

1) Influence of Crop Types:

The first experiment aimed at verifying if the model was able to correctly identify differences exclusively related to irrigation practices, or if it was merely classifying irrigation systems based on the crop type present in each field. There is a proven relationship between crop types and irrigation systems, where in most of the cases few prevalent systems of irrigation are present for each particular crop type, as shown in Table III. Different models of Time Series Forest were trained on each crop type to predict the irrigation systems separately. Their aggregated accuracy was then compared with a general Time Series Forest model trained without discriminating by crop types in order to detect which approach was more favourable. The experiment was performed to assess the accuracy of SM , ET a , and both the variables together, SM +ET a .

2) Importance of crop vegetative period: A second experiment was also designed using the same variables and the time-series forest model. This second experiment was used to check whether time-series classification models required to be manually cropped in advance or whether the model was able to independently select the period of more intensive irrigation. Time-series forest was run with only a part of the time-series, which was cropped to isolate only the vegetative period of crops, in which there was a greater intensity of irrigation. Cropping implied selecting the spring and summer period (from the 15th of May) for all the crops except for winter cereals, where the winter season (until the 15th of July) was used. The need for cropping time-series was evaluated through a comparison of the overall accuracies from the classification of irrigation systems using cropped time-series with the accuracies retrieved using entire time-series (e.g. ET a vs ET a,cropped , SM vs SM cropped , etc.). Another approach used to evaluate the need for cropping time-series was to visualize Temporal Importance Curves in order to verify if the Time Series Forest model trained with the complete time-series was able or not to independently select the most important part of the year for the classification of irrigation systems. These curves do not only show the most important period of the year for the classification task, but they also provide information on which of the extracted features is most useful (among the ones selected in this study: mean, standard deviation, and slope).

B. Importance of crop vegetative period

Table V also shows two additional results: first, it is possible to notice how combining ET a +SM leads to consistently higher classification accuracies than when using these two variables separately. All crop types show better accuracy when both variables are used for classification. The only exception to these results seems to be the classification of irrigation types for alfalfa, where accuracies remains low even when using ET a + SM , when comparing the two different irrigation systems present, flood and sprinkler. A possible explanation could be the multiple rapid vegetative growth cycles that characterize this crop. The rapidity may result in very similar irrigation practices between flood and sprinkler, since a shorter period of time available for irrigation reduces the variability of the two irrigation systems both in terms of water amount and intervals of time between consecutive irrigations.

A second result is that cropping does not lead to improvements in terms of accuracy. Removing parts of the timeseries actually degrades the final results, confirming that ML models for time-series are able to independently select and exploit the most valuable part for the time-series with no need for this preprocessing step. As an additional proof of the capability of the ML model to independently select the most interesting part of the time-series, Temporal

Importance Curves from the Time Series Forest were calculated. Fig. 4 shows the importance curves for the irrigation classification task for the different crop types and the corresponding variables (ET a and SM ) used to generate them. Temporal importance curves show that the model only focuses on a specific part of the year that corresponds to the irrigation and growing periods for each crop. It is noticeable that the feature importance curves for ET a and SM do not always select the same period and the same feature, suggesting that both variables provide complementary information for this classification task. Among the three feature importance curves for the ET a time-series, the 'mean' curve has a higher value for most of the crop types. This is also visually evident when looking at the time-series ET a under different irrigation systems: ET a shows a clearer difference in magnitude for fields irrigated with different irrigation systems, which has a direct relationship with the amount of water supplied to the field. Flood-irrigated fields tend to have the highest values of ET a , especially during the warmest period, when water is almost exclusively supplied by irrigation. This behavior was expected since flood irrigation is a traditional system that notably employs the largest quantities of water. Flood-irrigated fields show higher ET a for all the annual crops, from Fig. 4(a) to 4(d), especially for winter cereals, where flood-irrigated fields are compared to non-irrigated fields. For fields cultivated with double crops, alfalfa, and maize, flood is compared with sprinkler-irrigated fields and it is notable how flood-irrigated fields are still evapotranspirating more, even if the difference between these two systems is less evident. The second row corresponds to orchards, from Fig. 4(e) to 4(g): it is noticeable that when comparing non-irrigated and drip-irrigated fields (for the case of olive and vineyard fields), the model is still selecting the 'mean' as the most important among the feature importance curves since drip-irrigated fields shows higher values than non-irrigated fields. Only in the case of fruit and nut-bearing trees 'slope' is most important. The 'slope' curve produced from the ET a time-series is calculated as the first derivative of the time-series, and can be interpreted as the degree of change of the time-series during the season: Fig. 4(e) shows a clear different timing in irrigating with drip and flood, so there is a clear difference in temporal changes of ET a , while non-irrigated fields have almost constant values during the summer season, thus not showing temporal changes.

Regarding SM time-series, it is most evident how fields irrigated differently have a different range of variation of SM , especially during the warmest periods. For this reason, the 'std' feature importance curve is selected by the model as the most informative: For winter cereals, it is evident that the 'std' curve detects differences between the different irrigation systems during the crop growing period. Similarly, the 'std' feature importance curve is also showing the highest peaks during the summer period for the case of alfalfa, double crops, vineyard, and olive fields.

The only exception is present for the importance curves created from the SM time-series of the fruit and nut-bearing trees. In this case, the 'slope' curve presents the highest peak, which is during the month of May/June: this is also visible in the SM curves, where during May the fields irrigated by a drip system show marked differences with respect to flood-irrigated and non-irrigated fields.

C. Model and variable selection

After verifying the capability of the model to classify irrigation systems, we proceeded to investigate which variables are most suited for the classification of irrigation systems and which model performs better for this classification task. Fig. 5 summarizes the accuracies retrieved for the three different models trained using different sets of the input variables described. All models have higher accuracy when using L4 variables, with the highest accuracy in terms of average being reached when all the L4 variables are used together, followed by the combination of ET a and SM only. Despite performing best, the accuracy reached combining all the L4 variables is very close to the accuracy reached using ET a and SM together, but it shows a larger standard deviation, caused by the addition of K s and LAI which probably do not positively contribute to the classification. For this reason, we considered ET a and SM to be the best combination of variables. These two variables show good accuracy even when used separately, so it could be possible to use these single variables in the classification process, only losing a small amount of accuracy.

Nevertheless, ET a and SM together reach higher accuracies because their information complement each other:

while SM is capable of distinguishing better large wet surfaces caused by flood irrigation, ET a can detect higher plant evaporation from drip with respect to non-irrigated fields. Combining the two variables always brings a general improvement of the prediction, as also demonstrated in Table VII.

Regarding the different models used, the ResNET model consistently shows higher results. ResNET did not only out-perform the other models in terms of accuracy, but it showed consistently higher results for all the different metrics used to compare the different models, as presented in The final overall accuracy of ResNET in classifying different irrigation systems is comparable (and in some cases higher) to accuracies presented in literature for studies involving irrigation mapping at field level, such as Gao et al., 2018 [14] and Bazzi et al., 2020 [15]. This result shows that with the presented approach it is possible to keep a high overall accuracy even when adding complexity to the problem of irrigation classification. The only drawback of ResNET is that it requires a significantly higher computational cost for the model training than the other two traditional models. In case computational cost is an issue, a good trade-off between accuracy and computational time is offered by the Rocket model, which is less accurate than ResNET for this specific classification task but it is around one order of magnitude faster in training and suggested to be used as a default model for multivariate classification tasks given its remarkable results for large scale studies [START_REF] Ruiz | Benchmarking multivariate time series classification algorithms[END_REF]. Another advantage for Rocket is that its training time is linearly scalable with the size of the training set [START_REF] Dempster | Rocket: exceptionally fast and accurate time series classification using random convolutional kernels[END_REF], which could be of great value in case this classification approach is applied to large irrigated areas, where the number of training samples inevitably grows.

Finally, Fig. 6 shows the confusion matrix for ResNET with ET a and SM as inputs. The matrix is calculated on the final results, after the post-processing step, which included aggregating the model prediction at the field level and performing temporal post-processing for the three different years. Precision and recall values are also presented and indicate how all the values of the metrics are very close to each other: an additional indicator of the robustness of the classification, which is not imbalanced towards any particular irrigation system. The lowest metric is represented by the precision for the drip irrigation system, which appears to be the label that is most mis-classified by the ML model: as a matter of fact, in a few cases, the model appears to classify drip irrigation as flood or non-irrigated. Drip irrigation is sometimes confused with non-irrigated fields due to the low soil wet surface around the emitter, which minimizes losses through evaporation and runoff. Additionally, there is also a more marked misclassification in those irrigation districts with limited water allocation, where sustained deficit irrigation strategies are usually adopted, such as some areas of Segarra-Garrigues or Garrigues Sud. We have also realized that some recently planted fields of grapevines, almonds, and pistachios trees were classified as non-irrigated, but instead were drip-irrigated. This probably occurred due to its still low canopy vigor, evapotranspiration, and soil moisture values throughout the growing season. On the other hand, the confusion between drip and flood irrigation could be explained either due to a decrease in SM and ET a on dates between irrigation events or due to a higher ET a caused by crop cover between rows. In both cases, time-series between flood and drip-irrigated fields may look similar. In order to improve classification for these particular cases, the selection of more fields with these characteristics in further studies will help to obtain a more robust classification. showed some inaccuracies already from the comparison with the ground truth dataset. In particular, SIGPAC-DUN showed a tendency to misclassify modern irrigated fields as not irrigated or irrigated with the traditional system of flood. This suggested that the database probably reflects a picture that is not up-to-date, and it is also visible in this general comparison of the entire study area with the map produced by the ResNET model. A secondary cause for the discrepancy is also the limitation in the classification of the ResNET model, which reaches a final accuracy of around 90% when compared to the ground truth database. 1 Literature are administrative data taken from [START_REF] Dari | Exploiting high-resolution remote sensing soil moisture to estimate irrigation water amounts over a mediterranean region[END_REF] for "Canal d'Aragó i Catalunya" and from [START_REF] Rubió | El regadiu a lleida. resultats de diverses avaluacions a la zona regable dels canals d'urgell (lleida)[END_REF] for "Canals d'Urgell" and "Canal de Pinyana".

V. CONCLUSION

A key missing information for irrigation management and for hydrological studies over irrigated regions is the precise knowledge, at field level, of the different irrigation systems installed in each field, and the trends and changes of these systems over different years. This study has provided for the first time a method to classify irrigation systems (flood, sprinkler, drip) and not irrigated fields, using remotely sensed time-series at sub-field scale resolution. Key hydrological variables were used as inputs for the classification of irrigation systems through the use of different ML models.

Two main hydrological variables, ET a and SM at 20 m led to the best performance in the classification of irrigation systems when combined together, regardless of the ML model used for the classification. This result is indicative of the usefulness of these datasets in providing complementary information: SM directly detects large surface soil wetting, thus easily detecting large uses of water, as in flood and sprinkler irrigation. ET a is instead able to detect fields that are irrigated with drips, which does not create dramatic changes in the soil water content but keeps the plant at high ET a levels (close to potential evapotranspiration) with respect to non-irrigated fields.

An initial experiment was run in order to verify if the difference in crop types was interfering with the prediction of irrigation systems. Results showed that crop type does not interfere with the task of classifying irrigation systems, and similar performances are reached when using a general model for all the crops or multiple models specialized by crop types. These proved the feasibility of classifying systems of irrigation using only hydrological variables from remote sensing.

Among the three ML models tested, ResNET showed the best performance for all the metrics used for this classification task. ResNET is the only deep neural network model proposed for this study, and its architecture was shown to be more suitable for detecting more complex variations from the analyzed data. This characteristic is an advantage since irrigation practices strongly vary, both in timing and amount among different fields, even for fields that grow similar crops and employ the same irrigation system.

Finally, we compared a map of irrigation systems derived from the ResNET model using ET a and SM of the year 2020 against the irrigation systems map provided by the SIGPAC-DUN administrative catalog. Results at the district level showed a general agreement in the percentages of irrigation systems, even though the ResNET map appears to classify more fields with more efficient irrigation systems (drip and sprinkler irrigation). These differences are present since the ResNET map delivers a more updated depiction of irrigation systems compared to the map from SIGPAC-DUN, and captures the ongoing conversion from flood to more efficient irrigation systems.

This work represents the first study dedicated to the automatic detection at the field level of irrigation systems from satellite remotely sensed data. It shows how this is achievable with good accuracy when applied to semi-arid areas. Over semi-arid areas, the low cloud cover allows for high availability of thermal and optical satellite data, and irrigated and non-irrigated areas are easily distinguishable due to the marked difference in soil water availability from the surrounding dryland areas.

This study only focused on semi-arid regions since they are areas where water availability is a topic of increasing concern, and where optimization in the use and distribution of water for agriculture can bring the most noticeable improvements. One research avenue could be to apply the proposed methodology to more temperate areas. This application is expected to be challenging for two main reasons: first, in those areas crop phenology is more similar between irrigated and non-irrigated crops [START_REF] Pageot | Detection of irrigated and rainfed crops in temperate areas using sentinel-1 and sentinel-2 time series[END_REF], secondly, there will be larger gaps in the time-series due to more frequent cloud cover, that will reduce the amount of information usable by the model.

Other Mediterranean regions are expected to be well suited for applying this methodology, where transfer learning techniques could be explored. This approach will require a small amount of local training data, which will only be used to tailor the weights of the last layer of the pre-trained ResNET model in order to improve accuracies for the specific region of interest. Unsupervised learning could also be explored as a solution that avoids ground truth data collections, but which is often less performing than supervised methods.

Fig. 1 .

 1 Fig. 1. Depiction of the study area (black square in the top right figure) with eight of its biggest irrigation districts. The bottom figure depicts the spatial distribution of the fields sampled during the field campaign. The number above each circle represents the number of fields clustered by proximity, and the color goes from green to orange to indicate the increasing quantity of clustered fields.

Fig. 2 .

 2 Fig. 2. Confusion Matrix representing the discrepancies between the Ground Truth dataset collected during the field campaign and SIGPAC-DUN information.

  This article has been accepted for publication in IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing. This is the author's version which has not been fully e content may change prior to final publication. Citation information: DOI 10.1109/JSTARS.2022.3222884 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/

  was used to estimate the latter. The theoretical base of the SW model was provided by the Penman-Monteith energy combination equation, which has two parts: one for the soil surface and another for the plant surface. ET pot was computed with the SW model by setting a minimum stomatal resistance value of 100 sm -1 .SM was created from the disaggregation of low-resolution original data employing the DisPATCh algorithm.DisPATCh uses a semi-empirical soil evaporative efficiency model and a linearized relationship between Soil Evaporative Efficiency (SEE) and low resolution SM to perform the disaggregation[START_REF] Merlin | Disaggregation of smos soil moisture in southeastern australia[END_REF]. Differently from the classical version of the algorithm, a modification of DisPATCh for areas under high vegetation cover was added to the classic DisPATCh algorithm, as proposed by Ojha et al., 2021[START_REF] Ojha | Extending the spatio-temporal applicability of dispatch soil moisture downscaling algorithm: A study case using smap, modis and sentinel-3 data[END_REF]. For this study, SM was retrieved from the disaggregation of the original SMAP enhanced L2 SM product (L2_SM _P _E) gridded at 9 km. The Sentinel missions provided the high resolution optical and thermal data: N DV I maps at 20 m were extracted from the combination of bands 4 and 8A of the MSI instrument from Sentinel-2, while thermal maps were retrieved from the Sentinel-3 mission and sharpened at 20 m using Sentinel-2 reflectances bands. A Digital Elevation Map (DEM ) at 30 m from the Shuttle Radar Topography Mission was also used in order to account for topographic effects during the disaggregation process.

Fig. 3 .

 3 Fig. 3. LAI time-series from Sentinel-2, smoothed (with Hampel filter) and linearly interpolated to be used for the crop-type classification. A simple peak-detection algorithm was employed to differentiate the crop types. The bold blue lines represent the numerical average of all the time-series inside each category.

3 )

 3 Model and variable selection: After verifying the capability of the model to classify irrigation systems, a final experiment was designed to investigate which variables are most suited for the classification of irrigation systems and which model performs better for this classification task. The classification of irrigation systems was performed with both L2 and L4 variables, training the different models with each variable separately and with a combination of them. All the different ML models proposed in this study were used for the comparison: the two classic ML learning models (time-series Forest and Rocket) applied to both L4 and L2 variables and trained at field level (one

Fig. 4 .

 4 Fig. 4. Temporal Importance Curves from the Time Series Forest model trained over specific crop types and corresponding input time-series of ETa and SM .

Fig. 5 .

 5 Fig. 5. Accuracies of different approaches for the classification of irrigation systems (median and standard deviation of 10 different runs). Multiple ML models and input variables are employed. Input variables are grouped as (a) Level 2 variables, representing the initial satellite products, and (b) Level 4 variables, representing variables produced from an elaboration of satellite data with different models.

Fig. 6 .

 6 Fig. 6. Confusion Matrix at field level (grouped with majority voting of the pixels contained in the field) and post-processed. The matrix shows the results of one run of the ResNET model using ETa and SM as inputs. Each cell shows a percentage over the total True label and number of Fields. Precision and Recall are also shown on the side.

Fig. 7 .

 7 Fig. 7. Irrigation systems maps (a) as delivered by SIGPAC-DUN [68] and (b) as created by the ResNET model with ETa and SM time-series for 2020, with 3 zoomed-in areas (c) for a visual comparison, showing an increasing zoom level to verify spatial consistency of the product in different irrigation districts.

Fig. 8

 8 Fig. 8 shows a direct comparison between the different systems of irrigation as estimated by SIGPAC-DUN and predicted by the ResNET model used for this study. It is noticeable how the ResNET model consistently predicts a lower percentage of traditional irrigation by flooding, and in almost all cases (but for Canals d'Urgell and Algerri-Balaguer) a decrease in the percentage of non-irrigated fields. This discrepancy was expected, since SIGPAC-DUN

Fig. 8 .

 8 Fig. 8. Comparison of systems of irrigation estimated by SIGPAC-DUN (S) and predicted by ResNET (R) for the 8 irrigation districts considered for this study.

  This article has been accepted for publication in IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing. This is the author's version which has not been fully e content may change prior to final publication. Citation information: DOI 10.1109/JSTARS.2022.3222884

This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ et al., 2021

  This article has been accepted for publication in IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing. This is the author's version which has not been fully e content may change prior to final publication. Citation information: DOI 10.1109/JSTARS.2022.3222884

This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/

Table I summarizes all the samples collected during the field campaign divided by irrigation systems, as detected during the campaign, and crop types, as indicated in SIGPAC-DUN.

 I This article has been accepted for publication in IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing. This is the author's version which has not been fully e content may change prior to final publication. Citation information: DOI 10.1109/JSTARS.2022.3222884 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/

TABLE I INITIAL

 I CATALOG OF THE FIELDS COLLECTED DURING THE FIELD CAMPAIGN IN 2020. CROP TYPES ARE INITIALLY TAKEN FROM SIGPAC-DUN, WHILE IRRIGATION SYSTEMS COME FROM THE LABELLING ACTIVITY DURING THE FIELD CAMPAIGN.

	IRRIGATION SYSTEM	CROP TYPE		TOTAL
		FRUIT and NUT TREES 78
	DRIP	VINEARD	12	101
		OLIVE	11
		WINTER CEREALS	53
	SPRINKLER	MAIZE	11	71
		ALFALFA	7
		WINTER CEREALS	41
	FLOOD	MAIZE ALFALFA	13 10	82
		FRUIT and NUT TREES 18
		WINTER CEREALS	40
	NOT IRRIGATED	FRUIT and NUT TREES 13 VINEARD 7	77
		OLIVE	17
	Despite being very detailed, only the latest version of the dataset, corresponding to the year 2021, contains
	information about two important aspects of the fields: which secondary crop type (if present) was cultivated during

TABLE II OVERVIEW

 II OF REMOTE SENSING VARIABLES USED IN THIS STUDY.

		Variable name Spatial Resolution Details	References
		SMAP SM	33 Km	SMAP L3 soil moisture enhanched product.	[69]
	L2 Variables		(gridded at 9 Km)
		S3 LST	1 Km	Sentinel-3 daily Land Surface temperature.	[70]
		S2 NDVI	20 m	Sentinel-2 NDVI from Band 4 (VIS) and Band 8A (NIR).	-
		ETa	20 m	Obtained with the Priestley-Taylor Two-Source Energy Balance (TSEB-PT)	[47]
				using Copernicus-based inputs.
	L4 Variables	SM	20 m	DisPATCh algorithm to disaggregate SMAP SM using S3 LSt and S2 reflectances. [57]
		Ks	20 m	Ratio of ETa and potential evapotranspiration (ETpot).	[47]
		S2 LAI	20 m	Calculated mainly from Sentinel-2 reflectances using PROSAIL	[71]
				radiative transfer model and Neural Networks algorithm.

  This article has been accepted for publication in IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing. This is the author's version which has not been fully e content may change prior to final publication. Citation information: DOI 10.1109/JSTARS.2022.3222884

This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/

  This article has been accepted for publication in IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing. This is the author's version which has not been fully e content may change prior to final publication. Citation information: DOI 10.1109/JSTARS.2022.3222884

This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/

TABLE III SUMMARY

 III OF THE TOTAL NUMBER OF FIELDS COLLECTED DURING THE FIELD CAMPAIGN PERFORMED IN 2020. FIELDS ARE DIVIDED BY CROP TYPE AND IRRIGATION SYSTEM, TOTALS ARE SHOWN ON THE SIDE TOGETHER WITH THE SUM OF ALL THE 20M PIXELS CONTAINED IN THE FIELDS COLLECTED.

	IRRIGATION SYSTEM CROP TYPE	Number of Fields 2018 2019 2020	FIELD LEVEL	PIXEL LEVEL
		FRUIT and NUT TREES 78	78	78	234	24201
	DRIP	VINEYARD	12	12	12	36	4599
		OLIVE	11	11	11	33	3201
		MAIZE	8	8	8	24	10950
	SPRINKLER	DOUBLE CROPS	55	56	56	167	43849
		ALFALFA	7	7	7	21	3777
		WINTER CEREALS	9	9	9	27	444
		MAIZE	14	14	13	41	1322
	FLOOD	DOUBLE CROPS	32	33	33	98	5859
		ALFALFA	9	9	9	27	2733
		FRUIT and NUT TREES 18	18	18	54	1734
		WINTER CEREALS	40	36	40	116	27584
	NOT IRRIGATED	FRUIT and NUT TREES 13 VINEYARD 7	13 7	13 7	39 21	1578 867
		OLIVE	17	17	17	51	6231
	TOTAL		330	328	331	989	138929
	C. ML Models

280

Time Series Classification (TSC) is a specific machine learning class of algorithms that classifies data taking 281 into account their ordered structure. This class of algorithms is chosen in this study since timing is a key element 282

  Table IV summarizes the reason of the selection of these particular models. This article has been accepted for publication in IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing. This is the author's version which has not been fully e content may change prior to final publication. Citation information: DOI 10.1109/JSTARS.2022.3222884

284

This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/

TABLE IV OVERVIEW

 IV OF THE MODELS USED IN THIS STUDY.

  This article has been accepted for publication in IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing. This is the author's version which has not been fully e content may change prior to final publication. Citation information: DOI 10.1109/JSTARS.2022.3222884

This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ ResNET (Residual networks) is a Deep neural Network model which was designed initially for computer vision task and adapted to time-series classification

  This article has been accepted for publication in IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing. This is the author's version which has not been fully e content may change prior to final publication. Citation information: DOI 10.1109/JSTARS.2022.3222884

This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/

  This article has been accepted for publication in IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing. This is the author's version which has not been fully e content may change prior to final publication. Citation information: DOI 10.1109/JSTARS.2022.3222884 time-series per field), and ResNET, the DNN model, applied only to L4 variables at pixel level, since the low spatial resolution of the L2 variables did not provide a large enough dataset for training and testing of this model.Each model and each variable were trained and evaluated 10 times, changing each time the train and test datasets' distributions and the models' random initial weights. Median and standard deviation of the overall accuracies were Table V summarizes the accuracy retrieved from the simulations. The experiment was performed to assess the accuracy of SM , ET a , and both variables together, SM +ET a . Results from these initial experiments show that a general model (last column of TableV) has comparable accuracy with respect to the aggregated accuracy (second to the last column) of multiple models trained separately for each crop type (with a small difference of ∆ = 2.16%).Using a general model, trained on every crop type has a noticeable advantage of not requiring a crop type map, which makes the approach more versatile since it can be adapted to areas where crop types are not known or where there are different crop types from the ones analyzed in this study.

	used as a comparison metric.
	IV. RESULTS AND DISCUSSION
	A. Influence of Crop Types

This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/

TABLE V ACCURACIES

 V

	Variables	Winter	Maize	Double	Crop types Alfalfa Fruit & Nut	Olives	Vineyards	Aggregated	RESULTS General
		Cereals		Crops		Trees			Models		Model
	ETa-TSEB	94.00%	76.36%	90.00%	58.33%	92.59%	100.00% 85.56%	89.94% ± 2.23	86.11% ± 3.20
	ETa-TSEB cropped	93.20%	77.27%	88.18%	48.33%	92.96%	100.00% 84.44%	89.07% ± 1.70	-
	SM Dispatch	90.80%	79.09%	86.36%	58.33%	90.74%	99.23%	86.67%	88.02% ± 2.33	83.77% ± 1.85
	SM Dispatch cropped 89.60%	82.73%	85.23%	58.33%	90.37%	96.92%	82.22%	87.22% ± 1.81	-
	ETa+SM	93.20%	82.73% 92.95% 58.33% 93.70%	100.00% 88.89%	91.60% ± 2.11	89.44% ± 2.97
	ETa+SM cropped	93.20%	81.82%	92.27%	53.33%	94.26%	100.00% 86.67%	91.23% ± 1.63	-

OF THE CLASSIFICATION OF THE TIME-SERIES FOREST APPLIED TO DIFFERENT INPUT VARIABLES, WITH ANNUAL LENGTH OR CROPPED TO THEIR VEGETATIVE PERIODS. FINAL RESULTS ARE PRESENTED IN THE LAST TWO COLUMNS, BETWEEN THE AGGREGATED ACCURACIES OF THE MODELS TRAINED FOR EACH CROP TYPES AND A GENERAL MODEL TRAINED WITH ALL CROP TYPES TOGETHER. HIGHEST ACCURACY FOR EACH COLUMN IS PRESENTED IN BOLD.

  This article has been accepted for publication in IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing. This is the author's version which has not been fully e content may change prior to final publication. Citation information: DOI 10.1109/JSTARS.2022.3222884

This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/

  This article has been accepted for publication in IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing. This is the author's version which has not been fully e content may change prior to final publication. Citation information: DOI 10.1109/JSTARS.2022.3222884

This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/

  This article has been accepted for publication in IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing. This is the author's version which has not been fully e content may change prior to final publication. Citation information: DOI 10.1109/JSTARS.2022.3222884 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/

  Table VI for experiments performed with ET a and SM as input variables. article has been accepted for publication in IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing. This is the author's version which has not been fully e content may change prior to final publication. Citation information: DOI 10.1109/JSTARS.2022.3222884 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/

	This

TABLE VI CLASSIFICATION

 VI METRICS FOR THE THREE MODELS WHEN USED WITH ETa AND SM . MEAN AND STANDARD DEVIATION DERIVED FROM 10 DIFFERENT RUNS OBTAINED SHUFFLING THE TRAIN AND TEST DATASET.

			MODELS
	METRICS (%)	tsForest	ROCKET	ResNET
	Accuracy	85.29 ± 2.41 87.56 ± 2.95 90.10 ± 2.70
	Average Precision 85.43 ± 2.53 88.80 ± 3.12 90.33 ± 2.78
	Average Recall	84.76 ± 2.51 86.81 ± 3.17 90.02 ± 2.76

  This article has been accepted for publication in IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing. This is the author's version which has not been fully e content may change prior to final publication. Citation information: DOI 10.1109/JSTARS.2022.3222884

This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/

  This article has been accepted for publication in IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing. This is the author's version which has not been fully e content may change prior to final publication. Citation information: DOI 10.1109/JSTARS.2022.3222884

This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/

  This article has been accepted for publication in IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing. This is the author's version which has not been fully e

content may change prior to final publication. Citation information: DOI 10.1109/JSTARS.2022.3222884 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/

TABLE VII PERCENTAGES

 VII OF THE AREA UNDER DIFFERENT SYSTEMS OF IRRIGATION, DIVIDED PER IRRIGATION DISTRICT. THE PERCENTAGES IN LITERATURE ARE TAKEN FROM LITERATURE AND COMPARED WITH THE LATEST SIGPAC-DUN DATASET [68] AND THE MAP RETRIEVED WITH THE RESNET MODEL.

	Irrigation District	IRR systems	Literature 1 SIGPAC-DUN ResNET
		FLOOD	90 %	81 %	55 %
	Canals d'Urgell	DRIP	4 %	10 %	11 %
		SPRINKLER	2 %	5 %	24 %
		NOT IRRIGATED 0 %	4 %	11 %
		FLOOD	79 %	50 %	27 %
	Canal de Pinyana	DRIP	10 %	33 %	34 %
		SPRINKLER	10 %	13 %	38 %
		NOT IRRIGATED 0 %	3 %	2 %
		FLOOD	18 %	14 %	2 %
	Canal d'Aragó i Catalunya	DRIP	28 %	37 %	32 %
		SPRINKLER	54 %	46 %	65 %
		NOT IRRIGATED 0 %	2 %	1 %

  This article has been accepted for publication in IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing. This is the author's version which has not been fully e content may change prior to final publication. Citation information: DOI 10.1109/JSTARS.2022.3222884

This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/

  This article has been accepted for publication in IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing. This is the author's version which has not been fully e content may change prior to final publication. Citation information: DOI 10.1109/JSTARS.2022.3222884

This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/

This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/

This article has been accepted for publication in IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing. This is the author's version which has not been fully e content may change prior to final publication. Citation information: DOI 10.1109/JSTARS.2022.3222884 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/

ACKNOWLEDGEMENTS

The first author received the grant DIN2019-010652 from the Spanish Education Ministry (MICINN) and DI-2020-093 from the Catalan Agency of Research (AGAUR). We are grateful for the financial support given by the ACCWA project, funded by the European Commission Horizon 2020 Program for Research and Innovation, in the context of the Marie Skłodowska-Curie Research and Innovation Staff Exchange (RISE) action under the grant agreement No. 823965, and by the PRIMA ALTOS project (No. PCI2019-103649) of the Ministry of Science, Innovation and Universities of the Spanish government. The PRIMA IDEWA project is also acknowledged.