Integration of Remote Sensing and Field Observations in Evaluating DSSAT Model for Estimating Maize and Soybean Growth and Yield in Maryland, USA

: Crop models are useful for evaluating crop growth and yield at the ﬁeld and regional scales


Introduction
Agriculture is one of the major economic sectors in the United States of America (U.S.), contributing approximately $992 billion in revenue to the U.S. economy annually and about 6% of the gross domestic product (GDP) [1,2]. As an example, during the third quarter of The objectives of this study were to (1) derive LAI and crop growth stages using remote sensing observation, and (2) use remotely sensed data along with field observations to calibrate the DSSAT model for soybean and maize growth and yield in an experimental agricultural area. The calibrated model can be scaled up for regional yield prediction and climate change impact study.

Description of the Study and Experimental Area
Our study site was carried out at the United State Department of Agriculture, Agricultural Research Service, Beltsville Agricultural Research Center (BARC), Beltsville, Maryland ( Figure 1). This site is a Long-Term Agroecosystem Research (LTAR) network and production research site in Maryland and includes~2670 hectares of agricultural fields. The sizes of the fields vary from 1 hectare to more than 20 hectares [47]. No-till practices were used for all fields in this study. The major crops grown in these fields are maize, soybean, wheat, barley, rye, and cover crops. The dominant soil types range from sandy loams to loamy sands. Climatological conditions at BARC are characteristic of the coastal physiographic zone of the mid-Atlantic Region. The research center keeps a comprehensive record of all planting operations, such as pre-planting operations, planting dates, cultivar name, planting population and row spacing, irrigation and fertilizer applications, and harvest dates for each field using FarmLogic Systems. In this study, for maize, the cultivar used was a medium season cultivar, whereas for soybean, maturity group three and group four were used. The days to flowering and maturity for maturity group three ranged from 37 to 52 and from 93 to 128 days after planting, respectively. For maturity group 4, the days to flowering and maturity ranged from 43 to 50 and from 106 to 126 days after planting, respectively. There were no clear differences in days to flowering and maturity between the maturity groups. This may have been a result of the differences in fields and years used. For the maize cultivar, plants were no-till seeded between 38 cm and 76 cm rows at a density of between 69,000 and 74,000 plants ha −1 from 25 April to 10 June. Soybean crops were no-till and seeded between 19cm and 76cm rows at a density of between 383,000 and 445,000 plants ha −1 from 20 April to 28 June. Both cultivars were grown under rainfed conditions. All of the above planting and management systems were used in our DSSAT for the calibration and evaluation of our model.
Grain yield measurements were obtained from yield monitor data acquired from the production fields at the BARC located in Beltsville, Maryland ( Figure 1). For our research purposes, we set up criteria for selecting field plots based on field size and shape. Fields of a relatively large size (i.e., 5.0 ha or larger) allowed for the use of combines outfitted with yield monitors for the collection of harvest data (See [48] for more details). In addition, some fields had installed permanent and mobile PhenoCams that were used for phenology estimations, because in situ phenology observations were rarely collected (see Figure 2). The derived phenology for a few of these fields were used as a reference to derive our phenology from the remote sensing observations as described in Section 2.2. Ground LAI measurements were made in 2018 and 2020 using LI-COR LAI-2200 Plant Canopy AnalyzerPlant in one of these fields (#6-60).    Figure 1a).

Leaf Area Index
Leaf area index (LAI) was retrieved using a reference-based approach [49]. We used Sentinel-2 LAI as a reference to generate PlanetScope-based LAI at high spatial and temporal resolutions. The Sentinel-2 constellation includes 2A and 2B satellites which acquire optical imagery at 10-60 m spatial resolution collectively every 5 days [50]. The Plan-etScope constellation consists of over 200 CubeSats in low earth orbits that can capture earth imagery daily at a 3 m resolution [51]. In order to generate daily LAI at a high spatial resolution, we combined Sentinel-2 LAI and Planet Fusion surface reflectance data using a machine learning approach, as described below.
Planet Fusion [54] constitutes a comprehensive fusion methodology based on the Cubesat Enabled Spatio-Temporal Enhancement Method (CESTEM) algorithm [55] to intercalibrate, harmonize, and fuse multi-sensor data streams leveraging rigorously calibrated public mission (i.e., Sentinel-2, Landsat 8/9, Moderate Resolution Imaging Spectroradiometer (MODIS), Visible Infrared Imaging Radiometer Suite (VIIRS)) products in concert with the higher-spatial-and-temporal-resolution CubeSat images from the PlanetScope constellations. The Framework for Operational Radiometric Correction for Environmental Monitoring [56] generates a harmonized Sentinel-2 and Landsat 8/9 Bidirectional Reflectance Distribution Function (BRDF)-adjusted surface reflectance (SR) product that is used as the cross-calibration reference target during the radiometric harmonization step. Planet  Figure 1a).

