Calibration and Validation of CERES-Wheat ( Triticum Aestivum ) Model for Simulating Fertilizer Application Rates in Management Zones

Precision agriculture requires precise urea fertilizer application rates for site-specific applications to maximize crop yield across the management zones (MZs). A two years (2010-11 to 2011-12) field experimental study was conducted at Postgraduate Agricultural Research Station, University of Agriculture, Faisalabad, Pakistan to simulate urea fertilizer application rates for four MZs using CERES-Wheat (Triticum Aestivum) model. The model was calibrated using grain yield data of urea fertilizer application rate of 247 kg-urea/ha during growing season of 2010-11 in MZ 1. It was validated against two years independent yield data sets for all treatments ranging from no urea application to 247 kg-urea/ha application in each MZ. The model simulations were found to be acceptable for calibration as well as validation period, as the model evaluation indicators showed root mean square error of 314 kg/ha having its range from 77 to 566 kg/ha, model efficiency of 66% ranging from 24 to 98%, mean percent difference of -4.83%, ranging from -9.93 to 3.70%, against all observed grain yield data in four MZs. Scenario simulations revealed that urea fertilizer application rates of 221, 210 , 208 and 197 kg-urea/ha simulated maximum wheat grain yield of 3679, 3582, 3689, 3690 kg/ha, in MZs of 1, 2, 3 and 4, respectively. These simulated urea fertilizer application might be used to maximize wheat grain yield for each MZs within the field. Furthermore, field verification should be required by applying the simulated urea fertilizer application rates in each MZ.


Introduction
The food and fiber requirements of Pakistan are increasing because of its population growth rate of 2.05%, which has increased population to 176 million in 2011, it is likely to be doubled by 2050 (GOP, 2011).To meet demands of this growing population, there is tremendous pressure on land and water resources of the country (Ali et al., 2011).Whereas, agricultural production of the major food crops including wheat in Pakistan is lower when compared with other wheat producing countries of the world.The average wheat yield per hectare in Pakistan is 27 % lower than the world's average wheat yield per hectare (Arifullah et al., 2009).This situation needs attention of the scientist and engineers to find the factors responsible for low wheat productivity in the country.In addition to cultivar potential, several scientists have reported that topography, climatic conditions, management practices and spatial variability effects within the field are considered as the main factors affecting wheat productivity (Bakhsh et al., 2000a;Abbas et al., 2005;Ayoubia et al., 2007;Marques & Silva, 2008;Begue et al., 2010).
Hence, there is a need to control these factors by developing new strategies.One such strategy has been recognized as the precision farming.Precision farming practices are used to identify, analyze and manage the spatial and temporal variability effects within the field in order to optimize profitability, sustainability and environmental protection (Duffera et al., 2007).To make precision farming as a viable practice, there is a need to delineate site-specific MZs, each with similar values for likely crop response and nutrient concentrations.It has also been reported that site-specific MZs play an important role to manage the spatial variability effects by defining the homogenous sub-field areas within the entire field.The zoning strategy can avoid problems of over-application of inputs in areas where they are not needed and under-application of inputs where they are needed (Cahn et al., 1994;Ferguson et al., 2002).After delineating the site specific MZs, there is difficult to decide optimum/appropriate rate of inputs to maximize wheat yield.
Several researchers have used statistical methods such as decision making tools to understand the functional relationships of crop yield and landscape attributes and to determine optimum/appropriate rate of crop inputs (Chintala et al., 2012a(Chintala et al., , 2012b)).These methods have failed to explain yield variability effects as a function of landscape attributes because crop yield is not only a function of these factors but also a function of temporal variability (Pierice et al., 1995;Mallarino, 1996;Sudduth et al., 1996;Paz et al., 1999;Fraisse et al., 2001;Miao et al., 2006).The temporal variability also has effects on crop yield.Therefore, the process based crop growth models have been used as decision making tools in precision agriculture (Thorp et al., 2008;Awodele & Jegede, 2009).These models have the ability to consider variability in the landscape attributes as well as weather and climatic condition of the specific area.Crop growth models can be used to predict crop growth and development for the particular seasons on basis of the inputs applied in the season.Different researchers have reported that process oriented crop simulation models can be used as efficient and cost effective management tools for evaluation of the impact of different farming practices on the soil and landscape attributes to determine cause effect relationship among soil properties, environment, management practices and crop yield (Awodele & Jegede, 2009;Zhai et al., 2010;Hakojarvi et al., 2010).
Among the numerous crop growth models, the Decision Support Systems for Agrotechnology Transfer (DSSAT) model is the most widely used in agriculture systems.Different modules are included in DSSAT family which can simulate growth of 16 different crops such as maize, wheat, soybeans, rice, and others (Thorp et al., 2008).These modules can simulate crop yield data as a function of climatic factors, management practices, nitrogen and water availability on homogeneous area of the field, which can be helpful in decision making.This capability of DSSAT model makes it more suitable to predict response of the complex system affected by many factors such as crop growth and crop yield in wake of variable soil potential (Paz et al., 1999;Miao et al., 2006;Solaimani, 2009;).Many researchers have used DSSAT crop growth model to study the climate change impact, sustainability, crop management practices and precision agricultural research and have validated the model for different areas and crops (Fetcher et al., 1991;Alexanderrove & Hoogenboon, 2001;Paz et al., 2003;Thorp et al., 2008).Keeping in view, the above cited research wok, this study was designed to use of CERES-Wheat model, which is imbedded in DSSAT family to simulate urea fertilizer application rates for each MZ to maximize wheat yield within the field.

