Global Optimization of Cultivar Trait Parameters in the Simulation of Sugarcane Phenology Using Gaussian Process Emulation

The global optimization of parameters in process-based crop models is often considered computationally expensive. Gaussian process (GP) emulation is a widely used method for reducing the computational burden of the optimization process. Total above-ground biomass and cane dry weight of three Thai sugarcane cultivars (KK3, LK92-11 and 02-2-058) collected under rainfed and irrigated conditions were used to optimize cultivar-specific parameters in the Agricultural Production Systems sIMulator (APSIM)-Sugarcane crop model through a GP emulation. GP emulators were trained and validated to approximate APSIM-Sugarcane model and then used for optimizing the cultivar-specific parameters through the differential evolution algorithm. Resulting optimized parameters allowed to obtain simulations that quite well approximated the observed biomass and CDW (validation results between simulated and observed yields: R2 0.93–0.98; normalized root mean squared error: 5–22%; Willmott’s agreement index: 0.87–0.99). The best parametrization was obtained under the lowest water stressed conditions. Based on these results, we suggest that GP emulation can be efficiently implemented for the parameterization of computationally expensive simulators.


Introduction
Sugarcane is an important crop for sugar and bioenergy worldwide. As the fourthlargest sugar producer and the second-largest sugar exporter in the world, sugarcane has become the most important crop in Thailand's agriculture [1].
However recent evidence indicates that Thailand's sugarcane production is highly affected by climate change. For instance, sugarcane output in the crop year 2020/2021 was recorded as 66.7 million tonnes, which is down from 74.9 million tonnes in the last crop year primarily due to drought [2]. Thus, identification of suitable management strategies to cope with such temporal and spatial variability of sugarcane production has become important for the Thai sugarcane industry.
Obtaining accurate model predictions in process-based crop modelling requires cultivar trait parameterization. Parameter-optimizing techniques are frequently used in the parameterization process [6,16,17]. Seidel et al. [18] have reported that almost 50% of the 211 respondents to their survey on parametrization practices searched for the best-fit parameters using trial and error. However, Harrison et al. [19] indicated that such manual parameterization techniques allow only a small number of parameters to be calibrated. As a result, many of the parameter combinations remain uninvestigated.
Holzworth et al. [20] suggested that the use of a much more objective and reproducible parameterization and validation methodology is of great importance to the formulation of models that are to be applied to the growth of agricultural crops. Studies have recently focused on the use of complex statistical approaches for global optimization of parameters in crop models. Formal and informal Bayesian methods are common approaches that have been used for parametrization of crop models. For instance, Sheng et al. [21] used a Generalized Likelihood Uncertainty Estimation (GLUE) method and Differential Evolution Adaptive Metropolis (DREAM) method for an APSIM model of maize production. Sheng et al. [21] have concluded that similar performance is obtained with both GLUE and DREAM, but they recommend GLUE because it is easier to use. Sexton et al. [22] have used GLUE and a Markov Chain Monte Carlo (MCMC) model in an APSIM model of sugarcane production. Although both methods produced acceptable results, Sexton et al. [22] recommend MCMC because its statistical background is well documented. However, due to the relative advantages of diverse parametrizing methods and internal structural differences of crop models, evaluating the performances of various parameterization methods is further important in providing comprehensive guidance for the method selection.
The current study was focused to use a differential evolution (DE) algorithm for global optimization of cultivar trait parameters implemented in the APSIM-Sugarcane model. According to Georgioudakis and Plevris [23], DE has become one of the most popular optimization algorithms used in continuous optimization problems. For instance, the extensive survey conducted by Bilal et al. [24] reviewed nearly 283 research articles focused on the use of DE during the last 25 years.
Although there are a number of promising methods that could be used in optimization, Harrison et al. [19] have explored several parameterization techniques and have reported that the length of time required for convergence is a factor that limits obtaining accurate results. Moreover, according to Saltelli et al. [25], global optimization of parameters in process-based crop models is often considered computationally expensive. Therefore, carrying out the number of simulations required for optimization may not be feasible and may be extremely time-consuming. A versatile way to reduce the computational time is the use of a Gaussian process (GP) [26] for DE optimization. GP emulators are surrogate models, and they are computationally low expensive. An emulator of sufficient accuracy can be used as a substitute for the original simulator (APSIM), and parametrization can be based on the emulator. The feasibility of substituting APSIM-Sugarcane for sensitivity analysis based on GP emulation has previously been evaluated in several studies [27][28][29][30]. However, there has not yet been any study of the use of the GP emulators to optimize parameters in the APSIM-Sugarcane model nor an examination of the accuracy of the optimization method.
The aim of the present study was the utilization of a GP-based DE algorithm for global optimization of cultivar trait parameters implemented in the APSIM-Sugarcane model. Total above-ground biomass (Biomass) and cane dry weight (CDW) at 30-day intervals from 90 to 390 days after planting of three sugarcane cultivars from Thailand obtained under rainfed (Rf) and irrigated (Ir) conditions were used for the optimization of cultivar trait parameters. The influence of field experimental characteristics in the parameter optimization of APSIM-Sugarcane was examined.

Materials and Methods
The current study used observations of sugarcane biomass and CDW of three Thai sugarcane cultivars (Section 2.1) for global optimization of cultivar trait parameters implemented in the APSIM-Sugarcane model. APSIM-Sugarcane simulations were initially prepared (Section 2.2) and used for training and validating several emulators (Section 2.3), one for each experimental condition. Emulator accuracy was evaluated. Emulators were then used instead of APSIM-Sugarcane for DE-based global optimization (Section 2.4). Several optimized parameter ensembles for each cultivar were obtained from the above process and used in APSIM-Sugarcane to obtain simulated outputs (sugarcane biomass and CDW). Simulated outputs were then validated with observations for the selection of an optimal parameter ensemble for each Thai sugarcane cultivar (Section 2.5).

