AquaCrop Simulation of Winter Wheat under Different N Management Practices

AquaCrop is a well-known water-oriented crop model. The model has been often used to simulate various crops and the water balance in the field under different irrigation treatments, but studies that relate AquaCrop to fertilization are rare. In this study, the ability of this model to simulate yield and the water balance parameters was investigated in a wheat field under different nitrogen management practices. During the 2015–2016 and 2016–2017 growing seasons, meteorological data were provided from a nearby meteorological station, and the evolution of soil water content and final yields were recorded. The model showed a very good performance at simulating the soil water content evolution in the root zone. Notwithstanding its simplicity, AquaCrop based on a semi-quantitative approach for fertility performed well at the field level for the final yield estimation under different nitrogen treatments and field topography variation. Although the correlation coefficient between simulated and measured final yields was high, increased values of variations were observed in the various zones of this experimental field (−50% to +140%). The model appears to be an efficient tool for evaluating and improving the management practices at the field level. The experiments were conducted in Thessaly, which is the largest plain and the main agricultural area of Greece. Thessaly, however, has a strong negative water balance, which has led to a strong decrease in the level of the aquifer and, at the same time, to sea intrusion. There is also a significant risk of contamination of the groundwater aquifer due to increased use of agrochemicals. This analysis is particularly important for Thessaly due to the need for improvement of agricultural practices in this area, to decrease the pressure of agricultural activities on natural resources (soil, water) and reverse the consequences of current management.


Introduction
The winter and spring rainy seasons are the periods of wheat cultivation in Greece [1] with the possible water surplus drained or running off the farms. Some farmers apply supplementary irrigation in the spring that is usually not scheduled properly [2]. To reduce water losses and nitrogen leaching, a crop simulation model is needed that is able to combine the available field data, i.e., temperature, rain/irrigation, soil texture, and dry yield. Based on the concept of a direct link between crop water use and crop yield [3], the the in-season VRA approach significantly reduced N inputs without yield losses in winter wheat, thus offering higher revenue over fertilizer costs with the least uncertainty to the farmer [40]. However, VRA systems are effective when N does not escape from the rooting zone after fertilizer has been applied. Furthermore, N leaching raises environmental concerns in areas of intensive cultivation. European Union member countries are obliged to comply with the Water Directive (2000/60/CE) [48] that requires countermeasures for N water pollution. For these reasons (economic and environmental), synchronization of fertilization with the water cycle variability is the suitable strategy [49].
The objectives of this study were: (i) evaluation of the effectiveness of AquaCrop to estimate the soil water content evolution and water balance in the field and (ii) calibration and validation of the model to simulate winter wheat yield under various N treatments.