Objectives of the Study
The main objectives of this study were: 1) Calibration and validation of DSSAT model to predict wheat yield based on field experimental data; 2) To simulate management scenario based on MZs for determining urea fertilizer application rates to maximize wheat yield.

Study Area
Field experiments were conducted at the Postgraduate Agricultural Research Station (PARS) of the University of Agriculture, Faisalabad, Pakistan during wheat growing seasons of 2010-11 and 2011-12.The study area is located in Rachna Doab (land between river Ravi and Chenab) with coordinates of 73 o 1′E and 31 o 2′N (Figure 1).The climate of the study area touches two extremes with maximum daily summer temperature reaching 48 o C and winter minimum daily temperature of about 4.8 o C. The average normal precipitation at the study area is about 386 mm (ASP, 2010).Figure 1 shows the location of sampling points in the experimental field of 8 ha in size.The elevation survey of the field was conducted at the forty eight data points, following a regular grid of 24 × 67 m in size using Sokkia C330 Optical Surveying dumpy level.The highest elevation of 185.9 m, above mean sea level, occurred at the east corner of the field whereas the lowest elevation of 185.0 m occurred at the west corner of the field, showing slope in the direction from east towards west.A total of 48 soil samples were collected from top 30 cm of the soil at center of each grid of 24 × 67 m in size using augers prior to sowing of wheat (Figure 1).These soil samples were sent to the Soil Salinity and Water Testing Laboratory, Ayub Agricultural Research Institute, Faisalabad, for determining percent sand, silt, clay, soil EC, pH, soil nitrogen and phosphorus.

Management Practices and Yield Data
The field has previously been divided into four MZs based on the elevation, soil EC, soil nitrogen and %sand contents.The Management Zone Analyst (MZA) and Geographic Information System (GIS) packages were used to delineate the site-specific MZs with in the field.Each zone was of 2 ha in size having 12 experimental units.The size of each experimental unit was of 24 × 67 m (0.161 ha) (Farid et al., 2013a(Farid et al., , 2013b)).Wheat variety of AS-2002 was grown within each site-specific MZ during growing seasons of 2010-11 and 2011-12.Disc plough and tine cultivar were used for primary and secondary tillage operations, respectively.Seed drill was used for sowing wheat in the rows spaced at 15 cm.Weeds were controlled using herbicide.

