Assessing Spatio-Temporal Dynamics of Deep Percolation Using Crop Evapotranspiration Derived from Earth Observations through Google Earth Engine

Abstract

Excess irrigation may result in deep percolation and nitrate transport to groundwater. Furthermore, under Mediterranean climate conditions, heavy winter rains often result in high deep percolation, requiring the separate identification of the two sources of deep percolated water. An integrated methodology was developed to estimate the spatio-temporal dynamics of deep percolation, with the actual crop evapotranspiration (ET_{c act}) being derived from satellite images data and processed on the Google Earth Engine (GEE) platform. GEE allowed to extract time series of vegetation indices derived from Sentinel-2 enabling to define the actual crop coefficient (K_{c act}) curves based on the observed lengths of crop growth stages. The crop growth stage lengths were then used to feed the soil water balance model ISAREG, and the standard K_c values were derived from the literature; thus, allowing the estimation of irrigation water requirements and deep drainage for independent Homogeneous Units of Analysis (HUA) at the Irrigation Scheme. The HUA are defined according to crop, soil type, and irrigation system. The ISAREG model was previously validated for diverse crops at plot level showing a good accuracy using soil water measurements and farmers’ irrigation calendars. Results show that during the crop season, irrigation caused 11 ± 3% of the total deep percolation. When the hotspots associated with the irrigation events corresponded to soils with low suitability for irrigation, the cultivated crop had no influence. However, maize and spring vegetables stood out when the hotspots corresponded to soils with high suitability for irrigation. On average, during the off-season period, deep percolation averaged 54 ± 6% of the annual precipitation. The spatial aggregation into the Irrigation Scheme scale provided a method for earth-observation-based accounting of the irrigation water requirements, with interest for the water user’s association manager, and at the same time for the detection of water losses by deep percolation and of hotspots within the irrigation scheme.

1. Introduction

Given the climatic conditions of the Mediterranean regions, groundwater reserves play an important role in interannual regularization, minimizing the impact of climate variability, particularly droughts. However, the conservation of groundwater, in terms of quantity and quality, is a major environmental challenge [1]. The conservation of quantity requires the maintenance of its recharge, while the conservation of quality requires the existence of low concentrations of nitrates and other elements in the flow of water that reaches the aquifer. The challenge increases if we consider the additional effect of climate change, which could lead to a decrease in aquifer recharge, increased crop irrigation needs, and rising sea levels [2,3]. Due to the characteristics of the Mediterranean climate, where the hot season is also the driest, irrigation is mandatory for achieving high productivity and competitiveness in the agricultural sector. However, intensive land use in terms of irrigation and fertilizer consumption in areas overlying aquifers can lead to risks of contamination, as water, when in excess, or poorly distributed, can transport nitrates into the deepest layers [4]. This situation is particularly important in the Nitrate Vulnerable Areas (ZVN), where nitrate concentration in groundwater has reached/exceeded the limit imposed by the European community (50 mg NO₃⁻ L⁻¹). It is the deep percolation of water, resulting from precipitation and irrigation, that serves as a transport of nitrates to groundwater. Thus, the mismanagement of water and fertilizers plays an important role in modifying the natural recharge of aquifers and in the groundwater contamination [5].

Some studies present estimates of aquifer recharge as a percentage of the annual rainfall [6]. However, the deep percolation and recharge associated with the practice of irrigation has been subject of reduced attention. In irrigated areas, in the absence of precipitation, recharging during the summer months can be significant depending on the water balance in the crop rootzone, creating hotspots for potential aquifer contamination with nitrates from fertilizers, organic matter mineralization, etc.

The water balance in the crop rootzone, with the objective of calculating deep percolation, requires the determination of the crop water requirements which in turn relies on the correct estimation of their biophysical parameters, such as crop coefficients and the durations of phenological phases. Although there is information in the literature, as in the recent review by Pereira et al. [7,8], these authors state that these parameters may vary with several factors so regional characterization is advised.

Satellite remote sensing (RS) data have been widely used for phenological monitoring of vegetation in large geographic domains, using various observation methods [9,10,11,12]. Emerging techniques are based on time series analysis of remote detection multispectral imaging which can be used to estimate the crops’ phenological phases (e.g., beginning of the cycle, peak vegetative development, and the duration of phenological phases) [13,14,15]. As recently discussed by Mahlayeye et al. [16], there is a need to further study the cropping patterns at larger scales such as irrigation schemes and watersheds.

Remote sensing temporal data with high spatial and temporal resolution can be acquired, for example, from the Sentinel-2 satellite, which is an Earth Observation (EO) mission of the European Union, the Copernicus program, that systematically acquires high-resolution optical images. Together, Sentinel-2 satellites have high temporal resolution (5 days) and spatial resolutions of 10, 20, and 60 m. These images are made available free of charge by the European Space Agency (ESA) and can be downloaded on various platforms, such as the Google Earth Engine (GEE), which makes these images available in its library as a collection. It has two Sentinel-2 collections, the “COPERNICUS/S2” which corresponds to the reflectance at the top of the atmosphere and the “COPERNICUS/S2_SR” for the surface reflectance [17,18].

Vegetation indices (VI) correspond to an image treatment technique that is applied to extract information such as crop planting density and vegetation vigor (e.g., [19]). Based on these, the values of the actual crop coefficients (K_{c act}), planting dates, and the duration of the phenological phases can be estimated. They can also be used to identify the soil cropping patterns in a given period [20,21,22]. The Normalized Difference Vegetation Index (NDVI) is a widely used index to characterize crop phenology and it is calculated using the red and near infrared bands present in most sensors [23,24]. Several authors have proposed correlations between NDVI and K_{c act} based upon field data, establishing a linear relationship between the two variables [10,20,25,26,27], although this relation is affected by other factors including crop architecture. The approach proposed by Calera-Belmonte et al. [20] has been widely used both at the plot scale and at the regional level for obtaining the actual K_c, the length of the crop phenological phases, and the identification of crop types, with the aim of determining the irrigation water requirements [10,13]. The K_{c act} values are then used to calculate the actual crop evapotranspiration and subsequently perform the soil water balance [26,28,29,30].

