Development of irrigation schedule and management model for sustaining optimal crop production under agricultural drought

Agriculture is vulnerable to drought indicating that the increasing climate crisis requires the necessity of sustainable crop production. In this study, we developed the Irrigation Schedule and Management (ISM) model based on a simulation–optimization (Soil Water Atmosphere Plant-SWAP model with Genetic Algorithm-GA) framework. The ISM model finds an optimal combination of Irrigation Water Amount (IWA) and Irrigation Interval (II) by adjusting Water Stress (WS) responding to environmental conditions (weather, soils, crops and bottom boundary conditions) throughout growing periods. By conditioning the crop (WS) and water management (IWA and II) variables, ISM improves the sustainability of optimal crop productions under different climatic-land surface conditions. The Regional Agromet Center (RAC) site in Faisalabad (at Punjab, Pakistan) was selected to test the proposed ISM model for the field validation/synthetic numerical experiments with various crops (Wheat, Corn and Potato) and soils. We demonstrated that the ISM model that reflects the relationship between crop and water management variables improved the sustainability of crop productions and Water Productivity (WP) compared to those of the conventional irrigation method at the RAC site under various environment conditions. Additionally, the ISM-based long-term crop productions showed the variations along the yearly precipitation changes indicating that optimal combinations of the crop and water management variables are considerably influenced by environmental conditions. Although uncertainties exist, our proposed ISM model can contribute to the establishment of efficient irrigation schedule/management plans under agricultural drought.


Introduction
The frequency and intensity of natural disasters (e.g., floods, droughts, etc.) have been gradually increasing across the world as a result of climate changes (Lettenmaier et al. 1996;Aswathanarayana 2001;Gobiet et al. 2014;Dai et al. 2018). Mishra and Singh (2009) reported that global/local regions might be exposed to more severe droughts in the future . Drought causes widespread and severe economic losses in hydrology, agriculture, and ecosystem across various regions (FEMA 1995;Svoboda et al. 2002). Agriculture is highly vulnerable to drought due to the meteorological conditions. A large amount of available fresh water across the world (about 70% of the existing freshwater supplies) is utilized for agriculture (Evans and Sadler 2008). Therefore, an increase in water scarcity might be inevitable in the future indicating that the agricultural sector needs to be prepared to deal with drought.