Study Field
Data from the three field experiments (namely A, B1 and B2) conducted by Preecha et al. [6] in Khon Kaen (KK), northeast Thailand (16.48 • N 102.82 • E; 181 m elevation), was used for the study. Observed data of each experiment composed of sugarcane biomass and CDW of three sugarcane cultivars (namely, 02-2-058, KK3 and LK92-11) at specific dates (days after planting-DAP). Details of the experiments were indicated in Table 1. Mean values of the observed data from the field experiments are indicated in supplementary materials Table S1. More details about the experimental arrangement of the field experiment can be found in Preecha et al. [6]. According to the Köppen-Geiger system, KK has a tropical wet-dry climate (Aw) [31]. Figure 1 shows mean values of monthly rainfall, daily minimum temperatures and daily maximum temperatures and daily solar radiation during study period (2010/12-2012/12). At first view, in 2011, the weather seems wetter, albeit in the first three months and the last two months there was no precipitation (Figure 1). Temperatures do not differ much between 2011 and 2012. In 2012, it seems that the weather was drier and more irradiated. However, experiment A experienced the highest water stressed condition for sugarcane growth compared to B1 and B2. Because it was conducted under Rf conditions and a prolonged dry period was recorded in the early and late stages of sugarcane growth ( Figure 1: year 2011). There was less chance for water stress conditions manifest itself in experiment B1 because during the period rainfall was well distributed (Figure 1: year 2012) throughout the year compared to experiment A when irrigation was applied. However, B2 was conducted under Rf conditions in the same year with B1, but water stress conditions were lower than experiment A and higher than B1. Soil properties in the experimental field are shown in Table 2 [32]. 2012) throughout the year compared to experiment A when irrigation was applied. However, B2 was conducted under Rf conditions in the same year with B1, but water stress conditions were lower than experiment A and higher than B1. Soil properties in the experimental field are shown in Table 2.

APSIM Simulation
APSIM [33] is a comprehensive model developed to simulate biophysical processes in agricultural systems. APSIM contains interconnected biophysical and management models to simulate the processes that occur in systems comprising soil, crop, trees, pasture and livestock, and it has the flexibility to integrate non-biological farm resources. A more detailed description of the APSIM platform can be found on Holzworth et al. [33] and The APSIM Initiative [34].
The next generation of the APSIM [33] Sugarcane model [34] was used for this study. In the APSIM-Sugarcane model, crop dry weight accumulation is determined via the conversion of the intercepted radiation into biomass, which is based on radiation-use efficiency (rue). The area of the crop leaf canopy in which the radiation is intercepted is expanded as a function of temperature. The partitioning of biomass occurred among various components of the plant (sucrose, structural stem, leaf, cabbage and roots) based on the phenological stage of the crop: 1. sowing: from sowing to sprouting; 2. sprouting: from sprouting to emergence; 3. emergence: from emergence to the beginning of cane growth; 4. begin cane: from the beginning of cane growth to flowering;

APSIM Simulation
APSIM [33] is a comprehensive model developed to simulate biophysical processes in agricultural systems. APSIM contains interconnected biophysical and management models to simulate the processes that occur in systems comprising soil, crop, trees, pasture and livestock, and it has the flexibility to integrate non-biological farm resources. A more detailed description of the APSIM platform can be found on Holzworth et al. [33] and The APSIM Initiative [34].
The next generation of the APSIM [33] Sugarcane model [34] was used for this study. In the APSIM-Sugarcane model, crop dry weight accumulation is determined via the conversion of the intercepted radiation into biomass, which is based on radiation-use efficiency (rue). The area of the crop leaf canopy in which the radiation is intercepted is expanded as a function of temperature. The partitioning of biomass occurred among various components of the plant (sucrose, structural stem, leaf, cabbage and roots) based on the phenological stage of the crop:
emergence: from emergence to the beginning of cane growth; 4.
begin cane: from the beginning of cane growth to flowering; 5.
flowering: from flowering to the end of the crop; 6.
end of the crop: crop is not currently in the simulated system.
Several conditions, such as extremes of temperature, plant nitrogen deficits, or soil water shortage or excess, can limit growth during these stages. The APSIM-Sugarcane model controls the above-mentioned biophysical processes of sugarcane growth via several parameters. Such parameters could be categorized into four groups, namely plant-crop parameters, ratoon-crop parameters, cultivar-specific parameters, and soil-and-climate parameters. Plant_crop parameters govern the processes of sugarcane growth and partitioning, water usage, water and temperature stresses, frosting, nitrogen contents and nitrogen stresses. Same as Plant_crop parameters, Ratoon_crop parameters govern the above processes in ratoon stage of the sugarcane. When parametrizing the model for a new cultivar, cultivar-specific parameters are of special interest because they represent key traits of a particular cultivar. Table 3 provides a description of cultivar trait parameters used in the APSIM-Sugarcane model.
Although the APSIM-Sugarcane model does not contain rue and transp_eff_cf among the cultivar-specific parameters, they were included in the current study because Sexton et al. [29] have identified rue as a highly sensitive parameter for biomass yield estimation and transp_eff_cf as a highly sensitive parameter for biomass yield estimation under water stress conditions. Bandara et al. [27] have obtained similar results related to rue and transp_eff_cf for estimation of CDW. If the supply of soil water is not a growth-limiting factor, the APSIM-Sugarcane model controls dry matter assimilation via radiation interception and rue. However, water supply, vapor pressure deficit, and transp_eff_cf govern dry matter assimilation under conditions in which transpiration demands cannot be met by the soil water supply.
APSIM-Sugarcane was used to make simulations for each experiment using the environmental conditions indicated in Figure 1 and Table 2, default cultivar specific parameters and management conditions indicated in Table 1. Prepared simulations for; experiment A was denoted as A SIM , experiment B1 was denoted as B1 SIM and experiment B2 was denoted as B2 SIM . Biomass and CDW at different DAP (indicated in Table 1: Observed data) were selected as the outputs of each simulation. These simulations were used in building emulators (described in Section 2.3.2), emulator accuracy evaluation (described in Section 2.3.3) and checking the goodness of the optimization process (described in Section 2.5).
When preparing each simulation, default cultivar specific parameter values of APSIM-Sugarcane were selected. However, these parameters were replaced with random parameter ensembles to build emulators (Section 2.3.2) and also were replaced with optimized parameter ensembles in validation (Section 2.5). Parameters of APSIM-Sugarcane other than those listed in Table 3 remained at default at each simulation.