The availability of free satellite images with high spatial and temporal resolution makes the use of this technology very interesting for obtaining crop data to improve the estimation of actual ET_c. Furthermore, the Google Earth Engine (GEE) [31] which is a “cloud-based” platform for planetary-scale geospatial analysis that combines computing capabilities to the storage and processing of satellite imagery, enables the access to publicly available RS data, such as Landsat 4–8 or Sentinel 1–2, and its data processing and visualization. Thus, GEE allows performance of this type of study without the acquisition of high-performance computers, enabling a wider adoption of EO data on the estimation of ET_{c act}. This service contains several libraries with different types of data and allows the processing of this data on their own servers. Datasets can be acquired and manipulated with a JavaScript code editor platform [31]. GEE allows for rapid data processing, namely the automation of the extraction of temporal series of vegetation indices (e.g., NDVI), without the need for a large data storage capacity, or high-performance computers [31,32,33]. These characteristics of the platform made it possible to overcome important limitations to a greater adoption of high-resolution satellite images, as is the case of Sentinel-2, for crop data retrieval [31]. Vegetation indices such as NDVI can be extracted from satellite images on the platform, and they can be viewed in time series charts or transferred to external processing. GEE was successfully used to estimate crop evapotranspiration [34,35,36], crop productivity [37,38], to cartograph crop areas [39,40,41], for supporting irrigation scheduling [42,43], and water management [44].

Aiming at performing the spatial analysis of the deep percolation and detecting hotspots at an irrigation scheme level, it is necessary to integrate soil water balance modeling with spatially distributed information in a GIS environment. To perform this integration the following steps are required: (a) development of the GIS database, climatic database, soil database, and other databases (administrative limits, irrigation, topographic and land use information); (b) integration of the database and the customization of the GIS; (c) integration of the irrigation simulation model. Soil water balance modeling can be completed in a GIS environment, using a vector approach through the definition of homogeneous analysis units (soil × crop × irrigation system), or using data grids (raster) for climate, soil, and crops [45].

The question to be answered by the present study is: are there hotspot areas in the Caia Irrigation Scheme with high deep percolation due to irrigation practices? To answer this question, the following objectives were established: (1) estimate crop evapotranspiration using remote sensing images to improve the definition of the crop stages; (2) define homogeneous analysis units (soil × crop × irrigation system) in a GIS environment, for which the soil water balance is modeled; (3) estimate the deep percolation associated with the irrigation and with the precipitation along the crop season, and also off season, after crop harvesting; and (4) perform the spatial analysis of deep percolation and identify its potential causes.

2. Materials and Methods

2.1. Caia Irrigation Scheme Case Study

2.1.1. Localization and General Characteristics

The Caia Irrigation Scheme (CIS) is located in the Alentejo region, southern Portugal (Figure 1). The total irrigated area of the CIS, which is around 7000 ha, is included in the Elvas Vulnerable Zone to Nitrates. The CIS serves 887 farmers with water from the Caia Reservoir, located in the Caia River, which also provides water for domestic consumption.

Water is distributed to farms in open channel, with upstream control, and a fixed rotation scheme where each farmer receives a maximum amount of water per hectare.

2.1.2. Climate and Soil Characterization

According to the Köppen classification [46], the climate in the region of study is a Csa—temperate Mediterranean with hot and dry summers. Weather data used in the current study was collected from a nearby weather station (38°54′56″ N, 7°05′56″ W), located under reference site conditions. The long term (2003–2021) average annual precipitation is 519 mm, mostly occurring in winter, while the average annual reference crop evapotranspiration (FAO-PM, [47]), which represents the climatic demand, is 1404 mm. Figure 2 shows the monthly average values of the reference crop evapotranspiration (ET_o) and precipitation for the four studied years. It clearly shows a water shortage for the crops during the spring–summer seasons. When comparing the conditions of the studied years with the long-term data, it is shown that 2018 was a humid year with 835 mm of precipitation while 2019 was the driest one with 317 mm. The annual ET_o ranged from 1290 to 1375 mm during the studied years.

According to the World Reference Base of Soil Resource [48], the predominant soils in the irrigated perimeter are Fluvissols (44.9%), Luvissols (29.6%), and Calcissols (18.7%).

2.1.3. Main Crops and Irrigation Systems

In 2002, the main crops in CIS were maize (Zea mays L.) occupying 49% of the area, wheat (Triticum aestivum L.) with 17%, sunflower (Helianthus annuus L.) with 7%, tomato (Solanum lycopersicum L.) with 6%, and olive orchards (Olea europaea L.) (super-intensive and hedgerow) with 4% of the total area [49]. However, in the last few years (Figure 3) there has been a significant conversion from areas with annual crops to permanent crops such as olive orchards and nut-tree orchards, including walnut and almond.

These areas have been stabilizing in the last years (Figure 3). According to the Caia water user’s association, in 2020 the irrigated area was 8559.4 ha, because it includes areas that are outside the irrigation scheme but that receive water from the CIS. The main irrigation system is drip with 82% of the irrigated area, followed by sprinkler irrigation with 17.1%, mainly center pivots, while surface irrigation has a very small share of the irrigated area.

2.2. Data Collection and Georeference Database Building

A modeling tool (ISAREG model [50]) was used to support the estimation of the crops water consumption (ET_c) and of the deep percolation associated with both irrigation and precipitation along the studied years (off-season was included in the study). To feed the modeling tool, a georeferenced database with diverse information was built. Due to the diversity of data sources, they were edited and converted to the same spatial and temporal resolution.

Data used for the estimation of the deep percolation are summarized in

Table 1. Data sources used in the present study.

Data Set	Observations	Source
Soils	Soil water holding capacity and textural characteristics; land use capacity Map 1:25,000	Portuguese Soil Map (CSP) and Land Use Capacity (DGADR—Ministry of Agriculture)
Weather	Daily weather data (2002–2020): maximum and minimum temperatures (°C), maximum and minimum relative humidity (%), solar radiation (kJ·m−2·dia−1), wind speed (m·s−1), and precipitation (mm)	Meteorological station of Elvas (38°54′56″ N, 7°05′56″ W, 202 m a.s.l) (COTR)
Topography	Military Map of Portugal (1:25,000)	Instituto Geográfico do Exército
Hydrography and altimetry	DTM (Digital Terrain Model) for the slope Contour lines, elevation points, and water flow lines	Instituto Geográfico do Exército
Administrative limits	Caia Irrigation Scheme limits	Official Administrative Map of Portugal and Irrigation Scheme Maps of Portugal—Ministry of Environment
Nitrate Vulnerable Zone limits	Nitrates Vulnerable Zones limits Map (1:25,000)	Ministry of Agriculture
Crops	Cropping patterns, crop phenological stages	Sentinel-2 images (Level-2A), Google Earth Satellite Images (QMS)
	Crop coefficients	[7,8,51]
