Evaluation of Remote Sensing-Based Irrigation Water Accounting at River Basin District Management Scale

: The Water Framework Directive in Europe requires extending metering and water abstraction controls to accurately satisfy the necessary water resource requirements. However, in situ measurement instruments are inappropriate for large irrigation surface areas, considering the high investment and maintenance service costs. In this study, Remote Sensing-based Irrigation Water Accounting (RS-IWA) (previously evaluated for commercial plots, water user associations, and groundwater water management scales) was applied to over 11 Spanish river basin districts during the period of 2014–2018. Using the FAO56 methodology and incorporating remote sensing basal crop coefficient time series to simulate the Remote Sensing-based Soil Water Balance (RS-SWB), we were able to provide spatially and temporally distributed net irrigation requirements. The results were evaluated against the irrigation water demands estimated by the Hydrological Planning Offices and published in the River Basin Management Plans applying the same spatial (Agricultural Demand Units and Exploitation Systems) and temporal (annual and monthly) water management scales used by these public water managers, ultimately returning ranges of agreement (r 2 and d r ) (Willmott refined index) of 0.79 and 0.99, respectively. Thus, this paper presents an operational tool for providing updated spatio-temporal maps of RS-IWA over large and diverse irrigation surface areas, which is ready to serve as a complementary irrigation water monitoring and management tool.


Introduction
AQUASTAT, the global information system on water resources, quantifies and differentiates human water abstractions based on three use categories: municipal (11%), industrial (19%), and agricultural (70%) [1]. At the continental level, Africa and Asia have an agricultural proportion of 80%, while the water used in Europe for agriculture makes up 21%. However, the same balance changes drastically over Mediterranean regions, where 80% of freshwater is abstracted for agricultural purposes [2].
Considering "the need for action to protect Community waters in qualitative as well as in quantitative terms" and the fact that "Waters in the Community are under increasing pressure from the continuous growth in demand for sufficient quantities of good quality water for all purposes", the European Commission created DIRECTIVE 2000/60/EC [3], commonly known as the Water under the assumed criteria of homogeneity, simplicity, replicability, and reliability to estimate the irrigation water demands over large and diverse irrigation surface areas. Hence, the Landsat 8 and Sentinel-2 satellites were selected, as they offer adequate spectral, temporal, and spatial resolutions for crop monitoring and ET c estimations. Particularly, we used the reflectance-based basal crop coefficient approach, which is linearly related with the VI, as the K cb is later incorporated into a dual crop coefficient approach to derive the ET c adj under a soil water balance scenario, which ultimately provides spatially and temporally distributed maps of the net irrigation requirements (NIR). The RS-SWB was previously evaluated against different ground truths at different spatial and temporal scales, such as for commercial plots with volumetric meters and water user associations using irrigation water management data [33]. Evaluations using the national piezometric network after determining the groundwater withdrawals from the irrigation surface area greatly depend on these groundwater resources [34], as does using the annual irrigation water demands provided by the RBMPs at an RBD spatial scale in mainland Spain [35]. The next RS-IWA evaluation was performed using the same irrigation water management scales that the Hydrological Planning Offices take into consideration, including temporal (monthly and annual time scale) and ADU and ES spatial scales.

Study Area
The irrigated surface areas of 11  . With a total surface area of about 440,000 km 2 (87% of the total national surface), these 11 RBDs contain most Spanish irrigated surface areas. Hence, this extension includes great diversity and complexity in addressing diverse types of crops, irrigation systems, and agronomic management practices. The different scales of work represented in Figure 1 align with the irrigation water demands estimated by the Hydrological Planning Offices and published in the corresponding RBMPs, according to different temporal (annual and monthly) and spatial (ADU, ES, and (in the Duero RBD) Elemental Demand Unit, EDU) irrigation water management scales.
The most detailed analyses at both management scales were conducted in the Duero and Júcar RBDs and correspond to the ADU spatial management scale at a monthly and annual temporal frequency. We justify their choice over other RBDs by the diverse types of irrigated crops developed therein (herbaceous versus citrus and vineyards). In the Duero RBD, grain cereals stand out, with maize grown in about 19% of the irrigable area, alongside barley (≈17%) and wheat (≈16%), followed by forage crops such as alfalfa (≈6%) and industrial crops such as sunflower (≈5%). In the Júcar RBD, citrus occupies almost the half of the irrigated area (≈43%), followed by maize and rice (≈20%) and, finally, vineyards and horticultural crops, with about 9%. Moreover, both RBDs present different geographical conditions (eastern peninsular and Atlantic slopes versus western and Mediterranean slopes, respectively) and, consequently, different climatic characteristics according to the Köppen-Geiger characterization [36] (a temperate climate with a dry and temperate summer, Csb, in the Duero RBD, unlike the cold steppe climates, Bsk, and a temperate climate with a dry and warm summer, Csa, in the Júcar RBD). Finally, we also added transparency and accessibility to the published data (databases and shapefiles) using their webGIS viewers. In parallel, a second analysis including a larger surface area and an annual temporal scale was conducted in 69 ES spatial scales contained in the 11 RBDs mentioned above. Finally, the EDU spatial division already accounted for in the Duero RBD was also analysed at an annual temporal scale.