Emulation
The emulation concept has drawn increasing attention [35,36] as a way to reduce the computational burden of dynamical simulator runs (e.g., APSIM-Sugarcane). An emulator m, which is either a metamodel or a surrogate model, can be used to approximate in a statistical sense the underlying simulator model (m) to reduce the cost of training runs. Whenever a simulation run is needed at a point (e.g., q), a fast prediction from a metamodel m(q) can be used to replace the costly value m(q).

Gaussian Process-Based Emulation
A GP is one of the statistical processes commonly used for emulation. Here we provide a brief overview of how to use GP for emulator generation. More comprehensive details of the GP can be found in Williams and Rasmussen [26] and Kennedy and O'Hagan [37].
If we consider the crop simulator (APSIM-Sugarcane) as an unknown function Y = f (X), where Y is the model output and X is a vector of (p) parameters (X = [x 1 , x 2 , . . . , x p ]), then Y is a random function in a Bayesian framework. According to Kennedy and O'Hagan [37], GP can be considered to be a flexible and convenient class of distributions that can be used to represent prior understanding about f (X). A prior distribution can therefore be assumed to be a multi-variate normal function that is characterized by a linear additive mean (Equation (1)) and a covariance function (Equation (2)). The covariance function is specified to characterize the smoothness of the output.
where m(X) indicates the linear additive mean function of X, X is a vector of p parameters (X = [x 1 , x 2 , . . . , x p ]) and β 0 , β 1 , . . . , β p are unknown coefficients. cov ( f (x), f (x) ) is the covariance function of any pair of joint probability distribution ( f (x), f (x) ), r i is a scaling parameter determining how rough the function is with respect to the ith input and σ 2 is the overall variance of the mean function.
The training data used to build the posterior distribution were used to estimate the β hyper-parameters, σ 2 parameter, and roughness parameters (r i ).
The emulator's posterior distribution was generated by using training data that could be obtained from the runs of the simulator. The training data consisted of the model outputs (Y) generated from each simulator run and the corresponding vectors of the p-dimensional parameters (X = [x 1 , x 2 , . . . , x p ]).

Building Emulators
At first, 500-parameter ensembles were prepared to train and validate emulators. Each ensemble was composed of 14 parameters as listed in Table 3. Parameter ensembles were prepared by using 500 random numbers which are distributed uniformly between the maximum and minimum values of each parameter space (Table 3: Parameter space) using functions in R [38] software. Reasonable values for the parameter spaces were selected based on the descriptions of previous studies of Bandara et al. [27] and Sexton et al. [22,29].
Then, the parameter ensembles were used to run APSIM-Sugarcane under single experiments A SIM , B1 SIM and B2 SIM (described in Section 2.2) producing 500 simulations per single experiment. Required scripts were created using the jsonlite [39], DBI [40] and RSQLite [41] packages of R.
Three hundred outputs of biomass and CDW at the specified DAP were used to train GP emulators and remaining 200 outputs were used for testing the accuracy of the emulators (Section 2.3.3).
The GP implemented in "Gaussian Emulation Machine for Sensitivity Analysis (GEM-SA)" [42] software was used to build the emulators. Fifty emulators were built to approximate sugarcane biomass and CDW under each experiment (A, B1 and B2) at each DAP: (for Experiment A, biomass at 7 specific DAPs: 7 emulators; Experiment A CDW at 7 specific DAPs: 7 emulators; Experiment B1 biomass at 9 specific DAPs: 9 emulators; Experiment B1 CDW at 9 specific DAPs: 9 emulators; Experiment B2 biomass at 9 specific DAPs: 9 emulators; Experiment B2 CDW at 9 specific DAPs: 9 emulators were built).
All the emulators were evaluated (Section 2.3.3) and used in the global optimization (Section 2.4) as a surrogate for the simulator.

Emulator Accuracy Evaluation
The accuracy of the built emulators was evaluated by using coefficient of determination (R 2 ), leave-one-out cross-validated root-mean-squared standardized error (CV RMSSE ) and sigma-squared value (σ 2 ). R 2 was calculated between the APSIM-simulated values and the emulator-predicted values. While developing each emulator, the remaining 200 parameter ensembles that were not used to generate emulators (described in Section 2.3.2) were used in the GEM-SA for obtaining the emulator-predicted outputs. Additionally, corresponding APSIM simulated outputs of the 200 parameter ensembles mentioned in Section 2.3.2 were used here as simulator predicted outputs. Then, both the emulator and simulator predictions were graphed to determine the R 2 between them. The emulators with R 2 values close to 1.0 were considered to be highly accurate. R 2 of emulator accuracy is hereafter denoted by R 2 emu . Both CV RMSSE and σ 2 were computed internally by the GEM-SA while building emulators with 300 training data points. To calculate CV RMSSE in GEM-SA, leave-oneout-cross-validation was used. In briefly, it fits the emulator by leaving one data point from the training data and, the missing point is predicted from the fitted emulator. This is repeated for all combinations of training data to provide overall effectiveness of the emulator as CV RMSSE . According to Qin et al. [43] and Sexton [44], the emulator variance is considered to have accurately estimated the actual error variance if the CV RMSSE is close to 1; overestimation and underestimation are indicated by lower and higher values, respectively.
Equation (3) defines CV RMSSE as follows: where "y i is the true output for the ith training run,ŷ is the corresponding emulator approximation, s i is the standard deviation calculated with the ith training point removed and n is the number of runs" [45]. According to Petropoulos et al. [46], the σ 2 value effectively measures the nonlinearity of an emulator via indicting the emulator variance after the output standardization. For a linear model, the σ 2 value is close to 0, whereas moderately to highly nonlinear models have greater σ 2 values (without a defined cutoff value).