Land cover classes	Land use map (COS) Land cover classes, Corine Land Cover (CLC) 2018 (spatial resolution 20 m) Crop plots map	Ministry of Environment Copernicus Program (EEA, JRC)Caia Water Users Association
Farm/plot data	Farm/plot identification (WMS/WFS format)	Ministry of Agriculture
Groundwater	Water table depth	Water Resources National Information System—Ministry of Environment
Irrigation	Crop systems, irrigation calendars, irrigation systems, and soil moisture data	Caia Water Users Association
Irrigation systems	Irrigation systems efficiency	[52]

and the data sources used to collect the various types of information are also included.

2.2.1. Crop Characteristics Database

The data relative to crop characteristics included the identification of the crop type, spatial distribution of the diverse cropping systems, crop growth stages dates, and crop coefficients for the diverse crop stages. These data were collected from diverse sources as detailed in Table 1.

As previously stated, due to its capabilities the Google Earth Engine (GEE) platform was used in the current study.

• Satellite imagery processing and crop data collation

Aiming at decreasing the uncertainties in the characterization of the cropping systems and crop growing seasons, Sentinel-2 satellite imagery with high spatial (10 × 10 m) and temporal (5 days) resolution was used. All Sentinel-2 images available on GEE platform for the CIS region between October 2016 and October 2020 were used, with a cloud cover with less than 35%, which corresponds to a compromise between the number of images available to produce the VI time series and its quality. The Sentinel-2 information was then processed in the platform Google Earth Engine using the Java script code editor.

The flowchart presented in Figure 4 shows the global process of selection, validation, and pre-processing of information in order to build the database for crops.

2.

Crop identification and spatialization

Firstly, a layer containing the spatial distribution of the different types of crops was created in QGIS; it included the vector information of the cultivated fields. This layer was based on field information provided by the Caia water user’s association (WUA). In the following and based upon the land use map (COS) (Table 1) the polygons of permanent crops (olive orchards, orchards, and vineyards) were extracted, based on the assumption that these crops vary little from year to year within the studied years (Figure 3). For producing the NDVI [53] for each crop and each field the series of the spectral bands of the Sentinel-2 sensors B4 (red) and B8 (near infrared) and QA60 (cloud mask) were used from the image collection “COPERNICUS/S2_SR”, which corresponds to Level-2 data (orthorectified Bottom-Of-Atmosphere (BOA) reflectance) [18]. The GEE allowed the extraction of the average NDVI per pixel for the period March–May which allowed for the identification of the winter crops (Figure 5a) while the average NDVI for the period May–August allowed identifying the spring–summer crops (Figure 5b). During summer, the NDVI map allowed the precise identification of the irrigated areas. Figure 5c shows an example of the identification of the actual irrigated areas using the Sentinel-2 images. Images were previously processed to discard the pixels that correspond to small roads separating orchard blocks, roads, and edges as the soil conditions in these pixels are not the same as agricultural zones without vegetation; thus, avoiding its influence on the vegetation indices calculation. In addition, the crop NDVI temporal profile analysis allowed identifying the plot fields that did not match the typical crop NDVI pattern, which were removed from its respective crop area [14].

The cropping pattern in CIS includes more than 30 types of crops, and therefore crops were grouped according to their relative importance in terms of cropping area, crop cycle, and irrigation water requirements.

3.

Locally adjusted crop growing stages

Subsequently, the processed satellite images were used to obtain the observed planting dates and the diverse crop growth stage lengths, since the poor characterization of the crop cycle is one of the main sources of uncertainty in the estimation of crop water requirements (e.g., [54]).

The polygons corresponding to the fields with the different crop groups were imported into the GEE platform and the average values of NDVI for each polygon were extracted from the image collections “COPERNICUS/S2” and “COPERNICUS/S2_SR”.

In the following, the NDVI temporal profile for each crop was obtained from the spectral signature of the fields identified with a certain crop; the actual K_c curves (affected by diverse stresses, e.g., water and nutrients) were derived from those profiles using empirical equations [20,25,27]. Atypical NDVI profiles were identified and discarded.

The NDVI temporal profile for each crop allowed estimating the average sowing/planting dates and the lengths of the growing stages. Thus, allowing defining the crop growing stages required for constructing the K_c segmented curve [47] for each crop. These were assumed to be representative for the entire CIS.

4.

Crop coefficients along the growing seasons

The standard potential K_c values for the diverse crop growing stages (initial, mid-season, and end-season) were collected from the tabulated values recently updated by Pereira et al. [7,8] and Rallo et al. [51]. For each year, the K_c values for the initial stage were those provided by Allen et al. [47] but adjusted for soil wetted fraction (irrigation or precipitation). While the crop coefficients for the mid-season and end-season values were adjusted to the prevailing local climate conditions (RH_min and u₂) in each studied year using an algorithm based upon the methodologies described in Allen et al. [47].

The crop data used in the computations is presented in

Table A1. Crop stage duration and seeding/planting dates obtained from the literature.

Crop	Crop Stage Duration (Days)				Planting Date	Literature
	Initial	Development	Mid	Late
Garlic	70	30	47	27	01/Dec	[6,81]
Broccoli	35	45	40	15	01/Nov.	[47]
Chickpea	20	30	35	15	15/Mar
Winter cereals	30	140	40	20	15/Oct
Olive groves	30	90	60	90	01/Mar
Orchards	30	50	130	30	01/Apr
Vineyards	20	50	75	60	01/Mar
Almonds	32	65	140	41	23/Jan	[82]

and

Table A2. Actual crop coefficients.

	Kc act ini				Kc act mid				Kc act end
	2017	2018	2019	2020	2017	2018	2019	2020	2017	2018	2019	2020