Treatments Description
The year-wise urea fertilizer treatments were applied as given below.These fertilizer treatments were designed on the basis of recommended dose (120 kg-urea/ha) of urea fertilizer by Punjab Agriculture Department.These fertilizer treatments were designed to determine the urea fertilizer application rates for maximizing wheat grain yield.
Treatments for 2010-11 and 2011-12: Treatment 1 (T1) = 173 kg-urea/ha in single application with 1 st irrigation Treatment 2 (T2) = 123 kg-urea/ha in single application with 1 st irrigation Treatment 3 (T3) = 74 kg-urea/ha in single application with 1 st irrigation Treatment 4 (T4) = Control, no urea fertilizer was applied Treatment 5 (T5) = Variable rate of urea fertilizer application (Urea fertilizer application rates for variable treatments were determined based on the recommended dose of urea fertilizer application rate for the study area, by the Provincial Department of Agriculture, minus 50% of the soil nitrogen, considered as available to the crop during growing season).
Treatment 6 (T6) = 247 kg-urea/ha @ three split applications with 1 st , 2 nd , 3 rd irrigations At maturity, wheat was harvested manually and three yield samples of 1 m 2 in size per plot were collected with their position data using Global Positioning System receiver (Garmen, GPS60).Average of these three samples was considered as yield of the plot.A total of 48 yield samples were collected in each year.These yield samples were threshed manually and grains were dehydrated.The wheat grain and biological yields per unit area were determined.

CERES-Wheat Crop Growth Model
CERES-Wheat crop growth model was used to characterize the yield variability across the delineated MZs.The CERES-Wheat model has the ability to simulate growth, development and crop yield on homogenous area of the field either plot, field or regional scale (Paz et al., 1999;Gabrielle et al., 2004;Thorp et al., 2008).The model operates on a daily time steps and computes the state variable on each day of a year or growing seasons (Fraisse et al., 2001).The CERES-Wheat model requires input data such as management practices (sowing depth and dates, variety, row spacing, emergence date, plant population, irrigation and fertilizer application dates and amount), daily weather data (max and min temperature, rainfall and solar radiation) and soil data.In order to calibrate the model to specific local conditions, crop yield data is also required.
In the present study, the creation of input files for model simulation of each site specific MZs was accomplished by manually manipulating measured model parameters with MZA software.The daily weather data including maximum and minimum temperature, rainfall and solar radiation for the study area were obtained from agro-meteorological observatory of Department of Crop Physiology, University of Agriculture, Faisalabad, Pakistan.The crop management data was recorded throughout the growing seasons.In this study, it was assumed that weather and management data were not spatially variable for the study area.It was also assumed that crop genetic coefficients within single field and for same cultivar are spatially uniform (Thorp et al., 2008).The soil physical and chemical properties and crop yields were considered spatially variable across the each delineated MZ.Therefore, soil input files, file A (average measured data file) and file T (time series data file) were prepared and changed for each MZs involving the same variety to calibrate and validate the CERES-Wheat model site specifically (Pang et al., 1998;Fraisse et al., 2001;Asadi & Clemente, 2003).

Model Evaluation
The performance indicators such as Percent Difference (%D), Root Mean Square Error (RMSE), Coefficient of Determination (R 2 ) and Model Efficiency (EF) were used to estimate the model prediction capability as used by many researches (Nash & Sutcliffe, 1970;Hanson et al., 1999;Ahuja et al., 2000;Bakhsh et al., 2004).

Model Calibration
The CERES-Wheat model was calibrated against the measured data of grain and biological yields in MZ 1 for treatment T1 (173 kg-urea/ha) for growing season of 2010-11 and was validated using data of all treatments for growing seasons of 2010-11 and 2011-12.The calibration procedure involved adjusting the crop genetic coefficients for closer match between observed and predicted data.The best fit value of vernalization coefficient (P1V) was found to be 35 d for wheat cultivar (Wheat variety of AS-2002).It was reported that values of P1V vary from 10 to 60 d being independent of climate type and continent (Rezzouq et al., 2008).The photoperiod coefficient (P1D) was set to be 60% as reported by Saseendran et al. (2004) for semiarid climate of Eastern Colorado.The literature values of grain filling duration (P5) ranged from 332 to 690 GDDo (Saseendran et al., 2004;Nakayama et al., 2006;Yang et al., 2006;Rezzouq et al., 2008;Wang et al., 2012).The P5 value for this study was adjusted to be 650 GDDo.There is slight difference between adjusted value and the value used by Nakayama et al. (2006) in Northen China.However, the value of P5 was acceptable because study area, cultivar and climate were different (Wang et al., 2012).The kernel number coefficient (G1) was found to be 15 k/g which is closer to value reported by Yang et al. (2006) for semiarid climatic conditions in Northern China.
The values of kernel weight coefficient (G2), spike number coefficient (G3), Phonological interval (PHINT) (GDDo) were adjusted as 100 mg, 1.00 g and 85 GDDo, respectively.These values were also within the range as described by Godwin et al. (1989) and Nakayama et al. (2006).The calibrated values of crop genetic coefficients remained same for simulation in all four MZs.Fraisse et al. (2001) reported that same crop cultivar planted at all the monitoring sites had the same calibrated values for these genetic coefficients.The calibration process revealed that model predicted wheat grain and biological yields reasonably well showing %D of 1.54 and -2.99, respectively, between observed and simulated data.

