Integrating remotely sensed water stress factor with a crop growth model for winter wheat yield estimation in the North China Plain during 2008-2018

Wen Zhuo ,Shibo Fng,b,* ,Dong Wu ,Lei Wng ,Mengqin Li ,Jinsu Zhng ,Xinrn Go

a State Key Laboratory of Severe Weather,Chinese Academy of Meteorological Sciences,Beijing 100081,China

b Collaborative Innovation Centre on Forecast and Evaluation of Meteorological Disasters,Nanjing University of Information Science &Technology,Nanjing 210044,Jiangsu,China

c School of Geography and Earth Sciences,McMaster University,Hamilton,Ontario L8S 4L8,Canada

Keywords:WOFOST Evapotranspiration Drought Data assimilation Winter wheat yield

ABSTRACT Accurate estimation of regional-scale crop yield under drought conditions allows farmers and agricultural agencies to make well-informed decisions and guide agronomic management.However,few studies have focused on using the crop model data assimilation(CMDA)method for regional-scale winter wheat yield estimation under drought stress and partial-irrigation conditions.In this study,we developed a CMDA framework to integrate remotely sensed water stress factor(MOD16 ET PET-1)with the WOFOST model using an ensemble Kalman filter (EnKF) for winter wheat yield estimation at the regional scale in the North China Plain(NCP)during 2008-2018.According to our results,integration of MOD16 ET PET-1 with the WOFOST model produced more accurate estimates of regional winter wheat yield than open-loop simulation.The correlation coefficient of simulated yield with statistical yield increased for each year and error decreased in most years,with r ranging from 0.28 to 0.65 and RMSE ranging from 700.08 to 1966.12 kg ha-1.Yield estimation using the CMDA method was more suitable in drought years(r=0.47,RMSE=919.04 kg ha-1)than in normal years(r=0.30,RMSE=1215.51 kg ha-1).Our approach performed better in yield estimation under drought conditions than the conventional empirical correlation method using vegetation condition index(VCI).This research highlighted the potential of assimilating remotely sensed water stress factor,which can account for irrigation benefit,into crop model for improving the accuracy of winter wheat yield estimation at the regional scale especially under drought conditions,and this approach can be easily adapted to other regions and crops.

1.Introduction