Garlic	0.70	-	-	-	1.00	-	-	-	0.65	-		-
Tomato	0.50	-	0.50	0.50	1.10	-	1.04	1.05	0.90	-	0.84	0.85
Melon	0.40	-	-	-	0.99	-	-	-	0.79	-		-
Sweet peppers	0.50	0.5	0.5	0.5	1.09	1.03	1.05	1.05	1.04	0.98	1.00	1.00
Broccoli	0.40	0.4	0.40	0.40	1.04	0.94	0.96	1.00	1.04	0.94		1.00
Sunflower	0.35	0.35	0.35	0.35	1.19	1.13	1.15	1.15	0.34	0.3	0.3	0.30
Rapeseed	0.35	0.35	0.35	0.35	1.03	1.02	1.05	1.00	0.28	0.35	0.3	0.25
Chickpeas	-	0.4	0.40	0.40	-	0.98	0.99	1.00	-	0.35	0.34	0.35
Rice	1.05	1.05	1.05	1.05	1.20	1.2	1.20	1.20	1.09	1.05	1.05	1.05
Maize grain	0.30	0.3	0.30	0.30	1.20	1.2	1.19	1.20	0.36	0.3	0,39	0.30
Maize silage and sorghum	0.30	0.3	0.30	0.30	1.20	1.06	1.14	1.15	1.01	0.86	0.94	0.95
Winter cereals	0.30	0.3	0.30	0.30	1.09	1.08	1.11	1.06	0.34	0.33	0.36	0.31
Fodder crop multiple cuts	0.40	0.4	0.40	0.40	0.99	0.89	0.91	0.95	0.94	0.84	0.86	0.90
Pastures	0.75	0.75	0.75	0.75	0.75	0.75	0.75	0.75	0.75	0.75	0.75	0.75
Almond	0.40	0.4	0.40	0.40	0.92	0.82	0.84	0.85	0.67	0.57	0.59	0.60
Orchards	0.45	0.45	0.45	0.45	1.06	0.97	0.99	1.00	0.76	0.67	0.69	0.70
Vineyards	0.30	0.3	0.30	0.30	0.75	0.67	0.69	0.70	0.60	0.52	0.54	0.55

in Appendix A.

2.2.2. Soil Characteristics Database

Due to the large spatial variability of CIS soils, they were grouped into classes of land suitability for irrigation. These classes result from the integration of three characteristics: soil texture, slope, and permeability. Class I includes soils with excellent suitability for irrigation, with sandy loam or clay texture, good structure and permeability, and a maximum slope of 2%. Soils in class II present good suitability for irrigation, are generally deep soils, and moderately drained. Soils within class III present reasonable suitability for irrigation and are shallow soils with heavy texture, with slopes varying between 4% and 6%. The soil textures were extracted form Portuguese Soil Map (CSP). Slopes were obtained from the digital terrain model (MDT), built using information on contour lines, elevation points, and the hydrographic network, with the support of QGIS. The permeability classes as defined by ARS/USDA [55] were used.

2.2.3. Homogenous Unit of Analysis

Using the information from the previously referred databases, diverse homogenous units of analysis (HUA) were defined at the CIS and the water balance was independently modeled for these units. Each HUA corresponds to an area where the following three features influencing the water balance terms do not spatially change: land suitability for irrigation, crop, and irrigation system. The HUA were obtained in QGis by intersecting the layers containing the three features above.

2.3. Soil Water Balance Calculation

2.3.1. Modeling Strategies and Field Validation

The soil water balance (SWB) was modeled for each HUA using the ISAREG software. The ISAREG model calculates crop evapotranspiration according to the FAO methodology and the SWB follows a reservoir approach. It was used because it has very good accuracy when adequately calibrated and validated against observations of the soil water (e.g., [56]) and/or of crop ET and because it has been applied to diverse climatic and soil conditions. Examples of application include those to field (e.g., [57,58,59]), vegetable crops [60,61], pastures and grasslands [62,63], and orchards [64,65]. Furthermore, the model has been used for the assessment of climate change on irrigation water requirements (e.g., [66,67]) and for irrigation management at irrigation scheme [68] and basin level [69].

The flowchart presented in Figure 6 shows how the water balance is calculated by ISAREG. The SWB was separated in two components, one for the irrigation season, to account for the percolation created by irrigation, and the other for the off-season period, to account for the percolation originated by autumn–winter precipitation. The general water balance equation (Equation (1)) is applied to the crop rootzone:

DP = P + I − ET_{c act} − ΔS

(1)

where DP is deep percolation, P is precipitation, I is irrigation (0 for the off-season period), ET_{c act} is actual crop evapotranspiration, and ΔS is the storage variation. Runoff was considered = 0 due to the large scale of application. All terms are in mm.

The information relative to the crop characteristics, derived from EO, along with the farmers actual irrigation calendars were introduced in the modeling tool. In a few cases, the information relative to the irrigation schedules was not available; thus, to overcome the issue the WUA provided for information regarding the total water consumption per field. This information enabled to estimate the irrigation schedules per crop based upon the assumption that the farmers behave similarly in terms of irrigation depths and timing.

Several fields with the most representative crops within the irrigation district (

Table 2. Farmers’ field locations and characteristics used to feed and to validate the water balance model simulations, Caia Irrigation Scheme.

Field	Crop	Area (ha)	Irrigation System	Soil Texture	Type of Field Data	Years with Data	Data Provider
A	Maize FAO600	58.8	Pivot	Loamy sand	Irrigation amounts and frequency Soil moisture	2018, 2019, 2020	APAP
B	Maize FAO200	37.4	Pivot	Loam	Irrigation amounts and frequency	2017	ABCaia
C	Olive grove_1	26.8	Drip	Clay loam	Irrigation amounts and frequency Soil moisture	2018, 2019, 2020	APAP
D	Olive grove_2	396.3	Drip	Clay loam	Irrigation amounts and frequency	2017	ABCaia
E	Processing tomato_1	20.2	Drip	Silty loam	Irrigation amounts and frequency Soil moisture	2019, 2020	APAP
F	Processing tomato_2	28.7	Drip	Loam	Irrigation amounts and frequency	2017	ABCaia

Note: APAP—Portuguese association of precision agriculture; ABCaia—Caia water users association.

) were used for validating the soil water balance simulations by comparing the soil water storage for the entire rootzone (ASW) predicted by ISAREG with soil water storage obtained from soil moisture data obtained from continuous monitoring with capacitive probes provided by the Portuguese Association of Precision Agriculture (APAP).

A set of goodness-of-fit indicators was used to assess the performance of the ISAREG model in estimating the soil water storage (ASW, mm). These indicators are: (i) the regression coefficient (b₀) describing a linear regression forced through the origin [70] between measured (ASW_OBS) and model-simulated ASW (ASW_SIM) values; (ii) the percent bias (PBIAS, %), which measures the average tendency of the simulated data to be larger or smaller than their corresponding observations (Gupta et al., 1999); (iii) the coefficient of determination (R²) of the ordinary least squares regression; (iv) the root mean square error (RMSE, mm); (v) the normalized RMSE (NRMSE, %), defined as the ratio between RMSE and the mean of observations $\overline{\mathrm{O}}$; and (vi) the modeling efficiency (EF, non-dimensional), which is an indicator proposed by Nash and Sutcliffe [71] used to assess the magnitude of the mean square error (MSE = RMSE²) relative to the observed data variance [72].