Field Experiments and Data
The field experiments were carried out in a winter wheat field near the village of Nea Lefki (Larisa, central Greece) at the coordinates of 39 • 31 30 N, 22 • 30 52 E and an altitude of 125 m [50,51]. The temperate Mediterranean climate has cold humid winters but hot and arid summers (Bsk-Csa according to Köppen's climatic classification) [1]. Based on historical data, the mean minimum and maximum temperatures are 0.5 • C (January) and 33 • C (July), respectively, with an average relative humidity of 67%. The mean monthly rainfall varies from 15.5 mm (August) to 58.2 mm (November) with a mean annual total of 313 mm. During the November to May growing season, the mean total rainfall is 277 mm, and the average monthly temperature varies from 5 • C in January to 17 • C in May [51]. The rainfall was 255.0 mm and 283.8 mm during 2015-2016 and 2016-2017 growing seasons, respectively.
A soil sampling was performed with regard to field topography, at lower, middle, and upper zones of the field. A composite soil sample (0-20 and 20-40 cm depth, n = 3) was collected from 1 m 2 sub-plots in three landscape positions of each treatment strip (lower zone, middle zone, and upland zone). Soil samples were analyzed in the laboratory for soil texture as determined by physical fractionation (Bouyoukos 1951) [52], and for soil organic matter as determined by the Walkley-Black method of wet oxidation [53]. Carbonate content, as an estimate of inorganic C, was determined by using a Bernard calcimeter to measure the released carbon dioxide (CO 2 ) after addition of dilute hydrochloric acid (HCl) solution [53]. Total N was measured by the semimicro-Kjeldahl method [54]. Soil K was extracted with 1 N ammonium acetate at pH 7 [55] and measured by a Corning 410 flame photometer. Available P was determined as per Olsen and Sommers (1982) [56].
The soil in the field is classified as a Vertic Xerochrept of clay-to-clay loam texture, but with differences in the soil properties relevant to topography (Table 1).
Strong effects of soil erosion in the slopes and upland areas of the field have been reported by Stamatiadis et al. (2018) [50]. Those areas of the field had a higher carbonate and lower organic matter and total N content than the lowland plots as a result of the removal of surface horizons that revealed the parent material rich in calcium carbonate.
The topographic variation in the field is shown in Figures 1 and 2. The lower ground level on the slope by 1-2 m, relative to the adjacent field on the right (Figure 1, see the red arrow), appears to be the product of soil erosion. In contrast to the adjacent field, winter wheat was grown without rotation over recent decades. The soil properties in different sections of the field were analyzed in the laboratory.
The grown wheat cultivar was Simeto, an Italian variety that is popular among Greek farmers. This cultivar was selected based on its high yield, climatic adaptability, and performance stability.  Strong effects of soil erosion in the slopes and upland areas of the field have been reported by Stamatiadis et al., (2018) [50]. Those areas of the field had a higher carbonate and lower organic matter and total N content than the lowland plots as a result of the removal of surface horizons that revealed the parent material rich in calcium carbonate.
The topographic variation in the field is shown in Figures 1 and 2. The lower ground level on the slope by 1-2 m, relative to the adjacent field on the right (Figure 1, see the red arrow), appears to be the product of soil erosion. In contrast to the adjacent field, winter wheat was grown without rotation over recent decades. The soil properties in different sections of the field were analyzed in the laboratory.    Strong effects of soil erosion in the slopes and upland areas of the field have been reported by Stamatiadis et al., (2018) [50]. Those areas of the field had a higher carbonate and lower organic matter and total N content than the lowland plots as a result of the removal of surface horizons that revealed the parent material rich in calcium carbonate.
The topographic variation in the field is shown in Figures 1 and 2. The lower ground level on the slope by 1-2 m, relative to the adjacent field on the right (Figure 1, see the red arrow), appears to be the product of soil erosion. In contrast to the adjacent field, winter wheat was grown without rotation over recent decades. The soil properties in different sections of the field were analyzed in the laboratory.    The field was divided into four blocks, and three N treatments were randomly assigned within each block to follow a randomized complete block design ( Figure 3). Each treatment was a field strip of 7 m wide running the entire length of the field (200 m) to simulate fullscale field conditions. The three treatments during the 2015-2016 growing season consisted of preplant application only (100 kg N/ha), a farmer preplant and in-season uniform N application (100 + 112 kg N/ha), and a preplant with in-season VRA (100 + VRA kg N/ha). During the 2016-2017 growing season, the farmer rates were reduced so the respective treatments were the preplant (50 kg N/ha), the farmer application (50 + 100 kg N/ha), and the VRA treatment (50 + VRA kg N/ha). The field was divided into four blocks, and three N treatments were randomly assigned within each block to follow a randomized complete block design ( Figure 3). Each treatment was a field strip of 7 m wide running the entire length of the field (200 m) to simulate full-scale field conditions. The three treatments during the 2015-2016 growing season consisted of preplant application only (100 kg N/ha), a farmer preplant and inseason uniform N application (100 + 112 kg N/ha), and a preplant with in-season VRA (100 + VRA kg N/ha). During the 2016-2017 growing season, the farmer rates were reduced so the respective treatments were the preplant (50 kg N/ha), the farmer application (50 + 100 kg N/ha), and the VRA treatment (50 + VRA kg N/ha).  Figure 4a shows the average daily temperatures for the two growing seasons, as well as the monthly average of a thirty-year period from the CLIMWAT base [57]. The data allowed us to assess whether the temperature conditions of the two experimentation years are close to those of an average year. Figure 4b shows the daily rainfall during the two growing seasons. Average daily temperatures and daily rainfall are inputs for Aqua Crop. Figure 4c shows the monthly rainfall of the two growing seasons and the average monthly rainfall of a thirty-year period from the CLIMWAT database [57]. These data allowed us to assess whether the rainfall conditions of the two experimentation years are close to those of an average year.  Figure 4a shows the average daily temperatures for the two growing seasons, as well as the monthly average of a thirty-year period from the CLIMWAT base [57]. The data allowed us to assess whether the temperature conditions of the two experimentation years are close to those of an average year. Figure 4b shows the daily rainfall during the two growing seasons. Average daily temperatures and daily rainfall are inputs for Aqua Crop. Figure 4c shows the monthly rainfall of the two growing seasons and the average monthly rainfall of a thirty-year period from the CLIMWAT database [57]. These data allowed us to assess whether the rainfall conditions of the two experimentation years are close to those of an average year.

Soil Moisture Measurement Network
EM50 and Em50b data loggers and EC5 and 10HS soil moisture capacitance sensors (Decagon devices, 2365 NE Hopkins Ct, Pullman, WA 99163, USA) were used to dynamically monitor the soil water content during the growing season at different field positions. Based on topographic features, the experimental field was divided into a lower, a middle, and an upper zone. Each zone of the field was equipped with soil moisture sensors installed at a 30 cm depth in the treatment strips of the two central blocks in the fall of 2015.

Soil Moisture Measurement Network
EM50 and Em50b data loggers and EC5 and 10HS soil moisture capacitance sensors (Decagon devices, 2365 NE Hopkins Ct, Pullman, WA 99163, USA) were used to dynamically monitor the soil water content during the growing season at different field positions. Based on topographic features, the experimental field was divided into a lower, a middle, and an upper zone. Each zone of the field was equipped with soil moisture sensors installed at a 30 cm depth in the treatment strips of the two central blocks in the fall of 2015. One additional multi-sensor system (with 10HS sensors) was installed at 10, 30, 50, 70, and Hydrology 2022, 9, 56 7 of 20 90 cm depths in the lower and middle zone of block 2 to monitor the soil moisture of the 1 m profile. Their exact position is shown with the star symbols of Figure 2. A similar arrangement was followed in the fall of 2016, but with four multi-sensor systems to cover all topographic zones. To improve the accuracy of soil moisture sensors, a specific calibration was carried out for the soil conditions of the experimental field.
The sensors were installed in a horizontal position a few days after sowing and the soil water content was measured every 2 h. Disturbed and undisturbed soil samples were taken from the trenches for texture, gravimetric soil moisture, and bulk density determination. To minimize disturbance to the soil structure during 10HS-sensor installation, the soil was excavated layer by layer and was backfilled in its place. The soil was then thoroughly compacted to make good contact with the sensors and restored to its natural density. The good sensor-soil contact was checked with a ProCheck device (Decagon, USA). New grain was sown at the surface to restore plants that were removed. Sensors were installed at a 30 cm soil depth (EC5 sensors) because a previous analysis [58] showed that the measurement at this depth is a good approximation of the average moisture of the 0-100 cm layer.

AquaCrop Model Description
AquaCrop was designed as an evolution from the original FAO Irrigation Drainage Paper 33 [3], which represents the yield response to water as a linear, crop-specific function of the ratio of actual to potential evapotranspiration over a growing season. Steduto et al., 2009 [16] contains a very comprehensive description of the model's conceptual design.
AquaCrop simulates soil water balance and crop growth processes as a function of crop, soil, weather, and management input data on a daily time step [59]. In addition, AquaCrop simulates soil evaporation and crop transpiration explicitly as individual processes. Crop yield response is understood as a function of water consumption [60]. The calculation cycle contains a simulation of green canopy development (CC), crop transpiration (Tr), aboveground biomass production (B), and crop yield [60,61]. AquaCrop has four sub-model components: (i) the soil (water balance); (ii) the crop (development, growth, and yield); (iii) the atmosphere (temperature, rainfall, evapotranspiration (ET), and carbon dioxide (CO 2 ) concentration); and (iv) the management (major agronomic practices such as planting dates, fertilizer application, and irrigation if any) [35].
The AquaCrop model is a water-driven model because it simulates aboveground biomass production in exchange for water transpired by the crop according to the following relation [16,19]: where B n is the cumulative aboveground biomass produced after n days (g m −2 ); Tr i is the crop transpiration per day (mm day −1 ); ETo i is the daily reference evapotranspiration (mm day −1 ); i = 1 . . . n are the days for B n production; WP* is the normalized crop water productivity (g m −2 ). The crop canopy development and phenology are driven by temperature. The water transpired by the crop is determined by the canopy cover [19].
It also simulates daily water movement in and out of the soil profile considering surface runoff, infiltration, capillary rise, soil evaporation, and crop transpiration. As for water balance, the simulation starts with the drainage of the soil profile. Subsequently, water infiltrates into the soil profile (after the subtraction of surface runoff) or moves upward by capillary rise from a shallow groundwater table. Finally, the amount of water lost by soil evaporation and crop transpiration is calculated. In each of the described subroutines, the soil water content is updated at the end of the time step and at each soil layer, resulting in the final daily water content [61]. When calculating the soil-water balance, the amount of water stored in the root zone can be expressed as an equivalent water depth (Wr) or as root zone depletion [2]. This allows the calibration and validation of the model with soil moisture measurements taken in the field.
Unlike all other models, AquaCrop uses canopy cover (CC), not leaf area index (LAI). The use of CC, as opposed to LAI, is meant to introduce simplicity by reducing the overall aboveground growth into just a single growth function [62].
AquaCrop does not explicitly consider nutrient cycles or balances [46,60]. However, soil fertility stress is determined by its expected effect on crop biomass production, using a semi-quantitative assessment to establish the degree of stress resulting from various levels of nutrient deficiency. Under these conditions, the model is based on local calibration to relative biomass under different fertility and salinity conditions. This approach yields a ratio (Brel), calculated as the total dry aboveground biomass at the end of the growing season in a field with soil fertility stress (Bstress) divided by that without soil fertility stress (Bref) (see Equation (2) and [37]).
As shown in Equation (2), Brel ranges from 0%, meaning complete crop failure from nutrient deficiency, to 100%, indicating no nutrient stress. This characteristic of the model allows the user to simulate the combined effect of soil fertility and water stress that is a major strength of the model [46]. Using Bref, harvest index, and Brel, the final yield is estimated.  Table 2 present the conservative and the nonconservative parameters used for the simulations during the two growing seasons. The model was adjusted to crop yield (kg/ha) as measured in plot samples. Measurement of soil water content is a prerequisite for establishing a reliable water balance at the field scale. (ii) AquaCrop was evaluated to estimate the wheat yield under different treatments. The simulation was a two-step procedure: (a) simulations at field positions with 30 cm soil moisture sensors (blocks/repetitions 2 and 3) were used to determine soil water content evolution, yield, and average fertility per zone and treatment (Calibration of the model); (b) using these fertility levels and local soil data, crop production was estimated in other parts of the field with available yield data (blocks/repetitions 1 and 4) under different N treatments. The simulated yields were compared to the field yields (Validation of the model). For the simulation, we used the default in the AquaCrop wheat.cro file (Valenzano).

Simulation Procedure
For the two cultivation periods, the reference evapotranspiration was calculated by the model using as input the minimum, mean, and maximum daily temperature according to FAO 56 protocol [63]. For this evaluation, data from a private database (larisa.meteoclub.gr) were used. For comparative reasons, the mean reference annual evapotranspiration was presented using 30 year data from the FAO database concerning the region of Larissa (near the study area). The CROPWAT model was used for the 30 year average ETo. Meteorological data were provided by CLIMWAT [57].
The physical soil characteristics of all the plots in the upper, lower, and middle zone of the experimental field were determined. They were processed with Soil Water Hydraulic Properties Calculator [64] to calculate hydraulic parameters required by AquaCrop. These included volumetric soil water content at field capacity (FC), permanent wilting point (PWP), saturation (SAT), and saturated hydraulic conductivity (Ksat). Curve number was also adjusted to slope variations according to the Sharply-Williams method [65].

Meteorological and Soil Conditions
The reference evapotranspiration was close to the 30 year average in both years of the experiment ( Figure 5). This was not the case with rainfall ( Figure 4c) and the mean air temperature (Figure 4a). During March and May 2016, the rainfall events were close or higher than 40 mm (Figure 4b). Monthly precipitation was higher than the 30 year average during March 2016, and January and May 2016 and 2017. On the contrary, during December and April of both growing seasons, the monthly rainfall was much lower than the 30 year average. Soil texture was significantly affected by the topographic variation. The slope had coarser soil due to erosion that lowered the level of the ground (see photo in Figure 1) as a result of management practices and crop rotations.
The soil texture analysis indicated the high spatial variability of the experimental field ( Table 3). The lower zone had a predominantly clay texture. The upper zone also had a clay texture (C), but some positions were characterized as clay loam (CL) and two positions with increased sand (sandy clay (SC), in Blocks 1 and 3). Due to slope, the soil of the middle zone was eroded so that SC and sandy clay loam (SCL) were often present. Only the left of the experimental field (Block 4) had reduced erosion due to the lower slope, and class C was the main soil texture.
Wilting point and field capacity were determined using soil texture at the sampling points and Soil Water Hydraulic Properties Calculator [51]. Their values were higher in the lowland area due to the increased fine particles (clay and silt).

Simulation: Model Calibration and Validation
The model was adjusted to the measured plot yields where the 1 m soil moisture profile was installed ( Table 4). The water balance was calculated for the same plots (Table 5). Figure 6 presents the measured and simulated 1 m soil water content during the growing season. Their correlation coefficients are shown in Table 5 (n = 87). For the estimation of water content in the root zone, data of five soil moisture sensors per profile were used. The soil texture analysis indicated the high spatial variability of the experimental field ( Table 3). The lower zone had a predominantly clay texture. The upper zone also had a clay texture (C), but some positions were characterized as clay loam (CL) and two positions with increased sand (sandy clay (SC), in Blocks 1 and 3). Due to slope, the soil of the middle zone was eroded so that SC and sandy clay loam (SCL) were often present. Only the left of the experimental field (Block 4) had reduced erosion due to the lower slope, and class C was the main soil texture.
The position of the sampling refers to the 2016-17 experimental design. C: clay, CL: clay loam, SCL: silty clay loam, SC: silty clay Wilting point and field capacity were determined using soil texture at the sampling points and Soil Water Hydraulic Properties Calculator [51]. Their values were higher in the lowland area due to the increased fine particles (clay and silt).

Simulation: Model Calibration and Validation
The model was adjusted to the measured plot yields where the 1 m soil moisture profile was installed ( Table 4). The water balance was calculated for the same plots (Table  5). Figure 6 presents the measured and simulated 1 m soil water content during the growing season. Their correlation coefficients are shown in Table 5 (n = 87). For the estimation of water content in the root zone, data of five soil moisture sensors per profile were used.
Τhe yield simulation results during the calibration phase and the average yearly fertility using blocks 2 and 3 are presented in Table 6. The simulated yields are close to the    The yield simulation results during the calibration phase and the average yearly fertility using blocks 2 and 3 are presented in Table 6. The simulated yields are close to the measured ones. This indicates the high adaptation of the model and its ability to predict the field data. In the second step, the model estimated crop yield (blocks 1 and 4) using the average soil fertility by treatment and zone during the calibration procedure. The calculated yields for the rest of the plots without soil moisture sensors in the same blocks are shown in Table 7 (validation procedure). The relative variation (%) of simulated final yields in comparison to the measured ones is presented in the Table 8. The correlation between measured and simulated yields is presented in Figure 7 and Table 9. In this table, the same correlation for the upper, middle, and lower zone is also presented. Table 5. Simulation results (E, evaporation; T, transpiration; R, runoff; D, drainage) and Pearson correlation coefficients between measured and simulated mean daily soil water content for each 1 m profile (n = 87). calculated yields for the rest of the plots without soil moisture sensors in the same blocks are shown in Table 7 (validation procedure). The relative variation (%) of simulated final yields in comparison to the measured ones is presented in the Table 8. The correlation between measured and simulated yields is presented in Figure 7 and Table 9. In this table, the same correlation for the upper, middle, and lower zone is also presented. The statistical analysis of the data used for validation (Table 7) revealed a strong relationship between measured yield and fertility, simulated yield and fertility, and measured yield and simulated yield (Table 10).            The statistical analysis of the data used for validation (Table 7) revealed a strong relationship between measured yield and fertility, simulated yield and fertility, and measured yield and simulated yield (Table 10).

Treatment E (mm) T (mm) R (mm) D (mm) R
In the following Figure 8a-d, the fertility factor and the measured yields for the three zones and the three treatments are presented.  In the following Figure 8a-d, the fertility factor and the measured yields for the three zones and the three treatments are presented. Figure 8. Fertility factor in the three zones of the field (%) and for the three treatments (%), respectively (a,b). Measured yields for the three field zones (kg/ha) and for the three treatments (kg/ha), respectively (c,d).

Simulations Using the 100 cm Soil Moisture Profile
In most cases, the simulated soil water content was in agreement with the water content measured in the rootzone ( Figure 6). The success of the model in simulating soil water changes was also evident by the correlation coefficients between simulated and measured moisture values in each soil profile (Table 5). Water balance parameters, such as water requirement (actual evapotranspiration), leaching, and runoff, are needed for scheduling irrigation and fertilization so that N losses can be avoided (Table 5).

Simulations Using the 100 cm Soil Moisture Profile
In most cases, the simulated soil water content was in agreement with the water content measured in the rootzone ( Figure 6). The success of the model in simulating soil water changes was also evident by the correlation coefficients between simulated and measured moisture values in each soil profile (Table 5). Water balance parameters, such as water requirement (actual evapotranspiration), leaching, and runoff, are needed for scheduling irrigation and fertilization so that N losses can be avoided (Table 5).
Considering the field and simulated data in blocks 2 and 3, the use of VRA resulted in high fertility for the first growing season (100% in the lowland and 59% in the slope, Table 4). On the contrary, the lower N rates of the 2016-17 growing season resulted in lower soil fertility for the VRA operation (72% in lowland, 34% in upland, and 24% in the slope). Soil fertility in the slope of the farmer treatment was even lower (21%).
The actual soil fertility (expressed as a percentage of the maximum that could be achieved in Equation (2)) was used for the estimation of simulated yield. Measured yields were higher in 2016 than 2017 when comparing areas with similar conditions, i.e., lower and middle zones of VRA treatment (Table 3). This is probably due to the January 2017 frost (Figure 4a), the lower precipitation (Figure 4b,c), and the lower preplant fertilization during the second growing season. The first growing season (2015-2016) was favored by the high rainfall of March 2016 at a suitable time for plant growth. In the second year, the combination of unfavorable weather and lower preplant fertilization halted crop development, as it is seen in the fractioning of evapotranspiration. Where there was lower plant development, evaporation was higher than transpiration (Table 5).
It is also noted that the model responded to the different locations of the field. This is imprinted on yields with corresponding fertility adjustments in both growing seasons. Fertility was higher in the lower zone than in the middle and upper zone ( Table 4). The higher yields in the lower zone were due to sediment accumulation and deeper and more fertile soil with higher organic matter content and greater mineralization potential. In the middle and upper zones, the fertility was lower as conditions are less favorable (increased erosion, stony and shallow soil). When comparing treatments, the VRA yields were 280 kg/ha higher than those of the conventional farmer practice by taking into consideration the middle zone of 2016-17 (Table 4).
We also observed that drainage losses were reduced, so the local environment is safe from nitrate contamination. Only in the case of the middle zone 2015 and 2016, the leaching reached up to 53-57 mm and was probably caused by lower clay content that allowed the increased infiltration of water (Table 5).
There is some concern about surface runoff, so care must be taken to avoid fertilizer application prior to big storms. The warning concerns even slightly sloped clay soil of lowland areas of the field that may have considerable runoff (Table 5).

Simulations in All the Zones of the Field for Yield Estimation
During the calibration procedure, the average fertility by treatment and zone was estimated using data from blocks 2 and 3 ( Table 6). As expected, the fertility was higher in the lower zone in both growing seasons due to a deeper soil profile with fine soil texture. The fertility of the middle and upper zones was also higher in 2015-2016 (61-78%) than in 2016-17 (27-39%).
The values of simulated and measured yields during the calibration procedure (Table 6) demonstrate the ability of AquaCrop to be adapted to the field results.
Fertility stress, explaining the lower-than-potential observed yields, changed the ranking of several plots. Plots that were expected to give high yields due to their hydraulic properties and nitrogen treatment in reality appeared with a fertility much lower than the optimal. This is partially indicated by the high calcium content of the middle and upper zone [50].
The results show that yield variations were caused mostly by pedoclimatic conditions. The fertility parameter helped us exclude these external factors for a clear comparison of fertilization treatments. The comparison shows that the conventional farmer fertilization was ineffective. In Blocks 2 and 3, the fertility of farmer treatment (average value) was lower than that of the preplant in all zones of 2015-16 and in the middle zone of 2016-17 ( Table 6). The results indicate that large amounts of nitrogen were lost to the environment.
Variable-rate application performed better than the other treatments in almost all cases (Table 6). In the first growing season, the VRA yields of the lower zone were 12% and 15% higher than those of the preplant and farmer treatments, respectively (Blocks 2 and 3). In the second growing season, the VRA yields of the lower zone were 14% and 10% higher than those of the preplant and farmer treatments, respectively. The fact that the higher fertility levels were achieved with less fertilizer indicates that VRA was effective as a means of targeted nitrogen application.
Taking into account the average fertility by zone and treatment (calibration procedure in Blocks 2 and 3), the model was applied to the corresponding plots of Blocks 1 and 4 (validation procedure, Table 7). The relationship between the obtained simulated yields and the measured field yields is presented in Figure 7. Regression analysis showed that this relationship is strongly linear with R 2 = 0.95 (Table 9). In the same table, we can observe that the model slightly overestimated the yields for all field zones. Additionally, the relationship between simulated and measured yields was stronger in the lower (R 2 = 0.99) and upper (R 2 = 0.98) zones compared to the relationship in the middle zone (R 2 = 0.89), probably due to its higher spatial variability (  [26]: R 2 = 0.90. Although the correlation coefficients were high, the range of relative variations observed in the various field zones is important: +6.85% to +47.09% in the lower zone, −41.97% to +143.44% in the middle zone, −18.31% to +50.57% in the upper zone (Table 8). The yield simulations were more successful in the lower zone.
The fertility factor was significantly higher in the lower zone (p = 0.041) (Figure 8a) reflecting the more favorable conditions for crop growth. The fertility factor did not significantly differ between the three treatments (p > 0.05) (Figure 8b). The fertility level was lower in the second growing season possibly due to less preplant fertilizer applied. The early-season frost probably contributed to the lower yields in 2017.
The measured yields were also higher in the lower zone but did not differ significantly from those obtained in the other zones (p > 0.05, Figure 8c). The measured yields did not differ significantly between the three N treatments (p > 0.05, Figure 8d). However, the lower fertilizer rate of VRA (131 kg/ha) in comparison to the farmer practice (212 kg/ha) led to a clear advantage of this technology in terms of improved fertilizer management and protection of the environment. The VRA was more effective under the favorable conditions of the low zone (deeper and more fertile root zone) and in the first growing season. The preplant N treatment was more efficient than the farmer practice.

Conclusions
The model was focused on the water factor and successfully simulated the evolution of soil moisture during the growing season, the estimation of the various water balance parameters, and the expected yield. The model was useful for scheduling irrigation and controlling fertilizer losses through leaching and surface runoff.
The irrigation schedule, soil texture, underground water (if any), initial soil moisture conditions, adjustment of curve number according to the slope variation, and the management practices (use of mulch or not) were needed to estimate the evolution of soil moisture during the growing season. It is important to note that the soil water content profile was close to the measured profile in most simulations, which is a strong indication of a good estimate of the water balance at the field level.
AquaCrop allows the impact of climatic conditions to be dissociated from soil fertility. It is evident that the VRA treatment increased soil fertility with a reduced nitrogen rate (81 kg N/ha less than that of the farmer in 2015-2016) and was more effective in the lower zone.
Based on the analysis of the experimental data, the ability of AquaCrop to estimate wheat yield was demonstrated in a field with high soil and topographic variability. The analysis showed the potential of AquaCrop to be adapted to the field data (calibration procedure) and its ability to give a good estimation (high R 2 ) of the measured yield at the field level (validation procedure). However, simulated yields with significant deviations from the real yields could be locally observed.
The water-driven AquaCrop model was proven particularly useful in optimizing agricultural practices regarding the combined effect of water and nitrogen fertilization. Analyses of this type are particularly useful for areas such as Thessaly, which has environmental issues from intense agricultural activity. The significant falls in the aquifer level, salinization, but also contamination by agrochemicals are important problems that must be addressed with the restructuring of crop patterns and the improvement of agricultural practices. AquaCrop can also help analyze the combined action of water fertilization in the context of climate change.
The weakness of the model is the qualitative approach of soil fertility. The level of fertility should be determined by experiments in the study area and/or by historic yield with corresponding fertilization data.
The model appears to be an efficient tool for evaluating and improving the management practices. To catch up with market and environmental challenges, concerted action is needed that includes improved decision-making supported by scientific evidence. More field trials are required to verify our findings in different fields, crops, and growing seasons in the study region.