Grain Yield
The calibrated model was used to simulate wheat grain yield data against all other applied urea fertilizer treatments in each MZ for growing seasons of 2010-11 and 2011-12.The average measured wheat grain yield was found to be 3451 kg/ha, while the simulated grain yield was of 3286 kg/ha.During the validation process, the average %D was found to be -4.83(Table 1) for all urea fertilizer treatments, applied during growing seasons of 2010-11 and 2011-12 for all MZs.The average RMSE and EF were determined as 314 kg/ha and 0.66, respectively (Table 3).The coefficient of determination (R 2 ) was found as 0.73.The validated model explained more than 70% average grain yield variability across all applied urea fertilizer treatments, growing seasons and MZs (Figure 2).In general, the average %D varied from 14.43 to -17.73 against the applied urea fertilizer treatments and MZs for both growing seasons.Model provided better simulation results for growing season of 2011-12 than that of growing season of 2010-11 with %D of -2.71 and -6.94, RMSE of 273 kg/ha and 354 kg/ha and EF of 76% and 57%, respectively.The growing season's rainfalls were found to be 49.3 mm during 2010-11 and 23.8 mm during 2011-12.These differences in rainfalls for both the growing season might affect the model performance because model was calibrated for treatment T1 in growing seasons of 2010-11.Overall, the model performance indicators were found to be satisfactory and in a reasonable range (Irmak et al., 2001;Rezzouq et al., 2008) indicating adequate calibration and model inputs.

Grain Yield
The simulated and observed wheat grain yield matched well for all the urea fertilizer treatments in MZ 1.The %D between simulated and observed wheat grain yield was in the range of -10.81 to 6.08% for growing season of 2010-11 and -3.37 to 2.46 for growing season of 2011-12 (Table 3).The RMSE, MPD and EF for growing season of 2010-11 were 207.6 kg/ha, -2.7 and 0.85, respectively, and for 2011-12 were of 77.04 kg/ha, 0.06 and 0.98, respectively (Table 3).In MZ 2, the %D between observed and simulated wheat grain yield for all treatments were found in the range of -14.67 to 8.08% for both the growing seasons (Table 1).The RMSE was found to be 366.06 and 249.47 kg/ha for 2010-11 and 2011-12, respectively.The EF was found to be 40% in MZ 2 for growing season of 2010-11, which is lesser than EF of 74% for growing season of 2011-12 (Table 3).The MPD for 2010-11 was higher than that of growing seasons of 2011-12.The %D for each treatment in MZ 3 was observed as -10.27 to 13.58 for growing season of 2010-11 and -12.43 to 16.62% for growing season of 2011-12 (Table 1).The EF for both the seasons was found to be 79%.The RMSE and MPD were observed as 278.30kg/ha, 336.04 kg/ha and -3.65, 3.70 for both the growing seasons, respectively (Table 3).
The %D between observed and simulated grain yield was in the range of -18.88 to 12.63 for all the treatments for both growing seasons in MZ 4. It was observed that %D range in MZ 4 was higher than those of other three MZs (Table 1).The RMSE was also found to be higher than that of other three MZs, whereas, the EF of 21% in MZ 4 for 2010-11 is lower than the other MZs and than the growing seasons of 2011-12.It was also observed that %D, RMSE and MPD values were lower for MZ 1 than those of MZ 2 and MZ 3. The EF values were higher in MZ 1 than those of MZ 2 and MZ 3. The higher values of %D, RMSE and MPD and lower values of EF in MZ 1 due to the difference in soil and landscape attribute between MZs because the model was calibrated for MZ 1, which is at higher elevation and have higher %sand contents than that of MZ 4. It has been reported that the elevation and soil type affect the soil moisture availability and have impact crop yield (Paz et al., 1999;Bakhsh et al., 2001Bakhsh et al., , 2004)).Overall results indicated that the model performance indicators have closer agreement between observed and simulated wheat grain yield data for both growing seasons in all MZs (Irmak et al., 2001;Miao et al., 2006).Note.%D: Percent Difference.