2.3.2. Irrigation Water Requirements

Irrigation water requirements (IWR) were calculated by dividing the net irrigation water requirements calculated with ISAREG, by the average irrigation system efficiency (drip = 0.9; center pivot = 0.8; traveling gun = 0.75; paddy rice basins = 0.7). In the case of permanent crops irrigated by drip systems, the value was multiplied by the reduction factor 0.8 to account for the reduction in the irrigated area at the soil surface [73].

3. Results

3.1. NDVI Temporal Profiles and Crop Coefficient Curves

3.1.1. Permanent Crops

The NDVI profiles for the permanent crops, e.g., vineyards, orchards, and olive groves are shown in Figure 7. The NDVI profiles allowed distinguishing plots with young orchards from those in full production. For the orchards, the end of dormancy was identified to be by March (DOY 90), while leaf senescence occurred by October (DOY 290). The NDVI increase during spring may be associated with the interference of wet soil upon the reflectance [74]. As shown in Figure 7, the presence of clouds along the crops’ seasons affected the precise identification of the various crop growth stages. For the case of the olive groves, the NDVI profiles allowed to distinguish between hedgerow and intensive groves. For the hedgerow, the NDVI is high, generally ranging from 0.3 to 0.6, and is minimally influenced by the presence of weeds in the inter-row spacing, while for the intensive groves, due to its management and higher space between trees the weeds develop particularly during spring and autumn; thus, presenting higher NDVI values. Overall, for all permanent crops, the high NDVI values at the beginning of the year may be associated with the presence of weeds that grow in the inter-row spacing due to winter and spring precipitation events.

3.1.2. Winter Crops

As per the permanent crops, the lack of images for some dates due to the presence of clouds prevented the adequate characterization of the crop growth stages for the winter crops (garlic, cereals, pastures) (Figure 8); thus, data was complemented with that available in the literature (Table A1 in Appendix A).

3.1.3. Annual Spring Crops

For the spring crops, the results of the NDVI temporal profiles allowed for an accurate estimation of the planting/sowing dates as well as for the duration of the crop growth stages; thus, allowing improvement of the definition of the crop coefficient curves, and therefore increasing the accuracy of the water requirement computations. Figure 9 presents examples of the NDVI time series obtained for several spring crops along the 2016/2017 season. The first attempt for the definition of the FAO K_c curves was to adjust it to the mean values of the K_c NDVI after applying the Calera Belmonte et al. [20] relationship to the NDVI profiles; however, this procedure created, in some cases, a flattening of the K_c curve. The use of the average K_c curves showed to be not appropriate for the definition of the crop growth stages because farmers used different varieties and sowing/planting dates. Therefore, a different approach was used in the current study, consisting of adjusting the K_c FAO curves, particularly the duration of the crop phases, to the K_c-NDVI clouds of points. Figure 9 shows the results of the K_c curves obtained for various crops while Table A2, in Appendix A, presents the values for the initial, mid- and end-season actual K_c for all crops in all the considered years.

Abnormal behavior of NDVI occurred mainly during the initial stage resulting in a less accurate determination of the duration of this stage. This behavior occurs mainly for the lower NDVI values, due to instability associated with soil background influence during the initial stage. Similar findings were reported by Cao et al. [74] and in the recent review by Pôças et al. [10]. Other sources of uncertainty in the definition of the K_c-NDVI curves is the use of empirical formulas calibrated and validated for other regions [10,21,27]. Similar findings were reported by Vilar et al. [75] using Landsat-8 images and by Rolim et al. [14] using images from SPOT-5 TAKE-5.

Table 3. Example of the average lengths of the spring crop growth stages and of the seeding/planting dates estimated from NDVI temporal profiles for 2017.

Crop	Lengths of Crop Growth Stages (Days)				Seeding/Planting Date	Obs.
	Initial	Development	Mid	Late
Paddy rice	30	27	50	13	25/Apr.
Tomato	30	37	50	13	25/Apr.
Melon	35	27	33	27	30/Apr.
Bell pepper	30	60	30	13	02/May
Rapeseed	30	30	40	25	02/Apr.
Sunflower	25	32	40	30	20/May
Grain maize	20	27	53	37	25/Apr.
Silage maize	22	30	35	12	13/May
Fodder crop multi-cuts	30	117	20	---	16/Nov.	1st cut
	6	31	20	---		2nd cut
	3	16	13	---		3rd cut

presents the length of the crop growth stages and the seeding/planting dates for the spring crops, estimated from the NDVI temporal profiles. Results show the adequacy of using the NDVI temporal profiles for improving the characterization of the crop cycles.

For tomato crop, the use of plastic mulch was assumed. In the case of pepper and melon crops, plastic mulch was not considered since most farmers in the CIS do not use it. For the winter vegetables, the K_c values that were ingested in the modeling tool correspond to the most representative crop, i.e., the broccoli because they are the main crop in terms of cropped area. The same reasoning was applied to pepper, which were considered to represent the spring–summer vegetables, while almonds were selected for representing the group of dry fruits. For the pastures and the autumn–winter forages, the winter cereal data were used because they presented similar crop cycles. For the irrigated pastures, an average K_c value of 0.75 was considered for the entire irrigation season.

The cropping patterns in the CIS, for the four studied seasons, after EO data validation against the data provided by the WUA are presented in Figure 10. The maps produced in the GEE platform with the average NDVI per pixel, for each growing period, allowed to identify the areas with winter and spring crops. During summer, the NDVI map allowed the precise identification of the irrigated areas. As previously stated, during the entire study period, olive groves were the predominant crop in CIS spreading across the irrigated perimeter. Winter cereals also occupy an important area, but they are preferentially located on the soil classes of low suitability for irrigation. Maize and fodder crops are the most representative annual spring–summer crops.

3.2. Irrigation Systems

The spatial distribution of the irrigation systems in the CIS for the studied years of 2017–2020 is presented in Figure 11. This information was provided for 2020 by the water users association, while, for 2017–2019, it was estimated from the high-resolution satellite images. As depicted in Figure 11, drip irrigation systems were predominant, followed by the sprinkler irrigation system with center pivots while travelling guns and surface irrigation have little significance.

3.3. Land Suitability for Irrigation

As previously explained, the classes of land suitability for irrigation, which is one of the layers needed for the definition of the HUA, were obtained by overlaying information relative to soil use capacity, land slope, and soil permeability. Figure 12 shows the spatial distribution of the land suitability classes for irrigation within the CIP. Results show that most of the CIS area consists of class III soils, presenting moderate limitations for irrigation, while classes II and I, which present lower limitations for irrigation, are mainly located along the banks of the Caia River. Results also show that the class with the lowest suitability for irrigation (exc) occupies 15% of the total area.