Leaf Area Index
Leaf area index (LAI) was retrieved using a reference-based approach [49]. We used Sentinel-2 LAI as a reference to generate PlanetScope-based LAI at high spatial and temporal resolutions. The Sentinel-2 constellation includes 2A and 2B satellites which acquire optical imagery at 10-60 m spatial resolution collectively every 5 days [50]. The Plan-etScope constellation consists of over 200 CubeSats in low earth orbits that can capture earth imagery daily at a 3 m resolution [51]. In order to generate daily LAI at a high spatial resolution, we combined Sentinel-2 LAI and Planet Fusion surface reflectance data using a machine learning approach, as described below.
Planet Fusion [54] constitutes a comprehensive fusion methodology based on the Cubesat Enabled Spatio-Temporal Enhancement Method (CESTEM) algorithm [55] to intercalibrate, harmonize, and fuse multi-sensor data streams leveraging rigorously calibrated public mission (i.e., Sentinel-2, Landsat 8/9, Moderate Resolution Imaging Spectroradiometer (MODIS), Visible Infrared Imaging Radiometer Suite (VIIRS)) products in concert with the higher-spatial-and-temporal-resolution CubeSat images from the PlanetScope constellations. The Framework for Operational Radiometric Correction for Environmental Monitoring [56] generates a harmonized Sentinel-2 and Landsat 8/9 Bidirectional Reflectance Distribution Function (BRDF)-adjusted surface reflectance (SR) product that is used as the cross-calibration reference target during the radiometric harmonization step.
Planet Fusion imagery acquired on seven cloud-free reference Sentinel-2 dates were resampled and used to build a combined training dataset at 10 m resolution. A data-driven machine learning approach (Cubist regression tree) was used to generate an LAI model by training the Sentinel-2 LAI and Planet Fusion surface reflectance at 10 m resolution. The trained regression tree model was applied to Planet Fusion data at the original 3 m resolution. The daily PlanetScope-based LAI was averaged for each field used in this study to generate field-level LAI from 2017 to 2021.

Crop Growth Stages
Crop growth stages were estimated using PhenoCam observations and Planet Fusion data. For the fields that contain installed PhenoCams, crop growth stages were visually interpreted using daily PhenoCam photos. This includes six maize and four soybean site years. Because some crop growth stages (especially reproductive stages) are hard to identify and differentiate using PhenoCam photos, the identified maize growth stages were limited to planting, tasseling (VT), silking (R1), maturity (R6), and harvest. For soybean growth stages planting, beginning bloom (R1), beginning maturity (R7), full maturity (R8), and harvest were identified based on the PhenoCam and other available digital images. These fields served as reference sites for the rest of the fields and years.
For target fields that did not contain PhenoCams, the two-band enhanced vegetation index (EVI2) computed from Planet Fusion data was used to estimate crop growth stages. The 2-band EVI2 approximates the 3-band enhanced vegetation index (EVI). EVI2 can be computed using surface reflectance from red and near-infrared (NIR) bands as: We associated crop growth stages observed from PhenoCam reference sites with the EVI2 time series from the same field and year. Then, we compared the EVI2 time series from each target field to the reference sites using a dynamic time warping (DTW) technique [57]. The DTW measures the similarity between two non-linear time series. The reference site with the best match (smallest DTW distance) to the target site was selected. The dates of crop growth stages were retrieved using the warping path. For the fields without planting records, we used the within-season emergence (WISE) algorithm [47] to estimate crop emergence dates and then extrapolated planting dates utilizing the relationship between planting and emergence dates from reference sites.