Global Optimization
Global optimization seeks the best parameter ensemble which provides the best agreement between observed values and model-simulated outputs (global optima) in the presence of multiple local optima. Global optimization was conducted using the DE function implemented in the DEoptim [47] package of the R platform. Here, we provide a brief overview of the DE function; a detailed description can be found in Ardia et al. [47].
DE is among a class of genetic algorithms that have been inspired by biological processes. More details about genetic algorithms can be found in Mitchell [48]. Processes such as selection on a population, mutation, and crossover are used by genetic algorithms during optimization to minimize an objective function over the course of successive generations.
Like other evolutionary algorithms, DE evolves a population of parameter vectors to solve optimization problems. A single population is composed of a number of parameter vectors (NP). In many cases, NP should set at least 10 times the number of parameters used for optimization to avoid possible misconvergence in the optimization. The initial population is generated with values given by the user or with random values between lower and upper bounds that are defined by the user. The next successive generation of the population is then created using the parameter vectors of the current population by implementing differential mutation. To generate the first mutant parameter vector (v i ) of the next successive generation, three random parameter vectors of the existing population (x r0 , x r1 and x r2 ) are selected, and v i is generated by using v i where F is a differential weighting factor, effective values of which are typically between 0 and 1. After the first mutation, the remaining parameter vectors of successive generations are created by continuing the mutation with a crossover probability CR ∈ [0, 1]. The fraction of the parameter values that are copied from a mutant is controlled by the CR. During the process, all elements of the vector are created with respect to the defined lower and upper bounds of the parameters. The corresponding values of the objective function with respect to each trial vector are then determined. The previous vector in the population is replaced if the objective function value of a particular trial vector is equal to or lower than the previous vector; otherwise, the previous vector is retained.
The "DEoptim" function in the DEoptim package of R requires that several arguments be defined before the optimization process, including the objective function, the lower and upper bounds of the parameters and several control parameters (defined in Ardia et al. [47]). The root mean square error (RMSE) between the observed data and corresponding emulator predictions are defined as the objective function to be optimized (minimized) (Equation (4)). We ran the set of files produced by GEM-SA in R during each emulator generation to obtain emulator predictions. The R script, which can be modified to accomplish this task, can be found in Kennedy and Petropoulos [45]. The same lower and upper bounds used to generate emulators (Table 3: Parameter space) were defined as lower and upper bounds of the DEoptim function. The default control parameters described in Ardia et al. [47] were used for the optimization.
where O i is the value of the ith observation, Ep i is the emulator-predicted value of the ith observation, and n is the number of observations. The parameter vector that gave the lowest RMSE for the objective function after 1000 iterations from a population of 500 NP of the DEoptim function was obtained as the return result (a vector composed of 27 elements including sub levels of parameters [see Table 3: Code]). For a single cultivar, six optimized parameter ensembles (or vectors) were obtained from above process. Each of this parameter ensemble was optimized using biomass or CDW yield of the particular cultivar in the three field experiments (A, B1 and B2) [six different sets of observed data]. Emulator predictions given by those parameter ensembles were also recorded (Ep i ).

Validation-Step One
In the first step, optimized parameter ensembles were used to estimate the sugarcane yield under the same experimental condition (which was used to optimize the particular parameter ensemble) by using APSIM-Sugarcane. Then accuracy of each parameter ensemble was evaluated by comparing simulation predictions with observed data.
By running APSIM-Sugarcane using the parameter ensembles (described in Section 2.4), APSIM predictions (Sp i ) were obtained. Then, Ep i , Sp i , and O i were graphed together, and the normalized root mean square error (NRMSE), R 2 , and Willmott's agreement index (AI) were calculated to evaluate the accuracy of the optimized results. The R 2 was calculated from the linear regression between the observed and simulated values (Ep i and Sp i ). The NRMSE (Equation (5)) was equated to the RMSE divided by the output range and was reported as a percentage. The AI is a measure of non-parametric goodness of fit (Equation (6)), and the desired value is close to one.

Validation-Step Two
In the second step, optimized parameter ensembles were used to estimate the sugarcane yield under other experimental conditions (which were not used to optimize the particular parameter ensemble) by using APSIM-Sugarcane. Accuracy of each parameter ensemble was evaluated by calculating R 2 , NRMSE, and AI between the observed values and simulator predictions. Comparing all optimized parameter ensembles based on the validation results of step one and two, an optimal parameter ensemble was selected for each cultivar to parameterize APSIM-Sugarcane for estimating biomass and CDW. Validation process of step two was showed in Figure 2.

Emulator Accuracy
This section discusses the performances of the built emulators in terms of the calculated (Figure 3), the GEM-SA internally generated σ 2 (Figure 4), and the CVRMSSE values ( Figure 5).
The linear relationships between the APSIM-simulated values and the emulator-predicted values are indicated as scatter plots in Figure 3. Because it was difficult to show all the graphs, Figure 3 shows only the emulator performances of experiment B1. Results related to experiment A and B2 are indicated in Supplementary Materials: Figure S1 (results of experiment A), Figure S2 (results of experiment B2). Each scatter plot represents the results of two emulator simulations of biomass and CDW and the corresponding APSIM simulations.
With the exception of a few conditions, the calculated values all ranged be- As indicated in Figure 2, hereafter parameter ensembles optimized from the observed data of; experiment A biomass was denoted as P1, experiment A CDW was denoted as P2, experiment B1 biomass was denoted as P3, experiment B2 CDW was denoted as P4, experiment B2 biomass was denoted as P5 and experiment B2 CDW was denoted as P6 of each cultivar. Further calculated NRMSE, R 2 and AI of second validation step are hereafter denoted by NRMSE val , R 2 val and AI val , respectively. Results were depicted using box diagrams and dot plots under Section 3.2.2. Finally, the selected optimal parameter ensembles were compared.