Table 4. Average soil properties for each class of land suitability for irrigation.

Land Suitability for Irrigation	θFC (g·g−1)	θWP (g·g−1)	Z (cm)	Bd (g·cm−3)
I	34.72	17.02	120	1.22
II	30.84	13.46	80	1.23
III	29.12	12.52	65	1.25
exc	30.83	14.18	50	1.25

Note: θ_FC is the soil moisture content at field capacity, θ_WP is the soil moisture content at the permanent wilting point, Z is the soil depth, and B_d is the bulk density.

presents the average soil properties used to characterize the different land suitability classes that were inputted in the modeling tool. These properties were determined considering the representativeness of each soil type within each class.

3.4. Homogeneous Units of Analysis

The produced homogeneous units of analysis (HUA), obtained by intersecting in QGIS the layers containing information of the spatial distribution of crop types, the land suitability for irrigation classes and the irrigation systems, are depicted in Figure 13. The total numbers of HUA were 104, 88, 112, and 109, respectively, for 2016/2017, 2017/2018, 2018/2019, and 2019/2020 (Figure 13). In a certain area, the changes in the HUA from one year to the other are mainly due to changes in the cultivated crops. An example of the detailed description of the HUA can be found in

Table A3. Description of the homogeneous units of analysis for the 2019/2020 season.

ID	Soil Type	Irrigation System	Crop	ID	Soil Type	Irrigation System	Crop
0	Exc.	n.a.	n.a.	55	II	Drip	Chickpeas and peas
1	Exc.	Fallow	Fallow	56	II	Drip	Winter vegetables
2	Exc.	Drip	Almond, walnut, and pistachio	57	II	Drip	Summer vegetables
3	Exc.	Drip	Chickpeas and peas	58	II	Drip	Olive groves
4	Exc.	Drip	Winter vegetables	59	II	Drip	Other crops
5	Exc.	Drip	Summer vegetables	60	II	Drip	Orchards
6	Exc.	Drip	Olive groves	61	II	Drip	Tomato
7	Exc.	Drip	Other crops	62	II	Drip	Vineyards
8	Exc.	Drip	Pastures	63	II	Travelling gun	Chickpeas and peas
9	Exc.	Drip	Orchards	64	II	Center pivot	Winter cereals
10	Exc.	Drip	Tomato	65	II	Center pivot	Rapeseed
11	Exc.	Drip	Vineyards	66	II	Center pivot	Fodder crops
12	Exc.	Travelling gun	Chickpeas and peas	67	II	Center pivot	Sunflower
13	Exc.	Center pivot	Winter cereals	68	II	Center pivot	Chickpeas and peas
14	Exc.	Center pivot	Rapeseed	69	II	Center pivot	Winter vegetables
15	Exc.	Center pivot	Fodder crops with multiple cuts	70	II	Center pivot	Summer vegetables
16	Exc.	Center pivot	Fodder crops	71	II	Center pivot	Maize grain
17	Exc.	Center pivot	Sunflower	72	II	Center pivot	Maize silage and sorghum
18	Exc.	Center pivot	Chickpeas and peas	73	II	Center pivot	Pastures
19	Exc.	Center pivot	Winter vegetables	74	II	Flooded paddies	Rice
20	Exc.	Center pivot	Maize grain	75	II	Rainfed crops	Fodder crops
21	Exc.	Center pivot	Maize silage and sorghum	76	II	Rainfed crops	Maize grain
22	Exc.	Center pivot	Pastures	77	II	Rainfed crops	Olive groves
23	Exc.	Flooded paddies	Rice	78	III	n.a.	n.a.
24	Exc.	Rainfed crops	Fodder crops	79	III	Fallow	Fallow
25	Exc.	Rainfed crops	Almond, walnut, and pistachio	80	III	Drip	Almond, walnut, and pistachio
26	Exc.	Rainfed crops	Maize grain	81	III	Drip	Chickpeas and peas
27	Exc.	Rainfed crops	Olive groves	82	III	Drip	Winter vegetables
28	Exc.	Rainfed crops	Other crops	83	III	Drip	Summer vegetables
29	Exc.	Rainfed crops	Vineyards	84	III	Drip	Maize grain
30	I	n.a.	n.a.	85	III	Drip	Olive groves
31	I	Fallow	Fallow	86	III	Drip	Other crops
32	I	Drip	Almond, walnut, and pistachio	87	III	Center pivot	Pastures
33	I	Drip	Chickpeas and peas	88	III	Drip	Orchards
34	I	Drip	Summer vegetables	89	III	Drip	Tomato
35	I	Drip	Olive groves	90	III	Drip	Vineyards
36	I	Drip	Other crops	91	III	Travelling gun	Chickpeas and peas
37	I	Drip	Orchards	92	III	Center pivot	Winter cereals
38	I	Drip	Tomato	93	III	Center pivot	Rapeseed
39	I	Drip	Vineyards	94	III	Center pivot	Fodder crops with multiple cuts
40	I	Travelling gun	Chickpeas and peas	95	III	Center pivot	Fodder crops
41	I	Center pivot	Winter cereals	96	III	Center pivot	Sunflower
42	I	Center pivot	Fodder crops with multiple cuts	97	III	Center pivot	Chickpeas and peas
43	I	Center pivot	Fodder crops	98	III	Center pivot	Winter vegetables
44	I	Center pivot	Sunflower	99	III	Center pivot	Maize grain
45	I	Center pivot	Summer vegetables	100	III	Center pivot	Maize silage and sorghum
46	I	Center pivot	Maize grain	101	III	Center pivot	Pastures
47	I	Center pivot	Maize silage and sorghum	102	III	Flooded paddies	Rice
48	I	Center pivot	Pastures	103	III	Rainfed crops	Fodder crops
49	I	Flooded paddies	Rice	104	III	Rainfed crops	Almond, walnut, and pistachio
50	I	Rainfed crops	Fodder crops	105	III	Rainfed crops	Maize grain
51	I	Rainfed crops	Olive groves	106	III	Rainfed crops	Olive groves
52	II	n.a.	n.a.	107	III	Rainfed crops	Other crops
53	II	Fallow	Fallow	108	III	Rainfed crops	Vineyards
54	II	Drip	Almond, walnut, and pistachio