Biological Yield
The comparison between observed and simulated wheat biological yield in MZs (1-4) indicated that the model was simulated biological yield well for both growing seasons (Table 2).The %D between the observed and simulated values of biological yield for all urea fertilizer treatments was found in the range of ±16% within each MZ for both the growing seasons.The RMSE ranged from 873.651 to 347.81 kg/ha and MPD of -7.08 to 0.89 were observed for the 1-4 MZs during both the growing seasons.The EF in case of biological yield was observed more than 70% in each MZ for both the growing seasons except in MZ 4 during 2010-11.The EF in MZ 4 for 2010-11 was found to be 21% which is much lower than that of other MZs (Table 4).Overall results indicated that during model evaluation and validation, the model performance was satisfactory for each delineated MZ under the given set of conditions.The results also indicated that model can be used as decision making tool regarding the urea fertilizer application rates.Nasim et al. (2010) also reported that CERES-Wheat model can simulate crop growth as well, yield fairly well under the semiarid conditions of the study area.

Grain Yield Production Function for Each MZ
Figures 4a and 4b show the functional relationship between the amount of applied urea fertilizer (kg-urea/ha) and wheat grain yield (kg/ha) in MZ 1 for growing seasons of 2010-11 and 2011-12.The function starts relatively with a steep slope indicating the efficiently use of urea fertilizer to increase the grain yield.The slope starts diminishing as the applied urea fertilizer level increases.In fact, the slope goes to zero as the function reaches to maximum yield.It is clear that beyond this point, further increase in applied urea fertilizer decreased the wheat grain yield.
The maximum grain yield was predicted at average 221 kg-urea/ha in MZ 1 for both the growing seasons.Beyond this rate, the wheat grain yield decreased for increase in urea fertilizer application rates.Similarly, a wheat grain yield production functions were simulated for MZ 2 for both the growing seasons (Figures 4c and  4d).The analysis indicated that maximum grain yield was predicted at 210 kg-urea/ha of applied urea fertilizer levels.Any further increase in applied urea fertilizer level decreased the grain yield.The simulated grain yield production function in MZ 3 showed the similar trend to those of MZ 2 (Figures 4e and 4f).The maximum grain yield was predicted at 208 kg-urea/ha application rate for both the growing season.Figures 4g and 4h showed simulated grain yield production function in MZ 4 for both the growing seasons.The function showed the maximum grain yield at 197 kg-urea/ha application rates.These results showed the spatial distribution of urea fertilizer rates for each MZ across the field.The scenario simulation revealed that urea fertilizer rates of 221, 210, 208 and 197 kg-urea/ha simulated maximum wheat grain yields for MZs of 1, 2, 3 and 4, respectively.These simulated urea fertilizer application might be used to maximize wheat grain yield for each MZs within the field.Furthermore, field verification should be required by applying the simulated urea fertilizer application rates in each MZ.

Figure 2 .
Figure 2. Average simulated vs. observed wheat grain yield (kg/ha) across MZs, applied urea fertilizer treatments over the years

Figure 4 .
Figure 4. Grain yield production functions simulation across the MZs

Table 1 .
Observed and simulated wheat grain yield in four MZs for six urea fertilizer treatments.

Table 2 .
Observed and simulated biological yield in four MZs for six urea fertilizer treatments

Table 3 .
Model performance indicators for wheat grain yield in four MZs

Table 4 .
Model performance indicators for biological yield in four MZs Note.MZ: Management Zones, RMSE: Root Mean Square Error, MPD: Mean Percent Difference, EF: Model Efficiency.