Emulator Accuracy
This section discusses the performances of the built emulators in terms of the calculated R 2 emu (Figure 3), the GEM-SA internally generated σ 2 (Figure 4), and the CV RMSSE values ( Figure 5).

11, x FOR PEER REVIEW 12 of 24
predictions make by the emulator will have understated uncertainty. The new output values will then be further from the emulator approximations than the emulator expects. We therefore further analyzed σ 2 and CVRMSSE values to assess the performances of the emulators.    Agronomy 2021, 11, x FOR PEER REVIEW 12 predictions make by the emulator will have understated uncertainty. The new outpu ues will then be further from the emulator approximations than the emulator expect therefore further analyzed σ 2 and CVRMSSE values to assess the performances of the em tors.     predictions make by the emulator will have understated uncertainty. The new outpu ues will then be further from the emulator approximations than the emulator expect therefore further analyzed σ 2 and CVRMSSE values to assess the performances of the em tors.            The σ 2 values ranged between 0.06 and 1.64 for all emulators (Figure 4). We observed relatively high σ 2 values for the early stages of CDW (e.g., 96_DAP of experiment A and 99_DAP of experiment B1). This pattern was consistent with the observations of low values under those conditions. The computed values of CVRMSSE of the emulators were ranged between 0.74 and 1.01 ( Figure 5).
Petropoulos et al. [46] have reported that the σ 2 values for their emulators ranged between 0.13 and 1.6, Qin et al. [43] have reported that the σ 2 values for their emulators ranged between 0.6 and 2.1 and the parameters they used showed only moderate deviation from the linearity. Gunarathna et al. [30] have indicated that their emulators showed good to moderate linearity with σ 2 values that ranged from 0.10 to 1.43. It can hence be concluded that good linearity was shown by our emulators under all considered environmental and management condition.
In comparison with the CVRMSSE values obtained by Gunarathna et al. [30], Kennedy and Petropoulos [45] and Petropoulos et al. [46] the values we obtained were lower, and the fact that they were close to 1.0 in all the experiments suggested that the built emulators could represent the true model well.

Validation-Step One
This section evaluated the accuracy of the optimized results with respect to the observed data which was used to optimize the particular parameter ensemble. The parame- The linear relationships between the APSIM-simulated values and the emulatorpredicted values are indicated as scatter plots in Figure 3. Because it was difficult to show all the graphs, Figure 3 shows only the emulator performances of experiment B1. Results related to experiment A and B2 are indicated in Supplementary Materials: Figure S1 (results of experiment A), Figure S2 (results of experiment B2). Each scatter plot represents the results of two emulator simulations of biomass and CDW and the corresponding APSIM simulations.
With the exception of a few conditions, the calculated R 2 emu values all ranged between 0.9 and 1. The indication was that all the emulators could approximate the APSIM simulators successfully. However, the R 2 emu of CDW in 96_DAP of experiment A and 99_DAP of experiment B1 and B2 were only 0.43, 0.47 ( Figure 3) and 0.3, respectively. In these cases, results were noisy and highly variable. This noise reflects the fact that around 96_DAP and 99_DAP, the sugarcane plant in the APSIM model is at the "emergence" stage, where emergence is the beginning of cane growth. Although the CDW (sucrose weight and structural stem weight) is very low, it is highly variable (high coefficient of variation) at the emergence stage. As a result, the smoothness of the output of the GP, which is characterized by its covariance function (Equation (2)), can be high, therefore predictions make by the emulator will have understated uncertainty. The new output values will then be further from the emulator approximations than the emulator expects. We therefore further analyzed σ 2 and CV RMSSE values to assess the performances of the emulators.
The σ 2 values ranged between 0.06 and 1.64 for all emulators (Figure 4). We observed relatively high σ 2 values for the early stages of CDW (e.g., 96_DAP of experiment A and 99_DAP of experiment B1). This pattern was consistent with the observations of low R 2 emu values under those conditions. The computed values of CV RMSSE of the emulators were ranged between 0.74 and 1.01 ( Figure 5).
Petropoulos et al. [46] have reported that the σ 2 values for their emulators ranged between 0.13 and 1.6, Qin et al. [43] have reported that the σ 2 values for their emulators ranged between 0.6 and 2.1 and the parameters they used showed only moderate deviation from the linearity. Gunarathna et al. [30] have indicated that their emulators showed good to moderate linearity with σ 2 values that ranged from 0.10 to 1.43. It can hence be concluded that good linearity was shown by our emulators under all considered environmental and management condition.
In comparison with the CV RMSSE values obtained by Gunarathna et al. [30], Kennedy and Petropoulos [45] and Petropoulos et al. [46] the values we obtained were lower, and the fact that they were close to 1.0 in all the experiments suggested that the built emulators could represent the true model well.