Note: n.a.: not applicable since it concerns non-agricultural areas.

in the Appendix A. As previously pointed out, the soil water balance model was run independently for each HUA and the results were spatialized using the GIS platform.

3.5. Soil Water Balance Modeling Results

3.5.1. Model Validation for Soil Water Storage Predictions

The soil water balance simulations were validated using soil moisture time series obtained, for the study years, in three farms planted with maize, processing tomato, olive groves. Based upon the soil moisture measured within the rootzone with capacitance probes presenting sensors each 10 cm, the soil water storage (ASW) was calculated and compared with the simulated values, as shown in Figure 14. It can be observed that the model was able to adequately predict the behavior of the ASW along the seasons for the different crops and soils. In the example, maize was cropped in a loamy sand, tomato in a silty loam, and olives in a clay loam. Furthermore, the model was able to simulate well the impacts on ASW of the diverse irrigation schedules applied by the farmers along the different seasons.

The goodness-of-fit indicators relative to the ISAREG model test using farmers’ field observations of ASW are presented in

Table 5. “Goodness-of-fit” indicators obtained by comparing measured and simulated available soil water for the diverse crops during the 2018–2020 season parameters.

Crop	Year	Observations	b0	R2	PBIAS (%)	RMSE (mm)	NRMSE (%)	EF
Maize	2018	57	1.00	0.77	−0.1	3.1	2.5	0.76
	2019	71	1.00	0.82	0.6	6.8	5.9	0.80
	2020	78	1.01	0.84	1.2	5.8	5.0	0.83
Olive grove	2018	328	1.03	0.85	4.0	18.6	7.3	0.80
	2019	364	1.01	0.93	1.1	13.6	5.6	0.92
	2020	348	1.00	0.80	−0.2	17.5	6.2	0.77
Tomato for	2019	121	0.99	0.90	−1.2	5.5	3.4	0.88
processing	2020	99	0.99	0.87	−0.7	6.1	3.9	0.86

Note: b₀—coefficient of regression forced to the origin; R²—coefficient of determination; PBIAS—percent bias; RMSE—root mean square error; NRMSE—normalized root mean square error; EF—modeling efficiency.

. Results confirmed the good accuracy of the model in predicting the observed ASW after proper model parameterization. For maize crops, the model shows no tendency for under- or over-estimation of ASW with b₀ close to 1.00 and quite low PBIAS (0 < PBIAS < 1.5%). Results also show that the estimated ASW values were statistically close to the measured ones. The R² is high, indicating that the model could explain most of the observed variance. The errors of estimates are quite small: RMSE ranging from approximately 3 to 7 mm and the NRMSEs are lower than 6%. Finally, the EF values are higher than 0.76; thus, indicating that the variance of residuals is much lower than the variance of the measured ASW data. Similar results were obtained for maize with other applications of SWB models in Portugal (e.g., [76,77]).

Very good modeling results were obtained using the two data sets of ASW_OBS in tomato (Table 5) with b₀ close to 1.0 and low PBIAS values; thus, the model shows no trend for under- or over-estimation. R² values range from 0.87 to 0.90; thus, higher than those obtained for maize. Estimation errors were low, with RMSE close to 6 mm and NRMSE lower than 4%. EF results are high (EF > 0.86) confirming that the variance of residuals is much lower than the variance of measured data. These results are comparable with those in the literature for the Mediterranean (e.g., [78]).

For the olive grove, the model was run for the growing and non-growing seasons showing very good accuracy of the model along the years (Figure 14b). Results in Table 5 show just a slight tendency for the model to over-estimate the ASW_OBS by 2018 (b₀ = 1.03 and PBIAS = 4%) while in the following years the model showed no tendency to under- or over-estimate. Errors of estimates were also relatively low with RMSE < 19 mm and representing less than 8% of the mean of the ASW_OBS. The EF values are generally high. Few studies are available in literature relative to the use of SWB models for predicting ASW in olive orchards. The results of the current study are better than those reported using the WABOL model for an olive orchard in Spain [79].

Overall, results show that the parameterization of the ISAREG model was adequate and that the model was able to simulate the soil water storage along the seasons for the three studied crops. Thus, it can be assumed that the model may be further explored for estimating the soil water balance for the other crops present in the CIS.

3.5.2. Irrigation Water Requirements

The spatio-temporal distribution of crop irrigation water requirements (IWR), computed for each HUA (Figure 13), is presented in Figure 15, while Figure 16a,b shows the IWR per crop and per year, respectively.

The crop with the highest IWR was paddy rice due to the use of flooding irrigation, with lower application efficiencies than the other systems resulting from the large evaporation losses. Regarding the crops irrigated by pressurized systems, maize presents the highest IWR, with values above 750 mm during the 2017 and 2019 seasons, while lower values of IWR were estimated for 2018 and 2020. The higher IWR by 2017 and 2019 relates to the occurrence of less precipitation during the spring combined with higher ET_o (Figure 16). The relatively large amplitudes of the boxes associated with each crop type (Figure 16a) are related to crop management practices (e.g., plant density, varieties) as well as the cultivation in different soil classes, which influences the water balance results. The second highest IWR are those for fodder crops, followed by tomato for industry and vegetable crops. The high IWR of fodder crops are, in part, associated with the lower efficiency of the traveling gun irrigation system, as compared with the center pivot and drip irrigation (see Section 2.3). For the lowest IWR, there are the olive groves, averaging 275 mm, and legume crops averaging 180 mm. The IWR are also influenced by the soil depth that constraints the soil available water for the crops. In the case of spring–summer vegetable crops, IWR showed a reduced spatial variability, which may be due to the shallow rooting system, so the depth of the soil does not impose limitations on the root system and therefore does not affect the IWR. Differently, in the case of maize and olive groves, there was a greater spatial variation in the IWR. In the olive grove, this spatial variation was most noticeable in 2017, because it was a precipitation scarce year. Thus, in the cases that orchards were installed in shallower soils, the IWR were higher, and irrigation started by early March.

Figure 16b shows the IWR for each irrigation campaign, considering all crops together. The higher IWR in 2017 relates to both the available water from precipitation as well as the climatic demand (ET_o). On one hand, during the spring of 2017, precipitation was lower than the historical average for the months of April and May, when the spring crops are planted (Figure 2). On the other hand, the climatic demand (ET_o) was higher than the average for April, June, and September, increasing crop evapotranspiration. During the winter of 2019, precipitation was very low, leading to a small amount of water storage in the soil profile at the time of planting in the spring, therefore, the irrigation season began earlier, as in 2017, leading to higher irrigation requirements when compared with 2018 and 2020.