DSSAT Model Description
The Decision Support System for Agro-technology Transfer (DSSAT) is a process-based crop model, evolved from the International Benchmark Sites Network for Agrotechnology Transfer Project [58][59][60][61]. DSSAT version 4.7.5 was used for this study. In order to simulate crop growth, development, and yield, the model requires a defined minimum dataset on weather, soil profiles, and crop management. The minimum weather dataset required to run the model includes maximum and minimum temperatures, solar radiation, wind speed, rainfall, and relative humidity. The soil parameters needed for the model include soil depth, and physical and chemical characteristics (i.e., drained lower limit, drained upper limit, saturation, and saturated hydraulic conductivity). In this study, soil profiles were obtained from the Soil Survey Geographic Database (SSURGO) (https://www.nrcs. usda.gov/resources/data-and-reports/soil-survey-geographic-database-ssurgo, accessed  Supplementary Table S1).
The model also requires crop management information such as crop type, cultivar, planting date, row and plant spacing, fertilizer application, tillage practices, and information on organic amendments for the fields or region [58][59][60][61]. In this study, two of the required crop variables (phenology and LAI) were not measured in most fields but were derived from the remotely sensed observations and subsequently used for model calibration and validation (see Section 2.2). Default cultivars were selected for the initial coefficient calibration based on the common cultivars planted or grown in the study area (see Tables 1 and 2). Daily weather data were obtained from the USDA Soil Climate Analysis Network (SCAN) station located at BARC, nearby to our study area (Powder Mill (2049), https:// wcc.sc.egov.usda.gov/nwcc/site?sitenum=2049, accessed on 16 September 2021)). Weather variables used included daily maximum and minimum temperatures, solar radiation, wind speed, rainfall, and relative humidity. The data were quality controlled (i.e., missing values were filled with values from nearby stations or other data sources, and outliers were removed) before the DSSAT weather file format was created.

Model Calibration
To apply the DSSAT model in a particular location, one typically has to first run the model over a time period before the target application to calibrate key model parameters, such that it will reproduce, on average, the measured or reference growth and yield in most situations. One such key parameter set includes the cultivar coefficients which determine crop growth, development, and yield. To calibrate the model cultivar coefficients, crop information such as phenology, LAI, and yield are critical. This is because the calibration process entails fitting the model to measured or reference variables and where these variables are not measured, the model cannot be calibrated or adapted to the local agronomic conditions. Therefore, we need to estimate the set of parameter values that can accurately reproduce measured growth and yield. In our study, we used 17 treatments (site years) for maize (10 treatments as the calibration dataset, and seven treatments as the validation dataset). For soybean (maturity groups three and four), a total of 18 treatments (site years) were used (nine treatments as the calibration dataset, and nine treatments as the validation dataset) from the 2017 to 2020 period.  Table 2  For the calibration of cultivar coefficients, it is known that if water, nutrients or other stresses affect experiments or trials (the crop data used for calibration), the estimate of cultivar coefficients may be poor [58][59][60][61]. Thus, we carefully selected our calibration treatments to include a range of weather conditions, day lengths, and crops growing in a minimumstress or stress-free environment for the best results. Note that using data from multiple years/sites generates coefficients that are more robust. For this study, the treatments selected for calibration were based on the best selection criteria described above. We used the data described from Sections 2.1-2.3 as input to our model for calibration, including a combination of the field observations and remote sensing-derived parameters. We used all available crop management (seeding rates, fertilizer and irrigation applications, etc.), soil physical characteristics, and climate information (maximum and minimum temperature, precipitation, solar radiation, wind and relative humidity) data as model inputs [62].
We used the Generalized Likelihood Uncertainty Estimation (GLUE) tool in the DSSAT model for our initial calibration of cultivar coefficient for both maize and soybean cultivars. The GLUE procedures required two passes: the first pass estimates phenological parameters, and the second pass estimates growth and yield parameters. The remotely sensed LAI and phenological variables were key in our cultivar coefficients estimation, supporting model calibration and evaluation even though these variables were not measured in the field. The derived remote sensing and field variables were used for fitting the model to measured or reference variables. The cultivar coefficients were carefully calibrated with the field and remotely sensed observable crop variables such as phenology, LAI, and yield. We used more than 6000 runs to generate random sets of parameters based on prior distribution of the known genetic coefficients to generate our calibrated cultivar coefficients [63]. We then manually refined our results where necessary based on DSSAT calibration criteria in order to simulate the measured or reference results of the maize and soybean for our study [58,63]. The calibrated genetic cultivar coefficients for maize and soybean are listed in Tables 1 and 2.

Model Evaluation
The agreement between the observed and simulated crop parameters was analyzed using the coefficient of determination (R 2 ), which describes the proportion of the total variance explained by the model, as well as Willmott's statistics measures [64,65] which include the index of agreement (d), which measures the degree to which simulated values match the observed values, the root mean squared error (RMSE), which is the sum of the differences between simulated and the observed values, and the normalized root mean squared error (NRMSE), which is a measure of the relative difference (%) between modelsimulated and observed values (see the details of the index descriptions, equations and rating in [62]. We also used percentage difference (PD) for our analysis of the model, calculated as: where S i and O i are the simulated and observed or reference values of a given variable, respectively. Lower PD indicates better model performance. Figure 3 shows the map of the PlanetScope-based LAI on 1 August 2020, and the time series of LAI for the field 6-60 from 2018 to 2020. The remote-sensing-retrieved LAI at the field level was assessed using in situ LAI measurements collected in 2018 and 2020. Figure 3 shows that remote-sensing-retrieved LAI (field average) agrees well with field measurements (points) and follows a similar temporal pattern. The mean differences between remote-sensing-retrieved LAI and field measurements for the field 6-60 were −0.54 in 2018 and 0.27 in 2020, which is in the range of Sentinel-2 LAI retrieval accuracy for crops [53]. The daily PlanetScope-based LAI field averages were used as a parameter for the calibration and evaluation of our model. Note that the Planet team is developing an LAI data product (Planet Fusion LAI) using Planet Fusion surface reflectance. The Planet Fusion LAI can be used directly in this study once it becomes available. In this paper, we used the Sentinel-2 LAI as a reference to retrieve PlanetScope-based LAI, which is different from the Planet Fusion LAI by the Planet team.

Remote Sensing LAI and Phenology Results
The daily PlanetScope-based LAI field averages were used as a parameter for the calibration and evaluation of our model. Note that the Planet team is developing an LAI data product (Planet Fusion LAI) using Planet Fusion surface reflectance. The Planet Fusion LAI can be used directly in this study once it becomes available. In this paper, we used the Sentinel-2 LAI as a reference to retrieve PlanetScope-based LAI, which is different from the Planet Fusion LAI by the Planet team. Figure 3. PlanetScope-based LAI over BARC on 1 August 2020 (a) and LAI time series for the field 6-60 (old airport converted field, red circle in LAI map) from 2018 to 2020 (b). LAI ground measurements (points) from 2018 and 2020 (b) generally followed the LAI time series trend at the field level. In plot (a), the PlanetScope-based LAI map is scaled from 0 to 5 using a grayscale ramp.
The estimated crop growth stages were assessed using the leave-one-out cross-vali- Figure 3. PlanetScope-based LAI over BARC on 1 August 2020 (a) and LAI time series for the field 6-60 (old airport converted field, red circle in LAI map) from 2018 to 2020 (b). LAI ground measurements (points) from 2018 and 2020 (b) generally followed the LAI time series trend at the field level. In plot (a), the PlanetScope-based LAI map is scaled from 0 to 5 using a grayscale ramp.
The estimated crop growth stages were assessed using the leave-one-out cross-validation approach. We found that crop growth stages estimated from the Planet Fusion time series agreed well with the PhenoCam interpreted dates. The average absolute difference between the two was within one week from planting to maturity for both maize and soybean (see Supplementary Table S2). The harvest dates vary from year to year depending on management purposes; however, the remotely sensed harvest dates were not used in the DSSAT calibration process. The DSSAT model was set up to harvest at physiological maturity for both maize and soybean. Considering the uncertainty of photo interpretation, the accuracy of crop growth stages can satisfy the requirement for crop modeling. We used these derived growth stages for our model calibration and evaluation. Table 3 and Figure 4 present the results for flowering and physiological maturity dates, and for the final grain yield of maize for the calibration dataset for 10 site years at BARC. The results are within 9.6 to 12% NRMSE of observations, which indicates that the calibrated model was able to estimate phenology and grain yield reasonably well for most of the fields analyzed. The Willmott's statistics show a good agreement between the simulated and the observed yield (0.65 to 0.85). The results of the RMSE are reasonably low, especially for the flowering and grain yield of maize at BARC. These results demonstrated that remotely sensed derived parameters can be used to augment field observations, especially in datapoor regions. Additionally, these results are consistent with the conclusion of previous studies that integrating remote sensing with field observations as inputs to crop growth models provides an effective solution for assessing crop growth conditions and yield in data-poor regions [18][19][20]40].   For the soybean maturity groups 3 and 4, Tables 4 and 5 and Figures 5 and 6 show the results of phenology (days to flowering and maturity) and yield for the calibration dataset for 4 site years (note, only one year of yield was available for maturity group three). The NRMSE ranged from 5 and 7%, indicating an excellent result for the phenology. This indicates that the calibrated model was able to estimate phenology reasonably well. However, the results of the Willmott's statistics (d index) show a very good, calibrated model of soybean phenology for maturity group three but not for maturity group 4. Yield data for three out of the four fields planted in soybean maturity group three were not available; however, the field with available yield data showed an absolute deviation of 6.6%, which is very good. For the maturity group four, the model-simulated yield for the six calibrated fields did reasonably well (NRMSE of 16%).     Note that when the default model parameters were used, the results were very poor, as some fields could not reach maturity due to frost occurrence (results not shown). Calibration is critical to adapt DSSAT to local conditions. Agronomy 2023, 13, x FOR PEER REVIEW 15 of 23  Table 6 show the results of the flowering, physiological maturity, and the final grain yield of maize for the evaluation (validation) dataset for 7 site years at BARC. The days to flowering, maturity and final grain yield of maize were simulated with fairly good results (NRMSE < 20%) ( Table 6). These results indicate that the simulated parameters strongly agreed with the measured or reference values (Figure 7 and Table 6). The deviation of the simulated from the observed evaluated parameters (flowering, maturity, and grain yield) for soybean maturity group three ranged between 1 and 8% (NRMSE , Table 7), which indicates a very good, evaluated model. The model produced somewhat good results in some fields which is consistent with other evaluation studies [68]. The deviations from the observed values were larger for the maturity stage and grain yield of maize than for the flowering (Table 7). This is expected considering that it is harder to determine maturity based on the PhenoCam images which were used as a reference for the remotely-sensed derived phenology. Additionally, the yield monitor data show considerable variability of grain yields within the maize field ([48,69] see Figure 1a). However, our study used field averages for our calibration and evaluation (validation). For field 5-3A, the measured yield was unusually high (Table 6). Our investigation discovered that the field was heavily manured in the past, but this information was not adequately RMSE=368.97 NRMSE=11.1 D=0.38  Table 6 show the results of the flowering, physiological maturity, and the final grain yield of maize for the evaluation (validation) dataset for 7 site years at BARC. The days to flowering, maturity and final grain yield of maize were simulated with fairly good results (NRMSE < 20%) ( Table 6). These results indicate that the simulated parameters strongly agreed with the measured or reference values (Figure 7 and Table 6). The deviation of the simulated from the observed evaluated parameters (flowering, maturity, and grain yield) for soybean maturity group three ranged between 1 and 8% (NRMSE, Table 7), which indicates a very good, evaluated model. The model produced somewhat good results in some fields which is consistent with other evaluation studies [68]. The deviations from the observed values were larger for the maturity stage and grain yield of maize than for the flowering (Table 7). This is expected considering that it is harder to determine maturity based on the PhenoCam images which were used as a reference for the remotely-sensed derived phenology. Additionally, the yield monitor data show considerable variability of grain yields within the maize field ( [48,69] see Figure 1a). However, our study used field averages for our calibration and evaluation (validation). For field 5-3A, the measured yield was unusually high (Table 6). Our investigation discovered that the field was heavily manured in the past, but this information was not adequately available to be captured during our calibration process. However, we subsequently added organic amendment in the management input file. A total of 2000 kg/ha of forage grass residue material was added. This improved our results, and the percent deviation between observed to simulated decreased. In addition, these results will be improved when more ground measurements become available to be used as a reference for the remotely sensed derived parameters (phenology and LAI) which were used as parameters for the calibration and evaluation of our model.    Note, Obs = observed, Sim = simulated, dev = deviation, R7 = Physiological maturity, * observed yield data missing so not evaluated. D = Index of agreement, RMSE = root-mean-square error; NRMSE = normalized root-meansquare error. Ratings are as follows: Excellent, NRMSE < 10% of mean observed; Good, 10% < NRMSE < 20%; Satisfactory, 20% < NRMSE < 30%; Unsatisfactory, NRMSE > 30% (Threshold based on the recommendation by [66,67] Overall, the evaluation (validation) statistics show that our calibrated models were able to simulate the observed or referenced parameters adequately which is consistent to the previous evaluation studies of the model [70]. Therefore, these calibrated models can be used in assessing the effects of climate change and adaptation option in this region, considering that the cultivars used in this study are common for the region. Additionally, these results show that remotely sensed crop information can be used to assess crop growth conditions and yield forecasting, especially when the needed parameters are not measured or are in a data-poor region.

Figure 7 and
For the soybean maturity groups three and four, Tables 7 and 8 and Figures 8 and 9 show the results of phenology (days to flowering and maturity) and yield over 4 and 6 site years for maturity groups three and four, respectively. The NRMSE ranges from 2 to 8%, indicating an excellent result for the phenology. The yield results indicated an excellent model (NRMSE of 1%, see Table 7) for maturity group three. For maturity group four, the results indicated satisfactory results (NRMSE of 30%). This is expected, owing to the fact that the management information was not readily available for maturity group four. This created uncertainties in many of our input data. In most cases, the cultivar numbers were not indicated in the management data record and assumptions were made about maturity group four basedF on the information available from other fields for the same year. In general, our results were consistent with the other evaluation studies of the model.

Conclusions
This paper assessed the integration of remotely sensed crop information and field observations to evaluate the DSSAT model in simulating soybean and maize growth and yield. The study demonstrated that remotely sensed crop variables such as LAI and phenology can be used for fitting model cultivar coefficients for yield estimation. The results show that the model's comparison with the observed crop variables assessed produced a good result based on our limited range of data values in both the calibration and validation period. The model reasonably simulated measured soybean and maize phenology and produced a reasonable yield estimate, but not for all treatments. The evaluation of the model against the independent multi-year sites produced good results for phenology and somewhat good results for yield. Overall, the model performance with respect to phenology and grain yield was consistent with other evaluation studies of the model. These results show that remotely sensed crop information can be used to augment field observations for crop modeling, especially for data-poor regions or where no field records are available for modeling purposes.

Conclusions
This paper assessed the integration of remotely sensed crop information and field observations to evaluate the DSSAT model in simulating soybean and maize growth and yield. The study demonstrated that remotely sensed crop variables such as LAI and phenology can be used for fitting model cultivar coefficients for yield estimation. The results show that the model's comparison with the observed crop variables assessed produced a good result based on our limited range of data values in both the calibration and validation period. The model reasonably simulated measured soybean and maize phenology and produced a reasonable yield estimate, but not for all treatments. The evaluation of the model against the independent multi-year sites produced good results for phenology and somewhat good results for yield. Overall, the model performance with respect to phenology and grain yield was consistent with other evaluation studies of the model. These results show that remotely sensed crop information can be used to augment field observations for crop modeling, especially for data-poor regions or where no field records are available for modeling purposes.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/agronomy13061540/s1, Table S1: The cross comparison for maize showing the absolute difference (in days) between the remote sensing estimation and PhenoCaminterpreted growth stages for some fields at Beltsville Agricultural Research Center in Beltsville, Maryland; Table S2: Decision Support System for Agro-Technology (DSSAT) soil file for soybean fields.