Simulating water lateral inflow and its contribution to spatial variations of rainfed wheat yields

European Journal of Agronomy 137 (2022) 126515 1161-0301/© 2022 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/bync-nd/4.0/). Simulating water lateral inflow and its contribution to spatial variations of rainfed wheat yields Tom´ as R. Tenreiro a , b , * , Jakub Jeˇ r´ abek c , Jos´ e A. G´ omez a , David Zumr c , Gonzalo Martínez d , Margarita García-Vila b , Elías Fereres a , b a Institute for Sustainable Agriculture (CSIC), C´ ordoba 14004, Spain b Department of Agronomy, University of C´ ordoba, C´ ordoba 14014, Spain c Department of Landscape Water Conservation, Czech Technical University in Prague, Prague 16629, Czech Republic d Department of Applied Physics, University of C´ ordoba, C´ ordoba 14071, Spain ARTICLE INFO Keywords: Crop modeling Spatial modeling Water-balance Lateral flows HYDRUS AquaCrop Machine Learning Artificial Neural Network ABSTRACT Spatial variations of crop yields are commonly observed in typical rainfed systems worldwide. It is accepted that such variations are likely to be associated, among other factors, with water spatial variations due to lateral water flows occurring in fields with undulating topography. However, some of the main processes governing water spatial distribution such as lateral flow are not entirely considered by the most commonly adopted crop simulation models. This brings uncertainty to the process of yield simulation at field-scale, especially under waterlimited conditions. Although it is expected that lateral water movement determines spatial variations of crop yields, it is still unclear what is the net contribution of lateral water inflows (LIF) to spatial variations of rainfed yields in fields of undulating topography. In this sense, by combining field experimentation, simulation models (HYDRUS-1D and AquaCrop), and the use of artificial neural networks, we assessed the occurrence and magnitude of LIF, and their impact on wheat yields in Cordoba, Spain, over a 30-year period. Seasonal precipitation varied over 30 years from 212.8 to 759.5 mm, and cumulative LIF ranged from 30 to 125 mm. The ratio of seasonal cumulative LIF divided by seasonal precipitation varied from 10.7% to 38.9% over the 30 years. The net contribution of LIF to spatial variations of rainfed potential yields showed to be relevant but highly irregular among years. Despite the inter-annual variability, typical of Mediterranean conditions, the occurrence of LIF caused simulated wheat yields to vary +16% from up to downslope areas of the field. The net yield responses to LIF, in downslope areas were on average 383 kg grain yield (GY) ha −1 , and the LIF marginal water productivity reached 24.6 ( ±13.2) kg GY ha −1 mm −1 in years of maximum responsiveness. Decision makers are encouraged to take water spatial variations into account when adjusting management to different potential yielding zones within the same field. However, this process is expected to benefit from further advances in inseason weather forecasting that should be coupled with a methodological approach such as the one presented here. 1. Introduction Modern agriculture aims to optimize the efficiency and profitability of farming systems while sustaining the increase in food production needed for a growing population (Connor and Mínguez, 2012; Fischer and Connor, 2018; Kirkegaard and Hunt, 2010). In most typical rainfed Abbreviations: ANN, artificial neural network; CC, canopy cover, expressed in %; CP, capacitance probe; CUM. LIF, season cumulative lateral inflow (LIF), expressed in mm; CUM. P, season cumulative precipitation, expressed in mm; DEM, digital elevation model (raster); FAI, flow accumulation index (the absolute number of upslope cells flowing to each assigned cell of the DEM raster); GY, grain yield, expressed in Mg (dry mass) ha −1 ; LIF, lateral inflow, expressed in mm; LIF. MWP, LIF marginal water productivity (expressed in kg GY ha −1 mm −1 ); NFAI, normalized flow accumulation index; NP, neutron probe; NYRLIF, net yield response to LIF, expressed in Mg GY ha −1 ); P, daily precipitation, expressed in mm day −1 ; KSAT, saturated hydraulic conductivity, expressed in mm day −1 ; PF. LIF, postflowering LIF (the fraction of CUM. LIF taking place at post-flowering stages), expressed in %; Relative. T, mean relative crop transpiration (estimated as the season average of daily crop actual transpiration divided by potential transpiration), expressed in %; SWC, soil water content, expressed in mm. * Corresponding author at: Institute for Sustainable Agriculture (CSIC), C´ ordoba 14004, Spain. E-mail address: [email protected] (T.R. Tenreiro). Contents lists available at ScienceDirect European Journal of Agronomy journal homepage: www.elsevier.com/locate/eja https://doi.org/10.1016/j.eja.2022.126515 Received 5 July 2021; Received in revised form 11 March 2022; Accepted 11 April 2022 European Journal of Agronomy 137 (2022) 126515 2 systems worldwide, where spatial variations of crop yields are commonly observed (Bramley, 2009; Maestrini and Basso, 2018; Sadler and Russell, 1997; Sadras and Bongiovanni, 2004; Sida et al., 2021), there is an opportunity for increasing productivity of resource use by determining the management options to exploit the site-specific conditions within fields (Cassman, 1999; McBratney et al., 2005). Site-specific variations of crop yields caused by differences in water availability due to lateral inflow from up to downslope areas have been identified (Ciha, 1984; Batchelor et al., 2002; Halvorson and Doll, 1991; Rockstom and Valentin, 1997; Schmitter et al., 2015). However, the contribution of lateral inflow to spatial variations of yields has not been systematically explored. The intra-plot heterogeneity associated with lateral water flows has implications in input allocation, allowing for spatial variations in crop management in the context of precision agriculture (Ahuja et al., 2019; Nielsen et al., 1973; Sadler and Russell, 1997; Verhagen et al., 1995; Wallor et al., 2018; Ward et al., 2018). Precision agriculture would surely benefit from advances in the spatial simulation of water variations over fields. However, to model accurately rainfed yields in fields of undulating topography, we need simulation tools capable of forecasting spatial variations in water availability within a field for assessing their impact on crop performance. Over the last decades, there has been great expansion in modeling agricultural processes at the point scale (Jones et al., 2017; Spiertz, 2014), but insufficient efforts have been devoted to scale up water-related processes, which vary spatially, from point to field level (Ahuja et al., 2019; Wallor et al., 2018). Tenreiro et al. (2020) have recently reviewed some of the most widely adopted crop and hydrologic models and the main opportunities to simulate spatial water variations at crop field level, and concluded that the most promising strategies for scaling up are related to the incorporation of both surface and subsurface lateral flows when simulating crop performance. The incorporation of surface and subsurface lateral flows within simulation modeling requires innovative approaches and datasets, which should be spatially distributed and related to the geomorphological properties with implications for plant available water (Nielsen and Wendroth, 2003; Wallor et al., 2018). However, data collection requires field experimentation conducted at “real scales”, which are relatively expensive and difficult to replicate over long periods of time. Therefore, the combination of both experimentation and modeling is a valid strategy for making progress (Jones et al., 2017; Kamilaris et al., 2017; Toreti et al., 2018; Wolfert et al., 2017). Considering the typical inter-annual variation of the processes governing crop-water spatial relations (de Wit and van Keulen, 1987), the present study investigated the following question: “what is the net contribution of lateral water inflows to spatial variations of rainfed wheat yields in fields of undulating topography? Fig. 1 illustrates graphically the relevance of our research question. To address our research question, we developed a novel methodology to explore the linkage between lateral inflows (including both surface and subsurface flows) and yield variations in specific zones within a field. By combining field experimentation, simulation models and the use of artificial intelligence, we assessed the occurrence and magnitude of lateral inflows, and their impact on wheat yields in Cordoba, Spain, over a 30-year period. 2. Materials and methods 2.1. Experimental sites The experimental sites consisted of two nearby hydrologically independent fields, located in Cordoba, southern Spain (37.8 ∘ N, 4.8 ∘ W, mean altitude 170 m amsl.) of 42 and 36 ha, respectively. Two catchment areas within the selected fields, 9.5 and 6.2 ha respectively, were delineated from a flow direction raster obtained with the SAGA - Wang Liu algorithm (Wang and Liu, 2006), from a Digital Elevation Model (DEM) with 5 m spatial resolution collected with LiDAR (CNIG, 2019). The soils are formed on Miocene marls and have been classified as Typic Haploxerert (Soil Survey Staff, 1999), or Vertisol according to the FAO classification, characterized by clay texture with a shrinking-swelling nature, high bulk density (1.65–1.88 g cm −3 ) and of 1.2–1.6 m depth. The catchments have mean slopes of 2–6%, respectively N-S and W-E oriented, and with elevation varying from 140 to 195 m. Due to crop rotations, wheat was monitored in catchment one in 2019/2020 (Fig. 2-A and 2-C), and in catchment two in 2020/21 (Fig. 2-B and 2-D). 2.2. Sampling scheme and experimental design 2.2.1. Geomorphological properties and sampling points The spatial variation of soil geomorphological properties was characterized in catchment one by using an electromagnetic induction sensor (DUALEM-21S) to measure soil apparent electrical conductivity (ECa, dS m −1 ) within the top 0–50 and 0–90 cm soil layers, four days after a rainfall event of approximately 10 mm (McCutcheon et al., 2006) and before sowing. Soil samples (%Clay, %Sand, pH) were collected at 35 cm depth following a multistage sampling scheme that was based on two different ECa-based clusters. In catchment two, soil properties were averaged for the entire field using farm records (Table 1). Topographic attributes were computed with SAGA-GIS (version 2.3.2) from the DEM raster. Soil moisture sampling zones were delineated according to both elevation and the flow accumulation index (FAI). The flow accumulation index (FAI) is expressed as the absolute number of upslope cells flowing to each assigned cell of the DEM raster (Tarboton et al., 1991; Jenson and Domingue, 1988). Since it is dependent on field scale and input data spatial resolution, a normalization of the index (NFAI) was computed as follows: NFAI =1−FAIMAX −FAI FAIMAX (1) The NFAI was estimated with SAGA-GIS (Conrad et al., 2015) and it was represented as a raster with the same spatial resolution of the input DEM. Three sampling zones (N =3) were distributed along two converging flow pathways in 2019/20 (with 2–3 replicate samples of soil moisture per zone, Fig. 2-A and -C), and along the same flow channel in 2020/21 (with three replicates per zone, Fig. 2-B and 2-D). In catchment one (2019/20), the sampling points (2–3 per zone) were positioned according to the peak of the mean ECa histogram, which was done to select sites of maximum representativeness of soil properties within each zone. In catchment two (2020/21), where soil properties were averaged, the sampling points were simply positioned within each zone according to the flow direction and spaced in 1 m intervals. Fig. 1. Visual symptoms of early crop senescence apparently caused by spatial water variations in Cordoba, Spain. The higher elevation zones show yellowing patterns due to lower water availability, which limits crop yield (Sadras et al., 2016). Photo credits: T. R. Tenreiro. T.R. Tenreiro et al. European Journal of Agronomy 137 (2022) 126515 3 2.2.2. Rainfall and meteorological data Rainfall was monitored upslope with an autonomous rain gauge system (ECRN-100, ZENTRA Cloud ZL6, 16 cm collector diameter), with 10 min time resolution. The rain gauge was installed in an intermediate site, located 600–900 m from each catchment, at a spot above the highest point of each catchment (Appendix-A1). In 2019/20, manual pluviometers (TFA 47.1008, 12 cm collector diameter) were positioned at each observation point (N =7) to capture rainfall coefficient of variation. The coefficient of variation was computed in relation to the autonomous rain gauge measurements for two separate rainfall events before crop emergence in 2019/20. Weather data (Fig. 3) were obtained from a meteorological station nearby, located less than 5 km away from each field (Appendix-A1). These included global radiation, wind speed, air temperature and relative humidity, which were daily averaged over half hourly measurements, and used to compute daily mean ETo values according to FAO Penman-Monteith (Allen et al., 1998). 2.2.3. Soil water content and lateral inflow calculations Soil water content (SWC), expressed in mm, was monitored with multisensor capacitance probes (SENTEK-D&D 90 cm, Sentek Technologies Ltd., Australia), installed at each observation point (N =7 in Fig. 2. Maps of the experimental catchments: the elevation maps with contour lines of catchment one and two, respectively (A and B); the Normalized Flow Accumulation Index (NFAI) rasterized with 5 m spatial resolution for catchment one and two, respectively (C and D). Sampling zones (A, B, C) and sampling points (A1, A2, A3, B1, B2, B3, C1, C2, C3) are represented by solid black lines and purple circles, respectively. T.R. Tenreiro et al. European Journal of Agronomy 137 (2022) 126515 4 2019/20 and N =9 in 2020/21). Probes were respectively installed at day after sowing (DAS) 20 and 5 in 2019/20 and 2020/21. Each probe integrates nine sensors, distributed at every 10 cm, from 5 to 85 cm depth. The sensors provide near real-time data on soil temperature and SWC (30 min time resolution), which were accessed through the software IrriMAX Live (www.irrimaxlive.com/). An error assessment was conducted for each probe, based on discrete SWC measurements (N =20 date ×depth(z) per observation point) made with a neutron probe in 2019/20 (NP, Campbell Pacific Nuclear Scientific, Model 503). A relative error assessment of the capacitance probes (CP) was conducted by taking the NP measurements as control. The NP access tubes were installed 40 cm apart from the CP. One tube per observation point was considered sufficient to meet the requirements for precision and statistical power (Evett et al., 2009). NP measurements were taken at five different dates, and at 15, 30, 60 and 90 cm depth, and those values were correlated with the capacitance probe sensors located at 15 cm depth, the average of 25–35 cm, the average of 55–65 cm, and the last sensor located at 85 cm depth, respectively. The NP was calibrated with gravimetric measurements of SWC, obtained in the same soil type in a farm nearby (Appendix-A1). The NP calibration linear functions varied from depth to depth (Appendix-Table A) and were characterized by R 2 values of 0.96–0.98 for deep soil layers (30–90 cm) and 0.81–0.85 for the surface depth (0–30 cm). Multiple calibration functions were tested in the IrriMAX Live software in order to minimize error fluctuations per probe and depth. The ‘Sentek D&D cracking clays’ calibration function, available in IrriMAX Live software, was chosen as the most suitable option for our soil type (Paltineanu and Starr, 1997; RoTimi Ojo et al., 2015). IrriMAX Live values (i.e., from CP measurements) were corrected with the NP measurements as follows: SWCc=∑ z=90cm z=0cm (SWCrz⋅Ωz)(2) where SWC c represents the corrected SWC (expressed in mm), for the entire profile and SWCr z represents the SWC value (in mm) provided by the IrriMAX calibration function from raw input data, measured with the CP at depth z. Ω z is the correction factor (unitless) estimated for each probe ×depth(z) combination and it was computed as the mean ratio between the SWC measurements with the NP and the capacitance probes for depth z. The same calibration was used in 2020/21. Daily values of SWC were computed for each observation point as the simple average of all SWC c values registered within the same day. Lateral inflows (LIF) were calculated in daily time-steps for each observation point, from daily values of SWC c . Daily variations of SWC c were computed as the difference between SWC c in day n and day n−1. Observed lateral inflows were assumed to be the absolute difference between the daily SWC c variation and the daily rainfall registered by the rain gauge, which was computed as follows: SWCc(n)−SWCc(n−1)>P(n−1)⇒LIFn=SWCc(n)−SWCc(n−1)−P(n−1)(3) Every day SWC c varied by an amount greater than the registered Table 1 Main geomorphological attributes within each sampling zone (mean values and standard deviations). The standard deviations are presented in parentheses. Codes: ECa =soil apparent electrical conductivity, %Clay =percentage of clay content, %Sand =percentage of sand content, NFAI =“Normalized flow accumulation index” (Tarboton et al., 1991; Jenson and Domingue, 1988). The amount of chosen sampling points was supported by Chanzy et al. (1998). %Clay, %Sand and bulk density were averaged for the three zones in catchment-2 according to farm records. Parameter Units Catchment-1 Catchment-2 Zone A B C A B C Sampling points [N] 2 3 2 3 3 3 Elevation m (amsl.) 187 (2.4) 184 (5.1) 168 (2.3) 147 (0.5) 141 (0.3) 139 (0.3) ECa dS m −1 0.31 (0.08) 0.28 (0.09) 0.48 (0.06) – – – %Clay % 45 (3.4) 42 (3.3) 50 (3.4) 44 (2.8) %Sand % 18 (2.6) 22 (2.9) 15 (2.8) 22 (3.2) Bulk density g cm −3 1.78 (0.06) 1.81 (0.04) 1.74 (0.05) 1.66 (0.05) NFAI [0;1] 0.09 (0.021) 0.12 (0.082) 0.72 (0.218) 0.07 (0.005) 0.42 (0.002) 0.98 (0.011) Fig. 3. Daily mean values for (A) rainfall (dashed lines represent cumulative rainfall computed from the 1st of October, values shown by right vertical axis), (B) temperature and (C) reference evapotranspiration (ETo), respectively expressed in mm, degree Celsius and mm. Lines represent season time-series expressed in terms of days after sowing (DAS). Gray areas represent the monitoring time window of each observation year (i.e., DAS 20–157 in 2019/20 and DAS 5–158 in 2020/21). Daily mean temperature values are represented by the heavy solid lines. Maximum and minimum temperature values are respectively represented by a solid and a dashed skinny line. Each season values are shown in different colors (2019/20 in green and 2020/21 in brown). T.R. Tenreiro et al. European Journal of Agronomy 137 (2022) 126515 5 precipitation (P), a lateral inflow was assumed to take place of a magnitude equivalent to the difference between the increase in SWC and the rainfall amount. For the calculation of LIF, deep percolation and evapotranspiration (ET) were not considered for the following reasons. In our soil type, deep percolation approaches zero (Giraldez and Sposito, 1985). In the case of ET, the crop water extraction following a rainfall event is quite limited, as intercepted water evaporates first from wet canopies as shown by Tolk et al. (1995). Therefore, considering their very small magnitudes in our case, we did not attempt to measure or estimate either deep percolation or ET, as it would have added uncertainty to our LIF calculations. Our approach determines the minimum LIF quantity that actually takes place under field conditions. Therefore, if under our approach LIF is relevant for determining crop yield, then it must play an even more important role in actual yield variations within a field. 2.2.4. Crop data and field observations The wheat cultivars KIKO-NICK-R1 and Avispa-R1 were sown in Table 2 Crop data used for the parameterization of the HYDRUS and the AquaCrop models. CC MAX is maximum green canopy cover (Steduto et al., 2009), used for parameterization of the soil cover partitioning method, as described in Tenreiro et al. (2020). Mean sowing rate was 180 kg ha −1 for all trials. The standard deviations are presented in parentheses. Data Catchment-1 Catchment-2 Mean Parameter Zone A B C A B C - [units] Sowing date date 13-Dec 18-Nov 1-Dec Crop emergence DAS 10 8 9 Plant density plants m −2 200 (21.8) 280 (19.2) 225 (16.3) 250 (18.2) 225 (16.44) 225 (10.4) 236 (27.4) CC MAX % 80 (3.5) 82 (4.6) 84 (2.1) 90 (1.7) 95 (0.6) 92 (0.5) 87 (5.9) Root growth cm day −1 0.7 (0.1) 0.9 (0.1) 0.8 (0.1) Vegetative stage days 120 120 120 Anthesis stage days 10 14 12 Reproductive stage days 58 76 67 Senescence duration days 20 40 30 Crop maturity date 8-Jun 1-Jun 4-Jun Harvest date date 10-Jun 8-Jun 9-Jun Fig. 4. Sketch of the methodological design: Two field experiments were conducted, each sampling point was defined according to a spatial analysis (step 0). A daily calendar of lateral water inflow (LIF) was calculated based on field observations (step 1) and simulations of LIF were conducted for 30 years through a hydrologic modeling approach (step 2). LIF predictions were assimilated into the crop modeling stage (step 3) and the net yield responses were simulated (step 4). Dashed lines delineate different methodological stages, rounded parallelograms indicate experimental sites and data, solid line rectangles indicate different sub-steps, solid line circles represent simulation tools, losenge and arrows indicate conditional steps. More information on the simulation settings of HYDRUS-1D and AquaCrop is respectively provided by ˇ Simunek et al. (2018) and Steduto et al. (2009). Additional details related to the use of artificial neural networks (ANN) for hydrological modeling are provided by Maier and Dandy (2000). T.R. Tenreiro et al. European Journal of Agronomy 137 (2022) 126515 6 catchment one (2019/20 season) and two (2020/21 season), respectively. Seeding rates were 180 ( ±20) kg ha −1 . Catchment one was fertilized with two applications of calcium ammonium nitrate plus sulfur (62 +50 kg N ha −1 ) and catchment two was fertilized with ammonium sulfate plus urea (160 kg N ha −1 ), and diammonium phosphate (60 kg P ha −1 ). Crop nutrient status was controlled through foliar analysis conducted at flowering, with no critical deficiencies observed. Ground measurements of canopy cover (CC) were conducted in both trials, every 10–20 days, during the monitoring period. CC was measured at each observation point using a digital camera (Canon EOS 550D +EFS 18–135 mm CMOS APS-C 18.7 MP) at 1.5 m height and an image processing package (Patrignani and Ochsner, 2015). CC curves were used to parameterize crop growth related factors in the simulations with both the HYDRUS and the AquaCrop models (Table 2). Crop stages duration were obtained from field observations of phenological development and adjusted according to the CC curves obtained from satellite NDVI, as described in detail by Tenreiro et al. (2021). Both seeding rates and site-specific plant density values were registered and used for modeling parameterization as well. Rooting depth trend was inferred in each observation point using the SWC information (Table 2). Maximum rooting depth was considered to be equal to average soil depth (1.4 m). Grain yield was harvested by combine, using the ’New Holland’ Precision Land Manager (PLM) software which took as an input the shapefiles generated by the combine harvester monitor (Fendt PLI C 5275). Yield values were computed with R-studio (Lovelace et al., 2019), under a spatial resolution of 100 m 2 with the equation of Reitz and Kutzbach (1996). The accuracy of the yield data from the combine monitor was assessed by comparing manual samples taken at each observation point in catchment one (sampled areas of 0.9 m 2 ) against the combined harvest data. More information regarding the yield spatial assessment is provided in Tenreiro et al. (2021). Yield observations were obtained from yield maps of catchment two (2015–2020 harvests). Zonal means and standard deviations were estimated and values were plotted against simulated yield values. 2.3. Modeling approach To simulate water lateral flow and its contribution to spatial variations of crop yield, a multi-stage modeling approach was designed (Fig. 4). The modeling approach was based on field experimental data (step-0, Fig. 4) and it was divided in four subsequent steps (steps 1–4, Fig. 4). Field experimental data were collected according to a spatial analysis based on GIS data aimed at identifying observation points according to standard hydrological connectivity rules (i.e., field channel networks and flow accumulation index), topographic attributes (i.e., elevation, slope orientation) and other geomorphological properties (step-0, Fig. 4). The first step consisted on lateral inflow calculations from field observation data, in daily time-steps and for each observation point (step-1, Fig. 4). Then, a hydrologic routine simulated the lateral outflows (with HYDRUS-1D, ˇ Simunek et al., 2018), generated upslope in form of surface run-off, and then the occurrence of lateral water inflows (LIF) was predicted at each sampling point (step-2, Fig. 4). For each sampling point, daily calendars of LIF were determined according to a hydrological analysis that combined field measured data with both HYDRUS simulations and a machine learning (empirically based) approach (step-2, Fig. 4). An Artificial Neural Network (ANN) was used to simulate LIF over a period of 30 years (1990–2020). Maier and Dandy (2000) reviewed in detail the main applications of ANN models for the prediction and forecasting of hydrological variables. The ANN model architecture was delineated according to a trial-and-error procedure (step-2, Fig. 4). Then, the outputs of the hydrology-based modeling routine (step-2, Fig. 4) were used as inputs to the crop-based modeling stage (step-3, Fig. 4). The crop routine incorporated the calendars of daily LIF values as additional water supply (step-3 and -4, Fig. 4). The forecasted LIF over a 30-year period, were then introduced as additional water supply through the irrigation module in wheat yield simulations with the AquaCrop model (Steduto et al., 2009). Both the HYDRUS and the AquaCrop modeling routines functioned at one dimensional space. Every simulation run was performed at one-dimension, but we followed a ‘feedforward scheme’, which allowed us to simulate and incorporate LIF throughout consecutive modeling steps. 2.3.1. Simulating the hydrology with HYDRUS The HYDRUS-1D model (ˇ Simunek et al., 2018) was used to simulate water infiltration and estimate surface run-off and the spatial variation of soil hydraulic properties. The HYDRUS-1D is a physically-based model that solves Richards’ equation for transient water transport in variably saturated porous media, and incorporates processes such as soil evaporation, crop transpiration, root growth, and plant water uptake (ˇ Simunek et al., 2018, Tenreiro et al., 2020). The standard ’van Genuchten-Mualem model’ was used to represent the soil hydraulic characteristics (Mualem, 1976; van Genuchten, 1980). The soil profile was modeled in one dimension at each of the measured points in both experimental years. Two soil layers were considered: surface and sub-surface. The surface layer was set at the first 10 cm of the modeled profile, while the sub-surface layer was set down to 140 cm depth. The meteorological conditions governing evaporative demand were set as the upper boundary condition. A ‘free drainage’ condition was considered at the bottom boundary of the soil profile because well drainage conditions with lack of soil reduction symptoms were observed at the BC horizon below 140 cm depth, during a pit excavation conducted prior to this study, and no indications of a water table were found. A steady state of SWC based on the measured data was used to set initial conditions (i.e., 0.20 cm 3 cm −3 in the season 2019/20 and 0.27 cm 3 cm −3 in the season 2020/21, following an average 60 day long warm-up period to minimize the interference of soil-water altered conditions due to probes installation). The potential transpiration was estimated through the soil cover partitioning method as described in Tenreiro et al. (2020). The same canopy cover curves used in simulations with the AquaCrop model were used for this purpose. Root growth rates were assumed to be linearly constant (Table 2). The rooting depth distribution function was based on the model of Hoffman and van Genuchten (1983) and the water stress function of Feddes et al. (1978), adapted for wheat and described in Wesseling et al. (1991), was considered. Since the inverse optimization procedure is sensitive to the rooting depth distribution, this was adjusted by means of SWC (Zumr et al., 2006). HYDRUS-1D parameterization was optimized using the measured soil water content (section 2.2.3) and the best fitted saturated hydraulic conductivity (K SAT ) values estimated for each soil layer at each observation point. The parameters were optimized by the Marquardt-Levenberg Optimization Algorithm (ˇ Simunek et al., 2012). The initial estimation of parameters was based on field measurements of soil texture and bulk density (Table 1) using the Rosetta neural network predictor (Schaap et al., 2001). Initial K SAT ranged according to van Genuchtens’ α and n shape parameters, which were estimated through the Rosetta predictions (more information is provided in Table A.1 in Tenreiro et al., 2021). The minimum and maximum values of saturated and residual water content were estimated from the measured SWC data. Parameter intervals used for the optimization are given in Table 3. Due to spatial variability in measured SWC, these ranges were further adjusted for each sampling Table 3 The parameter intervals used for model optimization. Values were adjusted during the optimization process for each individual point (Fig. 2). Parameter units Min Max residual vol. water content (θ r ) cm 3 cm −3 0 0.2 saturated vol. water content (θ s ) cm 3 cm −3 0.3 0.5 van Genuchten’s parameter ( α ) cm −1 0.0001 0.1 van Genuchten’s parameter (n) – 1.2 1.5 saturated hydraulic conductivity (K SAT ) mm day −1 2 100 T.R. Tenreiro et al. European Journal of Agronomy 137 (2022) 126515 7 point separately. Several initial estimations were used to better explore the parameter spatial dimension and to avoid falling in local minimum values (ˇ Simunek et al., 2018). The objective function minimum was found typically after 7–20 iterations of the Marquardt-Levenberg Optimization Algorithm. However, due to the numerical instabilities during several optimization runs, multiple non-converging model runs needed to be performed for each sampling point (i.e., typically 5–50 runs). A daily calendar of surface run-off was estimated for each sampling point and experimental year. 2.3.2. Forecasting lateral inflows Daily calendars of lateral inflows (LIF n , expressed in mm day −1 ) at each sampling point were computed as a function of the run-off generated at that same point, as well as the run-off flowing from neighboring and upslope points within the same field. The daily values of simulated run-off with HYDRUS were used as input features into a static multilayer (feedforward) ANN model (de Vos and Rientjes, 2005). Different performance criteria were used to determine network inputs and model architecture (Maier and Dandy, 2000). Input features (N =6) were selected according to a principal component analysis (PCA) associating different potential predictors (N =12) with LIF observations at each sampling point (Table 4). Model architecture was delineated according to a trial-and-error procedure (Roadknight et al., 1997; Senthil-Kumar et al., 2005; Shukla et al., 1996). The ANN prediction accuracy and its computation-training speed were assessed with R-studio (Günther and Fritsch, 2010). The ANN processes multiple algebraic operations over several input features which are defined by a single column vector (X →). Each of the inputs is attenuated by a weight factor (w) that is linked to a transfer function (i.e., a logistic transformation of the data). The calculation scheme of an ANN using i input features (computed in daily time steps n) and k hidden layers can be simplified as following: CUM.LIFn=∑ k⎛ ⎝ xn1 ⋮ xni ⎞ ⎠⋅⎛ ⎝ ( α 11⋅w11)…( α 1i⋅w1i) ⋮ ( α j1⋅wj1)…( α ji⋅wji) ⎞ ⎠+ϵ(4) where α corresponds to the transfer function operator and the subscript j delineates the number of nodes in each layer (de Vos and Rientjes, 2005). The bias coefficient (ϵ) corresponds to the overall error of the network, i.e., the sum of each node residuals, and the dependent variable is expressed in daily cumulative terms (CUM. LIF n ) to reduce the effects caused by temporal deviations between observations and predictions of LIF. The relation between daily (n) cumulative LIF (CUM. LIF n ) and daily LIF (LIF n ), at each point-site, computed from the beginning of simulation until day n=i is expressed as following: CUM.LIFn=∑ n=i 0 LIFn(5) Daily LIF n was therefore estimated as the absolute difference between CUM. LIF n at day n and CUM. LIF n−1 at day n−1: LIFn=CUM.LIFn−CUM.LIFn−1(6) Catchment one data (2019/20) were used to train the ANN model while catchment two data (2020/21) were used for testing. For each trial, ANN-LIF predictions were plotted against the observations of LIF (section 2.2.3) and the Willmott index of agreement (d), the R 2 and the RMSE were used as statistical indicators of ANN performance. The bestperforming ANN was used to forecast CUM. LIF n time-series over a period of 30 years, from which the LIF n values were derived for each zone. CUM. LIF n values were smoothed by a ’monotonically increasing function’ to preserve a strictly increasing pattern over season. The best performing ANN was applied to forecast daily LIF in catchment two over the same period of 30 years. The model HYDRUS-1D was used to simulate run-off at each sampling point and for each growing season [1990–2020]. Catchment two properties (Table 1) and mean values for crop data (Table 2) were used as input features to the best performing ANN. Weather records (i.e., global radiation, wind speed, air temperature, relative humidity and rainfall) for Cordoba [1990–2020] were obtained from the same weather station (AppendixA1). 2.3.3. The crop model stage with AquaCrop The AquaCrop v6.1 model (Steduto et al., 2009, 2012) was used to simulate net yield responses to lateral inflows (NYR LIF ). NYR LIF was assumed to be the absolute difference in terms of simulated (water-- limited) yield, between these two different scenarios: 1) yield simulation without lateral water inflow, following the standard assumption in AquaCrop where water inflow takes only place due to vertical infiltration (Y 1 );. 2) yield simulation including lateral inflow as an additional water supply, which was set in the model through the irrigation module and computed as the LIF forecasted by the ANN (Y 2 ). NYR LIF was calculated as following: NYRLIF =Y2−Y1(7) where all terms are expressed in Mg ha −1 . Yield response to lateral inflow was also estimated in relative terms (NYR LIF−rel ) as following: NYRLIF−rel =NYRLIF YMean (8) where Y Mean represents the mean yield estimated in scenario 1. A sequence of 30 years LIF [1990–2020] and its impact on grain yield was simulated to derive probability distribution functions of both NYR LIF Table 4 Artificial neural network (ANN) potential predictors assessed through principal component analysis (PCA) and trial-and-error procedure (Roadknight et al., 1997). Main PCA outcomes are shown in Appendix-A3. More information regarding the PCA is provided in Supplementary material. id Input feature Estimation procedure/ description Contribution to LIF x1 Cumulative surface run-off [mm] at location (x, y) Simulated with HYDRUS-1D through multiple iterations optimized by means of the measured soil water content High x2 Day after run-off at location (x, y) Computed with R-studio as a function of x1 calendarization Low x3 NFAI [0;1] at location (x, y) Estimated with SAGA-GIS (2.3.2) and computed with Rstudio (Lovelace et al., 2019) High x4 Slope [%] at site i Estimated with SAGA-GIS (2.3.2) and computed with Rstudio (Lovelace et al., 2019) Medium x5 Cumulative surface run-off at upslope contributing points [mm] The daily median of x1 values, simulated with HYDRUS-1D for upslope hydrologically contributing points High x6 Canopy cover [%] According to Tenreiro et al. (2021) High x7 Surface saturated hydraulic conductivity [mm day −1 ] Optimized K SAT at location (x, y) with HYDRUS-1D for surface layer (0–10 cm depth) Medium x8 Day after run-off Computed with R-studio as a function of x5 calendarization Medium x9 Day after precipitation event Computed with R-studio as a function of x11 calendarization Low x10 Saturated hydraulic conductivity (mean) [mm day −1 ] The average between surface and sub-surface K SAT (optimized with HYDRUS-1D) High x11 Cumulative daily precipitation [mm] Field measured with a rain gauge (section 2.2.2) High x12 Boolean ‘satunsat′ parameter Boolean parameter set to define days of saturated vs. unsaturated conditions in the vadose zone Medium T.R. Tenreiro et al. European Journal of Agronomy 137 (2022) 126515 8 and NYR LIF−rel . From the obtained series of NYR LIF , the LIF marginal water productivity (LIF. MWP) was computed as following: LIF.MWP =NYRLIF ×1000 CUM.LIF (9) where for each marginal unit of water supplied as LIF (expressed in mm), the crop was assumed to respond with additional grain yield (Passioura and Angus, 2010). LIF. MWP was expressed in kg grain ha −1 mm −1 . The AquaCrop model was parameterized with field data (catchment two values, Table 1) and crop data (mean values, Table 2). Long-term simulations were conducted for catchment two because it showed lower standard deviations for zonal means (Table 1 and 2). The hydraulic conductivity values used in AquaCrop simulations (K SAT mean, expressed in mm day −1 ) were the mean values reported in Table A.1 in Tenreiro et al. (2021). Two soil horizons were also considered (surface above 30 cm depth and sub-surface from 30 to 140 cm depth). The initial curve number was set at a value of 84 (i.e., hydrologic group D). More information regarding the water balance approach that is followed by AquaCrop found in Tenreiro et al. (2020). Simulated canopy growth and grain yield were validated through field observations (Section 2.2.4). 2.4. Statistical analysis Differences among sampling zones were tested for significance, the null hypothesis was checked for the mean differences of both observed SWC and LIF with the non-parametric Tukey’s HSD (honestly significant difference) test because these variables were not normally distributed. Non-normality was checked with the Shapiro-Wilk test (Acutis et al., 2012). Within field spatial variations (among the three sampling zones) were assessed with standard coefficients of variation. The residuals of the ANN features were checked to be randomly distributed (Supplementary material). HYDRUS simulations of SWC were tested with the Nash-Sutcliffe model efficiency coefficient, the R 2 and the RMSE (Moriasi et al., 2007; Nash and Sutcliffe, 1970; Yang et al., 2014). AquaCrop simulations of grain yield were tested with the Willmott d index, the R 2 and the RMSE (Willmott, 1981). Spearman’s correlation analysis was used to explore relationships among several crop variables and both NYR LIF and LIF. MWP. Table 5 Mean soil water content (SWC), mean daily lateral inflow (LIF n ) and mean cumulative LIF (CUM. LIF n ) over season, for the two years of field data. All values are expressed in mm. The standard deviations are presented in parentheses. Mean values followed by a common letter are not significantly different according to the HSD-test conducted at the 5% level of significance (p−value <0.05). Zone Catchment-1 (2019/20) Catchment-2 (2020/21) SWC [mm] A 209.36 (39.6)c 289.19 (65.1)b B 237.36 (34.4)b 316.52 (61.2)a C 269.37 (63.1)a 310.87 (60.4)a LIF n [mm] A 0.08 (0.46)b 0.13 (0.8)b B 0.31 (1.52)ab 0.23 (1.6)b C 0.42 (2.09)a 0.71 (2.7)a CUM. LIF n [mm] A 6.98 (3.9)b 11.24 (11.2)c B 17.88 (15.7)a 23.94 (17.4)b C 19.76 (16.9)a 60.28 (42.1)a Fig. 5. Corrected soil water content (SWC c ) for each sampling zone (A-C), error bars are shown for daily SWC c values. Daily lateral inflows (LIF) per sampling zone (inferred from probes data). Values are expressed in mm. T.R. Tenreiro et al. European Journal of Agronomy 137 (2022) 126515 9 3. Results 3.1. Experimental data CPs overestimated SWC in comparison with the NP measurements. The mean correction factor (Ω z ) varied with both probe location and sensor depth (Appendix-A2). The amplitude of error variation was greatest at the surface layers (0–30 cm) but it tended to stabilize for deeper layers (Appendix-A2). The soil water content (SWC c ) varied spatially in both experimental catchments/years (Table 5 and Fig. 5). According to the HSD-test conducted at a 5% level of significance, the annual mean SWC values varied significantly among sampling zones (Table 5). Our sampling scheme did also capture significant differences among zones in LIF values (Table 5), both in daily (LIF n ) and cumulative terms (CUM. LIF n ). However, as observed for SWC, zone B showed a ‘transient profile’, where both the SWC and LIF values were not ‘systematically and significantly’ different from the adjacent sampling zones. For an accurate distinction among zones within our experimental sets, only zones A and C were consistently different in terms of SWC (and LIF) variations (Fig. 2 and Table 5). The magnitude of both SWC and LIF values also varied among catchments and/or years (Table 5 and Fig. 5). Catchment one showed lower mean values than catchment two (Table 5), partly explained by the differences in initial SWC (Fig. 5). There might be also differences in soil properties among sampling zones influencing SWC measurements in catchment two, because clay and sand content, and bulk density were averaged for all zones (Table 1). These differences can also be attributed to the performance of the CP readings according to the calibration functions used. The mean coefficient of variation of daily rainfall, measured for two events in 2019/20, was 9.5%, which was lower than the LIF coefficient of variation (13%) estimated from probes data among zones in the same season. 3.2. HYDRUS run-off simulation and ANN-LIF forecasting The HYDRUS-1D simulations of SWC showed a mean Nash-Sutcliffe model efficiency of 0.61 ( ±0.14) and 0.92 ( ±0.04) for each of the experimental years, respectively (Fig. 6). In addition, the model has also increased the R 2 (while reducing the RMSE), from the first to the second year (Fig. 6). A higher variation of performance indicators was also observed in catchment one than in two. The run-off coefficients (cumulative simulated run-off divided by the cumulative precipitation over the same period) varied from 4% to 13% in 2019/20 year and from 7% to 24% in 2020/21. The best-performing ANN had six input features, four hidden layers with three to four nodes per layer (Fig. 7-A and -B). The best-performing features were x1, x3, x5, x6, x10, x11 (Table 4). According to the PCA euclidean scores, these features were the most correlated with LIF, significantly contributing to the first and second components, which explained together 52.2% of LIF variance (Appendix-A3). The network was solved in 9962 steps (Fig. 7-D). The best performing ANN had a R 2 of 0.97 and 0.78, respectively for the training and the testing set, and the RMSE of predicted CUM. LIF values were 10.1 and 14.8 mm, respectively (Fig. 7-C). Observed CUM. LIF were scatter-plotted against simulated values and no biased error trend was observed for any of the subsets, as the residuals did not vary with the level of predicted LIF. Despite showing a solid forecasting capacity, the ANN exhibits a general trend to overpredict LIF in comparison to measured values (Fig. 7-C). Daily LIF calendars were predicted with the ANN (Fig. 7) over a 30year period (Fig. 8). While seasonal precipitation varied over 30 years from 212.8 to 759.5 mm, cumulative LIF ranged from 30 to 125 mm (Fig. 8). The ratio of seasonal cumulative LIF divided by seasonal precipitation varied from 10.7% to 38.9% over the 30 years. 3.3. Yield simulations with AquaCrop The AquaCrop simulation outcomes in terms of yield response to LIF were highly variable from year to year. Simulated plotted against observed yields [2015–20] showed RMSE, R 2 and Willmott d index respectively equal to 0.374 (Mg GY ha −1 ), 0.35, 0.76 (Appendix-A4). Differences among zones were only experimentally tested for the two systematically distinct zones (i.e., A and C according to Table 5). In this sense, the cumulative probability of NYR LIF for zone C can be also interpreted as the absolute difference between the two curves shown in Fig. 9-A. Since cumulative LIF in zone A was negligible (Appendix-A6), the NYR LIF was only estimated for the water-receiving zone C (Fig. 9-B). In this case, NYR LIF corresponds to Y 2 minus Y 1 , in zone C, or Y 2 in zone C minus Y 2 in zone A (Fig. 9-A). Fig. 9 shows the cumulative probability distribution curves for simulated yield responses in the two significantly distinct zones (A and C). Results are shown both in absolute and relative terms (Fig. 9-B and -D). Mean values of simulated NYR LIF and NYR LIF−rel were 383 kg ha −1 (Fig. 9-C) and 16.2%, respectively. Absolute values are expressed in terms of grain yield (GY) in dry mass (DM). Following crop yield simulations over a period of 30 years, the simulated seasonal NYR LIF varied from −0.14 to 2.08 Mg ha −1 , corresponding to a relative net contribution of LIF to grain yield (NYR LIF−rel ) that ranged from −3% up to +168%. NYR LIF was larger than 800 kg GY ha −1 in five out of 30 years (i.e., 1992/93, 1999/00, 2003/04, 2004/05, 2008/09). By contrast, in other five years (i.e., 1990/91, 1996/97, 2001/02, 2007/08 and 2019/20), LIF caused yield losses due to water excess (Fig. 9-C and Appendix-Table B). For the remaining years, the NYR LIF ranged from none to 670 kg GY ha −1 (e.g., 2014/15 as shown in Fig. 9-C), being below 265 kg GY ha −1 for at least 50% of the years (Fig. 9-B). According to the Spearman’s correlation matrix, shown in Appendix (Fig. A7), NYR LIF was negatively correlated with CUM. P and simulated yields. The lower the yields simulated in higher zones (i.e., zone A), the higher was NYR LIF at lower zones (i.e., zone C), which is caused by marginal water productivity gains (Appendix-A7 and -Table B). Fig. 6. HYDRUS-1D simulations performance according to the means of measured soil water content and the best fitted soil parameters. “n.a.” indicates non-applicable situations, which correspond to the sampling point A3 and C3 that were not considered in catchment one (Fig. 2). Both the Nash-Sutcliffe and the R 2 coefficients are indicated by the left axis, while the RMSE is indicated by the right axis. Mean R 2 values ranged from 0.65 ( ±0.13) in 2019/20 to 0.93 (±0.03) in 2020/21 and RMSE (cm 3 cm −3 ) from 0.04 ( ±0.01) to 0.03 (±0.01), by the same order. T.R. Tenreiro et al. European Journal of Agronomy 137 (2022) 126515 16 Fig. A2. The mean correction factor (Ω z ) plotted for each probe location and sensor depth in 2019/20, at the 95% confidence level. Ω z was computed as the mean ratio between the SWC measurements with the NP and the capacitance probes for depth z (Equation 2). Fig. A3. The principal component analysis (PCA) plot. Vectors are colored according to the contribution degree to LIF. PCA scores are equal to the module of each vector, indicating the weight associated with the combination of the two principal components (i.e., Dim 1 and Dim2, respectively explaining 34% and 18.2% of LIF variance). Variables are defined in Table 4. T.R. Tenreiro et al. European Journal of Agronomy 137 (2022) 126515 17 Fig. A4. Simulated vs. observed yields (obtained from historical yield maps). Units are expressed in Mg GY ha −1 . Circles and triangles represent yields in zone A and C, respectively. Simulated yields correspond to the Y 2 scenario. RMSE, R 2 and Willmott d index are respectively equal to 0.374 (Mg GY ha −1 ), 0.35, 0.76. The horizontal bars indicate the error associated with the process of yield mapping (i.e., 172–809 kg GY ha −1 ). Fig. A5. Cumulative precipitation (CUM. P) and HYDRUS run-off simulations over 30 years [1990–2020]. CUM. P is expressed in mm. Blue bars represent daily surface run-off values, expressed in mm day −1 . Fig. A6. Cumulative probability distribution curves of CUM. LIF and CUM. P over 30 years. Curves shown for zone A (A) and zone C (B), representing both cumulative precipitation and cumulative LIF, values expressed in mm day −1 . T.R. Tenreiro et al. European Journal of Agronomy 137 (2022) 126515 18 Fig. A7. Spearman’s correlation pairs-plot. Codes: PAW 0 =initial plant available water (at sowing date); CUM. P=season cumulative precipitation; CUM. LIF =season cumulative lateral inflow (LIF); CUM. ETo =season cumulative evapotranspiration (ETo); PF. LIF =postflowering LIF, corresponding to the fraction of CUM. LIF taking place at post-flowering stages; ’Relative.T′=mean relative crop transpiration (the simple average of the relative crop transpiration between zone A and C, expressed in %); Mean YP =mean yield potential (the simple average of simulated yield between zone A and C, expressed in Mg GY ha −1 ); LIF. MWP =LIF marginal water productivity (expressed in kg GY ha −1 mm −1 ); NYR. LIF =Net yield response to LIF (expressed in Mg GY ha −1 ); PAW 0 , CUM. P, CUM. LIF and CUM. ET0 are expressed in mm. ’Relative.T′is estimated at the season average of daily crop actual transpiration divided by daily potential transpiration. Simulation files are provided in Supplementary materials. Season is defined from sowing to harvesting date. Input values are synthesized in Appendix-Table B. Significant correlations at the 5% level of significance are highlighted with the symbol ’*’. Significance codes: ’***’ 0.1%, ’**’ 1%, ’*’ 5%. Table A1 The parameters of the neutron probe (NP) calibration functions (y=a[x std] − b). y is volumetric SWC, expressed in %, x is probe measured value (unitless) and ’std’ is the standard correction value. More information on the neutron probe calibration site is provided by Soriano et al. (2018). Depth a b std R 2 0–15 cm 25.416 0.264 7830 0.82 15–30 cm 21.974 7.381 7830 0.98 30–90 cm 27.210 16.897 7434 0.96 T.R. Tenreiro et al. European Journal of Agronomy 137 (2022) 126515 19 Appendix A. Supporting information Supplementary data associated with this article can be found in the online version at doi:10.1016/j.eja.2022.126515. References Abbate, P.E., Dardanelli, J.L., Cantarero, M.G., Maturano, M., Melchiori, R.J.M., Suero, E.E., 2004. Climatic and water availability effects on water-use efficiency in wheat. Crop Sci. 44, 474–483. Acutis, M., Scaglia, B., Confalonieri, R., 2012. Perfunctory analysis of variance in agronomy, and its consequences in experimental results interpretation. Eur. J. Agron. 43, 129–135. Ahuja, L.R., Ma, L., Anapalli, S.S., 2019. Biophysical system models advance agricultural research and technology: some examples and further research needs. In: Bridging Among Disciplines by Synthesizing Soil and Plant Processes, Advances in Agricultural Systems Modeling. American Society of Agronomy, Crop Science Society of America, and Soil Science Society of America, Inc., Madison, WI. Allen Jr., L.H., Kakani, V.G., Vu, J.C.V., Boote, K.J., 2011. Elevated CO 2 increases water use efficiency by sustaining photosynthesis of water-limited maize and sorghum. J. Plant Physiol. 168, 1909–1918. Allen, R.G., Pereira, L.S., Raes, D., Smith, M., et al., 1998. Crop evapotranspiration - guidelines for computing crop water requirements. In: Fao Irrigation and Drainage Paper 56, 300. Fao,, Rome, p. D05109. Balafoutis, A.T., Beck, B., Fountas, S., Tsiropoulos, Z., Vangeyte, J., van der Wal, T., SotoEmbodas, I., G´ omez-Barbero, M., Pedersen, S.M., 2017. Smart farming technologies - description, taxonomy and economic impact. In: Progress in Precision Agriculture. Springer International Publishing, Cham, pp. 21–77. Batchelor, W.D., Basso, B., Paz, J.O., 2002. Examples of strategies to analyze spatial and temporal yield variability using crop models. Eur. J. Agron. 18, 141–158. Bramley, R.G.V., 2009. Lessons from nearly 20 years of Precision Agriculture research, development, and adoption as a guide to its appropriate application. Crop Pasture Sci. Campbell, J.E., 1990. Dielectric properties and influence of conductivity in soils at one to fifty megahertz. Soil Sci. Soc. Am. J. 54, 332–341. Cassman, K.G., 1999. Ecological intensification of cereal production systems: yield potential, soil quality, and precision agriculture. Proc. Natl. Acad. Sci. USA 96, 5952–5959. Ciha, A.J., 1984. Slope position and grain yield of soft white winter wheat. Agron. J. 76, 193–196. Chanzy, A., Chadoeuf, J., Gaudu, J.C., Mohrath, D., Richard, G., Bruckler, L., 1998. Soil moisture monitoring at the field scale using automatic capacitance probes. Eur. J. Soil Sci. 49, 637–648. CNIG., 2019, Centro Nacional de Informaci´ on Geogr´ afica (CNIG).〈http://centrodedescar gas.cnig.es/CentroDescargas/index.jsp〉(accessed 8.20.19). Connor, D.J., Mínguez, M.I., 2012. Evolution not revolution of farming systems will best feed and green the world. Glob. Food Secur. 1, 106–113. Conrad, O., Bechtel, B., Bock, M., Dietrich, H., Fischer, E., Gerlitz, L., B¨ ohner, J., 2015. System for automated geoscientific analyses (SAGA) v. 2.1. 4. Geosci. Model Dev. Discuss. 8, 2. De Veaux, R.D., Ungar, L.H., 1994. Multicollinearity: a tale of two nonparametric regressions. In: Selecting Models from Data. Springer, New York, pp. 393–402. de Vos, N.J., Rientjes, T.H.M., 2005. Constraints of artificial neural networks for rainfallrunoff modelling: trade-offs in hydrological state representation and model evaluation. Hydrol. Earth Syst. Sci. https://doi.org/10.5194/hess-9-111-2005. de Wit, C.T., van Keulen, H., 1987. Modelling production of field crops and its requirements. Geoderma 40, 253–265. Evett, S.R., Schwartz, R.C., Tolk, J.A., Howell, T.A., 2009. Soil profile water content determination: spatiotemporal variability of electromagnetic and neutron probe sensors in access tubes. Vadose Zone J. 8, 926–941. Feddes, R.A., Kowalik, P.J., Zaradny, H., 1978.Water uptake by plant roots.Simulation of field water use and crop yield 16–30. Fischer, R.A., Connor, D.J., 2018. Issues for cropping and agricultural science in the next 20 years. Field Crops Res. 222, 121–142. Fischer, R.A., MorenoRamos, O.H., OrtizMonasterio, I., Sayre, K.D., 2019. Yield response to plant density, row spacing and raised beds in low latitude spring wheat with ample soil resources: an update. Field Crops Res. 232, 95–105. Florin, M.J., McBratney, A.B., Whelan, B.M., 2009. Quantification and comparison of wheat yield variation across space and time. Eur. J. Agron. 30, 212–219. Franz, T.E., Pokal, S., Gibson, J.P., Zhou, Y., Gholizadeh, H., Tenorio, F.A., Rudnick, D., Heeren, D., McCabe, M., Ziliani, M., Jin, Z., Guan, K., Pan, M., Gates, J., Wardlow, B., 2020. The role of topography, soil, and remotely sensed vegetation condition towards predicting crop yield. Field Crops Res. 252, 107788. Table B1 AquaCrop simulation outcomes. Codes: PAW 0 =Initial plant available water (at sowing date); CUM. ETo =Season cumulative ETo (from sowing to harvesting date); CUM. P=Season cumulative precipitation; CUM. LIF =Season cumulative lateral inflow (LIF); Post-flowering LIF =The fraction of season LIF taking place at postflowering stages (expressed in %); Relative T [Zone A] =Mean relative crop transpiration in zone A (or in the absence of LIF); Relative T [Zone C] =Mean relative crop transpiration in zone C (including LIF); Yield.1 [Zone A] =Yield in scenario 1 (in the absence of LIF); Yield.2 [Zone C] =Yield in scenario 2 (including LIF); NYR. LIF =Net yield response to LIF; NYR. LIF rel =Net yield response to LIF in relative terms; LIF. MWP =LIF marginal water productivity (expressed in kg GY ha −1 mm −1 ). Relative T is estimated as the season average of daily crop actual transpiration divided by potential transpiration. Simulation files are provided in Supplementary materials. Season PAW 0 CUM. ETo CUM. P CUM. LIF Post-flowering LIF Relative T [Zone A] Relative T [Zone C] Yield.1 [Zone A] Yield.2 [Zone C] NYR. LIF NYR. LIF rel LIF. MWP - [mm] [mm] [mm] [mm] [%] [%] [%] [Mg GY ha −1 ] [Mg GY ha −1 ] [Mg GY ha −1 ] [%] [kg GY ha −1 mm −1 ] 1990/91 157.5 434.4 419.1 46 0 84.2 83.4 3.68 3.58 -0.10 -3 – 1991/92 145.5 351.8 361.4 118 62 93.4 98.2 4.25 4.49 0.24 6 2.03 1992/93 135.4 356.8 224.3 48 25 87.0 90.0 1.24 3.32 2.08 168 43.33 1993/94 157.3 370.8 212.8 34 26 82.6 83.8 3.54 3.66 0.12 3 3.53 1994/95 151.2 391.8 223.8 37 32 78.3 83.4 3.77 3.87 0.10 3 2.70 1995/96 139.7 391.3 747.0 120 56 87.0 88.4 4.14 4.28 0.14 3 1.17 1996/97 162.4 387.5 759.5 122 0 72.5 70.4 3.36 3.34 -0.02 -1 – 1997/98 165.0 354.8 689.6 125 3 83.0 83.1 4.15 4.18 0.03 1 0.24 1998/99 117.9 389.6 301.6 34 41 88.6 90.1 2.56 2.87 0.31 12 9.12 1999/00 118.2 452.1 326.2 47 57 89.5 92.9 3.56 4.45 0.89 25 18.94 2000/01 124.0 350.7 523.9 115 48 85.7 89.6 4.22 4.51 0.29 7 2.52 2001/02 116.4 392.0 398.2 66 23 87.1 87.0 4.44 4.43 -0.01 0 – 2002/03 160.1 401.3 376.2 108 56 77.7 78.1 3.68 3.76 0.08 2 0.74 2003/04 121.6 352.1 310.0 120 0 91.1 97.4 3.83 4.81 0.98 26 8.17 2004/05 120.6 439.8 325.5 45 11 85.5 89.0 1.99 3.38 1.39 70 30.89 2005/06 154.3 394.9 286.4 44 20 91.6 95.2 4.37 4.58 0.21 5 4.77 2006/07 134.6 389.8 407.8 44 16 98.4 99.4 4.93 5.00 0.07 1 1.59 2007/08 133.1 424.9 359.8 85 40 94.9 93.2 4.89 4.75 -0.14 -3 – 2008/09 122.0 406.1 257.7 48 8 86.1 88.3 1.83 2.86 1.03 56 21.46 2009/10 145.9 397.1 215.8 36 19 82.1 88.0 3.86 4.28 0.42 11 11.67 2010/11 157.9 425.6 300.0 46 30 87.4 90.9 4.24 4.69 0.45 11 9.78 2011/12 182.3 446.6 266.9 39 23 88.9 92.1 4.01 4.44 0.43 11 11.03 2012/13 170.4 397.8 273.8 38 21 89.6 92.9 4.46 4.83 0.37 8 9.74 2013/14 172.5 453.6 299.8 44 18 83.5 87.0 3.89 4.38 0.49 13 11.14 2014/15 159.8 445.7 228.6 41 10 90.2 92.3 2.28 2.95 0.67 29 16.34 2015/16 173.2 396.7 228.9 30 10 85.8 88.1 4.10 4.51 0.41 10 13.67 2016/17 168.3 421.4 375.4 104 20 89.8 91.2 4.27 4.61 0.34 8 3.27 2017/18 176.9 416.0 297.1 46 0 84.0 85.6 4.97 5.10 0.13 3 2.83 2018/19 186.9 458.4 412.3 46 0 83.7 85.8 5.08 5.31 0.23 5 5.00 2019/20 161.8 453.2 535.2 112 0 88.6 86.5 5.35 5.21 -0.14 -3 – T.R. Tenreiro et al. European Journal of Agronomy 137 (2022) 126515 20 French, R.J., Schultz, J.E., 1984. Water use efficiency of wheat in a Mediterranean-type environment. I. The relation between yield, water use and climate. Aust 35, 743–764. García-Ruiz, J.M., 2010. The effects of land uses on soil erosion in Spain: a review. Catena 81, 1–11. Giraldez, J.V., Sposito, G., 1985. Infiltration in swelling soils. Water Resour. Res. 21, 33–44. Günther, F., Fritsch, S., 2010.Neuralnet: Training of neural networks.R J.2, 30. Halvorson, G.A., Doll, E.C., 1991. Topographic effects on spring wheat yields and water use. Soil Sci. Soc. Am. J. 55, 1680–1685. Hochman, Z., van Rees, H., Carberry, P.S., Hunt, J.R., McCown, R.L., Gartmann, A., Holzworth, D., van Rees, S., Dalgliesh, N.P., Long, W., Peake, A.S., Poulton, P.L., McClelland, T., 2009. Re-inventing model-based decision support with Australian dryland farmers. 4. Yield Prophet® helps farmers monitor and manage crops in a variable climate. Crop Pasture Sci. 60, 1057–1070. Hoffman, G.J., van Genuchten, M.T., 1983. Soil properties and efficient water use: water management for salinity control. Limitations to Efficient Water Use in Crop Production, pp. 73–85. Hsiao, T.C., 1993. Effects of drought and elevated CO2 on plant water use efficiency and productivity. Interacting Stresses on Plants in a Changing Climate. Springer, Berlin Heidelberg, pp. 435–465. Jenson, S.K., Domingue, J.O., 1988. Extracting topographic structure from digital elevation data for geographic information system analysis. Photogramm. Eng. Remote Sens. 54, 1593–1600. Jones, J.W., Antle, J.M., Basso, B., Boote, K.J., Conant, R.T., Foster, I., Godfray, H.C.J., Herrero, M., Howitt, R.E., Janssen, S., Keating, B.A., Munoz-Carpena, R., Porter, C. H., Rosenzweig, C., Wheeler, T.R., 2017. Toward a new generation of agricultural system data, models, and knowledge products: State of agricultural systems science. Agric. Syst. 155, 269–288. Kamilaris, A., Kartakoullis, A., Prenafeta-Boldú, F.X., 2017. A review on the practice of big data analysis in agriculture. Comput. Electron. Agric. 143, 23–37. Kempenaar, C., Lokhorst, C., Bleumer, E.J.B., Veerkamp, R.F., 2016.Big Data analysis for smart farming: results of TO2 project in theme food security. Kirkegaard, J.A., Hunt, J.R., 2010. Increasing productivity by matching farming system management and genotype in water-limited environments. J. Exp. Bot. 61 (15), 4129–4143. Kirkegaard, J.A., Lilley, J.M., Howe, G.N., Graham, J.M., 2007. Impact of subsoil water use on wheat yield. Aust 58, 303–315. Klaij, M.C., Vachaud, G., 1992. Seasonal water balance of a sandy soil in Niger cropped with pearl millet, based on profile moisture measurements. Agric. Water Manag. 21, 313–330. Kravchenko, A.N., Bullock, D.G., 2000. Correlation of corn and soybean grain yield with topography and soil properties. Agron. J. 92, 75–83. Lasanta, T., García-Ruiz, J.M., P´ erez-Rontom´ e, C., Sancho-Marc´ en, C., 2000. Runoff and sediment yield in a semi-arid environment: the effect of land management after farmland abandonment. Catena 38, 265–278. Lovelace, R., Nowosad, J., Muenchow, J., 2019. Geocomputation with R. CRC Press. Maestrini, B., Basso, B., 2018. Drivers of within-field spatial and temporal variability of crop yield across the US Midwest. Sci. Rep. 8, 14833. Maier, H.R., Dandy, G.C., 2000. Neural networks for the prediction and forecasting of water resources variables: a review of modelling issues and applications. Environ. Model. Softw. 15, 101–124. Maina, F.Z., Siirila-Woodburn, E.R., 2020. The role of subsurface flow on evapotranspiration: a global sensitivity analysis. Water Resour. Res 56. https://doi. org/10.1029/2019wr026612. McBratney, A., Whelan, B., Ancev, T., Bouma, J., 2005. Future directions of precision agriculture. Precis. Agric. 6, 7–23. McCutcheon, M.C., Farahani, H.J., Stednick, J.D., Buchleiter, G.W., Green, T.R., 2006. Effect of soil water on apparent soil electrical conductivity and texture relationships in a dryland field. Biosyst. Eng. 94, 19–32. Monzon, J.P., Calvi˜ no, P.A., Sadras, V.O., Zubiaurre, J.B., Andrade, F.H., 2018. Precision agriculture based on crop physiological principles improves whole-farm yield and profit: a case study. Eur. J. Agron. 99, 62–71. Moriasi, D.N., Arnold, J.G., Van Liew, M.W., Bingner, R.L., Harmel, R.D., Veith, T.L., 2007. Model evaluation guidelines for systematic quantification of accuracy in watershed simulations. Trans. ASABE 50, 885–900. Mualem, Y., 1976. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 12 (3), 513–522. Mwale, S.S., Azam-Ali, S.N., Sparkes, D.L., 2005. Can the PR1 capacitance probe replace the neutron probe for routine soil-water measurement? Soil Use Manag. 21 (3), 340–347. Nash, J.E., Sutcliffe, J.V., 1970. River flow forecasting through conceptual models part I - A discussion of principles. J. Hydrol. 10, 282–290. Nielsen, D., Biggar, J.W., Erh, K.T., 1973. Spatial variability of field-measured soil-water properties. Hilgardia 42, 215–259. Nielsen, D.R., Wendroth, O., 2003. Spatial and Temporal Statistics. Schweiz. Verl. Paltineanu, I.C., Starr, J.L., 1997. Real-time soil water dynamics using multisensor capacitance probes: Laboratory calibration. Soil Sci. Soc. Am. J. 61, 1576–1585. Passioura, J.B., Angus, J.F., 2010. Chapter 2Improving productivity of crops in waterlimited environments. In: Sparks, D.L. (Ed.), Advances in Agronomy. Academic Press, pp. 37–75. Patrignani, A., Ochsner, T.E., 2015. Canopeo: a powerful new tool for measuring fractional green canopy cover. Agron. J. 107, 2312–2320. Pearl, J., 2019. The seven tools of causal inference, with reflections on machine learning. Commun. ACM 62, 54–60. RattalinoEdreira, J.I., Guilpart, N., Sadras, V., Cassman, K.G., van Ittersum, M.K., Schils, R.L.M., Grassini, P., 2018. Water productivity of rainfed maize and wheat: a local to global perspective. Agric. . Meteor. 259, 364–373. Reitz, P., Kutzbach, H.D., 1996. Investigations on a particular yield mapping system for combine harvesters. Comput. Electron. Agric. 14, 137–150. Roadknight, C.M., Balls, G.R., Mills, G.E., Palmer-Brown, D., 1997. Modeling complex environmental data. IEEE Trans. Neural Netw. 8, 852–862. Rockstr¨ om, J., Valentin, C., 1997. Hillslope dynamics of on-farm generation of surface water flows: the case of rain-fed cultivation of pearl millet on sandy soil in the Sahel. Agric. Water Manag. RoTimi Ojo, E., Bullock, P.R., Fitzmaurice, J., 2015. Field performance of five soil moisture instruments in heavy clay soils. Soil Sci. Soc. Am. J. 79 (1), 20–29. Sadler, E.J., Russell, G., 1997. Modeling crop yield for site-specific management. In: Pierce, F.J., Sadler, E.J. (Eds.), The State of Site-Specific Management for Agriculture. Sadras, V., Bongiovanni, R., 2004. Use of Lorenz curves and Gini coefficients to assess yield inequality within paddocks. Field Crops Res. 90, 303–310. Sadras, V.O., Angus, J.F., 2006. Benchmarking water-use efficiency of rainfed wheat in dry environments. Aust 57, 847–856. Sadras, V.O., Reynolds, M.P., de la Vega, A.J., Petrie, P.R., Robinson, R., 2009. Phenotypic plasticity of yield and phenology in wheat, sunflower and grapevine. Field Crops Res. 110, 242–250. Sadras, V.O., Villalobos, F.J., Orgaz, F., Fereres, E., 2016. Effects of water stress on crop production. In: Villalobos, F.J., Fereres, E. (Eds.), Principles of Agronomy for Sustainable Agriculture. Springer International Publishing, Cham, pp. 189–204. Schaap, M.G., Leij, F.J., van Genuchten, M.T., 2001. rosetta: a computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions. J. Hydrol. Schmitter, P., Zwart, S.J., Danvi, A., Gbaguidi, F., 2015. Contributions of lateral flow and groundwater to the spatio-temporal variation of irrigated rice yields and water productivity in a West-African inland valley. Agric. Water Manag. 152, 286–298. Senthil Kumar, A.R., Sudheer, K.P., Jain, S.K., Agarwal, P.K., 2005. Rainfall-runoff modelling using artificial neural networks: comparison of network types. Hydrol. Process. Int. J. 19, 1277–1291. Shukla, M.B., Kok, R., Prasher, S.O., Clark, G., Lacroix, R., 1996. Use of artificial neural networks in transient drainage design. Trans. ASAE 39, 119–124. Sida, T.S., Chamberlin, J., Ayalew, H., Kosmowski, F., Craufurd, P., 2021. Implications of intra-plot heterogeneity for yield estimation accuracy: evidence from smallholder maize systems in Ethiopia. Field Crops Res. 267, 108147. Silva, J.V., Tenreiro, T.R., Sp¨ atjens, L., Anten, N.P.R., van Ittersum, M.K., Reidsma, P., 2020. Can big data explain yield variability and water productivity in intensive cropping systems? Field Crops Res. 255, 107828. ˇ Simunek, J., ˇ Sejna, M., Saito, M., Sakai, M., van Genuchten, M.T., 2018. HYDRUS-1D Version 4.17 - Manual. Department of Environmental Science, University of California, Riverside. ˇ Simunek, J., Van Genuchten, M.T., ˇ Sejna, M., 2012. HYDRUS: Model use, calibration, and validation. Trans. ASABE 55 (4), 1263–1274. Smith, M.J., 2020. Getting value from artificial intelligence in agriculture. Anim. Prod. Sci. 60, 46–54. Soil Survey Staff, 1999. A basic system of soil classification for making and interpreting soil surveys. Soil Taxonomy. 2nd ed., USDA Agr.Hbk.436, WA. Soriano, M.A., Cabezas, J.M., G´ omez, J.A., 2018. Soil water content and yield a vertisol in a rain-fed olive grove under four different soil management practices in a four year experiment. EGU General Assembly Conference Abstracts, p. 5390. Spiertz, H., 2014. Agricultural sciences in transition from 1800 to 2020: exploring knowledge and creating impact. Eur. J. Agron. 59, 96–106. Steduto, P., Hsiao, T.C., Fereres, E., 2007. On the conservative behavior of biomass water productivity. Irrig. Sci. 25, 189–207. Steduto, P., Hsiao, T.C., Fereres, E., Raes, D., 2012. Crop Yield Response to Water (Vol. 1028). Food and Agriculture Organization of the United Nations, Rome. Steduto, P., Hsiao, T.C., Raes, D., Fereres, E., 2009. AquaCrop–The FAO Crop model to simulate yield response to water: i. concepts and underlying principles. Agron. J. 101, 426–437. Tarboton, D.G., Bras, R.L., Rodriguez-Iturbe, I., 1991. On the extraction of channel networks from digital elevation data. Hydrol. Process. 5, 81–100. Tenreiro, T.R., García-Vila, M., G´ omez, J.A., Jimenez-Berni, J.A., Fereres, E., 2021. Using NDVI for the assessment of canopy cover in agricultural crops within modelling research. Comput. Electron. Agric. 182, 106038. Tenreiro, T.R., García-Vila, M., G´ omez, J.A., Jimenez-Berni, J.A., Fereres, E., 2020. Water modelling approaches and opportunities to simulate spatial water variations at crop field level. Agric. Water Manag. Tolk, J.A., Howell, T.A., Steiner, J.L., Krieg, D.R., Schneider, A.D., 1995. Role of transpiration suppression by evaporation of intercepted water in improving irrigation efficiency. Irrig. Sci. 16. Toreti, A., Maiorano, A., De Sanctis, G., Webber, H., Ruane, A.C., Fumagalli, D., Ceglar, A., Niemeyer, S., Zampieri, M., 2018. Using reanalysis in crop monitoring and forecasting systems. Agric. Syst. https://doi.org/10.1016/j.agsy.2018.07.001. Torralba, M.A., 2013. Evaluaci´ on de la erosion hídrica en parcelas experimentales en campos agrícolas de secano mediterraneo (Doctoral dissertation, Universidad Complutense de Madrid). van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44, 892–898. Verhagen, A., Booltink, H.W.G., Bouma, J., 1995. Site-specific management: balancing production and environmental requirements at farm level. Agric. Syst. 49, 369–384. Wallor, E., Kersebaum, K.-C., Ventrella, D., Bindi, M., Cammarano, D., Coucheney, E., Gaiser, T., Garofalo, P., Giglio, L., Giola, P., Hoffmann, M.P., Iocola, I., Lana, M., T.R. Tenreiro et al. European Journal of Agronomy 137 (2022) 126515 21 Lewan, E., Maharjan, G.R., Moriondo, M., Mula, L., Nendel, C., Pohankova, E., Roggero, P.P., Trnka, M., Trombi, G., 2018. The response of process-based agroecosystem models to within-field variability in site conditions. Field Crops Res. 228, 1–19. Wang, L., Liu, H., 2006. An efficient method for identifying and filling surface depressions in digital elevation models for hydrologic analysis and modelling. Int. J. Geogr. Inf. Sci. 20, 193–213. Ward, N.K., Maureira, F., St¨ ockle, C.O., Brooks, E.S., Painter, K.M., Yourek, M.A., Gasch, C.K., 2018. Simulating field-scale variability and precision management with a 3D hydrologic cropping systems model. Precis. Agric. 19, 293–313. Wesseling, J.G., Elbers, J.A., Kabat, P., Van Den Broek, B.J., 1991. SWATRE: Instructions for Input. Internal Note, Winand Staring Centre, Wageningen, The Netherlands. International Waterlogging and Salinity Research Institute, Lahore, Pakistan, p. 29. Willmott, C.J., 1981. On the validation of models. Phys. Geogr. 2, 184–194. Wolfert, S., Ge, L., Verdouw, C., Bogaardt, M.-J., 2017. Big data in smart farming - a review. Agric. Syst. 153, 69–80. Yang, J.M., Yang, J.Y., Liu, S., Hoogenboom, G., 2014. An evaluation of the statistical methods for testing the performance of crop models with observed data. Agric. Syst. 127, 81–89. Zumr, D., Dohnal, M., Hrnˇ cír^, M., Císlerov´ a, M., Vogel, T., Doleˇ zal, F., 2006. Simulation of soil water dynamics in structured heavy soils with respect to root water uptake. Biologia 61, S320–S323. T.R. Tenreiro et al.