Integration of Sentinel-2 Derived Phenology with SoilGrids and ERA-5 Meteorological Data for Operational Crop Yield Forecasting Under Large Scale Applications Emmanuel Lekakis 1, Ioannis Pourikas 2, Vaggelis Oikonomopoulos 1, Theano Mamouka 1 and Polymachi Simeonidou 1 1 AgroApps PC, Koritsas 34, 55133, Thessaloniki, Greece 2 Melissa Kikizas, Larissas - Farsalon Rd (5th km), 41110, Larissa, Greece Abstract Food and feed production must be increased or maintained in order to meet the demands of the earth's population. However, it is obvious that climate change will have a serious negative impact and threaten the productivity and sustainability of food production systems. Therefore, understanding and predicting the final crop production, with a view to adaptation and sustainability, is essential. The need for information on decision-making at all levels, from crop management to adaptation strategies, is constantly increasing and methods of providing such information are urgently needed in a relatively short period of time. Thus arises the need to use effective data such as satellite and meteorological, but also operational tools, to assess crop yields over local, regional, national, and global scales. In this work, an operational approach built on a fusion of satellite-derived vegetation indices, agro-meteorological indicators, and crop phenology is put to test and evaluated in terms of data-intensity, in predicting the yield of durum wheat, at a large-scale application. Keywords 1 Food security, yield prediction, AquaCrop, Operational application, Yield prediction 1. Introduction Understanding and predicting crop production outcomes and identifying changing patterns in major cereals such as wheat, under various climate scenarios and farm management practices geared towards adaptation and sustainability, is of the essence. Hence, the need for the development and efficient use of tools such as models to project food crop cultivation under various scenarios and time scales. Adaptation strategies are probably the only means by which food availability and stability can be maintained or increased to meet future food security needs [1, 2]. In this context, crop monitoring, in- season and post-season yield forecasting play a major role in anticipating supply anomalies, allow well- informed adaptive policy action and market adjustment, prevent food crises, market disruptions, and contribute to overall increased food security [3]. Model-based global estimates show that even incremental adaptation strategies could result in mean yield increases of ~7% at any level of warming [4-6]. This suggests that substantial opportunities may exist if more significant changes in cropping systems are implemented through yield modelling [2]. Process-oriented crop growth models have been extensively applied to simulate and predict production around the world. Among the existing crop models, AquaCrop is used for the estimation of crop yield, under multiple environments. It is a crop water productivity model developed by the Land and Water Division of FAO in 2009 [7, 8]. It simulates yield response to water of crops, crop growth stages, biomass accumulation and yield, among other Proceedings of HAICTA 2022, September 22–25, 2022, Athens, Greece EMAIL: mlekakis@agroapps.gr (A. 1); gpourikas@kikizas.gr (A. 2); voikonomopoulos@agroapps.gr (A.3); thmamouka@agroapps.gr (A. 4); msimeonidou@agroapps.gr (A. 5) ORCID: 0000-0002-1534-7686 (A. 1) ©️ 2022 Copyright for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0). CEUR Workshop Proceedings (CEUR-WS.org) 52 parameters. Models require a substantial number of input data, which creates a limitation in their usefulness for decision-making. Studies support that they are not suitable for predicting crop yields at a large scale because of the massive requirements of inputs and calibration data [9, 10]. However, with the introduction of freely available datasets, the advancement of computational ability and the development of massive data processing technology, the integration of data with different spatial and temporal resolutions into crop growth models has enabled operational and large-scale applications [11, 12]. Such freely available products are the satellite and agro-meteorological data provided by the Copernicus program, as well as Soilgrids [13] which produces maps of soil properties. The objective of this study is to put into test AquaCrop and evaluate the model in forecasting the yield of 184 durum wheat parcels, in two distinct regions, for two growing seasons. A second objective is to assess the model in terms of data intensity, while a third objective is to examine the accuracy of the model for operational large-scale applications by making full exploitation of freely available Copernicus data sets and soil databases. 2. Materials and Methods 2.1. Study Sites and Parcels’ Data Data from a total of 184 durum wheat parcels (Triticum turgidum subsp. durum) were used to evaluate the accuracy of AquaCrop. The parcels were located in the prefectures of Kozani (n=22), northern Greece, and Larissa (n=162), central Greece. The crops were established during the 2019/20 (n=70 Larissa) and 2020/21 (n=114 Larissa and Kozani) growing seasons, under contractual farming (Fig. 1). In these regions, wheat and barley are the main winter (i.e. sowing from October to December and harvest from mid-June to early July) crops grown in rotation with summer crops (e.g. cotton, maize, sunflower). Winter wheat is cultivated under rainfed conditions, but when water is available for irrigation, it is applied during flowering. Figure 1: Study sites and parcels The size of the parcels ranged from 0.1 to 10.5 ha. The geospatial data such as the location and perimeter of the fields, the final production which ranged from 500 to 6740 kgha -1, the sowing dates (3 Nov to 9 Dec 2019 and 25 Oct to 8 Dec 2020), the harvest dates (June 8 to 31) and the irrigation dates were provided by Melissa Kikizas SA. Rainfall during the 2019/2020 growing season amounted to 509 53 mm in Larissa, while in 2020/21 amounted to 450 mm in Kozani and Larissa, with significant variations in rain distribution between the two years and the two areas, especially during the flowering period (April-May). 2.2. Agro-Meteorological, Soil and Satellite Data Gridded agro-meteorological data for the 2019/20 and 2020/21 growing season at a deg resolution of 0.1 × 0.1 (approx. 12.5 km) and of 0.25 × 0.25 (approx. 25 km) were derived from the ERA5-Land and ERA5 reanalysis, respectively, generated by the European Centre for Medium-Range Weather Forecasts and freely distributed through the Copernicus Climate Data Store. They included hourly data consisting of 2-m air temperature, total precipitation, as well as all the necessary variables to obtain the reference evapotranspiration values with the Penman-Monteith equation (FAO56-PM). The soil physical properties for the 184 durum wheat parcels were obtained from the SoilGrids database. SoilGrids is a gridded multiple depth dataset at a 250 m spatial resolution, and it is available worldwide. These soil properties were used to calculate soil water constants (field capacity, permanent wilting point and saturated hydraulic conductivity) with the application of pedotransfer functions [14]. An average number of 25 and 22 cloud free multispectral high-resolution images from Sentinel-2A and 2B were acquired during 2019/20 and 2020/21 wheat growing seasons, respectively. The vegetation indices (NDVI, GreenWDRVI) were calculated based on atmospherically corrected Level-2A-Bottom- of-Atmosphere (BoA) reflectance data. The images were cropped to the polygons’ (parcels’) geometry. NDVI and GreenWDRVI average was assessed at pixels falling within each parcel. The equation to obtain the GreenWDRVI (Wide Dynamic Range Vegetation Index) is as follows [15, 16]: GreenWDRVI=(α∙NIR-Green)/(α∙NIR+Green)+(1-α)/(1+α) (1) where Green is the B3 band of Sentinel-2 MSI, NIR is the near-infrared B8 band of Sentinel-2 MSI and α is a Chl absorption coefficient, equal to 0.1. A 2nd order polynomial relating GreenWDRVI and Leaf Area Index [17] of wheat was used prior to determining the canopy cover: LAI=5.7⋅GreenWDRVI2+1.7⋅GreenWDRVI-0.08 (2) Lastly, the canopy cover of wheat was estimated throughout the growing season on each parcel, with the exponential equation [18]: CCRS (%)=94⋅[1-exp(-0.43∙LAI) ]0.52 (3) where CCRS is the canopy cover derived from remote sensing variables. According to the NDVI and reflected phenology, the wheat parcels were classified under 4 different growth dynamics patterns: 1. a high initial growth rate (n=53); 2. a moderate initial growth rate (n=68); 3. a low initial growth rate (n=41); and a very low initial growth rate (n=22). The average per parcel, NDVI time series graphs are displayed in Figure 2. 54 Figure 2: Reflected phenology of the four wheat growth rate patterns, expressed as the NDVI for the 2020/21 cropping season, versus time. The black line is the representative line of the growth rate. 2.3. AquaCrop Calibration AquaCrop simulates crop yield in four steps: Crop development, crop transpiration, biomass production and yield formation. It calculates the daily soil water balance and divides evapotranspiration into soil evaporation and crop transpiration. AquaCrop describes the foliage development of the crop by the canopy cover (CCAC) which is the fraction of soil surface covered by the green canopy. Transpiration is a function of canopy cover, while evaporation is proportional to the area of soil not covered by vegetation. The canopy cover is multiplied by the reference evapotranspiration (ET r) and the crop coefficient (Kc) to calculate potential crop transpiration. Actual transpiration (Ta) is calculated starting from the potential by accounting for water stress. Then, T a is used for the calculation of crop biomass though its multiplication with water productivity normalized for climate. By using a harvest index (HI), crop yield is obtained by the biomass [19]. Model parameters are grouped into two classes: Conservative and non-conservative. The canopy growth (CGC) and canopy decline (CDC) coefficients are considered two of the most important conservative parameters to calibrate the CCAC [20]. In this study, a number of non-conservative parameters were calibrated, in regards to local management information, such as sowing dates and plant densities, flowering date and duration, starting of senescence and maturity. These parameters are provided in Table 1. Three runs were performed per parcel containing different levels of information: 1. CGC and CDC retrieved from AquaCrop default values on wheat without considering any irrigation event (minimum data input); 2. CGC and CDC calibrated to the canopy cover derived from remote sensing data (CCRS) without considering any irrigation event (medium data input); 3. CGC and CDC calibrated to CCRS, including the irrigation events (one event: n = 45 parcels, two events: n = 12 parcels) applied in their vast majority during the 2020/21 growing season and only in Larissa region (full data input). 55 Table 1 User-specific parameters used in simulation Minimum, Medium Parameter Units Default & Full data input Soil surface covered by an individual seedling at (90%) cm2/plant 1.50 1.50 recover Number of plants per hectare hm-2 4500000 2500000 Maximum canopy cover (CCx) % 96 90 - 93 Calendar Days: from sowing to emergence d 13 13 Calendar Days: from sowing to maximum rooting depth d 93 93 Calendar Days: from sowing to start senescence d 158 178 Calendar Days: from sowing to maturity (length of crop d 197 221 cycle) Calendar Days: from sowing to flowering d 127 150 Length of the flowering stage (days) d 15 20 Maximum effective rooting depth m 1.5 0.3 Reference Harvest Index (HIo) % 48 42 Water productivity (WP) gm-2 15 17 The purpose of the selected simulations was to investigate the potential of AquaCrop to predict biomass and yield under a minimum, medium and a full input data scheme. In order for AquaCrop to simulate accurately canopy cover (CCAC) development, representative CCRS curves of the different wheat growth rates (Fig 2) were selected to calibrate the CGC and CDC. The calibrated values are provided in Table 2. Table 2 Calibrated conservative parameters of AquaCrop used in the simulation Medium & Full data Conservative Parameters Units Default input scheme High Initial Growth class 0.085 Moderate Initial Growth class 0.069 Canopy growth coefficient (CGC) fractiond-1 0.049 Low Initial Growth class 0.054 Very low Initial Growth class 0.035 High Initial Growth class Moderate Initial Growth class Canopy decline coefficient (CDC) fractiond-1 0.0718 0.0605 Low Initial Growth class Very low Initial Growth class 2.3.1. Performance Evaluation Metrics The ability of AquaCrop to predict the yield of the 184 parcels was assessed by adopting a number of statistical metrics. The Model Efficiency (ME), the coefficient of determination (R2), the root-mean- square error (RMSE), the normalized RMSE (nRMSE), the bias and the Willmott’s index of agreement (d) were selected as performance evaluation metrics. 3. Results and Discussion AquaCrop was evaluated for yield prediction under a minimum, medium and a full data input scheme. The calculated statistical metrics that were employed to evaluate the goodness of fit in these approaches, are summarized in Table 5. The reason to follow a minimum data requirements framework was to investigate whether the model has the potential to provide safe results on a simplified approach, 56 with regards to data intensiveness. A successful minimum or at least medium input simulation could encourage the regional or even national implementation of the model that could drive in-season or post- season adaptation strategies for food security. However, in the case of minimum data requirements scheme, the statistical metrics between measured and simulated values revealed that AquaCrop could not simulate wheat yield successfully. This is indicated by the low values of ME, bias, d and R2 of -0.4, -70.6, 0.485 and 0.05, respectively, as well as high absolute and relative magnitude of difference between simulated and measured final yield (RMSE = 1500 kgha -1 and nRMSE = 39.6%). The under- calibrated model significantly underestimated yield which can be attributed to limited rainfall and heat stress during flowering and grain filling stages in 2020/21, that resulted in predicting very low yields. High yield levels obtained by irrigation (>4500 kgha -1), during 2020/21, could not be achieved by the model running under rainfed mode, leading to large prediction inaccuracies. On the other hand, a medium data requirements scheme was applied, based on a calibrated simulation of canopy cover. The calibration of CCAC led to a significant improvement of the statistical metrics. Specifically, the tendency to underestimate was reduced from -70.6 to -38.6, R2 increased from 0.05 to 0.50, thus indicating that the calibrated CCAC model better explained the variance of observed yield values. Smaller estimation errors were received with nRMSE decreasing from 39.6% to 26.5%, RMSE decreasing from 1500 kgha-1 to 996 kgha-1 and a ME increasing significantly from -0.4 to 0.4. In spite of the improved accuracy obtained with the calibration of the CCAC, the results still discourage the large scale and limited data availability application of the model. However, by excluding the irrigated parcels from the statistical analysis (Table 4), the results obtained with CCAC calibration improved significantly and became comparable with those reported in previous studies. The RMSE values of 616 kgha-1 and ME of 0.80 are consistent with ranges reported by [21] (330 to 580 kgha-1) and model efficiencies of 0.78, when the researchers applied AquaCrop for wheat yield simulation. RMSE values of 580 kgha-1 (nRMSE = 11.9%) were also calculated by [22] applying AquaCrop for three years of wheat cultivation. A study conducted to assess the performance of AquaCrop to simulate spring wheat grain yield in Western Canada [23] produced R2 and RMSE values of 0.66 and 743 kgha- 1 , respectively. RMSE, nRMSE and R2 of 550 kgha-1, 8.77% and 0.82, respectively, were found by assimilation of winter wheat biomass retrieved from remote sensing data in AquaCrop, at Beijing, China [24]. AquaCrop was calibrated and evaluated for yield prediction in a field experiment of spring wheat for two growing seasons in Egypt, producing RMSE, d and R2 of 555 kgha-1, 0.93 and 0.84 respectively [25]. Although there was a rather extensive simplification to make the model less data intensive, the metrics of Table 4 clearly indicate the need for a careful calibration of the CCAC to reduce errors and bias when simulating grain yield. In fact, in this case, it is remarkable that the CCAC calibration simulates with great success a wide range of yields, from the very low (750 kgha -1) to high (6220 kgha-1) levels, obtained on rainfed parcels. The adjustment of the CGC and CDC shows a great potential for simulating yields of different wheat varieties with varying growth patterns, under multiyear and multi-environment conditions. AquaCrop can offer a balance between accuracy obtained with mechanistic models and robustness, and can be a valuable tool for yield prediction, particularly considering the fact that it requires a relatively small number of explicit data that can be readily available or easily interpreted to valuable information. 57 Table 3 Statistical evaluation of yield estimates with AquaCrop (n = 184) Minimum Medium Full Statistical metric Units data requirements data requirements data requirements Average estimates kgha-1 3080 3110 3731 Average measured kgha-1 3774 3774 3774 std estimates kgha-1 757 1067 1150 std measured kgha-1 1252 1252 1252 Value range estimates kgha-1 850 – 4670 750 – 4810 750 – 6370 Value range measured kgha-1 500 – 6740 500 – 6740 500 – 6740 ME - -0.4 0.4 0.8 RMSE kgha-1 1500 996 569 nRMSE % 39.6 26.5 15.1 bias kgha-1 -70.6 -38.6 -4.31 d - 0.485 0.827 0.942 R2 - 0.05 0.50 0.80 Table 4 Statistical evaluation of yield estimates with AquaCrop only for rainfed parcels (n = 127) Minimum Medium Statistical metric Units data requirements data requirements Average estimates kgha-1 320 343 Average measured kgha-1 345 345 std estimates kgha-1 600 1128 std measured kgha-1 1233 1233 Value range estimates kgha-1 1500 – 3200 750 – 6220 Value range measured kgha-1 1250 – 6700 1250 – 6700 ME - 0.2 0.8 RMSE kgha-1 1086 616 nRMSE % 31.5 17.9 bias kgha-1 -24.4 -1.3 d - 0.623 0.930 R2 - 0.27 0.75 4. Conclusions This study presents the large-scale simulation and assessment of winter wheat yield with AquaCrop, presenting different levels of input data intensity, based on predictors originating from open-access satellite data, SoilGrids, ERA5-Land and ERA5 climate reanalysis. The model was validated with 184 durum wheat yields in two distinct regions in Greece, for two growing seasons. Remote sensing has offered an unparalleled opportunity to collect data on crop phenology, vegetation development and canopy cover. Coupling this information with freely available datasets to feed crop growth models can significantly minimize time labor intensiveness, cost and delays in the period of data acquisition, making also large-scale applications viable for food security reasons. The results suggest that AquaCrop could support operational platforms for dynamic yield forecasting, operating at the administrative or regional unit scale. 58 5. References [1] D.K. Ray, J.S. Gerber, G.K. MacDonald, P.C. West. Climate variation explains a third of global crop yield variability. Nat. Commun. 6 (2015) 5989 [2] P.N. Kephe, K.K. Ayisi, B.M. Petja, Challenges and opportunities in crop simulation modelling under seasonal and projected climate change scenarios for crop production in South Africa. Agric & Food Secur. 10, 10 (2021) URL: https://doi.org/10.1186/s40066-020-00283-5. [3] D. Dimov, J.H. Uhl, F. Löw, G.N. Seboka, Sugarcane yield estimation through remote sensing time series and phenology metrics. Smart Agricultural Technology. (2022) 2, 100046, ISSN 2772- 3755. URL: https://doi.org/10.1016/j.atech.2022.100046. [4] A. Challinor, J. Watson, D. Lobell, S.M. Howden, D.R. Smith, N. Chhetri, A meta-analysis of crop yield under climate change and adaptation. Nature Clim Change. 4 (2014) 287–291. URL: https://doi.org/10.1038/nclimate2153 [5] C.H. Porter, C. Villalobos, D. Holzworth, R. Nelson, J.W. White, I.N. Athanasiadis, M. Zhang, Harmonization and translation of crop modeling data to ensure interoperability. Environ Model Softw. (2014a) 62, 495–508. [6] J.R. Porter, L. Xie, A.J. Challinor, K. Cochrane, S.M. Howden, M.M. Iqbal, M.I. Travasso, Food security and food production systems (2014b) [7] P., Steduto, T.C., Hsiao, D., Raes, E., Fereres, AquaCrop-the FAO crop model to simulate yield response to water: I. Concepts and underlying principles. Agron. J. (2009) 101, 426–437. [8] D., Raes, P., Steduto, T.C., Hsiao, E., Fereres, AquaCrop-the FAO crop model to simulate yield response to water: II. Main algorithms and software description. Agron. J. (2009) 101, 438–447. [9] Basso, B., Liu, L., 2018. Seasonal crop yield forecast: methods, applications, and accuracies. Adv. Agron., 154, 201–255, URL: https://doi.org/10.1016/bs.agron.2018.11.002. [10] Y., Li, K., Guan, A., Yu, B., Peng, L., Zhao, B., Li, J., Peng, 2019. Toward building a transparent statistical model for improving crop yield prediction: modeling rainfed corn in the U.S. Field Crops Res. (2019) 234, 55–65, URL: https://doi.org/10.1016/j.fcr.2019.02.005. [11] J.S., Bojanowski, S., Sikora, J.P., Musiał, E., Woźniak, K., Dabrowska-Zielińska, P., Slesiński, T., Milewski, A., Łaczyński, Integration of Sentinel-3 and MODIS Vegetation Indices with ERA-5 Agro-Meteorological Indicators for Operational Crop Yield Forecasting. Remote Sens. (2022) 14, 1238. URL: https://doi.org/10.3390/rs14051238. [12] A., Tewes, H., Hoffmann, M., Nolte, G., Krauss, F., Schäfer, C., Kerkhoff, T., Gaiser, How Do Methods Assimilating Sentinel-2-Derived LAI Combined with Two Different Sources of Soil Input Data Affect the Crop Model-Based Estimation of Wheat Biomass at Sub-Field Level? Remote Sens. (2020) 12, 925. URL: https://doi.org/10.3390/rs12060925. [13] L., Poggio, L.M., de Sousa, N.H., Batjes, G.B.M., Heuvelink, B., Kempen, E., Ribeiro, D., Rossiter, SoilGrids 2.0: producing soil information for the globe with quantified spatial uncertainty. SOIL (2021) 7, 217–240, URL: https://doi.org/10.5194/soil-7-217-2021. [14] H., Vereecken, Estimating the unsaturated hydraulic conductivity from theoretical models using simple soil properties. Geoderma (1995) 65(1-2), 81-92. doi:10.1016/0016-7061(95)92543-X. [15] A., Gitelson, Wide Dynamic Range Vegetation Index for remote quantification of biophysical characteristics of vegetation. Journal of Plant Physiology (2004) 161: 165–173. [16] Y., Peng, A.A., Gitelson, Application of chlorophyll-related vegetation indices for remote estimation of maize productivity. Agricultural and Forest Meteorology (2011) 151(9), 1267-1276, ISSN 0168-1923, URL: https://doi.org/10.1016/j.agrformet.2011.05.005. [17] A.L., Nguy-Robertson, Y., Peng, A.A., Gitelson, T.J., Arkebauer, A., Pimstein, I., Herrmann, A., Karnieli, D.C., Rundquist, D.J., Bonfil, Estimating green LAI in four crops: Potential of determining optimal spectral bands for a universal algorithm. Agric For Meteorol, (2014) 192– 193: 140-148, URL: https://doi.org/10.1016/j.agrformet.2014.03.004. [18] D.C., Nielsen, J.J., Miceli-Garcia, D.J., Lyon, Canopy Cover and Leaf Area Index Relationships for Wheat, Triticale and Corn. Agron. J. (2012) 104, 1569. URL: https://doi.org/10.2134/agronj2012.0107n [19] A., Dalla Marta, G.B., Chirico, S., Falanga Bolognesi, M., Mancini, G., D’Urso, S., Orlandini, C., De Michele, F., Altobelli, Integrating Sentinel-2 Imagery with AquaCrop for Dynamic Assessment 59 of Tomato Water Requirements in Southern Italy. Agronomy (2019) 9, 404. URL: https://doi.org/10.3390/agronomy9070404. [20] J., Toumi, S., Er-Raki, J., Ezzahar, S., Khabba, L., Jarlan, A., Chehbouni, Performance assessment of AquaCrop model for estimating evapotranspiration, soil water content and grain yield of winter wheat in Tensift Al Haouz (Morocco): Application to irrigation management, Agric Water Manag, (2016) 163, 219-235, ISSN 0378-3774, URL: https://doi.org/10.1016/j.agwat.2015.09.007. [21] H., Xing, X., Xu, Z., Li, Y., Chen, H., Feng, G., Yang, Z., Chen, Global sensitivity analysis of the AquaCrop model for winter wheat under different water treatments based on the extended Fourier amplitude sensitivity test. Journal of Integrative Agriculture (2017) 16(11), 2444–2458. doi:10.1016/S2095-3119(16)61626-X. [22] M.A., Iqbal, Y., Shen, R., Stricevic, H., Pei, H., Sun, E., Amiri, A., Penas, S., del Rio, Evaluation of the FAO AquaCrop model for winter wheat on the North China Plain under deficit irrigation from field experiment to regional yield simulation. Agric. Water Manage. (2014) 135, 61–72. doi:10.1016/j.agwat.2013.12.012. [23] M.S., Mkhabela, P.R., Bullock, Performance of the FAO AquaCrop model for wheat grain yield and soil moisture simulation in Western Canada. Agric. Water Manage. (2012) 110, 16–24. URL: https://doi.org/10.1016/j.agwat.2012.03.009. [24] X., Jin, L., Kumar, Z., Li, X., Xu, G., Yang, J., Wang, Estimation of Winter Wheat Biomass and Yield by Combining the AquaCrop Model and Field Hyperspectral Data. Remote Sens. (2016) 8, 972. URL: https://doi.org/10.3390/rs8120972. [25] A.M.S., Kheir, H.M., Alkharabsheh, M.F., Seleiman, A.M., Al-Saif, K.A., Ammar, A., Attia, M.G., Zoghdan, M.M.A., Shabana, H., Aboelsoud, C., Schillaci, Calibration and Validation of AQUACROP and APSIM Models to Optimize Wheat Yield and Water Saving in Arid Regions. Land (2021) 10, 1375. URL: https://doi.org/10.3390/ land10121375 60