Research Articles

Damage evaluation of soybean chilling injury based on Google Earth Engine (GEE) and crop modelling

  • CAO Juan 1 ,
  • ZHANG Zhao , 1, * ,
  • ZHANG Liangliang 1 ,
  • LUO Yuchuan 1 ,
  • LI Ziyue 1 ,
  • TAO Fulu 2, 3
  • 1. State Key Laboratory of Earth Surface Processes and Resource Ecology/MEM&MoE Key Laboratory of Environmental Change and Natural Hazards, Faculty of Geographical Science, Beijing Normal University, Beijing 100875, China
  • 2. Key Laboratory of Land Surface Pattern and Simulation, Institute of Geographic Sciences and Natural Resources Research, CAS, Beijing 100101, China
  • 3. College of Resources and Environment, University of Chinese Academy of Sciences, Beijing 100049, China
Zhang Zhao, Professor, E-mail:

Cao Juan, PhD Candidate, specialized in agricultural disaters and agricultural insurance

Received date: 2020-03-10

  Accepted date: 2020-05-12

  Online published: 2020-10-25

Supported by

National Natural Science Foundation of China(41977405)

National Natural Science Foundation of China(41571493)

National Natural Science Foundation of China(31561143003)

National Natural Science Foundation of China(31761143006)

National Key Research & Development Program of China(2017YFA0604703)

National Key Research & Development Program of China(2019YFA0607401)


Copyright reserved © 2020. Office of Journal of Geographical Sciences All articles published represent the opinions of the authors, and do not reflect the official policy of the Chinese Medical Association or the Editorial Board, unless this is clearly specified.


Frequent chilling injury has serious impacts on national food security and in northeastern China heavily affects grain yields. Timely and accurate measures are desirable for assessing associated large-scale impacts and are prerequisites to disaster reduction. Therefore, we propose a novel means to efficiently assess the impacts of chilling injury on soybean. Specific chilling injury events were diagnosed in 1989, 1995, 2003, 2009, and 2018 in Oroqen community. In total, 512 combinations scenarios were established using the localized CROPGRO-Soybean model. Furthermore, we determined the maximum wide dynamic vegetation index (WDRVI) and corresponding date of critical windows of the early and late growing seasons using the GEE (Google Earth Engine) platform, then constructed 1600 cold vulnerability models on CDD (Cold Degree Days), the simulated LAI (Leaf Area Index) and yields from the CROPGRO-Soybean model. Finally, we calculated pixel yields losses according to the corresponding vulnerability models. The findings show that simulated historical yield losses in 1989, 1995, 2003 and 2009 were measured at 9.6%, 29.8%, 50.5%, and 15.7%, respectively, closely (all errors are within one standard deviation) reflecting actual losses (6.4%, 39.2%, 47.7%, and 13.2%, respectively). The above proposed method was applied to evaluate the yield loss for 2018 at the pixel scale. Specifically, a sentinel-2A image was used for 10-m high precision yield mapping, and the estimated losses were found to characterize the actual yield losses from 2018 cold events. The results highlight that the proposed method can efficiently and accurately assess the effects of chilling injury on soybean crops.

Cite this article

CAO Juan , ZHANG Zhao , ZHANG Liangliang , LUO Yuchuan , LI Ziyue , TAO Fulu . Damage evaluation of soybean chilling injury based on Google Earth Engine (GEE) and crop modelling[J]. Journal of Geographical Sciences, 2020 , 30(8) : 1249 -1265 . DOI: 10.1007/s11442-020-1780-1

1 Introduction