Drought is an extreme climate phenomenon that severely impairs agricultural production [1,2].The North China Plain(NCP),which accounts for more than 60%of China’s wheat production (https://www.stats.gov.cn/),is usually affected by drought,and drought intensity,frequency and extent have increased in recent decades[3,4].Winter wheat is readily damaged by drought owing to its poor drought resistance[5],and more than 70%of the winter wheat area in the NCP is irrigated to ensure stable yield because only approximately 30%of the annual precipitation in this area occurs during the winter wheat growing season [6,7].For these reasons,accurate estimation of winter wheat yield under drought stress and partial irrigation conditions at large regional scale is crucial for optimizing agricultural water management and ensuring food security [8].

Existing regional-scale crop yield estimation methods(e.g.,statistical method,crop model method,or data assimilation method)have limitations to some extent.The conventional statistical method is usually applied by fitting a linear or nonlinear relationship to drought indicators and census yield data,and over 70 drought indices have been developed based on meteorological data(e.g.,rainfall,air temperature,air humidity or solar radiation) or remote sensing data(e.g.,visual light or thermal infrared or microwave radiation) [7,9-12].However,this approach offers poor universality for diverse crops and geographical regions,potential interactions between crop and environment are always overlooked,and it is hard to extend irrigation information to a large regional scale [13-17].Strictly speaking,machine learning (ML)and deep learning (DL) are also statistical methods.Recently,they have provided state-of-the-art results in regional crop yield or agricultural parameter estimation[18-21].Researchers usually use ML and DL models for crop yield estimation by integrating multisource data including remote sensing data,climate variables,and soil properties.However,the limited process-based interpretation of ML and DL models weakens the interpretability and traceability of crop yields [18].

Crop growth model(CGM)approaches such as WOFOST,DSSAT,APSIM,STICS,and AquaCrop,are powerful tools for simulating the soil-plant-atmosphere continuum [1,22-26].The superiority of CGM lies in interpreting crop growth mechanistically and characterizing the various crop and soil variables dynamically at the site scale with a daily time step [27].However,the accuracy of CGM depends strongly on local weather,soil,crop,and management strategies.Owing to the difficulties in regionalization of the model input parameters,such as crop growth parameters,soil characteristics,and irrigation information,crop models often use parameters calibrated in a few locations to represent a large spatial extent.Thus,large errors may be introduced when CGM is applied over a large area,as cultivars,management practice,and environments have high spatial heterogeneity [15,28].

The data assimilation method has been gaining increasing recognition as an effective approach to expanding the site-scale crop model to the regional scale by integrating remotely sensed crop or soil features with a crop model[29],permitting more accurate dynamic simulation of crop growth at the regional scale[15,30-34].The general workflow of a crop model data assimilation (CMDA) framework can be summarized as follows: remotely sensed crop parameters are assimilated into the crop growth model using data assimilation algorithms,and the crop model simulation accuracy of output variables (such as crop yield,leaf area index(LAI),and evapotranspiration(ET))can be improved by optimizing the relevant variables of the crop model.The most commonly used assimilation variable of CMDA framework is LAI[8,14-16],which is an important comprehensive parameter for crop growth monitoring,and it can reflect the stress and irrigation information to some extent.However,remotely sensed LAI is usually retrieved by vegetation index,which shows a time-lagged response to drought [35,36] and thus may not suitable for crop yield estimation under drought conditions.Consequently,the variables that are sensitive to water stress should be considered in the CMDA scheme for crop yield estimation under drought stress and partial irrigation conditions.

In recent years,soil water balance-associated variables(such as soil moisture (SM) and evapotranspiration (ET)) have been assimilated into crop models for crop yield estimation.Hu et al.[26]assimilated LAI and SM into the Soil Water Atmosphere Plant(SWAP) model for improving sugarcane growth simulation under diverse water stress conditions using three data assimilation approaches at the site scale,and the ensemble Kalman filter(EnKF)method most accurately estimated SM,LAI development and sugarcane yield.Zhuo et al.[33] assimilated remotely sensed soil moisture time series retrieved from Sentinel-1 and Sentinel-2 into the WOFOST model for regional-scale winter wheat yield estimation.The WOFOST model simulation process under water-limited mode was optimized and the accuracy of winter wheat yield estimation was also improved accordingly.Vazifedoust et al.[37]obtained the significant improvement in accuracy of winter wheat yield estimation in a small district in Iran by assimilating LAI and relative ET into the SWAP model using a constant-gain Kalman filter.

To our knowledge,no studies have used the CMDA method to assimilate ET-based water stress factor,which can contain crop irrigation information,into a crop model for crop yield estimation under drought stress and partial-irrigation conditions.Application of the CMDA method for large regional-scale crop yield estimation under drought conditions is rare,especially in China,although some encouraging results have been achieved at the site or small regional scales.

Therefore,the objective of this study was to estimate winter wheat yield at the regional scale in the NCP,which contains 350 counties,and to improve yield estimation accuracy by integrating remotely sensed water stress factor with the WOFOST model.Specifically,we conducted this research through: (a) evaluating the MOD16 ET PET-1for water stress diagnose at both site and regional scales,(b) assimilating remotely sensed water stress factor (MOD16 ET PET-1) into the WOFOST model for winter wheat yield estimation at the regional scale in the NCP during 2008-2018 using the EnKF method,and (c) comparing the yield estimation abilities of the CMDA and an empirical statistical approach using vegetation condition index (VCI) under drought conditions.The results are expected to strengthen our understanding for crop yield estimation under drought stress and partial irrigation conditions in the NCP.

2.Materials and methods

2.1.Study area

The North China Plain is located in northern China with an area of 4 × 105km2(Fig.1).It covers Beijing,Tianjin,Shandong,and most areas of Hebei,Henan,Anhui,and Jiangsu provinces.The climate of the NCP is a typical temperate monsoon climate with a mean annual precipitation of 472.7-889.2 mm and a mean temperature of 12.8-14.9 °C [27].The NCP is dominated by a typical double cropping system of rotational winter wheat and summer maize cultivation.Winter wheat in this area is usually planted from late September to early October and harvested in late May to June.In this study,we focused on the growing season of winter wheat from 2008 to 2018.The winter wheat region of the NCP comes from Huang’s research [15] that the spatial resolution is 30 m,we resampled the pixels to 500 m and selected winter wheat pixels with 85% purity.Although most winter wheat area on the NCP is generally irrigated and fertilized,drought still frequently occurs during the winter wheat growth period and affects wheat yields.

Fig.1.Study area.Green represents the winter wheat region and the blue triangles represent weather stations.

2.2.Datasets

The Moderate Resolution Imaging Spectroradiometer (MODIS)8-day ET product (MOD16A2) and 16-day normalized difference vegetation index(NDVI)product(MOD13A1)(https://ladsweb.modaps.eosdis.nasa.gov/) with 500-m spatial resolution were collected from January to June of 2008 to 2018.The MOD16A2 ET product was used to generate water stress factor (ET PET-1),and the MOD13A1 NDVI product was used for Vegetation Condition Index (VCI) calculation Eq.(1) [38,39],

where NDVIi，jrepresents the NDVI value for pixel i at time j and NDVIi，maxand NDVIi，minrepresent the long time series maximum and minimum NDVI for pixel i.

The Standardized Precipitation Evapotranspiration Index (SPEI)(https://spei.csic.es/database.html) was used to perform spatial comparison with ET PET-1and VCI.Monthly SPEI-6 data with 0.5° were collected from 2008 to 2018 during the winter wheat growing season.

Half-hourly eddy covariance data from the flux tower in the Yucheng station during 2008-2010 were compared with WOFOST simulated ET and MODIS ET data.Because the observed data from the flux tower is latent heat flux,the method of Huang et al.[40] was used to convert half-hourly latent heat flux to daily ET.

The input data for the WOFOST model include weather,crop,soil and management parameters.The WOFOST weather parameters comprise six elements (irradiation,early morning vapor pressure,maximum temperature,minimum temperature,wind speed and precipitation).Daily weather data for 2008-2018 with spatial resolution of 0.1°were obtained from the China Regional Surface Meteorological Elements Dataset produced by the National Tibetan Plateau Data Center (TPDC,https://data.tpdc.ac.cn/zh-hans/data/8028b944-daaa-4511-8769-965612652c49/?q=)and preprocessed to the WOFOST weather input format.The Chinese soil database(https://www.soil.csdb.cn) was used to derive soil moisture content at the wilting point(SMW),in saturated soil(SM0),and at field capacity (SMFCF).Some crop parameters,including the day of emergence (IDEM),the day of flowering,the day of maturity,and cumulative temperature from emergence to anthesis and from anthesis to maturity (TSUM1/TSUM2),were collected and calculated from field measurements and agrometeorological stations.Other parameters were calibrated using field-measured data from agrometeorological stations,or set according to previous studies,or set as default values.Huang et al.[15]provided details of parameterization of the WOFOST model for winter wheat.

Official government statistics on winter wheat yields were obtained at a county level from the 2008-2018 statistical yearbook of each province in the NCP.Relative error was used to evaluate the accuracy of winter wheat yield prediction with and without data assimilation.The relative error was calculated by Eq.(2):

where Yieldsimrepresents the WOFOST simulated winter wheat yield and Yieldstsrepresents official statistical winter wheat yield.

The 11 years(2008-2018)were divided into two parts:drought years(2009,2011 and 2014)and normal years(2008,2010,2012-2013,and 2015-2018),based on the drought-affected areas(statistic from the China National Bureau of Statistics).Winter wheat in the NCP is irrigated and water stress may not occur frequently,especially in normal years.However,winter wheat may suffer from water stress in drought years,which may lead to differences in crop growth conditions.Therefore,classifying years by drought was expected to lead to better estimates of the effect of water stress on crop yield.

2.3.WOFOST model

The WOFOST model was employed as the base model for daily winter wheat growth simulation in this study[41].The model is driven by a set of meteorological,crop,soil and management parameters.The major processes are phenological development,CO2assimilation,transpiration,respiration,partitioning of assimilates among various organs,and dry matter formation.A detailed description of the WOFOST model can be found at (https://www.wur.nl/en/Research-Results/Research-Institutes/Environmental-Research/Facilities-Tools/Software-models-and-databases/WOFOST.htm).There are three modes in the WOFOST model: a potential mode,a water-limited mode,and a nutrient-limited mode.The waterlimited mode was used,and irrigation was added to the model following Wang[42].In the water-limited mode,water stress factor(WSF)and actual daily CH2O assimilation rate(GASS)are defined in Eqs(3)and(4),respectively:

where Tarepresents actual transpiration,Tprepresents potential transpiration,and PGASS represents potential daily CH2O assimilation rate.However,the MODIS ET product only provides ET and potential evapotranspiration(PET)data.Because separation of transpiration flux from remotely sensed ET flux is difficult,we followed Vazifedoust’s strategy[37]of using ET PET-1as the remotely sensed WSF.

In this study,we assumed that the dominant winter wheat cultivar was planted in the NCP and that changes in winter wheat cultivar,soil properties,and management during 2008-2012 and 2013-2018 were negligible.Accordingly,the WOFOST model in 2008-2012 and 2013-2018 was calibrated based on the dominant winter wheat cultivar.The irrigation time in the model was fixed at April 1(Day of Year(DOY)91),May 1(DOY 121),and May 20(DOY 140).The irrigation quantity was determined empirically on each irrigation date.If the soil moisture content was lower than 90%of the SMFCF,SM was set equal to SMFCF;otherwise,no irrigation was applied on the date.Besides,three parameters (emergence date (IDEM),cumulative temperature from emergence to anthesis(TSUM1),and cumulative temperature from anthesis to maturity(TSUM2)) were regionalized with a spatial resolution of 500 m.First,winter wheat emergence date,flowering date,and maturity date images (Fig.S1A-C) were generated for each year based on agrometeorological station measurements using Thiessen polygons.Second,regionalized emergence,flowering and maturity date were used to calculate TSUM1 and TSUM2(Fig.S1E,F)based on the daily temperature data.

2.4.Data assimilation algorithm

The EnKF method was used to integrate WSF retrieved from MOD16 with the WOFOST model.The implementation of EnKF was based on our previous studies[8,33].The core part of the EnKF is the calculation of the Kalman gain matrix,

where y is the observation vector,ε and ν are Gaussian random error vectors with a mean of zero,and H is the observation operator for y and can be taken as an identity matrix in this study.A represents a linear state-transition model that links xtand xt-1,and in the CMDA system it represents the crop model.The forecast of xtat t=k is Gaussian with meanand error covariance,calculated as follows:

where f and a are indices of the prior and posterior estimates,respectively,I is the identity matrix,and K is the Kalman gain matrix,defined as

where Rt=kis the error covariance of the observation ensemble.For solving the Kalman gain,Houtekamer and Mitchell suggest calculatingdirectly from the ensemble members [43],rather than calculating each element of Eq.(10):

where Neis the number of ensemble members,n is a running index of ensemble member,andrepresents the ensemble mean calculated as Eq.(12).

We adopted the strategy of Huang et al.[8] incorporating an inflation factor E to solve the problem of ‘‘filter divergence” and enlarge K,and E is calculated by Eq.(14):

where r is a random value between 0 and 1,160 represents the total number of days from Jan.1 to maturity,and k represents the day number (from 1 to 160).

The ensemble number of the data assimilation scheme was set as 50,and WOFOST ensembles were generated by perturbing model input parameters by introducing 10% uncertainty using a Gaussian distribution based on previous studies [31,33],and the uncertainty of the MODIS ET product was accordingly set as 10%.The remotely sensed water stress factor was assimilated into the WOFOST model at the same time period (from April 1 to June 5)throughout the study area,because the LAI is generally greater than 2 m2m-2during this period.The data assimilation scheme was applied only when remotely sensed observations were available.

3.Results

3.1.Evaluation of ET PET-1 for water stress diagnosis at both regional and site scale

The spatial distribution of drought indicators during the winter wheat growing season of 2009,a typical drought year,is shown in Fig.2.From the distribution of SPEI we can see that meteorological drought occurs mainly in spring (from January to April),and that drought intensity decreased in May and June,which is generally consistent with the Yearbook of Meteorological Disasters in China.Compared to the meteorological drought index (SPEI),remotely sensed drought indicators (VCI and ET PET-1) monitor drought at a higher temporal and spatial resolution,and only one image per month was selected for analysis.VCI and ET PET-1showed similar drought change trends during the winter wheat growing season,but their performance differed in some local area.In western Shandong and northern Hebei,no water stress was observed in February and moderate drought in June by ET PET-1,whereas the opposite was observed with the VCI.Soil moisture at the weather-station level is also presented as an indicator of drought conditions.By this indicator,water deficiency occurred mainly in the western and central regions of the NCP.Overall,the drought indicators showed generally similar drought change trends and described drought well aside from local differences.

A comparison of the WOFOST simulated evaporation and transpiration is shown in Fig.3A.There were large differences betweenand ETa ETp-1during wintering stage (before March),and smaller differences after jointing stage(after April).was generally equal to ETa ETp-1when winter wheat LAI was greater than 2 m2m-2.Validation results are shown in Fig.3B and C: the coefficient of determination(R2)was 0.28 throughout the growing season (Fig.3B) and the correlation increased,with R2equal to 0.81,when LAI exceeded 2 m2m-2(Fig.3C).This finding confirms the reliability of using ET PET-1to characterize water stress in winter wheat,and the MODIS 16 ET PET-1data can then be assimilated into the WOFOST model for winter wheat yield estimation.

As shown in Fig.4,the WOFOST simulated ET showed close agreement with flux tower-measured ET,and the 8-day MODIS ET data generally captured the temporal variation of winter wheat ET.In general,these three ET data showed a good response to precipitation and irrigation,such that ET increased sharply when precipitation and irrigation occurred.Similar results were found in other years (Figs.S2 and S3).

3.2.Assimilation of MOD16 ET PET-1 into the WOFOST model

We assimilated MOD16 ET PET-1data into the WOFOST model using the EnKF algorithm during 2008-2018.Given that official government statistics of winter wheat yield are compiled at a county level,we aggregated the WOFOST simulated yield pixels with or without assimilation into the county scale using the ArcGIS Zonal Statistics toolbox.Fig.5 compares the mapped winter wheat yield simulated by the WOFOST model with and without data assimilation and official statistics at the county level.Overall,the WOFOST simulated yield with or without assimilation captured some of the spatial variation of winter wheat yield,and the simulated yields were generally between 5000 and 7000 kg ha-1,although inconsistent in some years,especially 2018.

Fig.2.Spatial distribution of four drought indicators in the North China Plain during winter wheat growing season in 2009.SPEI,standardized precipitation evapotranspiration index;VCI,vegetation condition index;ET PET-1,evapotranspiration divided by potential evapotranspiration;SM,soil moisture.

Fig.6 shows the validation results of the WOFOST simulated yield without and with EnKF assimilation at the county level using official statistics for each year.Table 1 shows detailed statistical validation results.These results indicate that although the WOFOST simulated yield without assimilation captured some of the spatial variation of winter wheat yield(Fig.5),it had a low correlation coefficient and large error with r ranging from 0.14 to 0.39 and RMSE ranging from 827.93 to 1742.87 kg ha-1.Compared with the WOFOST simulated yield without assimilation,the results with assimilation improved the accuracy of the winter wheat yield estimation,with the correlation coefficient increasing for each year and error decreasing in most years with r ranging from 0.28 to 0.65 and RMSE range from 700.08 to 1966.12 kg ha-1.Thus,assimilating MOD16 ET PET-1data into the WOFOST model reduced the uncertainty of model simulation and further increased the accuracy of yield estimation.

Table 1 Accuracy of simulated yield without and with assimilation at the county level for each year and three kinds of year combinations.

Validation results were generated for three categories:all years,drought years,and normal years(Fig.7).Higher r and lower RMSE were obtained when data assimilation was performed using remotely sensed water stress factor for all three categories.In particular,the accuracy of winter wheat yield estimates increased more markedly in drought years than in normal years.Fig.8 shows the distribution and relative error of winter wheat yield.Similarly,the winter wheat yield distribution clearly shows that the WOFOST simulated yield with assimilation is closer to the statistical results than that without assimilation in drought years (Fig.8B).Relative error (RE),was reduced in most years by integrating remotely sensed water stress factor with the WOFOST model,and the degree of reduction of RE was more pronounced in drought years(Fig.8D).

3.3.Comparison with yield estimation method using remotely sensed drought index in a typical drought year

Comparison of regional yield estimates using our approach and the conventional empirical correlation method was performed.Previous studies[2,44-46]have shown that the jointing,flowering,and grain-filling stages are highly sensitive to heat and drought stress.In this study,we chose VCI at DOY 105 and DOY 145,which correspond generally to the two sensitive growth stages,as the indicators for constructing the correlation with winter wheat yield.The statistical correlation model was first trained using fieldmeasured yield data from agrometeorological stations,and then the county-scale winter wheat yield was estimated.Fig.9 presents the correlations among winter wheat yield estimates at the county level in 2009,a typical drought year in the NCP.The county-scale yields include the WOFOST simulated yield without assimilation(YieldWOF)and with assimilation(YieldEnKF),yield estimation using vegetation condition index at the jointing-flowering stage(YieldVCI105)and the grain-filling stage(YieldVCI145),and official statistical yield (YieldSta).YieldEnKFperformed best among these county-scale yield estimations with a higher correlation coefficient(r=0.65)and a lower RMSE(RMSE=700.08 kg ha-1)in comparison with YieldSta.This finding further demonstrated that our approach integrating remotely sensed ET PET-1with the WOFOST model has higher potential for winter wheat yield estimation under drought conditions than the conventional correlation method.Comparison results for other years are provided in the Supplementary file(Figs.S4-S13).

Fig.3.Comparison of winter wheat evaporation and transpiration.(A)Comparison of evaporation and transpiration during the winter wheat growing season at the site scale,and (B) comparison between ETa ETp-1 and Ta during the full growing season and (C) when LAI >2 m2 m-2 (a sample assimilation unit in Hebei province).Ta,actual transpiration;Tp,potential transpiration;Ea,actual evaporation;Ep,potential evaporation;Eta,actual evapotranspiration;ETp,potential evapotranspiration;LAI,leaf area index.

Fig.4.Comparison of three evapotranspiration data sets during winter wheat growing season (a sample assimilation unit in Yucheng station).ETTower,eddy flux tower measured evapotranspiration;ETWOFOST,WOFOST simulated actual evapotranspiration;ETMODIS,MODIS evapotranspiration.

Fig.5.Comparison of the spatial patterns of winter wheat yield simulated by the WOFOST model with and without data assimilation and official statistics.(A) Data from 2008 to 2013.(B) Data from 2014 to 2018.

4.Discussion

Quantitative evaluation of crop yield under drought stress conditions in the NCP is desirable,given that drought is one of the major natural hazards in this area [7].However,previous studies did not well account for irrigation information when estimating crop yield at a large regional scale.Taking advantage of remotely sensed ET data in reflecting irrigation information over a large spatial extent,and the superiority of crop models for describing and modeling the timing and amount of water stress and crop water requirement [47],we integrated remotely sensed water stress factor(MOD16 ET PET-1)and the WOFOST model using an EnKF algorithm for winter wheat yield estimation at the regional scale in the NCP during 2008-2018.

In the WOFOST model,winter wheat suffers water stress if root water uptake cannot meet potential transpiration (Tp),and leaves begin to close stomata to reduce water loss [48].Relative transpiration,defined as the ratio of actual and potential transpiration(Ta),is the indicator of crop water stress in the WOFOST model,and it is essentially different from ETa ETp-1.In our study,was generally equal to ETa ETp-1when LAI exceeded 2 m2m-2with R2and slope equal to 0.81 and 0.75,respectively(Fig.3C).This might be because winter wheat transpiration is generally low and soil evaporation accounts for a large proportion of ET when ground cover is low(LAI<2 m2m-2).Our result is consistent with those of Feddes et al.[49]and Vazifedoust et al.[37],and confirms the reliability of coupling MOD16 ET product with the WOFOST model with ETa ETp-1as the assimilation variable.

Fig.6.Comparison of the WOFOST simulated winter wheat yield without (A) and with (B) assimilation and official statistical yield at the county level during 2008-2018.

Fig.7.Validation of the WOFOST simulated winter wheat yield at the county level without assimilation for three kinds of year combinations:(A)all years,(B)drought years,and (C) normal years;and with EnKF assimilation for three kinds of year combinations: (D) all years,(E) drought years,and (F) normal years using official statistical yield.

Assimilating remotely sensed water stress factor into the WOFOST model generally produced more accurate winter wheat yield estimates at the regional scale in each year than open-loop simulations (Fig.6B),which is consistent with that of Vazifedoust et al.[37].This demonstrated that the WOFOST model could take advantage of the benefits of the MOD16 ET product,which is a comprehensive characterization of crop transpiration and soil evaporation after precipitation,irrigation and fertilization.This finding means that remotely sensed water stress factor can partially compensate for the deficiency of the WOFOST model in considering irrigation,especially at the large regional scale.Validation results for three kinds of year combinations showed that our approach is more suitable in drought years (with higher r and lower RMSE) than in normal years (Figs.7,8),which has similar conclusion with Hu et al[26],who obtained better yield estimates under water-stressed than under well-watered conditions.Overall,our results indicated that the MOD16 ET product has the potential for optimizing the water stress information in the WOFOST model and further improving the accuracy of winter wheat yield estimation.

Fig.8.Winter wheat yield distribution for (A) all years,(B) drought years,and (C) normal years;and (D) relative error of WOFOST simulated yield without and with data assimilation during 2008-2018 and three kinds of year combinations.

Agricultural drought is a complex natural phenomenon that is related to many factors including weather systems,topographic conditions,crop types,and soil characteristics[7,46,50].Many previous studies[27,38,51,52]have constructed drought indices using various datasets(such as of meteorological elements,hydrological variables or spectral characteristics of vegetation) for quantitative evaluation of drought.VCI is a widely used remotely sensed drought index,and many successful results have been achieved in crop yield estimation using VCI [53-55].In this study,comparison among county-level yield estimates e.g.,YieldWOF,YieldEnKF,YieldVCI105,YieldVCI145and YieldSta(Fig.9) showed that our approach offers substantial advantages over the conventional correlation method using drought indices.This may benefit from the advantages of crop model that it is an effective means of daily crop growth simulation and permit explaining the growth and development of crops.For example,when crops suffer from drought,the crop model will make changes in water balance and vegetation physiology to respond to water stress in crop growth.However,remotely sensed drought indices often reflect instantaneous growth state of crops,and cannot accurately characterize the information of vegetation suffering from water stress [46].Specifically,NDVI,which was used for calculating VCI in this study,usually shows a time-lagged response to drought [35,36].In addition,the assimilation variable used in this study is generated from ET,which represents mass and energy exchange between vegetation and the atmosphere.Previous studies [9,56] indicated that ET responds faster to water stress (stomatal sensitivity) than do vegetation indices.

EnKF has been an effective approach in CMDA framework for crop yield estimation [8,30,31,33].However,the EnKF algorithm has the defect of filter divergence [31,57,58],which usually tends to reject observations in the late period of crop growth.In the present study,an inflation factor was constructed followed the strategy of Huang et al.[8] to enlarge the variance of the forecast ensemble and further reduce the effect of filter divergence.Although promising results have been achieved in this study that winter wheat yield estimation accuracy was essentially improved by integration of remotely sensed water stress factor with the WOFOST model using the EnKF algorithm,especially in drought years (Figs.6B,7),EnKF has its own limitation: it assumes that both observation and model errors are Gaussian,whereas the actual situation is much more complex.In future research,the Particle Filter (PF) algorithm,which assumes non-Gaussian errors,may be used in place of EnKF in our CMDA scheme for improving crop yield estimation [59,60].

In the NCP,the high heterogeneity of the winter wheat cultivation area poses a great challenge for winter wheat yield estimation by the CMDA method[29].In this study,we used remotely sensed ET PET-1data with 500 m spatial resolution,a size may greater than those of some fields,although we chose only winter wheat pixels with purity>85%.Pixels with coarse spatial resolution pixels contain mixed ground information,which may cause more uncertainty when integrated into the crop model.Meanwhile,remotely sensed pixels with high heterogeneity may cause greater scale mismatch between remotely sensed variables and crop model simulations.Therefore,in future studies,remote sensing data with high spatial resolution should be considered for integration in the CMDA framework to reduce the scale disparity.As the accessible of remote sensing data with a finer spatial resolution,promising results may be achieved by our approach using water stress factor data with 10-30 m spatial resolution in the NCP even for crop regions with fragmented fields.

In this study,the WOFOST model was calibrated under the assumptions that the dominant winter wheat cultivar was planted in the NCP and that winter wheat characteristics,soil properties,and management measures did not vary much throughout the study area.Moreover,only the IDEM,TSUM1,TSUM2,SMW,SM0 and SMFCF of the NCP were regionalized,and remaining parameters were calibrated using the dominant winter wheat cultivar.However,the choice of cultivar,soil characteristics,and management conditions vary within the NCP,possibly leading to large differences in crop growth and even yield[8].In future studies,studyarea partitioning,which accounts for spatial differences in crop cultivar,soil properties,and management measures by dividing a large area into subregions,should be considered for more accurate yield estimation [61].In view of the deficiency of the WOFOST model in irrigation [30,42],the SWAP model,which is an agrohydrological model for water flow,heat and solute transport as well as crop growth [62],could be incorporated for its superiority in characterizing water balance.Previous studies [26,37,40] have focused on assimilating multiple variables into the SWAP model for crop yield estimation and successful results have been achieved.

Fig.9.Correlation matrix of winter wheat yield at the county level in 2009 estimated by the WOFOST model without assimilation (YieldWOF),the WOFOST model with assimilation(YieldEnKF),vegetation condition index at the jointing-flowering stage(YieldVCI105),vegetation condition index at the grain-filling stage(YieldVCI145)and official statistical yield (YieldSta).

5.Conclusions

We assimilated remotely sensed water stress factor (ET PET-1)into the WOFOST model using the EnKF method for winter wheat yield estimation in the NCP during 2008-2018.The results illustrated that MOD16 ET PET-1product could characterize water stress at regional and site scales and could be further used to optimize the WOFOST model for winter wheat yield estimation.Yield estimation results showed that our approach improved winter wheat yield estimates and their spatial patterns at the regional scale over a relatively long term.In particular,assimilating ET PET-1into the WOFOST model led to higher accuracy (higher r and lower RMSE)of yield estimation in drought years than in normal years.The results also confirmed that our approach is superior to the conventional empirical correlation method using remotely sensed drought indices.Our study demonstrated the superiority of integrating remotely sensed water stress factor with a crop model for winter wheat yield estimation at the regional scale under drought-stress conditions and has great potential for application to other crops and regions.

CRediT authorship contribution statement

Wen Zhuo:Formal analysis,Visualization,Writing -original draft.Shibo Fang:Conceptualization,Funding acquisition,Writing-review&editing.Dong Wu:Data curation.Lei Wang:Visualization.Mengqian Li:Validation.Jiansu Zhang:Validation.Xinran Gao:Methodology,Formal analysis.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This research was supported by FengYun Research Plan (FYAPP-2021.0301),National Key Research and Development Program of China(2019YFC1510205),and National Natural Science Foundation of China (42075193).

Appendix A.Supplementary data

Supplementary data for this article can be found online at https://doi.org/10.1016/j.cj.2022.04.004.