3
Soil moisture is a key factor in crop production because plant roots directly absorb soil moisture for crop growth (Morita and Abe 1996). Deficient soil moisture contents disturb normal root growth (e.g., reduction of the number of roots, total root length, limiting of crop growth and yield) while excessive soil moisture conditions cause notable stimulation (Pardales et al. 2002;Sugimoto et al. 2000;Lv et al. 2011). A few studies reported that frequent irrigations with low water amount can improve crop productions (Jin et al. 2014). Appropriate irrigation supplies to wheat crop during the critical growth stage contribute to the improvement of photosynthetic rate, dry matter production, and transportation indicating that crop productions can be significantly increased compared to those under the rain-fed condition (Boughdiri et al. 2014;M'hamed et al. 2015;Saeedipour and Moradi 2011). Dilsher and Khalaf (2014) reported that crop productions were increased by 37.5% and 84.5% with irrigation applied done once/twice were increased compared with that of non-irrigation. Especially, timely irrigation significantly contributes to the improvement of the water use efficiency (WUE) of crops (Sun et al. 2006). Several studies on water-saving concepts (irrigation methods for crops) have been conducted for the efficient use of water in agriculture (Seckler 1996;Molden 1997;Ines et al. 2002). These watersaving concepts contributed to the improvement of water productivity in irrigation (Bos and Nugteren 1990;Bos 1997), but the efficiency indicators only adjusted the water flow ratio for the reuse of water in the irrigation system.
Crop growth highly depends on both climatic-land surface conditions (solar radiation, temperature, rainfall, soil physical characteristics, etc.) and artificial crop/water (Water Stress-WS, Irrigation Water Amount-IWA and Irrigation Interval-II) management (Olesen and Bindi 2002). Considering that artificial crop/water management are only adjustable for outdoor cultivation, irrigation is the key factor that can control crop growth, development and yield (Gornall et al. 2002;Prasad et al. 2008). Thus, the requirement of timely efficient irrigation schedule and management is inevitable for the sustainability of crop production against rapidly deteriorating drought in future. Water and crop management during crop growing seasons can improve the performance of the irrigation system (Ines et al. 2006), while Shin and Jung (2014) suggested the drought index-based irrigation water management approach that can adjust irrigation water based on drought condition. Basically, crop status (WS) under limited water use is highly related to irrigation water (IWA and II) management (Ines et al. 2006;Lee and Shin 2016;Domínguez et al. 2017). In other word, the degree of WS affecting by regional climatic-land surface characteristics significantly influences on determining IWA and II in crop cultivation and water efficiency under irrigation water deficit. Various research works have been done for optimal irrigation decisions, but few studies with respect to timely irrigation schedule and management that can consider the optimal combination (considering both crop and water management) of IWA and II by adjusting WS responding to environmental factors (weather, soil and crop) throughout crop growing seasons have been conducted for the sustainability of crop production. Considering that food security is becoming a global issue due to climate crisis, rapidly increasing world population, continental-/local-disputes, etc., an optimal combination of crop and water management is highly required for the sustainability of optimal (maximum) crop production and water efficiency.
The main goal of this study is to improve the sustainability of optimal crop production by considering timely efficient irrigation management under agricultural drought. The objectives of this research were as follows: (1) to develop an Irrigation Schedule and Management (ISM) model that can determine an optimal combination of IWA and II by adjusting WS responding to environmental factors (weather, soil and crop) based on a simulation-optimization framework and (2) to provide an efficient irrigation management plan for the sustainability of optimal crop production under different climatic-land surface conditions.

Concept of the ISM model
In this study, we developed the ISM model for improving the sustainability of optimal crop production by considering timely efficient irrigation schedule/management plans under irrigation water deficit. The ISM model adapting the simulation (Soil Water Atmosphere Plant-SWAP, Kroes et al. 1999;van Dam et al. 1997)-optimization (Genetic Algorithm-GA, Goldenberg 1989Holland 1975) framework is capable of sustaining optimal crop productions at field regions by determining the optimal combination of IWA and II by adjusting WS responding to environment conditions (weather, soil and crop) in Fig. 1. The ISM model incorporated with GA randomly generates the WS values and determines the IWA and II values based on the generated WS. The physically-based hydrological model of SWAP incorporated within ISM estimates the daily root zone (0-30 cm) soil moisture dynamics. Here, SWAP requires the input data of weather data (solar radiation, max. and min. temperature, wind speed, precipitation, humidity), soil hydraulic properties (α, n, θ res , θ sat , K sat ), and crop information for simulation. Also, WOFOST (WOrld FOod STudies) (van Dam et al. 1997;van Dam 2000) adapted within the SWAP model can estimate crop production (kg ha −1 ) for various crops at field-scales under different hydro-climatic conditions. To obtain the soil hydraulic properties, the ISM model was combined with the near-surface soil moisture data assimilation scheme (Ines and Mohanty 2008) based on Inverse Modeling (IM). The data assimilation scheme extracts the effective soil hydraulic properties that can represent Remotely Sensed (RS)/measured (in situ) soil moisture values by minimizing the differences between the measured and simulated soil moisture contents. The SWAP model also calculates the actual/potential evaporation (E a /E p ) and transpiration (T a /T p ) based on the Penman-Monteith equation. Basically, the ISM model determines the IWA and II variables based on the WS ([1−T a /T p ] × 100) levels. Considering that the WS levels are allowable between 5 and 40% in fields before irrigation (Berkoff and Huppert 1987;Bandaragoda 1998;Polinova et al. 2019;Bonfante et al. 2019;Tolomio and Casa 2020), the ISM model randomly generates the WS values in the range of 10 (wet) to 40% (dry) based on GA. Then, GA derives the IWA and II values based on the generated WS until GA finds the best combination (solution, k = {IWA i=1,…,I , II, WS}) of IWA, II and WS using the objective function (Z) adapting the maximization function (N∈M) in Eq. (1) along with the given GA generations (1 to G). Mathematically, the best combination is determined by obtaining the better Crop Production (CP) in the current generation (G = g) compared to those of the previous generations (G = 1,…,g−1) based on the objective function. Here, the Irrigation Schedule Date (ISD) value was initially determined based on the WS (e.g., [1−T a /T p ] × 100 > 40%) value across the daily time-scale. When the dry days are continued since the ISD is initially determined with the water scarcity (e.g., WS > 40%) for the individual irrigation events, ISM determines II within the constrained ranges (14-60 days used in this study) to trigger irrigation (flood irrigation method used in this study) for removing drought damages to crops. GA usually finds a single solution from the search space indicating that GAs have limitations that can be trapped at the local regions (minima). For this reason, we integrated an ensemble scheme with GA, called the Ensemble-based Genetic Algorithm (EGA), to overcome GA limitations for better searching solutions and determining uncertainty ranges of searched solutions. Here, Z is the objective function, e is the running index of ensembles (E = 10), E is the total number of ensembles, CP is the crop production, g is the running index of EGA, G is the total number of generations in the GA process, N is the individual population, and M is the total number of populations in EGA.

Soil-water-atmosphere-plant (SWAP)
The ISM model adapts the one-dimensional (1-D) SWAP model to simulate soil water flow through the soil profile by using the Richards' equation (Eq. 2) with the input variables of weather, soil, and crop information. The soil hydraulic functions can be described using the analytical expression (van Genuchten 1980) based on the relationship of the soil water content (θ), pressure head (h), and unsaturated hydraulic conductivity (K): where θ is the water content in soil (cm 3 cm −3 ), K is the hydraulic conductivity (cm d −1 ), h is the pressure head (-cm), z is the soil depth (cm), t is the time domain (d), C is the differential water capacity (cm −1 ), and S(h) is the actual soil moisture rate extracted by the crop root (cm 3 cm −3 d −1 ), as shown in Eq. (3): In the equation, T pot is the potential transpiration (cm d −1 ), Z r is the rooting depth (cm), and α w is the reduction factor by the pressure head (h), and it accounts for water deficit/oxygen stress (Feddes et al. 1978). The finite difference approach incorporated in Richards' Eq. (2) (Belmans et al. 1983) can simulate the soil water flow using the soil hydraulic properties and various management scenarios. The soil hydraulic functions of θ(h) and K(h) are derived based on the analytical expressions of van Genuchten (1980) and Mualem (1976) as below, where S e is the relative saturation (−); θ res is the water content (cm 3 cm −3 ) of soil in residual, and θ sat is the water content (cm 3 cm −3 ) of soil in saturation, α (cm −1 ), n(−), m(−), λ(−) are the shape parameters for the soil hydraulic functions, K sat is the hydraulic conductivity (cm d −1 ) of soil in saturation, and m = 1 − 1 n . Various top (weather variables) and bottom boundary (free-drainage and shallow ground water table depths) conditions need to be set manually to simulate the soil moisture dynamics (van Dam et al. 1997). SWAP is comprised of three crop modules: (1) a simple module to consider weather conditions, soil texture, and crop type, (2) a detailed module (WOFOST), and (3) the WOFOST module for simulating grass growth. Moreover, SWAP combines the water management modules to consider irrigation and drainage processes for irrigation events and schedules (van Dam et al. 1997;van Dam 2000). The SWAP model calculates potential evapotranspiration (ET p ) based on the Penman-Monteith equation. Soil evaporation (E p ) and potential transpiration (T p ) are calculated based on the leaf area and the soil cover fraction. ET p is converted into actual ET (ET a ) by adjusting soil evaporation and potential transpiration using empirical relationships as the soil becomes dry. The SWAP model performs well with various meteorological and environmental criteria; detailed information is included in the references (Wesseling and Kroes 1998;Sarwar et al. 2000;Droogers et al. 2000;Singh et al. 2006).

Ensemble-based genetic algorithm (EGA)
Genetic algorithms (GAs) are one of the various optimization algorithms that solve complex problems (Holland 1975;Goldberg 1989). The performance of GAs is highly affected by initial random generator seeds (e.g., idum: − 3000, − 2000, − 1000, etc.). Also, GAs usually have limitations that can be trapped at local regions (minima) at the unknown search space. To minimize the above-mentioned drawbacks of GAs, the Ensemble-based Genetic Algorithm (EGA) was developed for better-searching solutions and assessing uncertainties of searched solutions (Shin and Mohanty 2013). EGA derives the fittest chromosomes (k = {IWA i=1,…,I , II, WS}) in the individual population via the total generation (G). The derived chromosomes have new genetic values via the GA operators, which comprise the selection, crossover and mutation processes and search for more solution spaces. Moreover, EGA restarts from other search spaces to find solutions when the reproduced chromosomes converge to one region before all generations are completed. Using the restarting function, EGA can derive new chromosomes via the creep/jump mutation processes (Ines and Honda 2005). Furthermore, EGA remembers the best (elite) chromosomes in the previous generation (g − 1) and delivers them to the next generation (g) (Ines and Mohanty 2008). To search more unknown search spaces, a random resampling (ensemble) algorithm (IBM Programmers' Guide) (Efron 1982) was integrated with EGA. The constrained ranges of solutions by GA are shown in Table 1.

Numerical Experiments
The field validation and synthetic numerical experiments were conducted to test the performance of ISM model for improving the sustainability of crop productions. The numerical experiments were as follows: (1) Case Study 1: field validation experiments and (2) Case Study 2: synthetic numerical experiments, respectively. The Regional Agromet Center (RAC) site located in the district of Faisalabad at Punjab (Pakistan) under the semi-arid/arid climate condition was selected for conducting the Case Studies (Fig. 2). The wheat crop was cultivated during the winter periods (Nov. to Apr.) of 2007-2008 and 2008-2009 with sandy loam soil (predominant). The soil hydraulic properties (α: 0.018, n: 1.260, θ res : 0.044, θ sat : 0.340, K sat : 145.0) at 0-40 cm depth measured by Ahmad (2002) and Time Domain Reflectometry (TDR)-based in situ root zone (0-30 cm) soil moisture measurements (Hanif 2008) were used for the field validation of near-surface soil moisture data assimilation scheme that simulates the soil moisture dynamics. The in situ soil moisture data were measured at the multiple soil depths of 5, 10, 20, and 30 cm at the RAC site during the winter seasons (2007-2008 and 2008-2009). The averaged root zone soil moisture measurements in the multiple soil depths were used for the field validation experiments. The weather data (e.g., precipitation, maximum and minimum temperature, wind speed, humidity and solar radiation) were collected from the Pakistan Meteorological Department for the RAC site in Faisalabad (Pakistan). Also, the Field-scale Crop and Irrigation Management Information-FCIMI (Hanif 2008) such as crop growth period, sowing/harvesting dates, irrigation amounts (mm), irrigation type (flood), number of irrigations and measured crop production (kg ha −1 ) for wheat are shown in Table 2.
In Case Study 1, the near-surface soil moisture data assimilation scheme was calibrated for extracting the effective soil hydraulic properties (α, n, θ res , θ sat , K sat ) from the TDRbased in situ soil moisture measurements (Hanif 2008) during the winter season (2007)(2008). The soil hydraulic functions of θ(h) and K(h) based on the extracted (calibrated) soil hydraulic parameters were compared to those of the measured soil hydraulic parameters (Ahmad 2002). Based on the calibrated soil hydraulic parameters (2007)(2008), the root zone (0-30 cm) soil moisture dynamics were predicted and compared with the TDR-based measurements during 2008-2009 for validation of the near-surface soil moisture data assimilation. Then, the crop productions estimated by the SWAP (only)/ISM (integrated with SWAP) models were compared with the measured crop productions (Hanif 2008) during the winter seasons for calibration (2007-2008) and validation (2008-2009). In Case Study 2, the synthetic numerical experiments were conducted to analyze the impacts of various soils/crops (corn, potato and wheat) on crop productions. Various referencebased soil hydraulic properties (Clay Loam-CL, Loam-L, Sandy Clay Loam-SCL, Silty Clay Loam-SiCL, Silt Loam-SiL, and Sandy Loam-SL) were used from the UNSODA databases (Leij et al. 1999). Figure 3 shows the soil hydraulic functions of θ(h) and K(h) derived from the UNSODA database. The Bottom Boundary Condition (BBC) of free-drainage  [2007][2008]. We assumed that the soil profile is homogeneous at the soil depth of 0-200 cm from the soil surface (Fig. 4), and the soil profile was discretized into the 33 soil layers (1-10th layer: the intervals of 1 cm, 11-20th layer: the intervals of 5 cm, 21-32nd layer: the intervals of 10 cm, and 33rd layer: the interval of 20 cm). The RAC site has the ground water table depth of 800 ~ 1000 cm from the soil surface (Irrigation Department, http:// irrig ation. punjab. gov. pk/). For this reason, the BBC was set as the free-drainage condition under the assumption that the initial condition is in equilibrium with the presence of BBC under the modeling condition. Additionally, the Water Productivity (WP = CP/IWA, kg ha −1 mm −1 , French and Schultz 1984;Bessembinder et al. 2005;Passioura 2006) values were estimated for evaluating irrigation water use efficiency. To determine uncertainties of the model outputs, we set the number of ensemble (E = 10) for individual populations in EGA (Table 2). Also, the Pearson's correlation (R) and root-meansquare error (RMSE) were used to evaluate the model performance as shown in Eqs. (6) and (7), where obs t : the observed soil moisture with time (t),

Case study 1: Field validation experiment
In this study, we developed the ISM model that can improve the sustainability of crop productions with the optimal combination of IWA and II determined by adjusting WS responding to the environmental conditions (weather, soils and crops). The effective soil hydraulic properties as the input variables to SWAP were derived using the near-surface soil moisture data assimilation scheme. Note that the actual irrigation events (individually 75 mm for 1st: Dec. 11th in 2007 and 2nd: Jan. 8th, and 3 rd : Mar. 2nd, 4th: Mar. 25th in 2008) at the field site were included in the data assimilation process. Overall, the RAC site has the less precipitation amounts during the crop growing seasons (Nov. to Apr.) in the years of calibration (2007)(2008) and validation (2008-2009) while the relatively higher precipitation amounts and frequencies were generated in May to June. Figure 5a shows the comparisons of measured and nearsurface data assimilation (ensemble-averaged)-based soil moisture dynamics for calibration at the RAC site during the winter period (2007)(2008). The calibrated soil moisture dynamics (R: 0.908 and RMSE: 0.014) matched well with the measurements. The good agreement between the soil hydraulic functions of θ(h) and K(h) based on the measured (Ahmad 2002)/calibrated soil hydraulic properties support the robustness of near-surface soil moisture data assimilation scheme in Figs. 5b, c. For validation (2008)(2009), the soil moisture dynamics based on the calibrated soil hydraulic properties (2007)(2008) were estimated using the SWAP model in a forward mode and compared with the TDR measurements in Fig. 6. Although uncertainties were slightly increased compared to the calibration results (Fig. 5a), the estimated root zone soil moisture dynamics (R:  Near-surface soil moisture data assimilation (Ines and Mohanty 2008) 1 3 shown in Table 2) matched well with the measured data (3212 kg ha −1 ) during the period of 2007-2008 while the estimated crop production (3260 kg ha −1 ) for the period of 2008-2009 was slightly underestimated compared to the measured wheat production (3707 kg ha −1 ) with increased uncertainties. Although uncertainties slightly increased in the period of 2008-2009, these results support the good performance of SWAP for simulating crop productions at field-scales.
With the multiple WS values from 10 to 40%, we confirmed that the ISM-derived irrigation required (dry) days were increased, respectively. In general, the conventional irrigation method empirically determines the IWA and II by considering the increasing WS due to agricultural drought. Compared to the conventional method, we confirmed that the ISM model determines the ISD (indicating irrigation events) right after the generations of WS throughout the crop growing period under the fixed WS (10, 20, 30 and 40%) conditions. These results supported that the ISM model can respond immediately to changes of the environmental conditions (weather, soils and crops) throughout the growing period.
In Table 3, the ISM-based wheat crop productions were estimated and compared with the measurements for the years of 2007-2008 and 2008-2009. The measured wheat crop productions of 3212/3707 kg ha −1 were shown in Table 2 with the precipitation amounts (108 mm for 2007-2008 and 96 mm for 2008-2009) during the crop growing periods, respectively. Overall, the ISM-based crop productions (4289 kg ha −1 for 2007-2008 and 6297 kg ha −1 for 2008-2009) were significantly improved, while the estimated IWA values slightly increase compared to the measured data (FCIMI in Table 2) by the conventional irrigation method (manually operated by farm mangers). It is inferred that the increased IWA might be due to the prompt responses of ISM model to the generated WS at the field site. The decreased II values (6.1 and 4.4 days for 2007-2008 and 2008-2009) compared to those (about monthly irrigation intervals) of the conventional method (Table 2) support this inference. Also, the WP values of ISM for the periods of 2007-2008 (12.2 kg ha −1 mm −1 ) and 2008-2009 (22.7 kg ha −1 mm −1 ) were higher than those (11.8 and 13.6 kg ha −1 mm −1 ) of the conventional method (Table 2), respectively. Especially, the ISM-based WP value for 2008-2009 was significantly increased indicating that the ISM model contributed to both the optimal crop production and irrigation water use efficiency. These results supported that the ability that can consider the relationship between crop (WS) and water management (IWA and II) variables allows the ISM model to respond immediately to environmental changes. Figure 9 presents the rain-fed/ISM-based crop productions and IWA with the combinations of different soil types (CL, L, SCL, SiCL, SiL, and SL) and crops (wheat, corn and potato) at the RAC site during the crop growing period (2007)(2008) under the synthetic conditions. The crop productions under the rain-fed condition were significantly affected by different soil types. Physical characteristics of the L (38.3 cm/day), SiL (30.5 cm/day) and SL (41.6 cm/ day) soils in Fig. 3 weakens the crop productions compared to those of the other soils (1.8 ~ 9.78 cm/day for CL, SiCL, and SiL) indicating that crops with the predominant L, SiL and SL soils having the relatively high infiltration at field regions could be highly vulnerable to water deficit under agricultural drought. Overall, the ISM-based crop productions for the whole soils and crops were highly improved, especially for the L, SCL and SL soils. We also confirmed that the ISM-based IWA values were proportional to the increased crop productions, but the estimated crop productions using the ISM model were still fluctuated by the physical soil characteristics. These trends for wheat, corn and potato were similarly shown in Fig. 9a-c, but the differences ([rain-fed CP/ISM-based CP] × 100%) between the rain-fed and ISM-based CP values showed variations for the combination of soils and crops in Fig. 9. It is inferred that alternative crop selections by considering physical soil characteristics at local regions can improve to the sustainability of crop productions. Additionally, the SL soil having the highest infiltration rate in Fig. 3 was extremely vulnerable to water deficit condition, but the crop productions with the optimal IWA and II were dramatically increased with comparing to those of the other soils. The soil pore space for the SL soil is restively larger than the other soils indicating that air circulation and diffusion through large pore space actively occur in the soil matrix and are favorable to crop growth (Goorahoo et al. 2001;Du et al. 2018). The similar trends were shown for the L and SCL soils, respectively. Although the crop productions have variations based on the physical characteristics of soils, the ISM model sustained the relatively maximum levels of crop productions for the whole soils and crops.

Case study 2: Synthetic numerical experiment
In Fig. 10, we tested the impacts of BBCs such as free-drainage, SGW − 200 cm, SGW − 150 cm and SGW − 100 cm on estimating the wheat crop production at the RAC site during the crop growing period (2007 − 2008) under the synthetic condition, respectively. The crop production was only 4 kg ha −1 for free-drainage under the rain-fed condition, while SGW − 200 cm showed the crop production of 923 kg ha −1 . When the SGW table depths of − 150 and   . 11 The comparisons of long-term wheat crop productions and crop (WS)/water (IWA and II) management variables estimated by the ISM model with the weather data from 1989-1990 to 2008-2009 model can be expanded to large regions by combining RS soil moisture products such as advanced microwave scanning radiometer for Advanced Microwave Scanning Radiometer for EOS-AMSR-E (Njoku 2008), Soil Moisture Ocean Salinity-SMOS (Kerr et al. 2001), Soil Moisture Active and Passive-SMAP (Entekhabi et al. 2010), etc. Thus, these findings showed the potential application of ISM model in large spatial domains by combining the RS soil moisture data and near-surface soil moisture data assimilation schemes. We estimated the long-term wheat crop productions from 1989-1990 to 2008-2009 (18 years) at the RAC site with the long-term weather conditions in Fig. 11. Note that the years of 2001-2002 and 2005-2006 were excluded because of the weather data missing. Also, the same sowing/harvesting dates in the winter season of 2007-2008 were used for the whole years. The sowing/harvesting dates need to be adjusted based on the weather conditions, but these are not considered for the long-term crop production estimations owing to the lack of data. Overall, the precipitation amounts highly influence on the yearly crop productions compared to the relatively monotonous changes of yearly (growing period) averaged solar radiation (15,561.4/1412.5 for average/SD, KJ/m 2 ) and temperature (18.5/3.1 for average/SD, °C) values. Especially, the ISM-derived IWA, II and WS values showed the yearly variations and were inversely correlated with the occurrence days of effective precipitation (more than 5 mm). These findings indicated that optimal combinations of the crop and water management variables are considerably affected by yearly environmental conditions. Although uncertainties exist in these results, the ISM model improved the sustainability of crop productions with the optimal IWA, II and WS variables responding to the environmental conditions. Thus, these results showed the possibility that the ISM model can be used to establish longterm crop and irrigation management plans combining with climate change scenarios in future.

Conclusions
In this study, the ISM model was developed to improve the sustainability of optimal crop productions by providing timely efficient irrigation schedules and management under agricultural drought. The ISM model finds the optimal combination of IWA and II by adjusting WS responding to the environmental conditions based on the simulation-optimization scheme. The field validation and synthetic numerical experiments were conducted to test the proposed ISM model at the field region. The numerical experiments were conducted as follows: (1) Case Study 1: field validation experiments and (2) Case Study 2: synthetic numerical experiments, respectively. The RAC site in Faisalabad (at Punjab, Pakistan) was selected to conduct the Case Studies.
In Case Study 1, the ISM model improved both the sustainability of crop productions and water productivity during the periods of 2007-2008 and 2008-2009 compared to those of the conventional irrigation method. Especially, the estimated WP value was significantly increased for [2008][2009] indicating that the ISM model contributed to both the sustainability of optimal crop production and irrigation water use efficiency. These findings demonstrated that the ability that can consider the relationship between crop (WS) and water management (IWA and II) variables allows the ISM model to respond immediately to environmental changes.
In Case Study 2, the ISM model was tested and improved the crop productions under various environment conditions (weather, soils, crops and BBCs). The L, SCL and SL soils having the high infiltration rate were extremely vulnerable to water deficit, but air circulation and diffusion through large pore space for these soils (L, SCL and SL) are favorable to crop growth with the appropriate crop and water management plans suggested by ISM. Also, the CP values were improved with the presence of SGW becoming closer to the soil surface in the soil profile, because crop root can access easily to SGW due to the root moisture extraction force and capillary phenomena under the rain-fed condition. However, the ISM-based crop productions usually reached to the maximum levels with the presences of BBC (free-drainage and SGW − 200 ~ − 100 cm). Additionally, the long-term wheat crop productions from 1989-1990 to 2008-2009 were estimated and showed the variations along the yearly precipitation changes. Especially, the yearly crop productions were inversely correlated with the occurrence days of effective precipitation (more than 5 mm). These results supported that optimal combinations of the crop and water management variables are highly influenced by environmental conditions.
In this study, the applicability of newly developed ISM model was tested to improve the crop productions for various environment conditions under different hydro-climatic conditions. Although the field validation experiments were conducted at one region (the RAC site) in Faisalabad (Pakistan), our findings demonstrated the robustness of the ISM model for improving the sustainability of optimal (maximum) crop productions at field regions. Thus, our proposed ISM model can be used for establishing efficient irrigation schedule and management plans under agricultural drought.
Funding Korea Ministry of Environment, The SS (Surface Soil conservation, Yongchul Shin, management) projects, Yongchul Shin, 2019002820002, Yongchul Shin Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.