Soybean (Glycine max L.) is the world's most common oil-bearing crop and an important protein source for human populations, playing an essential role in global economic development (Tomar et al., 2009; Hellal and Abdelhamid, 2013; Das and Parida, 2014). Roughly two-thirds of the world's population does not consume enough protein. To balance population growth with protein needs, many countries have recently tried to increase soybean yields and planting areas (Li et al., 2019). Planting areas for soybeans in China are the most widespread after those for corn, rice, wheat, and potatoes. Soybean crops are also the most economically efficient and import intensive crops. Therefore, the soybean is an important basic and strategic material and occupies an extremely important position in China's national economy. The No.1 document of 2019 issued by Central Committee of the CPC describes the “implementation of a soybean and dairy revitalization plan” as one of the ten key objectives and calls to expand soybean planting areas; promote the demonstration and popularization of new varieties, technologies, and forms of mechanization; stabilize grain production; and improve subsidies for soybean producers. Northeastern China is one of the country's main soybean producing areas. Planting areas and yields in northeastern China are among the highest in the country, accounting for 1/3 of national soybean planting areas and yields and thus playing an important role in food production in China (Liu et al., 2008; Yang et al., 2020).
In terms of global warming, however, warming trends of the crop growing season have not been obvious. The crop planting belt has shifted north while varieties ripening later have increasingly extended their reach in recent years (Liu et al., 2016a; Fang et al., 2019). As a result, chilling injury is the most frequent agricultural problem experienced in northern China, resulting in the most serious yield losses (Han et al., 2010; Hao et al., 2010). From 1951 to 1980, eight severe cold periods occurred, three of which (1969, 1972, and 1976) resulted in yield losses of roughly 5.78 billion kg (Zhang et al., 2012). Soybeans grown in the northeast at high latitudes suffered more severe and frequent cases of chilling injury than those in other areas (Zhou et al., 2017). Due to the insufficient supply of heat and considerable temperature variations in northeastern China, soybean can suffer a series of problems throughout the growing period, i.e., the delay of the growth period, dysplasia, empty shells, and chaff grains, reducing soybean yield and quality and economic losses. Therefore, it is essential to monitor chilling injury occurrence and evaluate yield losses of soybean crops in northeastern China (Liu et al., 2012).
Several recent studies have monitored chilling injury and assessed yield losses mainly using crop modelling and remote sensing monitoring. The calibrated crop model can also be used to simulate chilling injury scenarios at the field level. While these methods are easily applied at the field level, inferring to larger areas has limited validity due to a need of extensive and detailed information (e.g., variety characteristics, climate variables, soil conditions, and field management and cultivation methods). While these data are easily available at the field level, errors are inevitable when using field-level data for regional estimations due to high levels of spatial heterogeneity (Sinclair et al., 1993). Remote sensing monitoring, the constructed vegetation index and the temperature index mainly focus on distinguishing occurrences and levels of chilling injury. In addition, chilling injury indices and ground-based yield measurements collected in years with natural disasters are used to establish empirical relationships. However, while remote sensing monitoring can explain over 80% of yield variations among individual fields (Sinclair et al., 1993; Liu et al., 2016b), it is difficult to extrapolate these equations to new locations or years due to their limited spatial generalizability. In addition, while numerous studies have focused on impacts of chilling injury on corn and rice crops in northeastern China (Liu et al., 2016b; Chen, 2017), few have explored the impacts of chilling injury on soybeans in the same region. The above limitations highlight the need to develop more general means to quickly and accurately evaluate the effects of chilling injury on soybean crops.
Yield loss assessments, coupling crop models and remote sensing monitoring have been major issues of focus in recent years. The Google Earth Engine (GEE), with its powerful data storage and cloud computing capabilities, allows for the free sharing of image data (Gorelick et al., 2017) and the CROPGRO-Soybean model can dynamically and intuitively simulate soybean growth, development, and yields. These new tools allow for soybean growth monitoring and disaster loss assessment at high spatiotemporal resolutions (Boote et al., 2003). Using the GEE platform and crop model, we develop a fast, accurate, and cross-regional method for evaluating the impacts of chilling injury on soybean yields in northeastern China. We first optimize parameters of a localized soybean model, simulate cold injury scenarios to construct soybean chilling injury vulnerability models, and then spatialize yield loss mapping at a high resolution. Finally, yield losses at the county level are used to evaluate the performance of the proposed technique. Our aim is to provide a new method of soybean yield estimation and disaster loss assessment.

2 Study area and data sources

2.1 Study area

Oroqen Autonomous Banner, lying between 48°50′-51°25′N and 121°55′-126°10′E, is located between the Da Hinggan Mountains and Hulun Beir Grassland. It covers an area of approximately 6×104 km2 with roughly 2.1×103 km2 of cultivated land. The local terrain, which includes plateaus, mountains, and hills, forms a typical plateau landscape. The study area is characterized by a sub-humid warm temperate continental monsoon climate with four distinct seasons. Average annual temperatures range from -2.7℃ to 0.8℃ and increase from west to east. Temperatures can reach as low as -30℃ in January and as high as 37.5℃ in July with average temperatures ranging from 17.9℃ to 19.8℃. Approximately 95 d in the year are frost-free. Annual precipitation reaches 459.3-493.4 mm. Moreover, wind speeds are relatively low at 1.8 to 2.9 m/s. The Oroqen area serves as a major grain producer in the Inner Mongolian Autonomous Region with an average annual production level of 200 million kg. Soybean is the most common crop grown in this region (Figure 1).
Figure 1 Overview of the study area

2.2 Dataset and preprocessing