Validation-Step One
This section evaluated the accuracy of the optimized results with respect to the observed data which was used to optimize the particular parameter ensemble. The parameter ensembles that gave the lowest RMSE values for the objective function (described in Section 2.4) are listed in Supplementary Materials: Table S2 (results for cultivar KK3), Table  S3 (results for cultivar LK92-11), and Table S4 (results for cultivar 02-2-058). Figure 6 compares those parameter ensembles based on simulated (Emulator and APSIM) sugarcane yields (Biomass and CDW) and observed sugarcane yields for cultivar KK3. Because of the difficulty of showing all the graphs, this manuscript presents only few examples to represent all conditions. We evaluated the parameter ensembles ( Figure 6) by comparing the R 2 opt and NRMSE opt percentages and the AI opt values between sugarcane yields simulations obtained with the optimized parameter ensembles and observed sugarcane yields.
The calculated R 2 opt and AI opt values and the NRMSE opt percentages fell in the ranges 0.95-1, 0.97-1 and 1-11.32%, respectively. The indication was that all parameter ensembles obtained from the optimization could probably be used to approximate observed values using APSIM.
As expected, the results of the APSIM using optimized parameters were less accurate than the emulator results ( Figure 6). This is because parameters were first obtained from the global optimization based on emulators and then used with the APSIM-Sugarcane simulations to obtain predictions for accuracy evaluation. However, as indicated by the above results, the APSIM using optimized parameters could still be used to simulate observed values in all cases with acceptable accuracy. The reason is that the built emulators indicated higher accuracy when approximating the simulator (Section 3.1).
Moreover, a GP-based optimization method is very efficient in terms of computational time. Because during our study we could observe that single simulation of APSIM-Sugarcane requires a CPU time of 1.46 ± 0.006 s (CPU had a Quad-core with 1.60 GHz clock speed and 6 Mb L3 cache with 7.86 GB usable RAM) while emulators require 0.02 ± 0.005 s to provide similar output (calculated by using proc.time() function of R). Therefore, emulators with sufficient accuracy can be used to reduce the computational burden of the process which often requires large number of simulator runs. For instance, Sexston et al. [22] used 30,000 simulator runs in calibration of varietal parameters in APSIM-Sugarcane and Sheng et al. [21] used 50,000 simulator runs in calibration of varietal parameters in APSIM-Maize. However, several simulator runs are required to build the emulators and for their accuracy evaluation. However, the number of simulator runs required to build an accurate emulator is less than the number of simulator runs required in case of only simulator is used for optimization. In our study, with 300 simulator runs, we could build sufficiently accurate emulators (additional 200 simulator runs were used for accuracy evaluation). Then emulators used for the optimizations with 1000 iterations each composed with 500 NP of emulator runs. More details about relative advantages of GP emulation could be find in Oakley and O'Hagan [49] and Kennedy et al. [50]. Our results therefore confirmed that the optimization method with the emulator could be used for improving the efficiency of model development.

Validation-Step Two
This section discusses the accuracy of the parameter ensembles obtained from optimization for use in cultivar trait parameterization of the APSIM-Sugarcane model under different environment and management conditions (see Figure 2). Comparisons between the observed sugarcane yields and the APSIM-predicted sugarcane yields are shown in Figures 7 and 8.
Calculated values ranged between 0.925 and 1.0 for all conditions (Figure 7). Based on the values, all parameter vectors indicated best performances for estimating observed yields. However, we observed large variations of the NRMSEval and AIval values among cultivars and parameter ensembles (Figure 8).

Validation-Step Two
This section discusses the accuracy of the parameter ensembles obtained from optimization for use in cultivar trait parameterization of the APSIM-Sugarcane model under different environment and management conditions (see Figure 2). Comparisons between the observed sugarcane yields and the APSIM-predicted sugarcane yields are shown in Figures 7 and 8. 11 and 02-2-058. All NRMSEval percentages were dispersed closer to 0 than other parameter ensembles (Figure 8). Moreover, those parameter ensembles corresponded to the lowest AIval values, and the data points were dispersed closer to 1 than the data points of the other parameter ensembles of the respective cultivars ( Figure 8). values calculated between observed and simulated yields for six different cases (A_Biomass, A_CDW, B1_Biomass, B1_CDW, B2_Biomass, and B2_CDW) of a single cultivar (see Figure 2 for more details). The median is indicated by thick black lines, the interquartile range (IQR) is indicated by the boxes, 1.5 times the IQR is indicated by the whiskers, and the outliers beyond 1.5 times the IQR are indicated by points in black color. val values calculated between observed and simulated yields for six different cases (A_Biomass, A_CDW, B1_Biomass, B1_CDW, B2_Biomass, and B2_CDW) of a single cultivar (see Figure 2 for more details). The median is indicated by thick black lines, the interquartile range (IQR) is indicated by the boxes, 1.5 times the IQR is indicated by the whiskers, and the outliers beyond 1.5 times the IQR are indicated by points in black color. Optimized parameters under high water stress conditions resulted in the lowest performances when estimating the sugarcane yields under the lowest water stress conditions. This pattern can be clearly observed in Figure 8 (Table 4). They could therefore be selected as the best parameter ensemble to parametrize the APSIM-Sugarcane model for cultivars KK3, LK92-11 and 02-2-058. All NRMSE val percentages were dispersed closer to 0 than other parameter ensembles (Figure 8). Moreover, those parameter ensembles corresponded to the lowest AI val values, and the data points were dispersed closer to 1 than the data points of the other parameter ensembles of the respective cultivars ( Figure 8). Optimized parameters under high water stress conditions resulted in the lowest performances when estimating the sugarcane yields under the lowest water stress conditions. This pattern can be clearly observed in Figure 8 (NRMSE val %), where optimized parameter ensembles (P1 and P2) for experiment A resulted in comparatively high NRMSE val percentages when estimating the observed sugarcane yield of experiment B1. To explain this pattern, we simulated the "soil water deficit factor for photosynthesis (swdef_photo)" in the APSIM for each experiment. For instance, the swdef_photo and the observed sugarcane biomass yields of experiments A, B1 and B2 of cultivar KK3 are indicated in Figure 9. We observed that water stress conditions were most apparent in experiment A (swdef_photo near 0), than experiment B1 and B2. Due to that the lowest sugarcane yield was also observed in experiment A. Therefore, when the parameters were optimized in experiment A, the parameters were estimated to result in lower yields than in experiment B1 and B2.
Use of parameters estimated to result in low yields under severe water stress conditions may lead to poor performances under low water stress conditions. However, parameters optimized under the lowest water stress conditions were found to provide better estimates of sugarcane yields under the highest water stress conditions. This pattern is apparent in the NRMSEval% in Figure 8. The optimized parameter ensembles for experiment B1 (P3 and P4) resulted in comparatively low NRMSEval values in the estimates of the observed sugarcane yields for experiment A. Because in APSIM-Sugarcane, parameter rue is previously identified as a highly sensitive parameter for the estimation of biomass and CDW of sugarcane [27][28][29][30]. Therefore, when optimizing parameters under low water stress conditions (B1), the estimated value of rue could be increased to estimate higher yields. When those optimized parameter ensembles are used for simulations under severe water stress conditions (A), APSIM-Sugarcane will reduce the rue because in APSIM-Sugarcane, rue tends to be reduced whenever a soil water shortage condition is met [9,51].