3.5.3. Deep Percolation

The results relative to the deep percolation associated with the irrigation events for the set of selected crops is shown in Figure 16c. The null values correspond to the dry years of 2017 and 2019. Results show that the median deep percolation value is lower than the average value (right skewed) which results from a longer tail in the high deep percolation values relative to the low ones. The relatively large amplitude of the boxes is due to the different classes of soil suitability for irrigation and, in the case of legume crops, due to the different irrigation systems used (center pivot, travelling gun, and drip irrigation) and respective application efficiencies. For the same crop, soil class, and year, the HUA with a traveling gun system presented the highest deep percolation values, 32% higher than for those using drip irrigation systems.

Considering all the years, crops and soil classes, deep percolation associated with irrigation averages 11 ± 3% of the total deep percolation that occurred during the year.

Crops associated with higher values of deep percolation resulting from irrigation are maize and vegetable crops, and at the same time, these crops present the largest variability of the results when all years and soil classes are considered. Farmers tend to apply more water to these crops since they are the ones presenting the highest economic returns. The opposite is found for sunflower and fodder crops. The irrigation schemes for olive groves are characterized by a period of low or null application of water, during summer dormancy, but prior to this period, farmers tend to irrigate more and even beyond soil field capacity and therefore some deep percolation occurs [67,80].

Regarding land suitability for irrigation, the classes with higher amounts of deep percolation associated with the irrigation events are III and excluded, due to the shallower depth of the soil profile and lower water retention capacity (Figure 16d).

The spatial variability of deep percolation that occurs during the periods when the crop is in the field is shown in Figure 17a. Overall, during the crop season, the sum of irrigation and precipitation are responsible for 29 ± 13% of the deep percolation, while the precipitation that occurred during the spring and autumn (off-season), typical of the Mediterranean climate was responsible for the higher share of total deep percolation. The spatialization of the deep percolation calculated for each HAU shows that there are some hotspots for deep percolation during the studied years. In this study, it was considered that deep percolation hotspots were associated with HUA with deep percolation higher than 200 mm during the crop season. The 2016/2017 season presented a low hotspot area corresponding to 10% of the total irrigated area, although 26% of the irrigated area presents deep percolation higher than 150 mm during the crop-growing season. Differently during the 2017/2018 season, the deep percolation hotspot area increased to 21%. Results show that 55% of these hotspot areas corresponding HUAs include soils with land suitability to irrigation from classes exc and III, mainly due to shallow depth and lower water retention capacity. In these HUA the crop is not a differentiating factor. However, 26 and 18% of the HUAs include classes II and I, classified as good and excellent, respectively.

In this case, crops with higher irrigation water requirements and higher economic return as maize and spring vegetables are predominant. The 2018/2019 season presents a low deep percolation hotspot area of 6%, like that for 2016/2017. During 2019/2020, 10% of the area is a deep percolation hotspot, from which 56% present suitability classes exc and III and all types of crops, while for the 44% of hotspots on soil classes I and II the crops are maize and pastures. On average, irrigation was responsible for 26 ± 8% of deep percolation during the crop seasons.

Figure 17b presents the deep percolation that occurred during the off-season, after crop harvesting; thus, resulting from precipitation only. The temporal pattern shows the inter-annual variation in precipitation, which in the case of the humid year of 2018 resulted in almost the entire CIS area being a hotspot. Spatially, the areas with lower deep percolation correspond to land suitability classes very good and excellent, i.e., classes I and II. On average, during 2018, deep percolation represented 54 ± 6% of the average annual value. High deep percolation during autumn–winter is favorable since it recharges the Caia aquifer. However, the Caia Irrigation Scheme is within a designated vulnerable zone, overlaying a nitrate-contaminated aquifer. Therefore, if a large amount of N remains in the soil profile after the harvest of the spring–summer crops, it will be leached to the aquifer with the first autumn precipitation. The soil water balance model predicted an amount of deep percolation associated to precipitation of about 23 ± 6% of the total precipitation, which agrees with other studies performed in the area [6].

Deep percolation in the referred hotspots leads to low efficiency of irrigation. Furthermore, it affects the aquifer quality, since drainage water is responsible for the transport of nitrates resulting from fertilization and crop residue mineralization.

4. Conclusions

An integrated methodology was developed to estimate the spatio-temporal dynamics of deep percolation with the farm information and crop coefficient curves being derived from high-resolution satellite images processed on the Google Engine platform.

The identification and characterization of crop cycles, fundamental for the quantification of crop irrigation water requirements, was based on the treatment of several high-resolution Sentinel-2 temporal images on the Google Earth Engine platform. In a first analysis of the NDVI temporal profiles extracted for each plot, in the GEE platform, the atypical profiles for each crop were identified and discarded. This method also allowed improvement of the identification of the crop growth stages and the planting dates at the different plots. This information was used to feed the soil water balance model, ISAREG, used to compute deep percolation below the root zone, for each homogenous unit of analysis.

The use of the proposed integrated methodology showed a good accuracy when it was validated using observed soil water storage data.

Results from the spatio-temporal analysis of the computed deep percolation along the years allowed to answer the question raised at the beginning of the present study which was: are there hotspot areas in the Caia Irrigation Scheme with high deep percolation due to irrigation practices? The results show that deep percolation during the crop season was quite variable among years, with average values ranging between 77.3 and 119.2 mm. However, these values present a very high spatial variability, ranging from 0 to more than 200 mm. Deep percolation was higher during the 2017/2018 campaign mainly resulting from precipitation due to the occurrence of a high amount of rainfall in a low number of events, which in parallel led to reduced irrigation water requirements. Deep percolation associated with irrigation averaged 11 ± 3% of total deep percolation. On one hand, the deep percolation hotspots associated with the irrigation events correspond to soils with low land suitability for irrigation, independently of the crop, probably due to the shallower depth of the soil profile and median low water retention capacity. On the other hand, maize and spring vegetables stand out when the hotspots correspond to land classes I and II.

Overall, during the crop season, irrigation plus precipitation was responsible for 29 ± 13% of the deep percolation, while the precipitation that occurred during the spring and autumn, typical of the Mediterranean climate, was responsible for the higher share. On average, during the off-season period, deep percolation corresponded to 54 ± 6% of the average annual precipitation. Thus, aiming at controlling the nitrates leaching in this vulnerable area, farmers need to have appropriate support to manage the irrigation scheduling according to the soil, crop, and weather conditions as well as to manage fertilization.