HidroMORE for Large Irrigation Areas
The RS-SWB (described in the following Section 2.3) was run using the latest version of the HidroMORE operational application to implement parallel processing with multithreading for running simulations over large areas [37]. Previously applied in different areas such as research [33][34][35][38][39][40], doctoral theses [41][42][43], and several research projects (SIRIUS EU project GA nº 262902, http://siriusgmes.es/; DIANA EU project GA nº 703109, https://diana-h2020.eu/en/; CAPRA national Chilean project http://www.inia.cl/proyecto/502145/), HidroMORE spatially and temporally distributes the FAO56 model and integrates multispectral remote sensing data to provide the crop parameter evolution of the linear relationships between VI and both K cb and the fraction of soil surface covered by vegetation (f c ) [41]. Hence, using a pixel-by-pixel spatial scale and a daily temporal step (extended to the desired time range, such as x weeks, x months, x years, etc.), maps of the soil water balance variables are obtained as outputs from the model (i.e., ET o , P, ET c adj , and NIR), thereby producing an irrigation thematic cartography time series (Scheme 1). The set of satellite images used determines the total surface simulated, while the smallest spatial unit of work is the pixel size from the input satellite images. To implement the RS-SWB, HidroMORE uses the following input data: (a) VI time series, (b) a map of irrigated crops, (c) a map of soil types according to their hydrological properties, and (d) agroclimatic databases (daily P and ET o records).

HidroMORE for Large Irrigation Areas
The RS-SWB (described in the following Section 2.3) was run using the latest version of the HidroMORE operational application to implement parallel processing with multithreading for running simulations over large areas [37]. Previously applied in different areas such as research [33][34][35][38][39][40], doctoral theses [41][42][43], and several research projects (SIRIUS EU project GA nº 262902, http://sirius-gmes.es/; DIANA EU project GA nº 703109, https://diana-h2020.eu/en/; CAPRA national Chilean project http://www.inia.cl/proyecto/502145/), HidroMORE spatially and temporally distributes the FAO56 model and integrates multispectral remote sensing data to provide the crop parameter evolution of the linear relationships between VI and both Kcb and the fraction of soil surface covered by vegetation (fc) [41]. Hence, using a pixel-by-pixel spatial scale and a daily temporal step (extended to the desired time range, such as x weeks, x months, x years, etc.), maps of the soil water balance variables are obtained as outputs from the model (i.e., ETo, P, ETc adj, and NIR), thereby producing an irrigation thematic cartography time series (Scheme 1). The set of satellite images used determines the total surface simulated, while the smallest spatial unit of work is the pixel size from the input satellite images. To implement the RS-SWB, HidroMORE uses the following input data: (a) VI time series, (b) a map of irrigated crops, (c) a map of soil types according to their hydrological properties, and (d) agroclimatic databases (daily P and ETo records).

VI Temporal Series
The normalized difference vegetation index (NDVI) was the VI selected to monitor crop evolution (see the justification in Section 2.3). While for the years of 2014 and 2015, the Landsat 8 satellite was used to process 938 scenes, during 2016-2018, the Sentinel-2 A and B satellites were selected to process 13,694 granules. Images were atmospherically corrected using an absolute normalization process [44], according to the selection of two different invariant surfaces to determine the NDVI linear regression values [45]. In addition, the inter-calibration of both satellites was achieved using coincident days to ensure the compatibility of the satellites with the NDVI time series over all the years. Moreover, as the selection of VI images is critical for correct irrigated crop monitoring, two criteria were followed: (i) selecting initial and final dates that covered the entire vegetative crop period and (ii) selecting images with the absence of clouds/shadows or images in Scheme 1. Overview framework for Remote Sensing-based Irrigation Water Accounting (RS-IWA).

VI Temporal Series
The normalized difference vegetation index (NDVI) was the VI selected to monitor crop evolution (see the justification in Section 2.3). While for the years of 2014 and 2015, the Landsat 8 satellite was used to process 938 scenes, during 2016-2018, the Sentinel-2 A and B satellites were selected to process 13,694 granules. Images were atmospherically corrected using an absolute normalization process [44], according to the selection of two different invariant surfaces to determine the NDVI linear regression values [45]. In addition, the inter-calibration of both satellites was achieved using coincident days to ensure the compatibility of the satellites with the NDVI time series over all the years. Moreover, as the selection of VI images is critical for correct irrigated crop monitoring, two criteria were followed: (i) selecting initial and final dates that covered the entire vegetative crop period and (ii) selecting Remote Sens. 2020, 12, 3187 7 of 28 images with the absence of clouds/shadows or images in which the presence of clouds/shadows could be correctly eliminated through masks and classification algorithms.
For the large and diverse irrigated surfaces monitored, the entire year was selected as the temporal date range. The final NDVI image selection after automatic cloud/shadow mask application was carried out through human supervision using the evolution of the VI time charts combined with colour RGB images loaded onto the SPIDERwebGIS ® platform, which is specialized in handling such time series. The cloud/shadow classification algorithms varied according to the reference satellite used: (a) the Quality Assessment band for Landsat 8 (CFmask algorithm) and (b) the Sen2Cor algorithm for Sentinel-2 (http://step.esa.int/main/third-party-plugins-2/sen2cor/). The reliability of these algorithms was determined from different validation tests. The accuracy of CFmask in detecting clouds and shadows was around 90% [46], while for Sen2Cor, the precision results were around 80%, both for the high probability cloud classes and their shadows [47]. Despite these high values of cloud/shadow identification, difficulties remained in the correct classification of fine and high clouds [47], so time series monitoring was used to identify discontinuities and thus remove useless images.

Map of Irrigated Crops
The RS-SWB model applied by HidroMORE requires an annual map with the spatial distribution of the different irrigated crop classes to apply the methodology. Supported by official cartographic sources that facilitate the spatial delimitation assessment of irrigated areas (Common Agricultural Policy and ADU official cartography), these maps were constructed using a supervised multitemporal classification based on the VI thresholds, which require phenological cycle knowledge of the crops to be classified [48,49]. The final irrigated crop classes produced were spring, summer, spring-summer (double harvests), autumn-winter, and woody crops such as vineyards, olives, citrus, and fruit trees. Moreover, since crops with similar phenological and ground cover patterns have similar VI evolution and crop water requirements [48], the collected cycle was then used for the selection of threshold values for each type of irrigated crop class by VI time series at different phenological phases; these values are freely available on the project description website (see the Introduction section) and in the final report [50]. Hence, using these crop group classes and their VI profiles, we used the decision-tree approach on certain dates for crop class discrimination. The adopted classification facilitated the discrimination of the wide range of cultivated crops, as the crops were grouped according to the period within the year in which their phenological development takes place. Finally, the national irrigated surface area accounts according to previously described remote sensing approaches were proven and included by the Ministry for the Ecological Transition in an official report about the actual state of RBMPs at the national level [28].
HidroMORE uses specific information about each crop class to determine the soil water balance. The information collected in tables refers to the following parameters collected in FAO 56: the evapotranspiration depletion factor (p), the linear relationships of VI-K cb and IV-f c , the maximum and minimum K c values, the maximum and minimum root depth, the maximum and minimum VI, the water stress coefficient (K s ), and the irrigation period and irrigation dose per irrigation event.

Soil Hydrological Maps
The RS-SWB model applied by HidroMORE requires an annual spatially distributed map of soil hydrological properties. There are two sources for this map at the national scale [51,52]. However, none of these maps are digital. Hence, we used the map and database provided by the "European Soil Database Derived data" [53], which is mainly based on the "European Soil Database v2.0" [54] (E 1/1,000,000).
This resource distributes open-source maps, including the total available water (TAW) and the maximum soil depth for root development, among other edaphological parameters.
These inputs are the most important for the necessary soil parameters to determine the cumulative net water depletion by roots in combination with the wilting point (θ WP ) and field capacity (θ FC ). The latter parameters are related with the provided TAW, so they are obtained from such information Remote Sens. 2020, 12, 3187 8 of 28 regarding the soil type knowing the TAW values. In addition, the soil evaporation from the exposed soil was modelled by applying the description of Torres and Calera [55] for the soil evaporation reduction coefficient and its correction factor under high evaporative rates, considering a depth of the surface soil layer of 0.1 m, as FAO56 recommends when the depth of the surface soil layer is unknown [10]. All this information is collected in tables that HidroMORE uses to run the RS-SWB.

Agrometeorological Databases
The RS-SWB model used with HidroMORE requires daily P and ET o data recorded by agrometeorological stations to spatially distribute the local information using the spatial interpolation (inverse distance weighting) of this local data. SIAR (the agrometeorological stations of the national network from the Agroclimatic Information System for Irrigation (http://www.siar.es/) of the Spanish Ministry of Agriculture, Food, and Fisheries) was used in combination with two regional agrometeorological networks (https://ruralcat.gencat.cat/web/guest/agrometeo.estacions and https://www.larioja.org/agricultura/es/informacion-agroclimatica/red-estaciones-agroclimaticassiar). The annual average number of agrometeorological stations used was 536 ( Figure 1a).

Remote Sensing-Based Irrigation Water Accounting for Large Areas
Remote Sensing-based Irrigation Water Accounting is the approach used to obtain the daily net irrigation requirements per pixel unit after applying a Remote Sensing-based Soil Water Balance to then spatially and temporarily aggregate these requirements to the desired work scale. The RS-SWB uses the well-known FAO56 methodology [10]. As expressed in Equation (1), this approach estimates the daily cumulative depth of evapotranspiration (depletion) from the root zone (D r ) that the crop extracts through adjusted crop evapotranspiration (ET c adj ) in irrigation areas to subsequently apply a net irrigation depth (I) for the scenarios in which the precipitation (P) does not fully satisfy the crop water requirements. In parallel, as expressed in Equation (2) by the dual crop coefficient methodology [11], the ET c adj is estimated for the scenarios in which crops do not completely cover the ground surface (woody crops, horticultural crops, or herbaceous ones in the first stages) or require high irrigation frequency events (e.g., horticultural crops). Hence, the described crop development scenarios are common for large and diverse irrigation surface areas such as the present study area.
Avoiding the use of tabulated values for the basal crop coefficient (K cb ) or the fraction of soil surface covered by vegetation (f c ) in large and diverse irrigation surface areas, the RS-SWB uses remote sensing capabilities to monitor crop evolution and obtain these crop parameters through dense VI time series [56,57] empirically derived from linear relationships (VI-K cb act and VI-f c act ). Besides focusing on operational applications, independent crop linear relationships are recommended [12,[57][58][59]. For these reasons, and considering the ongoing debate on the use of independent VI linear relationships to estimate biophysical parameters, we selected the NDVI as the operative VI.
In this sense, Equation (3) expresses the linear relationship of the NDVI-K cb selected [60]. Initially developed for vineyards, the NDVI-K cb does not strongly differ from other published and evaluated linear relationships for maize [61,62] or sorghum [63] and was positively evaluated at different locations, such as vineyards in Chile [26] and Australia [64]; for different vegetation covers, such as savahan dehesa [65]; and in comparison with different VIs, such as the soil-adjusted vegetation index (SAVI) [66], in apple orchards [25]. In parallel, Equation (4) express the NDVI-f c linear relationship [67]. Originally developed for herbaceous crops such as maize and wheat, the NDVI-f c was applied to dry herbaceous and vineyards crops and as part of the RS-SWB to evaluate its reliability in soil water content estimations against ground truth data provided by ground sensors [40]. Moreover, the NDVI-f c was used in combination with other linear relationships to determine the spectral response evolution of Remote Sens. 2020, 12, 3187 9 of 28 maize and wheat, including the critical green-up step [68], and was evaluated to have good behaviour in comparison to ground measurements in vineyards [69].
Finally, at a daily frequency and pixel-based spatial unit scale, the RS-SWB inverts Equation (1) to obtain Equation (5). In this way, the NIR (mm/day) is estimated, producing thematic cartography for the temporally and spatially distributed balance applied. In this equation, the capillary rise (CR i ) and superficial runoff (RO) were discarded. On the one hand, FAO56 does not recommend the use of CR when the water table is below 1 m from the soil surface [10] (a minority situation at the scale of the Iberian Peninsula) [70,71]. On the other hand, the RO was not computed by the HidroMORE software, as it was originally conceived for flat irrigation areas and works only with vertical fluxes, not lateral ones.
The methodological approach described above was positively evaluated at different irrigation management scales. At an annual temporal scale, along with three consecutive campaigns, the RS-IWA was compared with the ground truth data from commercial plots cultivated with wheat, barley, and maize and at the water user association management scale with a wide range of herbaceous crops and woody crops (mainly vineyards) [33]. Furthermore, the same approach was evaluated at a monthly temporal scale over four consecutive years and over a large irrigation surface area (about 100,000 ha, including herbaceous and woody crops) that is largely dependent on groundwater resources. The monthly groundwater withdrawal time series estimated by the RS-IWA were positively compared with strategically selected piezometric measurements from a piezometric network managed by the public authorities [34].

RBMP Irrigation Water Accounting Estimation (RBMP-IWA)
RBMPs are focused (among other objectives) on satisfying water demands, thereby rationalizing their use in harmony with the environment and other natural resources [31]. Thus, the knowledge of irrigation water demands is crucial for water demand management. However, due to the large extent of the territory that needs to be monitored and the high costs of the installation and maintenance of volumetric counters [72][73][74], the river basin authorities in Spain commonly determine water demands through spatially based indirect measurement approaches at the ADU, which is required by the national legal framework [30][31][32] as the minimum common spatial unit for irrigation water management in all Spanish RBMPs. The Hydrological Planning Offices are the public authorities in charge of this task, and large-scale efforts are made to combine different data sources to obtain the best estimations under the available resources. The 11 RBMPs consulted [75][76][77][78][79][80][81][82][83][84][85] followed similar approaches to irrigation water demand accounting, consisting of several sequential steps.
Firstly, the public water managers estimate the irrigable surface areas using sources such as the late national agricultural census (which dates from 2009), the surface related to water rights registration in the national private water registry (the variable dates for private water manager acquisitions), the information provided by the Common Agricultural Policy GIS (usually for the period 2010-2014), or, in some cases that do not cover the whole territory, RBD internal reports (based on field visits, regional surveys, or remote sensing techniques). In the second step, the proportion of crops or crop groups and their spatial distribution that covers the estimated irrigable surface areas are estimated by the public water managers. Again, the national agricultural census, the Common Agricultural Policy GIS, and regional surveys are the most common sources consulted. Thirdly, once the irrigable surface area and the proportion and distribution of crops or crop groups are estimated, the public water managers estimate the net irrigation water demand by multiplying the irrigable surface area proportion of each crop or crop group by their crop water requirements. The latter parameter relies on different data sources, such as the national recommendations provided in the ministerial law [30], which offer data ranges for parameters related to RBD spatial segmentation and crop or crop group division (i.e., 2900-7000 m 3 /ha·year in the Ebro RBD or 3400-6400 m 3 /ha·year in the Duero RBD for maize and sorghum, or 500-3350 and 1000-2900 m 3 /ha·year for the previous RBD in the case of winter cereals). However, previous recommendations are commonly complemented by and locally adapted with RBD internal reports based on local Irrigation Advisory Services that estimate the crop water demands in pilot areas. The fourth and final step consists of converting the estimated net irrigation water demands into gross irrigation water demands to account for the water that should be diverted to different irrigable surface areas. In this step, which is primarily based on national recommendations [30] and internal RBD reports, the Hydrological Planning Offices determine the irrigation efficiencies for the three irrigation water supply levels of application, distribution, and transport (i.e., 0.85-0.90 for open channel efficiency in the case of transport and distribution levels or 0.75-0.80 for the sprinkler application level) or a global efficiency coefficient that groups all levels.
For the previous approach, each Hydrological Planning Office fixes a base year for the irrigable surface areas and its probabilistic spatial crop distribution estimation over each ADU and then applies the crop water requirements and the irrigation efficiencies, thereby obtaining the irrigation water demand over large and diverse areas for this particular base year (in Spain, for the second hydrological planning cycle, the year of 2012 is commonly used). These values are then extrapolated to the next hydrological planning periods (2015-2021, 2021-2027, and so on) such that a unique yearly irrigation water demand value per ADU (or twelve values in the case of a monthly time scale) is used to cover all the years (or months) during the whole hydrological planning period (i.e., 2015-2021). The Hydrological Planning Offices extrapolate these values from the base year by considering future increments or decrements in the irrigable surface areas' growth and/or planned investments for improving the irrigation efficiency of water user associations or irrigation water supply networks. As these changes do not quickly occur over such large and diverse areas, the differences between the estimated values for the base year and the following hydrological planning period are insignificant and only relevant for a few ADUs.
Note that the RBMPs group ADUs with similar characteristics into greater spatial scales called ESs, thereby allowing irrigation water management over greater scales. Besides the Duero RBD, a spatial unit smaller than the ADU is taken into consideration. Named the EDU, this unit is a subdivision of the ADU with common characteristics, such as the same water supply source.
Hence, considering the data available for the 11 RBMPs [75][76][77][78][79][80][81][82][83][84][85] at an annual temporal scale, the RBMPs publish the irrigation water demands of the ADU and ES and even, in some cases, at a monthly temporal scale for the ADUs. This public information was used at two spatial scales (ES and ADU) and, consequently, at two temporal scales (annual and monthly) in this research to compare the estimated RBMP-IWA with the results produced by RS-IWA.

RS-IWA Evaluation Approach in Comparison with RBMP-IWA at River Basin Spatial Scales
The daily NIR values estimated by RS-SWB at the pixel spatial scale were then spatially and temporally aggregated (RS-IWA) into the required water management scales that the RBMPs publish: monthly and annual temporal scales and EDU, ADU, and ES spatial scales (Section 2.4). Thus, by employing GIS spatial statistics (using Quantum GIS version 3.10 La Coruña), we obtained spatially and temporally distributed spatial statistics of the RS-IWA that were compared with RBMP-IWA (Scheme 1). The official cartography was directly downloaded from the RBD GIS viewers or RBD websites.
For The statistical values used for the irrigation water accounting comparisons are the minimum (Min), maximum (Max), average (Av), standard deviation (σ), regression coefficient (r 2 ), and Willmott refined index (d r ) [86], which are commonly used in hydrological model evaluations. Both the absolute and relative values of the Root Mean Square Error (RMSE and RE RMSE ) and the Mean Absolute Error (MAE and RE MAE ) are also presented. However, although the RMSE is commonly used by researchers, the distribution and range of the comparison values in the present research work suggested the MAE to be the more representative of the two indexes [87].

Annual RS-IWA Evaluation for Exploitation Systems in the Duero and Júcar RBDs
The daily NIR values obtained by the RS-SWB at a pixel-based spatial scale were then annually aggregated by their summation to compute the spatial statistics (described in Section 2.5) against the different ES irrigation water management spatial scales for the whole RBD in Júcar (9) and Duero (13) (Figure 1). Table 1 and Figure 2 show the irrigation water accounting comparison. While the Júcar RBD provides the net annual irrigation water demands in the RBMP, the Duero RBD provides the annual gross irrigation water demands; thus, to provide a coherent comparison between the published data and the model results, the global irrigation efficiency coefficient of 0.62, published by Duero RBMP, was used to convert the gross irrigation water demands into the net irrigation water demands. Furthermore, the irrigable surface areas estimated by the RBMP were compared with the irrigated surface areas estimated by the RS-SWB.

Annual RS-IWA Evaluation Regarding Agricultural Demand Units in Duero and Júcar RBDs
The daily NIR obtained by the RS-SWB at a pixel-based spatial scale were then annually aggregated by their summation to compute the spatial statistics (described in Section 2.5) against the different ADUs' irrigation water management spatial scales for the whole RBD in Júcar (98) and Duero (290) (Figure 1b,c). Table 2 and Figure 3a,b show the irrigation water accounting comparisons. The 5 year scale annual average RS-IWA comparison with the official values for the 2015-2021 period showed a strong regression coefficient r 2 for both the Duero and Júcar RBDs (0.95 and 0.99, respectively), as well as a strong Willmott refined index d r (0.88 and 0.93, respectively). In that sense, the RE MAE provided values of 20% and 15% for the Duero and Júcar RBD, respectively. Case by case, the majority of the comparison values for the ESs in the Duero RBD showed a slight underestimation. However, this was not the case for the Júcar RBD, except for the Vinalopó-Alacantí ES. The latter case showed a larger underestimation, which, under our criteria, was based on the following: (a) the comparison between the irrigable surface area (RBMP) and irrigated surface area (RS-SWB) presented a larger negative relative error (−56%) compared to the other Júcar ESs, revealing many surface areas not considered by remote sensing techniques in this ES, and (b) the Júcar RBMP estimated a large irrigable surface area for successive annual crops (double harvest groups) with consequently high irrigation water demand estimations by the Hydrological Planning Offices. However, in this study, the crop classification for this area did not reveal an important surface for double harvest crop classes but did reveal an important surface for spring and summer classes (i.e., there was less irrigation water demand than that for the double harvest group). Even so, this underestimation does not invalidate the good r 2 and d r agreement reached in this RS-IWA evaluation due to its low statistical weight compared to the other ESs.

Annual RS-IWA Evaluation Regarding Agricultural Demand Units in Duero and Júcar RBDs
The daily NIR obtained by the RS-SWB at a pixel-based spatial scale were then annually aggregated by their summation to compute the spatial statistics (described in Section 2.5) against the different ADUs' irrigation water management spatial scales for the whole RBD in Júcar (98) and Duero (290) (Figure 1b,c). Table 2 and Figure 3a,b show the irrigation water accounting comparisons. While the Júcar RBD provides the net annual irrigation water demand in the RBMP, the Duero RBD provides the annual gross irrigation water demand. Thus, one irrigation efficiency coefficient for each ADU in the Duero RBD (provided by the RBMP) was compared against the RS-IWA values. The 5 year scale annual average RS-IWA comparison with the official values for the 2015-2021 period again shows a strong regression coefficient r 2 of 0.96 for both the Duero and Júcar RBDs, as well as a strong Willmott refined index d r (0.90 and 0.88, respectively). The RE MAE shows values of 26% and 28% for the Duero and Júcar RBDs, respectively. At this spatial scale, a common underestimation was observed for both RBDs. Hence, we highlight several ADUs in the Júcar RBD comparison in Figure 3a. The first group (black squares in Figure 3a) is represented by four ADUs, all located in the Vinalopó-Alacantí ES (082073A, 082075A, 082076A, and 092001A), that demand 84% of the irrigation water resources in this ES and were underestimated by the RS-IWA according to previous observations (commented in Section 3.1). The second group (red circles in Figure 3a) contains six ADUs (082034C, 082034D, 082054B, 082054C, 082054D, and 082054E) in which rice is the majority crop cultivated. In these ADUs, the RS-IWA underestimated the irrigation water demands by about 50% for the RBMP-IWA, except for the ADU 082054B, for which a 25% lower irrigation water demand was estimated by the RS-IWA. The difference between these two subgroups relies on their rice surface irrigation areas because, for the latter ADU, the proportion of rice is just 20%, compared to the other five ADUs, where the rice proportion is between 50 and 100%.   Figure 3a) is represented by four ADUs, all located in the Vinalopó-Alacantí ES (082073A, 082075A, 082076A, and 092001A), that demand 84% of the irrigation water resources in this ES and were underestimated by the RS-IWA according to previous observations (commented in Section 3.1). The second group (red circles in Figure 3a) contains six ADUs (082034C, 082034D, 082054B, 082054C, 082054D, and 082054E) in which rice is the majority crop cultivated. In these ADUs, the RS-IWA underestimated the irrigation water demands by about 50% for the RBMP-IWA, except for the ADU 082054B, for which a 25% lower irrigation water demand was estimated by the RS-IWA. The difference between these two subgroups relies on For the annual RS-IWA comparison, the model collects annual variability as observed previously in the ES spatial scale. However, in this case, we highlight the minimum values-i.e., the years with less irrigation water demands determined by the RS-IWA for both RBDs in 2016 and 2018.
Besides the annual irrigation water demands estimated by the Duero RBMP, annual monitoring measurement reports are available along three consecutive hydrological year (October to September for 2015-2016, 2016-2017, and 2017-2018) for the water consumed at the EDU spatial scale. Choosing the EDUs with a maximum bias of 20% between the irrigable and irrigated surface areas to use comparable data, we ultimately selected 49 EDUs. Figure 3c shows the irrigation water accounting comparison between the recorded measurements and the annual RS-IWA* (recomputed at annual temporal scales considering the hydrological year). The statistical analysis reveals results similar to previous findings, with a regression coefficient of r 2 = 0.94, a Willmott refined index of d r = 0.81, and an RE MAE of 24% (MAE = 2 hm 3 ).

Monthly RS-IWA Evaluation of Agricultural Demand Units in the Júcar and Duero RBDs
The daily NIR obtained for the RS-SWB at a pixel-based spatial scale were then aggregated monthly by their summations to compute the spatial statistics (described in Section 2.5) for the different ADU irrigation water management spatial scales for the whole RBDs in Júcar (98) and Duero (290) (Figure 4). While the Júcar RBD provides the annual net irrigation water demand in the RBMP, the Duero RBD provides the annual gross irrigation water demand; thus, the irrigation efficiency coefficient for each ADU in the Duero RBD was compared against the RS-IWA values. The r 2 and d r monthly time series are also represented in the figure.
Both time series for the 5 year scale monthly average RS-IWA and the Monthly RS-IWA follow similar monthly temporal patterns in their monthly estimated RBMP-IWA. In that sense, the maximum values are observed in spring or summer, while the minimum values occur in autumn and winter. However, for both RS-IWA time series, underestimations are common for all the years, especially in summer, while for some spring and autumn months, higher values were estimated. These values reflect time series decoupling for autumn in both RBDs, which is clearer in the Duero RBD, where higher monthly RS-IWA values are constant over the 5 years; therefore, the 5 year scale monthly average RS-IWA reflects this overestimation. Conversely, in the Júcar RBD, slight overestimations occur in spring.

Monthly RS-IWA Evaluation of Agricultural Demand Units in the Júcar and Duero RBDs
The daily NIR obtained for the RS-SWB at a pixel-based spatial scale were then aggregated monthly by their summations to compute the spatial statistics (described in Section 2.5) for the different ADU irrigation water management spatial scales for the whole RBDs in Júcar (98) and Duero (290) (Figure 4). While the Júcar RBD provides the annual net irrigation water demand in the RBMP, the Duero RBD provides the annual gross irrigation water demand; thus, the irrigation efficiency coefficient for each ADU in the Duero RBD was compared against the RS-IWA values. The r 2 and dr monthly time series are also represented in the figure.  For the monthly evolution of the statistical index, we observed strong results for the r 2 and d r (about 0.9 and 0.8, respectively) in months where more irrigation water resources are consumed-i.e., during the central months of irrigation campaigns in both RBDs, mainly April to September. However, these values strongly decrease in the winter months, when the minimum water resources are consumed for the monthly estimated RBMP-IWA, yielding r 2 values close to 0.0, especially in the Duero RBD, which also has d r monthly values close to −1.0. In the Júcar RBD, during the winter months, the statistical values are slightly higher than those in the Duero RBD. These results are related to the minimum (or absence, as in the Duero RBD) of irrigation water demand estimations from the RBMP. Hence, the comparison shows a remarkably similar numerical data range, albeit with a very low dispersion. Thus, while statistically speaking, r 2 values close to zero show no correlation, mathematically speaking (and according to the r 2 formulation), an output of values close to zero makes sense.
In addition to Figure 4, which shows the global situation in the RBDs according to the monthly irrigation water accounting and considering all the ADUs, Figure 5 shows the same irrigation water accounting time series but for the specific ADUs demanding the highest irrigation water volumes and belonging to different ESs. We included three ADUs per RBD: one of them principally depends on groundwater resources (Mancha Oriental and Bombeo Medina del Campo), while the other two depend greatly on surface water resources. Moreover, the monthly P and ET o values from the closest agrometeorological stations using the SIAR network were plotted. Including these climatic variables helps to explain the monthly variability of the RS-IWA. Remote Sens. 2020, 12, x FOR PEER REVIEW 16 of 28

Annual RS-IWA Evaluation for the ES Spatial Scale for the 11 RBDs in Mainland Spain
The daily NIR obtained by the RS-SWB at a pixel-based spatial scale were then annually aggregated by their summation to compute the spatial statistics (described in Section 2.5) for the Again, the maximum and minimum tendency values agree for both RS-IWA estimations in the monthly estimated RBMP-IWA. However, on a smaller scale and in response to climatic conditions, the RS-SWB shows variations in the estimation of the irrigation water demands when comparing the same months in different years (for example, May). Seasonal decoupling that does not manifest in other years is also observable for specific years (particularly in spring and autumn). In that sense, decoupling is also observable in the spring and autumn months. The highest decoupling was found every year for both of the ADUs with groundwater resource dependency, Mancha Oriental (Júcar RBD) and Bombeo Medina del Campo (Duero RBD). While the former ADU always reaches its monthly RS-IWA maximum values in spring, like the monthly estimated RBMP-IWA, in the latter case, these values are distributed in the spring and autumn seasons. Moreover, the RS-IWA distributes its maximum values over the summer months, while the RBMP-IWA estimates its maximum values in June.
For the statistical indices, we found strong r 2 agreements for ADUs with a superficial origin of water resources (between 0.89 and 0.94), which are much lower for ADUs with groundwater resource origins (0.66 in Medina del Campo and 0.49 for Mancha Oriental). In the same way, we observed stronger d r agreements for ADUs with superficial water resources (between 0.86 and 0.91) and the weakest ones for ADUs with groundwater resources (0.72 and 0.63, respectively).

Annual RS-IWA Evaluation for the ES Spatial Scale for the 11 RBDs in Mainland Spain
The daily NIR obtained by the RS-SWB at a pixel-based spatial scale were then annually aggregated by their summation to compute the spatial statistics (described in Section 2.5) for the different ES irrigation water management spatial scales of the 11 RBDs under study (four in GDN, one in SEG, two in GyB, one in TOP, five in CMA, four in CIC, five in TAJ, 13 in DUE, eight in GDQ, 17 in EBR, and nine in JCR). Table 3 and Figure 6 show the irrigation water accounting comparison. In the eight RBDs (TAJ, GDN, SEG, CMA, GyB, DUE, GDQ, and JCR), the global RBD efficiency coefficient was available, so the annual net irrigation water demands were computed. However, for the other RBDs, we did not find this coefficient, so we used the average value of these eight RBDs (a range of values between 0.41 and 0.81, with an average of 0.68) to convert the annual gross irrigation water demand estimations into their net values. Units: Av, Min, Max, σ, RMSE, and MAE are in hm 3 /year; relative errors are in %; and dr and r 2 are dimensionless.    (Table 3), which is clearly influenced by the annual RE MAE of 62% returned for 2018 in comparison with the four preceding years, where the RE MAE ranged from 30 to 45%. For the annual variability returned by the annual RS-IWA, the year of 2015 has the highest irrigation water demands, even larger than the annual estimated RBMP-IWA, while 2018 has the lowest. In terms of irrigation and irrigated surface area comparisons (Figure 6b), the statistical indices for each year separately, as well as the 5 year average, show good agreement, with an r 2 in the range of 0.95 to 0.98 and a d r of 0.86 to 0.88. In that sense, the RE MAE ranges from 27 to 34%, so any year could be highlighted as the best one. Figure 6 shows five ESs highlighted with red circles showing high under-or overestimations of the 5 year scale annual average RS-IWA, in contrast with published data on the annual estimated RBMP-IWA. Hence, in the four ESs (the Segura ES in SEG RBD, Regulación General ES in GDQ RBD, Sierra Nevada ES in CMA RBD, and Bajo Ebro ES in EBR RBD), the RS-SWB strongly underestimates the annual RS-IWA, with one high overestimation (Oriental ES in TAJ RBD). These four ESs also underestimate the irrigated surface area. These findings support the idea that the irrigated surface area determines good agreement for irrigation water accounting. However, a deeper analysis reveals important ES characteristics that deserve analysis.
Underestimation of 620 and 460 hm 3 /year was detected in the two ESs selected (Bajo Ebro (EBR RBD) and Regulación General (GDQ RBD), respectively), where the rice crop surfaces are presented as about 30,000 and 40,000 ha. As mentioned previously for the six ADUs from the Júcar RBD with rice crop surfaces (Section 3.2), a new and highlighted underestimation was revealed for such crops. A different situation was revealed in the Sierra Nevada ES (CMA RBD), where horticultural crops under greenhouse conditions were the predominant crops (the El Ejido zone, Almería). Hence, with a surface of about 27,000 ha (60% of the irrigable surface area estimated by the RBMP), an underestimation of 230 hm 3 /year was obtained. However, the NIR of this crop group cannot be simulated by the RS-SWB, as the spectral response from the satellite for greenhouses is linked to the plastic cover, not to crop development. Moreover, the other three ESs that comprise the CMA RBD also feature greenhouses, although not as densely, so underestimation is not as great as that for the aforementioned ES. The last selected underestimation is linked with horticultural crops in successive harvesting during the year, as is the case for the Segura ES (SEG RBD). The RBMP estimated around 50% of the irrigable surface area dedicated to this goal, while the supervised classification did not reveal such a figure, instead indicating mainly spring and summer crops. Therefore, the underestimation of about 450 hm 3 was related to an incorrect classification and, consequently, poor crop RS-SWB parametrization, while the RBMP estimates irrigation water demands considering several successive crop harvests in the same year.
On the other hand, we observed an overestimation (about 280 hm 3 ) for the Oriental ES (GDN RBD). Although the irrigated surface area returned by the RS-SWB was slightly higher than the irrigable surface area estimated by the RBMP, we believe that the annual estimated RBMP-IWA is restricted because in the ES, the irrigated agriculture is groundwater-dependent, and the aquifer is classified as overexploited.

Discussion
The present research evaluated the RS-IWA's reliability over large and diverse irrigation surface areas. Using the FAO56 approach and remote sensing images, an RS-SWB was applied to estimate the NIR at a daily temporal scale and a pixel-based spatial scale. The spatially distributed thematic cartography generated was then aggregated into the same spatial and temporal irrigation water management scales that public managers use for water management expressed in the RBMPs. Therefore, the RS-IWA was evaluated over three spatial scales (EDU, ADU, and ES) and two temporal scales (annual and monthly) using the official irrigation water demands published in the RBMPs.
The spatial adaptability of the NIR thematic cartography to a desired spatial management scale depends on the original satellite images' spatial resolutions. In that sense, for 2014 and 2015, the present research used 30 m-scale satellite Landsat 8 data, while for 2016-2018, we used 10 m-scale data from the twin satellites of Sentinel-2 (A and B). In parallel, the temporal adaptability to the desired temporal management scale depends on the product's temporal frequency generation. Thus, the RS-SWB operates with a daily time step using dense satellite VI time series to allow adequate crop evolution monitoring. Both irrigation management working scales are appropriate for irrigation water monitoring with different plot shapes and sizes [88], including Sentinel-2 products for precision agriculture [89]. Moreover, the temporal and spatial scales produced agree with different and actual remote sensing operative applications, such as ALEXI/DisALEXI, which generates daily evapotranspiration products up to a 30 m spatial scale [14]; MINARET, which annually estimates daily unstressed crop evapotranspiration at a 30 m spatial scale along the Guadalquivir RBD for water management purposes; and the daily TOP-SIMS maps in California on basal crop evapotranspiration, among other products, at a 30 m spatial scale.
To date, ET c has been commonly used to evaluate models over irrigation surface areas [90][91][92][93][94], as has the Crop Water Requirements approach [17,24,95]. However, increasingly, studies use the irrigation water volume, both the estimation and value recorded at the field plot level, to evaluate models over irrigation surface areas [26,96,97]. In this vein, the present research used NIR as a variable to evaluate the RS-IWA based on different criteria. On the one hand, there is no distributed instrument network for in situ ET c measurements [72] that can supply the required ground truths for the different irrigation management spatial scales used by water authorities. On the other hand, we can use the same magnitudes and products that water authorities publish in their RBMPs after their estimation by the Hydrological Planning Offices based on the best resources and tools available to these water managers. Thus, published data (without being strictly ground truths, since many data are estimations) play a legal role, whose numerical values are used as a reference for irrigation water management during each hydrological planning cycle. Therefore, the presented approach provides irrigation water accounting assisted by remote sensing, and the variables used to validate them represent a form of novelty. Recent work shows that similar procedures have already used data from water authorities and remote sensing to assess irrigation water requirements over large areas and provide recommendations such as for several River Basins located in Jordan [98] or the Tarim River Basin in China [99].
The RS-IWA evaluation was satisfactory at different temporal and spatial scales. The agreement indexes r 2 and d r returned positive values within the range of 0.79-0.99, except for the winter season, when the agreement indexes were closer to zero when evaluated at the monthly time scale (nevertheless, these findings are not statistically representative, as already explained in Section 3.3). Furthermore, the RE MAE maintained ranges of 11-40% in 28 of the 31 conducted evaluations, which is in line with the similar range of about 15-40% for typical errors when ET c is computed through the VI [22]. We used these ranges like qualitative indicators because their variables are not the same. Finally, we also highlighted the EDU spatial scale for the Duero RBD because the annual RS-IWA* evaluation compared to the ground truth volumetric measurements returned good results (r 2 = 0.94, d r = 0.81, and RE MAE = 24%), which were in line with similar findings produced (good agreements and RMSE of about 15%) by the same model over three consecutive years at the plot and water management user scales in a different location [33].
The presented operational approach incorporates annual and monthly climatic dynamics, mainly, the annual/monthly RS-IWA product, which corresponds to the precipitation recorded at different spatial scales and is in line with other research's crop water requirement estimations [100,101]. Conversely, the 5 year scale annual/monthly average RS-IWA product was used to aggregate these climatic dynamics to incorporate successive study years. Thus, according to the seasonal and annual climatological reports (http://www.aemet.es/es/serviciosclimaticos/vigilancia_clima/resumenes) by the Spanish Meteorological Agency along with the present research work, two years were considered wet (2014 and 2016, with 680 and 682 mm, respectively, of annual precipitation), one year was very wet (2018, with 808 mm), and two more years were classified as dry years (2015 and 2017, with 500 and 474 mm, respectively). This annual variability was incorporated into the RS-SWB (Table 3), returning lower annual RS-IWA values for wet years and higher ones for dry years, which suggests that the annual RS-IWA is more accurate for the RBMP-IWA when a dry year occurs. In parallel, the behaviour of the monthly RS-IWA is similar, as expected based on previous findings. Therefore, the RS-SWB model incorporates monthly variability, as shown in Figure 5, from particular ADUs in the Duero and Júcar RBDs, where the lower irrigation water demand estimations respond to higher precipitation (i.e., June 2018) and vice versa, including the decoupling of the monthly RS-IWA time series for the RBMP-IWA. Such monthly temporal mismatching between remote sensing and official water managers' approaches have been observed in similar studies in the case of the Doukkala region in Morocco [100] or the Haihe Plain in China [102].
For the RS-IWA spatial scale aggregation following the deep evaluation conducted over the Júcar and Duero RBDs, the results did not show strong evidence for which spatial scale provides better results. The RE MAE at the EDU and ADU spatial scales varied from 24 to 43%, providing slightly better results in the ES spatial scale (with a range of 11-37%). However, when increasing the study area to 69 ESs over the 11 RBDs, good agreement indexes were found, but the RE MAE for the product of the 5 year scale annual average RS-IWA was about 41% (because the year 2018 returned a 61% RE MAE value, in contrast with the average value of 38%, considering the previous four study years, which is in agreement with the worst RS-IWA results being found for wet years). At this spatial scale, as noted in Section 3.4, the underestimations that trigger a higher RE MAE are linked with ESs containing a high proportion of rice crops, greenhouses, or successive horticultural crops each year.
Moreover, such improvements are already available. The spatial locations of rice crops are correctly defined by the Common Agricultural Policy GIS, so actual public water methodologies from the RBMP methodology can be used to estimate the irrigation water demands. However, considering the remote sensing benefits of monitoring crop developments in near real time, we could use the operational approach of González-Dugo et al. [103] (applied using MINARET) to obtain rice crop extensions followed by the application of other vegetation indices, such as the SAVI, to estimate unstressed crop evapotranspiration (using certain adapted crop parameters), among other approaches.
In parallel, the horticultural crops under successive harvesting located in the Segura ES from the SEG RBD can be correctly monitored by remote sensing using the high temporal resolution provided by the twin satellites of Sentinel-2 (A and B). In this sense, recent research applying Random Forest classification approaches returned good results for carrot, pumpkin, and onion classifications [104]. However, we must bear in mind the diverse crop types under this crop group and the different strategies of farming using different water management approaches [105]. Consequently, aggregating all these elements into one group class will require further research to maintain an easy and operational methodology. Finally, the classification of greenhouses could be achieved using many recent advances by applying Sentinel-2 [106,107]. However, the RS-IWA under greenhouses remains a limitation for the methodology presented here because optical satellites are not able to monitor such crops.

Conclusions
Presently, operational applications in irrigated areas based on remote sensing techniques are available. Thus, enlarged remote sensing applications over large and diverse irrigation surface areas would offer public and private water managers a complementary tool for obtaining information and satisfy the European requirements at the Water Framework Directive level. In this sense, the good agreement reached between the RS-IWA evaluation and the RBMP-IWA data estimated by the Hydrological Planning Offices and published in the RBMP satisfies the initial objective of presenting a simple, homogeneous, and reproducible operational tool for large and diverse irrigation surface areas. Moreover, the thematic cartography is characterized by having great adaptability to the spatial and temporal irrigation water management scales used by water managers, so NIR maps are applicable at any management level.
Consequently, several RS-IWA products were introduced. First, the 5 year scale annual/monthly average RS-IWA, which aggregates the temporal behaviour of agroclimatic and agricultural strategies over five study years, is suggested for use when medium-/long-term accounting is required, such as for the six-hydrological-year cycle. Secondly, the annual/monthly RS-IWA, which incorporates variable climatological conditions and allows one to monitor the actual values of NIR, can be used to improve annual or monthly irrigation water resource monitoring.