Comparison of Optimized Parameters
This section discusses the differences between the three selected parameter ensembles (listed in Table 4) for varietal parameterization (Cultivar KK3, LK92-11 and 02-2-058) of the APSIM-Sugarcane.
Cultivars KK3 and 02-2-058 evidenced a higher rue parameter than cultivar LK92-11. This difference is apparent in the best estimated parameter ensembles listed in Table 4. The RUE3 and RUE4 values were lower for cultivar LK92-11 than for cultivars KK3 and 02-2-058. This difference reflects the fact that rue is a highly sensitive parameter for estimation of sugarcane biomass [29]. However, as explained previously, the value of rue may be lower when there is a soil water shortage. The lower values of rue (RUE3: rue of growth stage 3 [from emergence to the beginning of cane growth] and RUE4: rue of growth stage However, parameters optimized under the lowest water stress conditions were found to provide better estimates of sugarcane yields under the highest water stress conditions. This pattern is apparent in the NRMSE val % in Figure 8. The optimized parameter ensembles for experiment B1 (P3 and P4) resulted in comparatively low NRMSE val values in the estimates of the observed sugarcane yields for experiment A. Because in APSIM-Sugarcane, parameter rue is previously identified as a highly sensitive parameter for the estimation of biomass and CDW of sugarcane [27][28][29][30]. Therefore, when optimizing parameters under low water stress conditions (B1), the estimated value of rue could be increased to estimate higher yields. When those optimized parameter ensembles are used for simulations under severe water stress conditions (A), APSIM-Sugarcane will reduce the rue because in APSIM-Sugarcane, rue tends to be reduced whenever a soil water shortage condition is met [9,51].

Comparison of Optimized Parameters
This section discusses the differences between the three selected parameter ensembles (listed in Table 4) for varietal parameterization (Cultivar KK3, LK92-11 and 02-2-058) of the APSIM-Sugarcane.
Cultivars KK3 and 02-2-058 evidenced a higher rue parameter than cultivar LK92-11. This difference is apparent in the best estimated parameter ensembles listed in Table 4. The RUE3 and RUE4 values were lower for cultivar LK92-11 than for cultivars KK3 and 02-2-058. This difference reflects the fact that rue is a highly sensitive parameter for estimation of sugarcane biomass [29]. However, as explained previously, the value of rue may be lower when there is a soil water shortage. The lower values of rue (RUE3: rue of growth stage 3 [from emergence to the beginning of cane growth] and RUE4: rue of growth stage 4 [from emergence to the beginning of cane growth]) for cultivar LK92-11 therefore reflected the fact that the selected parameter for the cultivar was optimized under Rf (water stressed) conditions (Experiment B2). However, we observed a higher value for the RUE5 of cultivar LK92-11 (Table 4) compared to other cultivars. This difference could reflect the fact that in APSIM, growth stage 5 (from flowering to the end of the crop) is currently deactivated because of a lack of good physiological information on which to base predictions [9]. The rue of growth stage 5 (RUE5) thus has no influence on the determination of biomass yield.
Cultivars KK3 and 02-2-058 evidenced higher transp_eff_cf values for growth stage 4 (TEC4) than cultivar LK92-11. To explain this difference, we simulated the growth stages of APSIM-Sugarcane using the selected parameter ensembles (Table 4) under corresponding experimental conditions. Figure 10 shows the results of the simulated and observed yields of each cultivar under the same conditions. Growth stages of 1, 2, 5 and 6 were not included in the simulation in Figure 10 because the observed values were reported for DAP 99-390. It is obvious that transp_eff_cf of growth stages 1 (TEC1), 2 (TEC2), 5 (TEC5) and 6 (TEC6) would have less influence on the estimation of sugarcane biomass. However, transp_eff_cf has previously been identified as a highly sensitive parameter for estimating biomass yields under water stress conditions [29]. Even though the parameter ensemble of LK92-11 was estimated under Rf conditions, growth stage 3 ( Figure 10: DAP 99-128) was not affected by the water stress ( Figure 9: swdef_photo of DAP 99-128). In this case, assimilation dry matter was governed mainly by the RUE3. Similar values for TEC3 were thus observed for each cultivar. However, the TEC4 of cultivar LK92-11 was lower because during growth stage 4, water stress conditions were reported (Figure 9: swdef_photo of DAP 128-390). In this case, assimilation of dry matter in APSIM is governed by transp_eff_cf, water supply, and the vapor pressure deficit. As a result, a lower value for TEC 4 was reported in accord with the lower yield observed for LK92-11 compared to other cultivars ( Figure 10).
However, we observed better drought resistant characteristic of cultivar LK92-11 than the other cultivars. Because the parameter ensemble obtained for cultivar LK92-11 (listed in Table 4) underestimated the observed yields under high water stress conditions (results of estimates: Exp. A_ Biomass = NRMSE val : 18.14%, Exp. A_ CDW = NRMSE val : 21.61% [ Figure 8: cultivar LK92-11, NRMSE val % of P5]). However, parameter ensembles of the remaining cultivars indicated satisfactory results under the same conditions (Figure 8, cultivars KK3 and 02-2-058: P3). Our results were similar to the results obtained by Preecha et al. [6] using the same field experiments and are consistent with Peerasak [52], who reported that LK92-11 was less sensitive to water shortage than KK3. Cha-um et al. [53] have also indicated that LK92-11 is tolerant to water deficit.
Cultivar KK3 evidenced rapid sugar accumulation compared to the other cultivars under the same environmental and management conditions. This pattern was observed when simulating each cultivar under each experimental condition (A, B1 and B2) in APSIM-Sugarcane to obtain sucrose weight. This pattern is caused by the parameters that govern sucrose accumulation (cane_fraction, stress_factor_stalk, sucrose_fraction_stalk, sucrose_delay, min_sstem_sucrose_redn, and min_sstem_sucrose). The optimized KK3 parameters facilitated early sucrose accumulation compared to other cultivars (Table 4). For instance, the highest value for the sucrose_fraction_stalk and the lowest value for the min_sstem_sucrose were observed for KK3 (Table 4). Moreover, Khonghinta et al. [54] have conducted field experiments to compare KK3 with several other sugarcane cultivars in KK and have also indicated that KK3 accumulates sugar rapidly. y 2021, 11, x FOR PEER REVIEW 21 of 24 Figure 10. Observed biomass weight (gm −2 ) of selected parameter ensembles of each cultivar under optimized experimental conditions. Cultivars KK3 and 02-2-058: observed biomass yield in experiment B1 and cultivar LK92-11: observed biomass yield in experiment B2. Growth stages were obtained from the APSIM-Sugarcane simulations using the corresponding parameter ensemble and APSIM simulation files.
Cultivar KK3 evidenced rapid sugar accumulation compared to the other cultivars under the same environmental and management conditions. This pattern was observed when simulating each cultivar under each experimental condition (A, B1 and B2) in APSIM-Sugarcane to obtain sucrose weight. This pattern is caused by the parameters that govern sucrose accumulation (cane_fraction, stress_factor_stalk, sucrose_fraction_stalk, su-crose_delay, min_sstem_sucrose_redn, and min_sstem_sucrose). The optimized KK3 parameters facilitated early sucrose accumulation compared to other cultivars (Table 4). For instance, the highest value for the sucrose_fraction_stalk and the lowest value for the min_sstem_sucrose were observed for KK3 (Table 4). Moreover, Khonghinta et al. [54] have conducted field experiments to compare KK3 with several other sugarcane cultivars in KK and have also indicated that KK3 accumulates sugar rapidly.
All cultivars evidenced similar values (p < 0.05) for leaf area index (LAI) and phenological stages when simulated under similar environmental conditions in APSIM-Sugarcane. Although we observed higher values for leaf development parameters (green_leaf_no, average leaf_size, average tiller leaf size) of cultivar LK92-11 than of cultivars KK3 and 02-2-058 (Table 4), we observed similar LAI values because cultivar LK92-11 evidenced a lower rue than the other cultivars. In APSIM-Sugarcane, leaf area growth can be limited by the amount of biomass partitioned to the leaf. The amount of biomass that is produced from the conversion of intercepted radiation is largely governed by the rue.
We did not observe much difference between the parameters that control phenological development based on thermal time (tt_emerge_to_begcane, tt_begcane_to_flowering, and tt_flowering_to_crop_end) in APSIM-Sugarcane. Moreover, in APSIM-Sugarcane tt_flowering_to_crop_end is currently deactivated because of the absence of a good physiological basis for prediction. All cultivars evidenced similar values (p < 0.05) for leaf area index (LAI) and phenological stages when simulated under similar environmental conditions in APSIM-Sugarcane. Although we observed higher values for leaf development parameters (green_leaf_no, average leaf_size, average tiller leaf size) of cultivar LK92-11 than of cultivars KK3 and 02-2-058 (Table 4), we observed similar LAI values because cultivar LK92-11 evidenced a lower rue than the other cultivars. In APSIM-Sugarcane, leaf area growth can be limited by the amount of biomass partitioned to the leaf. The amount of biomass that is produced from the conversion of intercepted radiation is largely governed by the rue.
We did not observe much difference between the parameters that control phenological development based on thermal time (tt_emerge_to_begcane, tt_begcane_to_flowering, and tt_flowering_to_crop_end) in APSIM-Sugarcane. Moreover, in APSIM-Sugarcane tt_flowering_ to_crop_end is currently deactivated because of the absence of a good physiological basis for prediction.

Conclusions
This study was aimed on using the GP-based emulators to optimize cultivar trait parameters in the APSIM-Sugarcane model. According to the obtained values for R 2 , σ 2 and CV RMSSE values, the emulators we built for the optimization showed satisfactory results. The indication is that these emulators can approximate the original simulator (APSIM-Sugarcane) successfully. Via the GP-based emulator optimization, we could obtain acceptable parameter ensembles for parametrization of Thai cultivars KK3, LK92-11 and 02-2-058 by using the APSIM-Sugarcane model. The optimized parameters evidenced satisfactory results during the validation under the environmental and management conditions found in KK, Thailand. Based on our validation results, we suggest that GP emulation can be efficiently implemented for parameterization of computationally expensive simulators. Future studies will be needed to reach more robust conclusions concerning the use of emulation for parameter optimization with APSIM-Sugarcane.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/agronomy11071379/s1, Figure S1: Relationship between APSIM-simulated values and emulator-predicted values of biomass and CDW in experiment A at each reporting frequency at days after planting (DAP): 96, 117, 147, 173, 244, 293 and 388; Figure S2: Relationship between APSIM-simulated values and emulator-predicted values of biomass and CDW in experiment B2 at each reporting frequency at days after planting (DAP): 99, 128, 185, 238, 267, 299, 329, 360 and 390; Table S1: Summary of the data collected during field experiments; Table S2: Best estimated parameter ensembles obtained for biomass and CDW of cultivar KK3 based on the emulators of experiments A, B1 and B2; Table S3: Best estimated parameter ensembles obtained for the biomass and CDW of cultivar LK92-11 from the emulators of experiments A, B1 and B2; Table S4: Best estimated parameter ensembles obtained for the biomass and CDW of cultivar 02-2-058 from the emulators of experiments A, B1 and B2.