2.2.1 Satellite data
Remote sensing has been widely used for crop growth monitoring and yield prediction, as it rapidly, macroscopically, and dynamically captures a wide range of data on the earth's surface through various bands (Chen et al., 2016). The Landsat 5 and Sentinel-2A surface reflectance products were first preprocessed on the GEE platform, i.e., radiation, geometric, and atmospheric correction. We then selected all available images visualizing critical soybean growth stages (i.e., May-September) using JavaScript API, including Landsat 5 images for 1984 to 2012 with a 30-m resolution and 16-d repeat and Sentinel-2A images for 2017 to 2019 with a 10-m resolution and 5-d repeat. Their quality control (QA) bands were used to remove effects of clouds using codes for Landsat 5 and Sentinel-2A images available at and, respectively. After the above processing, we obtained the spectral curve and calculated the vegetation index with good fit to the ground-based curve measures.
2.2.2 Meteorological and soil data
Daily climate data for 1981 to 2018, including those for average, maximum and minimum temperature (℃); rainfall (mm); and sunshine hours (h), were obtained from the National Meteorological Information Centre ( Solar radiation levels were calculated according to previous research (Angstrom et al., 1924). Soil data used for the crop model are available from the Chinese soil science database ( and include data on bulk density; organic carbon; total nitrogen; pH; and percentages of silt, sand, and clay. Soil albedo, runoff potential, and drainage rate values were also estimated (Jones et al., 2003). The soil fertility factor (SLPF) was set as 1.0, and the initial soil moisture content level was set as the field moisture capacity.
Daily climate data for 1981 to 2018, including those for average, maximum and minimum temperature (℃); rainfall (mm); and sunshine hours (h), were obtained from the National Meteorological Information Centre ( Solar radiation levels were calculated according to previous research (Angstrom et al., 1924). Soil data used for the crop model are available from the Chinese soil science database ( and include data on bulk density; organic carbon; total nitrogen; pH; and percentages of silt, sand, and clay. Soil albedo, runoff potential, and drainage rate values were also estimated (Jones et al., 2003). The soil fertility factor (SLPF) was set as 1.0, and the initial soil moisture content level was set as the field moisture capacity.
2.2.3 Land use and field measurement data
Six ground-based field measures for 2014 to 2017 for model calibration and localization and a soybean planting area map of Oroqen with a 10-m resolution were obtained from a local insurance company. As shown in Tables 1 and 2, ground-based field data used mainly include location data, basic information on soybean crops (e.g., variety maturity, yield, cultivar, and development period), and field management data (tillage, weeding, irrigation, fertilization, etc.). Recorded Oroqen soybean yields at the county level were obtained from the Department of Plantation Management ( and the Statistical Yearbook for Inner Mongolia for 1981 to 2017 ( River and administrative divisions were derived from the Resource and Environment Data Cloud Platform ( Finally, projections of all data were uniformly converted to WGS-1984 form.
Table 1 The observation information of key growth period, planting density and yield of crop fields
Year Crop field Wuchagou Xingsheng Chuntinge Chunlin Naimuhe Xiaokumo
2014 Sowing date 16 May. 18 May. 1 May. 18 May. 28 May. 11 May.
Flowering date 14 Jul. 12 Jul. 30 Jun. 6 Jul. 18 Jul. 30 Jul.
Maturity date 26 Sept. 28 Sept. 18 Sept. 28 Sept. 28 Sept. 24 Sept.
Plant density (plants/m2) 39.69 30.42 47.73 36.75 37.98 38.47
Yield(kg/ha) 1160 1350 1650 1150 610 1050
2015 Sowing date 16 May. 16 May. 28 May. 8 May. 17 May. 18 May.
Flowering date 16 Jul. 10 Jul. 2 Jul. 12 Jul. 16 Jul. 8 Jul.
Maturity date 16 Sept. 14 Sept. 28 Sept. 28 Sept. 28 Sept. 8 Sept.
Plant density (plants/m2) 47.94 35.6 38.18 26.24 25.84 28.58
Yield(kg/ha) 980 1100 1410 1030 510 860
2016 Sowing date 6 Jun. 16 May. 6 May. 22 May. 18 May. 18 May.
Flowering date 20 Jul. 14 Jul. 28 Jul. 28 Jul. 2 Jul. 2 Jul.
Maturity date 30 Sept. 26 Sept. 18 Sept. 30 Sept. 18 Sept. 14 Sept.
Plant density (plants/m2) 32.97 30.55 40.04 42.9 40.8 35.82
Yield(kg/ha) 1280 1450 1850 1370 670 1110
2017 Sowing date 24 May. 23 May. 21 May. 18 May. 24 May. 2 Jun.
Flowering date 8 Jul. 12 Jul. 22 Jul. 8 Jul. 1 Jul. 14 Jul.
Maturity date 28 Sept. 24 Sept. 22 Sept. 18 Sept. 18 Sept. 24 Sept.
Plant density (plants/m2) 43.5 35.76 27.56 36.9 40.6 52.36
Yield(kg/ha) 1350 1500 1240 1500 750 1200
Table 2 The basic information of six crop fields
Crop field Wuchagou Xingsheng Chuntinge Chunlin Naimuhe Xiaokumo
Longitude (°) 124.42 124.49 124.49 124.42 124.11 124.17
Latitude (°) 50.10 49.97 49.83 49.83 49.4 49.7
Elevation (m) 486 393 380 365 381 447
Cultivar Heihe 38 Heilong 35 Hefeng 50 Hefeng 25 Hefeng 39 Heihe 18
Maturity Medium Medium Medium Mid-late Medium Medium
Cultivation way Drilling Drilling Drilling Ridge tillage Drilling Drilling

3 Method

3.1 Determination of historical cold years

Chilling injury assessment first involves determining historical cold years. According to previous research (Pang et al., 1983), two conditions should be met: (1) total Growing Degree Days (GDD), i.e., the mean daily temperature of above 10℃ for the growth period in a given year is less than 100℃ that of the historical mean soybean growth stage, and (2) lower soybean yields in the current year than in the preceding normal year (Pang et al., 1983). In this research, we only focus on identifying historical cold years and do not distinguish the degree of chilling injury. Years with GDD with temperatures below the historical mean of -100℃ include 1989, 1991, 1995, 2003, 2009, and 2018 (Figure 2), which is consistent with previous research (Zhang, 2017). However, no soybean yield losses occurred in Oroqen in 1991. Therefore, we only use 1989, 1995, 2003, 2009, and 2018 as historical cold years. Meanwhile, the preceding years within the study period (i.e., 1988, 1994, 2002, 2008, and 2017) are selected as normal years for our comparative analysis.
Figure 2 The cold years in Northeast China

3.2 Localization of the CROPGRO-Soybean model

The CROPGRO-Soybean model available through DSSAT 4.7 has been widely used to simulate soybean growth and development around the world (Wang et al., 1995; Zhu et al., 2008). In the CROPGRO-Soybean model, 15 genetic parameters can determine a soybean cultivar (Table 3). In this study, ground-based field measurements for 2014 to 2017 were used to adjust optimal parameters and then for model localization. First, genetic parameters were initialized, and the number of iterations was set as 6000. Soybean variety parameters were calculated based on the GLUE (Generalized Uncertainty Estimation) module of the DSSAT model. Specifically, field measurements for 2014 to 2015 were used for model calibration, and those for 2016 to 2017 were used for model testing. The genetic parameters of six crop fields were then successively simulated. Finally, Root Mean Square Error (RMSE) and Relative Root Mean Square Error (RRMSE) tests were used to estimate deviation between the simulated and observed variables (i.e., flowering date, maturity date, and yield).
Table 3 CROPGRO-Soybean model parameters and the genetic coefficients of crop fields
Coefficient Definition Default value Wuchagou Xingsheng Chuntinge Chunlin Naimuhe Xiaokumo
CSDL Critical short day length below which reproductive development progresses with no day length effects (h) 12.15 14.03 12.23 13.06 14.58 11.9 14.03
PPSEN Slope of the relative response of development to photoperiod with time (1/h) 0.2 0.235 0.294 0.287 0.229 0.146 0.304
EM-FL Time between plant emergence and flower appearance (day) 21 13.09 18.33 15.46 16.53 26.14 22.74
FL-SH Time between first flower and first pod (day) 6 6 6 6 6 6 6
FL-SD Time between first flower and first seed (day) 12 19.82 21.06 18.42 18.39 11.26 12.06
SD-PM Time between first seed and physiological maturity (day) 26 37.56 24.5 36.63 31.27 22.25 34.33
FL-LF Time between first flower and end of leaf expansion (day) 20 20 20 20 20 20 20
LFMAX Maximum leaf photosynthesis rate at 30℃, 350 vpm CO2, and high light (mg CO2/m2 s) 1.03 1.023 1.052 1.034 1.011 1.196 1.257
SLAVR Specific leaf area of cultivar under standard growth conditions (cm2/g) 385 311.2 303.8 301 317.6 337.9 301.1
SIZLF Maximum size of full leaf (three leaflets (cm2) 137 138.1 145.1 138 141.2 188.5 217.6
XFRT Maximum fraction of daily growth that is partitioned to seed+shell (-) 1 1 1 1 1 1 1
WTPSD Maximum weight per seed (g) 0.155 0.162 0.157 0.161 0.181 0.195 0.186
SFDUR Seed filling duration for pod cohort at standard growth conditions (day) 22 25.42 24.88 24.22 25.36 21.62 21.94
SDPDV Average seed per pod under standard growing conditions (numbers per pod) 2.2 2.415 2.276 2.09 2.24 2.397 1.794
PODUR Time required for cultivar to reach final pod load under optimal conditions (day) 13 13 13 13 13 13 13

3.3 Yield loss estimation

Yield loss assessment involves five steps: (1) simulating different combinations of cold and local field management scenarios based on the localized CROPGRO-Soybean model and then recording the simulated yield, daily Leaf Area Index (LAI), and the corresponding date (Day of year, DOY); (2) calculating GDD of each simulation scenario and then subtracting the mean number of GDD occurring from 1981 to 2018, namely, as Cold Degree Days (CDD); (3) constructing a vulnerability model of chilling injury from the simulated yield, CDD, and any combination of two LAIs for the critical growth period; (4) retrieving the actual LAI and acquisition date (DOY) of soybean planting grids on the GEE platform and calculating the actual number of CDD based on meteorological data; and (5) calculating the yield reduction rate on a pixel-by-pixel basis using the best available observation dates for each pixel.
3.3.1 Simulation of cold injury scenarios
A calibrated model can better simulate the impacts of climatic factors on crop growth and development and assess yield loss (Pang, 2015). In this study, a series of combined parameters of different sowing dates and planting densities was first generated. According to local field management records, sowing dates are set at equal intervals (5 d) from May 1st to June 6th and include May 1st, 6th, 11th, 16th, 22th, and 27th and June 1st and 6th. Planting densities were randomly generated at between 19.42 and 52.36 plants/m2, including values of 19.42, 22.77, 31.82, 35.86, 38.87, 41.06, 47.43, and 52.36 plants/m2. We also set up normal weather scenarios and seven cold injury scenarios: (1) three cold injury scenarios with daily temperature reductions of 1-3℃ across growth stages and (2) four other cold injury scenarios set randomly at a minimum temperature of 0℃ for 5 consecutive days from seedling emergence to the flowering stage, from the pod bearing stage to the pod bearing stage, and from the pod bearing stage to the mature stage. For each soybean field, we have 512 (8*8*8) simulation scenarios. Finally, the daily LAI and corresponding simulated yield of all soybean fields in the growth stage were recorded for analysis.
3.3.2 Calculating the chilling injury index
Effective accumulated temperatures exceeding 10℃ during the soybean growth periods of the simulation scenarios were calculated. Namely, the number of GDD and CDD was calculated by subtracting the mean number of GDD occurring in growth periods from 1981 to 2018. This CDD index was then used to determine the impact of chilling injury on final yields (Ma et al., 2003). Since the sowing dates and weather conditions are different for each simulation, the CDD values calculated from each simulation are also different. The CDD calculation formula is as follows:
where CDD is the effective accumulated temperature anomaly obtained from each simulation, D is the average daily temperature exceeding 10℃ on day i of the simulation growth period, n is the length of each simulation growth period, and $\overline{GDD}$ is the average effective accumulated temperature value of larger than 10℃ during the same period from 1980 to 2018.
3.3.3 Constructing vulnerability models of chilling injury
We expanded our analysis to 40 d before and after the soybean flowering date to create early and late windows. The simulated yields of 512 combination scenarios for six crop fields (6*512) were used as the dependent variable, and any two LAIs from the early and late windows and CDD were used as independent variables to construct multiple regression equations. A total of 1600 (40*40) equations were established.
$Yield={{\beta }_{0}}+{{\beta }_{1}}\times LA{{I}_{1,d}}+{{\beta }_{\text{2}}}\times LA{{I}_{2,d}}+{{\beta }_{3}}\times CDD$
where Yield represents the output of each simulation and LAI1,d and LAI2,d denote the LAI d days from the early window and late window, respectively. CDD represents the cool injury index for each simulation. β0, β1, β2, and β3 represent the regression coefficients of the constant, the LAI of the early window, the LAI of the late window, and CDD, respectively. Because the dates of usable imagery may differ by location and year and particularly for optical imagery for which cloud cover can vary by pixel, we compute regressions for an arbitrary number of combinations of image dates. Finally, we store the resulting coefficients for further analysis.
3.3.4 Estimations of yield loss in cold years
Previous studies have shown that the Wide Dynamic Range Vegetation Index (WDRVI) is sensitive to variations in the LAI and is closely related to the LAI (Anthony et al., 2012). Therefore, we first use the GEE platform to extract the maximum WDRVI and DOY of the early and late windows on a pixel-by-pixel basis and convert the WDRVI into the actual LAI of soybean planting grids (Lobell et al., 2015). We then apply regression equations to actual measurements obtained from satellite data (i.e., the actual LAI calculated from the WDRVI and DOY). This can generally be done on a pixel-by-pixel basis using the maximum WDRVI and DOY of the early and late windows for each pixel. The final step is to estimate yield and reduction rate at the pixel level and then generate yield estimations and yield reduction maps with a high resolution (i.e., 10 or 30 m). Due to a lack of other field measurements, we calculate the county-level yield and yield reduction rate from the mean of the pixel yield and reduction rate and use the recorded county-level yield to evaluate the performance of our proposed technique when applied to the Oroqen area. The yield reduction rate is expressed as “(the annual yield of a normal year - the annual yield of a cold year)/the annual yield of a cold year ×100%.”
$WDRVI=\frac{\alpha {{\rho }_{NIR}}-{{\rho }_{red}}}{\alpha {{\rho }_{NIR}}+{{\rho }_{red}}}$
$LAI=-\left\{ \left[ \ln \text{(}1.79-WDRVI\text{)}-0.532 \right]/0.3 \right\}$
where ρNIR and ρred represent reflectance at near infrared and red wavelengths, respectively. α is set to 1 to ensure linearity between WDRVI and LAI (Jiang et al., 2019).

4 Results and analysis

4.1 The performance of the CROPGRO-Soybean model

Comparisons of the observed and simulated variables for soybean crops for 2014 to 2017 show that CROPGRO-Soybean model effectively simulate the actual ADAT (Time to anthesis as days after planting) and MDAT (Time to maturity as days after planting) (Figure 3). The RMSE values for average flowing and maturity dates are less than 5 d and 6 d, respectively. The RRMSE values are below 9.3% and 7.9%, respectively. The RMSE for yields is 93.3 kg/ha and the RRMSE for yields is 8.0%. Overall, the simulation of the maturity date is superior to that of the yield and flowing date, but all simulation results falls within a reasonable range (RRMSE≤10%), indicating that the calibrated CROPGRO-Soybean model can reflect soybean growth and yields under different weather and management scenarios.
Figure 3 Comparison between simulated and observed variables of soybean

4.2 The results of vulnerability models

Figure 4 illustrates the simulated profiles of soybean LAIs determined from 512 simulations. The LAI values vary considerably from 2 to 7, which is largely due to the fact that soybean growth processes are variable under different weather conditions and farming management systems. Note that the early and late windows (denoted by shaded rectangles in the figure) basically cover the critical growth period of each simulation. These simulation results thus show a broad range of variability in training the regression model given in Eq. (1).
Figure 4 Daily leaf area index outputs from 512 simulations of the CROPGRO-Soybean model (the shaded rectangles indicate the “early-season” and “late-season” windows used for image observations)
Intuitively, the capacity for a linear regression model based on the LAI to explain simulated variability depends on the specific times at which observations are made (Figure 5), but this was generally found to be high (R2>0.8), especially for LAIs measured roughly DOY 191 and 210 from the early window and on all dates within the late window (R2>0.95). This is mainly the case because the vegetation index is closer to the maximum value during this period, and thus LAIs can capture more soybean crop canopy foliage and chlorophyll content and final yield variations. In addition, minimal training errors do not guarantee good data performance, as considerable input data (e.g., weather and soil) include errors in reality, and crop models do not perfectly capture all relevant processes of crop growth and development (Lobell et al., 2015). Nonetheless, it is found that the high R2 indicates any two LAI observations and CDD during the critical growth period should in many cases be sufficient to obtain accurate yield estimations.
Figure 5 The coefficient of determination (R2) for the regression models

4.3 Estimation of yield loss

4.3.1 Actual LAI based on the GEE platform
Figure 6 shows the two maximum LAIs (6a and 6b) and the corresponding DOY image (6c and 6d) for the early and late windows are based on the GEE platform. The actual maximum LAIs are similar to the simulated LAIs and range from 6 to 7.2, indicating that the seven simulated cold injury scenarios can accurately describe the relationship between local chilling injury and yields. While the specific timing at which satellite observations of soybean planting grids are taken is not consistent, limitations in spatial and temporal discontinuity affecting remote sensing images can be overcome by combining two LAI observations collected during the growth season. Thus, our method can creatively extrapolate crop model simulations from the field level to large-scale yield estimations and disaster monitoring at a high resolution.
Figure 6 The maximum LAI and its specific observation dates of early and late growing season windows obtained from Sentinel-2 in 2018 (a and c represent the maximum LAI and its specific observation dates of early growing season windows; b and d for late growing season windows, respectively)
4.3.2 Yield simulation under different cold injury scenarios
We simulated seven scenarios to explore the effects of different cold injury events on soybean yields. The simulated yields varied widely under different cold injury scenarios as illustrated in Figure 7. During growth stages (i.e., scenarios 1, 2 and 3), the more temperatures decreased, the more soybean yields decreased. During the soybean growing season, the daily minimum temperature decreased by 1℃, GDD decreased by roughly 100℃, and yields ranged from 910.5 to 1242.5 kg/ha. The daily mean temperature decreased by 2℃, GDD decreased by roughly 200℃, and yields ranged from 571.8 to 838.0 kg/ha. The daily mean temperature decreased by 3℃, GDD decreased by roughly 300℃, and yield ranged from 273.8 to 435.7 kg/ha. Soybean plants are thermophilic crops, especially after reaching the flower bud stage. Thus, chilling injury inhibits crop growth in the following stages, leading to filling stage obstruction and then (sometimes extreme) yield reductions. Additionally, scenarios 4, 5, 6 and 7 show that soybean yield responses to decreased temperatures in different stages vary with the strongest impacts observed from the pod bearing stage to the filling stage. This is the case because soybeans are more sensitive to temperatures in later stages, as chilling injury in the filling stage affects fruit formation, increases abortive grain levels, and delays maturation, resulting in (sometimes extreme) yield reduction, which is in line with previous research (Sang et al., 2013).
Figure 7 Estimated yields by calibrated CROPGRO-Soybean model under different cold injury scenarios (scenarios S1, S2 and S3 are set as reducing 3℃, 2℃ and 1℃ at the whole growth stages, respectively. Scenario N4 is set as actual weather. Scenarios S5, S6, S7 and S8 are set randomly as minimum temperature of 0℃ for 5 consecutive days during the periods from seedling emergence to flowering stage, flowering stage to pod bearing stage, and pod bearing stage to filling stage, and filling stage to maturity stage, respectively.)
4.3.3 Estimating the rate of yield reduction at the county level in historical cold years
According to our analysis of historical meteorological data for Oroqen, five years (1989,1995, 2003, 2009 and 2018) were set as historical cold years while the preceding years (1988, 1994, 2002, 2008 and 2017) were set as normal years. As no Sentinel-2A images for before 2015 were available, Landsat 5 images with a 30-m spatial resolution were used to generate yield maps. For 2017 and 2018, Sentinel-2A images were used to generate yield maps with a 10-m resolution. Finally, we calculated the yield reduction rate for five cold years via yield mapping. However, due to a lack of field measurements, the county-level yield reduction rate was calculated from the mean of all pixels. As shown in Figure 8 (please note that the actual recorded yield for 2018 is missing), the change in the historical yield reduction rate is similar to the distribution of GDD reduction (Figure 2) for soybean growth stages. For example, GDD was the lowest in 2003, causing the most serious chilling injury and the highest yield reduction rate for this year (actual of 47.7% vs. simulated of 50.5%); the lowest reduction rate occurred in 1989 when the simulated reduction rate was 9.6% while the actual reduction rate was 6.4%. While the simulated yield reduction rate of 1995 is slightly higher than the actual yield reduction rate, other simulated reduction rates for 1989, 2003, and 2009 are slightly lower than actual records. Errors of the simulations of historical cold years fall within 1SD, indicating that the proposed method can accurately evaluate soybean yield losses caused by chilling injury.
Figure 8 The comparison between actual and simulated yield losses
4.3.4 Yield spatial distributions in 2017 and 2018
Our evaluation results for historical cold years indicate that our method can accurately simulate the yield loss caused by chilling injury in Oroqen. We thus used this technique to generate a yield spatial distribution map for the last cold year (2018) and its preceding year (2017, defined as the normal year) (Figure 9). The county-level yield was obtained from the average of all pixels. Estimated county-level yields for 2017 and 2018 are 2004.6 kg/ha and 1540 kg/ha, respectively. As illustrated in Figure 9, the estimated yield for cold years (i.e., 2018) is significantly lower than that for normal years (i.e., 2017), which is consistent with actual conditions. According to our magnification figures (Figures 9a-9d), the yield losses of different fields vary in the midst of chilling injury. Figures 9a and 9b show that field yields were significantly lower in 2018 than in 2017 after chilling injury, that some areas (Figure 9c) maintained similar yields in both years, and that for a few fields (Figure 9d) yields slightly increased from 2017 levels. Our field investigations of Oroqen show that such trends are attributable to differences in field management measures, irrigation conditions, and topographic features.
Figure 9 Spatial distribution of estimated yields in normal year (2017) and cold year (2018)
4.3.5 Yield reduction rate assessment for 2018
Spatially, our results show that yields in Oroqen decreased significantly from 2017 to 2018 (Figure 10). Slight yield losses occurred throughout the study area with the exception of a few areas showing significant yield losses (i.e., the northeastern section of the study area) and others showing no changes (i.e., the southern section of the study area). Oroqen's yield reduction rate reached 23.1% in 2018. The reduction rate caused by chilling injury at the town level is shown in Table 4, illustrating that roughly 98% of the area suffered slight yield reductions except in a town not affected by the disaster, which shows a reduction rate of -29.94%. Thus, chilling injury had serious impacts on soybean crops in Oroqen in 2018.
Figure 10 The spatial distribution of yield losses in 2018 relative to 2017
Table 4 The estimated yield losses of soybean at town level
Town Yield loss (%) Town Yield loss (%) Town Yield loss (%)
Karichu 21.23 Wobei 29.77 Xiangyang 24.35
Lingnan 14.64 Zhaxi 19.54 Wuerqi 31.14
Neerkeqi 22.61 Wuchagou 15.86 Chaoyangliehu 21.19
Qingsonggou 24.52 Naimuhe 25.15 Shiliudongfang 32.28
Tuanjie 19.38 Maweishan 18.34 Nuominghe 6.54
Yuchang 8.76 Kuweidi 22.58 Xiaoerhong 18.08
Xinxing 9.55 Yuejing 28.83 Longtou 29.25
Oukenhe 18.24 Kuilehe 5.29 Xinfeng 20.63
Xingsheng 28.04 Dongsheng 21.52 Woluohe 30.79
Ershili -29.94 Xiaokumo 17.63 Tiedong 30.05
Wulubutie 9.53 Doushigou 21.30 Dakumo 21.58
Hongqi 22.65 Chaoyanggou 22.67 Ergenghe 7.97
Maanshan 18.52 Xinfa 26.92 Dongsheng 21.52
Chunlin 26.02 Chuntinge 16.09 Maojiapu 18.77

5 Conclusions and discussion

By combining the GEE platform with crop modelling, yield losses caused by chilling injury to Oroqen soybeans were evaluated at a 10-m resolution, and a high-resolution yield map was produced. We found that our calibrated crop model can simulate crop growth and yield variations under different climate scenarios and field management measures and that yield losses caused by continuous soybean cooling (i.e., 1-3℃) over several growth stages are more serious than those caused by a single growth period (i.e., four randomly selected growth periods with a minimum temperature of 0℃ lasting 5 consecutive days). In comparing simulated and recorded county-level soybean yields, we conclude that our novel technique can effectively simulate soybean yields. Moreover, county-level yield reduction rates are basically consistent with actual conditions observed in historical cold years, and spatial variations of yield losses can be effectively characterized at the pixel level. Compared to the traditional crop model and remote sensing monitoring, this novel method is more efficient and accurate and can be applied for large-scale, high-resolution field yield estimations and disaster assessments. Therefore, the present study provides a new framework for soybean yield and disaster loss assessment.
Nonetheless, due to a lack of data on yield losses at the field level, it is necessary integrate all yields of the pixel level with those of the county level to evaluate the model's performance. More ground-based yield measurements should thus be collected to further verify the capacities of the proposed technique. To produce yield loss maps with a 10-m resolution, we used Sentinel-2A images. However, as the Sentinel series was launched in 2015, for disaster years preceding 2015, only Landsat 5 series with a 30-m resolution can be used for map analysis. Therefore, future studies should focus on combining multi-source remote sensing data and multi-crop models for yield prediction and loss assessments.
Angstrom A, 1924. Solar and terrestrial radiation. Quarterly Journal of the Royal Meteorological Society, 50(210):121-126.


Anthony N R, Anatoly G, Yi P et al., 2012. Green leaf area index estimation in maize and soybean: Combining vegetation indices to achieve maximal sensitivity. Agronomy Journal, 104(5):1336-1337.


Boote K J, Jones J W, Batchelor W D et al., 2003. Genetic coefficients in the CROPGRO-soybean model. Agronomy Journal, 95(1):32-51.


Chen D, 2017. Study on monitoring and evaluating the chilling injury of rice in northeast China by remote sensing and crop model[D]. Nanjing: Nanjing University of Information Technology. (in Chinese)

Chen Z X, Ren J Q, Tang H J et al., 2016. Progress and perspectives on agricultural remote sensing research and applications in China. Journal of Remote Sensing, 20(5):748-767. (in Chinese)

Das A, Parida S K, 2014. Advances in biotechnological applications in three important food legumes. Plant Biotechnology Reports, 8(2):83-99.


Fang H, Zhang Y, Wei S, et al. 2019. Validation of global moderate resolution leaf area index (LAI) products over croplands in northeastern China. Remote Sensing of Environment, 233:11377.

Gorelick N, Hancher M, Dixon M et al., 2017. Google Earth Engine: Planetary-scale geospatial analysis for everyone. Remote Sensing of Environment, 202(3):18-27.


Han R Q, Li W J, Ai W X et al., 2010. The climatic variability and influence of first frost dates in Northern China. Acta Geographica Sinica, 65(5):525-532. (in Chinese)


Hao T Y, Wang S G, Shang K Z et al., 2010. Research progress and outlook of low temperature chilling injury in Northeast China. Meteorological and Environmental Research, 11(3):85-91.

Hellal F A, Abdelhamid T, 2013. Nutrient management practices for enhancing soybean (glycine max l.) production. Acta Biologica Colombiana, 18(2):3-14.

Jiang H, Hu H, Zhong R et al., 2019. A deep learning approach to conflating heterogeneous geospatial data for corn yield estimation: A case study of the US Corn Belt at the county level. Global Change Biology, 26(3):1754-1766.


Jones J W, Hoogenboom G, Porter C H et al., 2003. The DSSAT cropping system model. European Journal of Agronomy, 18(34):235-265.


Li W B, Song C X, Chang X C et al., 2019. Drought resistance evaluation of 20 soybean varieties under drought stress. Journal of Northeast Agricultural University, 50(4):1-10. (in Chinese)

Liu X, Jin J, Wang G et al., 2008. Soybean yield physiology and development of high-yielding practices in Northeast China. Field Crops Research, 105(3):157-171.


Liu D, Yu C L, Du C Y, 2016a. Dynamic monitoring rice delayed type chilling damage based on remote sensing in Northeast China. Transactions of the Chinese Society of Agricultural Engineering, 32(15):157-164. (in Chinese)

Liu X F, Ma S Q, Zhao H Q et al., 2016b. Evaluation on maize chilling injury based on WOFOST model in Hetao irrigation region in Inner Mongolia. Chinese Journal of Agrometeorology, 37(3):352-360.

Liu X F, Zhang Z, Shuai J B et al., 2012. Effect of chilling injury on rice yield in Heilongjiang Province. Acta Geographica Sinica, 67(9):1223-1232. (in Chinese)


Lobell D B, Thau D, Seifert C et al., 2015. A scalable satellite-based crop yield mapper. Remote Sensing of Environment, 164(11):324-333.


Ma S Q, Xi Z X, Wang Q, 2003. Risk evaluation of cold damage to corn in Northeast China. Journal of Natural Disasters, 12(3):137-141. (in Chinese)

Pang T H, Fang Z S, Zhao H K, et al., 1983. Crop Low Temperature Damage and Its Defense. Beijing: Agriculture Press. (in Chinese)

Pang Z K, 2015. Integration of remote sensing and crop growth model for regional low temperature impact monitoring early warning and yield estimation[D]. Hangzhou: Zhejiang University, 2015. (in Chinese)

Sang S P, 2013. Study on measures to cope with low temperature damage in different growth stages of soybean, Soybean Science & Technology, 5(1):53-54. (in Chinese)

Sinclair T R, Rawlins S L, 1993. Inter-seasonal variation in soybean and maize yields under global environmental change. Agronomy Journal, 85(2):406-409.


Tomar I S, Badaya A K, Vani D K et al., 2009. Impact of front line demonstrations on soybean by adoption of improved production technology. Soybean Research, 7:106-110.

Wang X L, Wang F T, Qiu G W, 1995. Application of optimizing theory to determining genetic parameters involved in CERES-Soybean model. Journal of Applied Meteorology, a01:49-54. (in Chinese)

Yang J, Ye W, Wang X et al., 2020. An Improved Method for the Identification of Soybean Resistance to Phytophthora sojae applied to germplasm resources from the Huanghuaihai and Dongbei regions of China. Plant Disease, 104(2):408-413.


Zhang J P, Wang C Y, Zhao Y X et al., 2012. Impact evaluation of low temperature to yields of maize in Northeast China based on crop growth model. Acta Ecologica Sinica, 32(13):4132-4138. (in Chinese)


Zhang S L, 2017. Analysis on the characteristics of chilling injury in Oroqen. Agricultural Engineering and Energy, 48(4):181-182. (in Chinese)

Zhou L W, 2017. Spatiotemporal characteristics of maize low temperature and cold damage and its effect on maize yield in Heilongjiang Province[D]. Harbin: Northeast Agricultural University. (in Chinese)

Zhu D W, Jin Z Q, 2008. Impacts of changes in both climate and its variability on food production in Northeast China. Acta Agronomica Sinica, 34(9):1588-1597. (in Chinese)