ZHANG Sha,YANG Shan-shan,WANG Jing-wen,WU Xi-fang,Malak HENCHIRI,Tehseen JAVED, ,ZHANG Jia-hua#,BAI Yun#

1 Research Center for Space Information and Big Earth Data,College of Computer Science and Technology,Qingdao University,Qingdao 266071,P.R.China

2 Aerospace Information Research Institute,Chinese Academy of Sciences,Beijing 100094,P.R.China

3 School of Surveying and Land Information Engineering,Henan Polytechnic University,Jiaozuo 454000,P.R.China

4 Department of Environmental Sciences,Kohat University of Science and Technology,Kohat 26000,Pakistan

Abstract

Accurate estimation of regional winter wheat yields is essential for understanding the food production status and ensuring national food security.However,using the existing remote sensing-based crop yield models to accurately reproduce the inter-annual and spatial variations in winter wheat yields remains challenging due to the limited ability to acquire irrigation information in water-limited regions.Thus,we proposed a new approach to approximating irrigations of winter wheat over the North China Plain (NCP),where irrigation occurs extensively during the winter wheat growing season.This approach used irrigation pattern parameters (IPPs) to define the irrigation frequency and timing.Then,they were incorporated into a newly-developed process-based and remote sensing-driven crop yield model for winter wheat (PRYM–Wheat),to improve the regional estimates of winter wheat over the NCP.The IPPs were determined using statistical yield data of reference years (2010–2015) over the NCP.Our findings showed that PRYM–Wheat with the optimal IPPs could improve the regional estimate of winter wheat yield,with an increase and decrease in the correlation coefficient (R) and root mean square error (RMSE) of 0.15 (about 37%) and 0.90 t ha–1 (about 41%),respectively.The data in validation years (2001–2009 and 2016–2019) were used to validate PRYM–Wheat.In addition,our findings also showed R (RMSE) of 0.80 (0.62 t ha–1) on a site level,0.61 (0.91 t ha–1) for Hebei Province on a county level,0.73 (0.97 t ha–1) for Henan Province on a county level,and 0.55 (0.75 t ha–1) for Shandong Province on a city level.Overall,PRYM–Wheat can offer a stable and robust approach to estimating regional winter wheat yield across multiple years,providing a scientific basis for ensuring regional food security.

Keywords: approximating irrigations,process-based model,remote sensing,winter wheat yield,North China Plain

1.Introduction

It is crucial to accurately estimate regional crop yields to evaluate national food security,make suitable food trade strategies,and develop sustainable agriculture (Lobelletal.2003; Huangetal.2015a; Lietal.2017; Chenetal.2018b; Wang Yetal.2019).China is the world’s largest wheat-producing country,and more than half of its production comes from the North China Plain (NCP) (Renetal.2008; Fangetal.2018; Zhengetal.2020).Therefore,it is essential to accurately estimate the regional winter wheat yield over the NCP to ensure food security and facilitate the development of national agricultural policies for China (Franchetal.2015; Huangetal.2019).

The model approach is effective in estimating regionalscale crop yields.Crop growth models (CGMs) are widely used to simulate crop yields (Liuetal.2012,2015; Ittersumetal.2013; Lietal.2014; Abi Saabetal.2021).Calibrated CGMs can effectively estimate yields for a certain crop in a certain region (Liuetal.2012,2015; Lietal.2014; Abi Saabetal.2021).However,this approach is dependent on many parameters,such as precise weather data,soil property data,variety data,and field management data; it thus has limited application over large-scale regions with scarce data (Becker-Reshefetal.2010; Wartetal.2013; Zhangetal.2013; Amarasinghaetal.2015).In addition,the calibration of CGMs involves large uncertainties over a large regional scale (Lietal.2014).Thus,local calibrations are required to further improve the estimations by CGMs.

Alternatively,remote sensing (RS)-driven crop yield models (i.e.,RS crop yield models) are another model that can be used to estimate crop yields.These models are primarily driven by satellite-retrievedin-situcropland surface parameters.RS data reflect temporally and spatially continuous crop physiological and ecological parameters,such as normalized difference vegetation index (NDVI),leaf area index (LAI),and phenological information (Lobelletal.2013; Changetal.2016; Lietal.2017; Xieetal.2018; Wuetal.2019).Thus,compared with CGMs,the RS crop yield models can be easily adapted to regional areas because they are less sensitive to meteorological data errors and field management information (Curneletal.2011; Maetal.2013; Huangetal.2015b).RS crop yield models include empirical models,light use efficiency (LUE) models,and process models based on RS (process-based RS models).Empirical models present a statistical relationship between yields and relevant factors,such as RS vegetation indices (VIs) and meteorological variables (Renetal.2008; Sonetal.2013; Franchetal.2015).Thus,these models are short in generalization ability over time and space (Cheng and Meng 2015; Chenetal.2016; Chaoetal.2019).In comparison,LUE models are typically straightforward to implement as they do not require numerous inputs,but their key parameter characterizing photosynthesis of crops,i.e.,maximum light use efficiency (εmax),is generally formulated empirically (Xiaoetal.2004; Yuanetal.2006; Zhuetal.2006).

Process-based RS models for regional crop yield estimation generally incorporate a mechanism-based method to simulate various plant physiological processes,such as photosynthesis,autotrophic respiration,and transpiration (Juetal.2010; Wangetal.2011; Yaoetal.2015,2016; Wangetal.2020; Zhangetal.2021).Compared to empirical models and LUE models,the process-based RS models can better characterize ecosystem biophysical traits.However,the existing process-based RS models consider crops to have only one source of available water,i.e.,precipitation (Chenetal.1999; Huangetal.2018),resulting in underestimations for irrigated croplands in water-limited regions (Zhangetal.2018a).Irrigation primarily affects the water status of winter wheat over the NCP.During the winter wheat growing season,the amount of precipitation meets only 25–40% of the water demand of winter wheat over the NCP.Thus,supplementary irrigation is crucial for maintaining winter wheat growth and yield (Shenetal.2013; Fangetal.2017; Wang and Li 2018; Wangetal.2019).However,commonly-used RS VIs can only partly characterize the water status of crops with hysteresis (Yuanetal.2015; Baietal.2017).In this regard,introducing soil moisture (SM) data retrieved from the microwave bands into the process-based RS models can potentially represent the effect of irrigation.However,the currently available SM data have limited applications due to coarse spatial resolution,a considerable amount of missing data,and insufficient representation of SM in deep layers (Zhu 2019; Zhuetal.2019).Some studies have proposed methods to retrieve high-resolution SM data,but their methods were generally locally applicable rather than over a large region (Ouaadietal.2021; Saidetal.2021; Guoetal.2022).Therefore,using only RS spectral VIs or SM data to represent crop water status may introduce additional uncertainty in estimating regional winter wheat yields.

In the existing studies,accurate information about irrigation amount and frequency was generally obtained using RS data-driven approaches,climate data,and observed data (Zhangetal.2007; Chenetal.2018a; Ouaadietal.2021).For example,Zhangetal.(2007) used the surface energy balance system to compute land surface evapotranspiration and derive the irrigation amount for winter wheat during the 2010–2011 growing season.Chenetal.(2018a) used the RS-derived timeseries greenness index to detect irrigation frequency and timing in four winter wheat growing seasons.Thus,irrigation information,namely,irrigation amount,frequency,and timing,retrieved following the approaches reported by Zhangetal.(2007) and Chenetal.(2018a) are only useful for specific growing seasons and regions.Thus,this hampers their applications over a large region and multiple years.Some other studies have reported a stable pattern of irrigation for winter wheat over the NCP,i.e.,irrigation usually occurs at certain growing stages,such as returning green,jointing,flowering,and so on (Wangetal.2013; Yangetal.2014).Moreover,RS timeseries data were used to estimate these growing stages of winter wheat (Renetal.2019; Wuetal.2019; Zhangetal.2019).Thus,a set of parameters can be used to depict the irrigation patterns for winter wheat over the NCP,facilitating more accurate estimates of winter wheat yield on a regional scale.

In this study,we propose a new method to characterize the irrigation pattern of winter wheat and incorporate this method into a process-based RS model to improve the regional estimates of winter wheat yields across multiple years and over the NCP.The objectives of this paper include: (1) proposing a new approach in which a set of irrigation pattern parameters (IPPs) define the irrigation frequency and timing of winter wheat to approximate irrigations of winter wheat over the NCP,(2) developing a process-based and RS-driven crop yield model for winter wheat (PRYM–Wheat) that incorporates the IPPs,and (3) evaluating the PRYM–Wheat model to estimate winter wheat yield over the NCP across multiple years.

2.Material and methods

2.1.Study area

The study area,NCP,is mainly located between 32.0–42.0°N and 111.0–121.0°E.It comprises Beijing Municipality,Tianjin Municipality,Hebei Province,Shandong Province,and Henan Province (Fig.1).The NCP accounts for more than 50% of total wheat production in China given that it is one of the most productive and intensively cultivated agricultural regions in the country (Wangetal.2012; Fangetal.2017).The study area generally experiences a temperate monsoonal climate and has an average annual precipitation of 500–600 mm,with 70% of rainfall concentrated between late June and September (Zhangetal.2017; Wuetal.2019).The main food crops in this region are winter wheat and summer maize.During the winter wheat growing season,irrigation is an effective approach to ensure winter wheat growth,as the amount of rainfall in this season usually does not meet the demand for water for wheat growth.

Fig.1 Spatial distribution of winter wheat in the study area (take winter wheat distribution in the year 2010 as an example,referring to Zhang et al. (2018b)),agro-meteorological sites,and collected observed winter wheat yield sites over the North China Plain (NCP).

2.2.Methodology

Configuration of PRYM–WheatThe configuration of PRYM–Wheat is shown in Fig.2.The model PRYM–Wheat runs on a daily step and consists of three submodules: a water balance module,a photosynthesis module,and a yield simulation module.The three-layer soil water balance module (Appendix A) developed by Baietal.(2018) is introduced into the PRYM–Wheat model.We improved the irrigation strategy in the three-layer soil water balance module proposed by Baietal.(2018) and incorporated it into the soil water balance module.The photosynthesis module and yield simulation module are based on previous studies,and their details are presented in Appendices B–D.

Fig.2 Configuration of the developed PRYM–Wheat.NDVI,normalized difference vegetation index; LAI,leaf area index; LULC,land use and land cover; ET,evapotranspiration; GPP,gross primary productivity; NPP,net primary productivity.

An overview of the water balance moduleThe threelayer soil water balance module introduced by Baietal.(2018) consists of eight sub-processes of water balance,i.e.,snow melting,canopy interception,saturated soil evaporation,canopy transpiration,moist soil evaporation,percolation of soil water above the field capacity,soil water input,and irrigation approximation.With respect to irrigation approximation,the module assumes that the available soil water content (SWC) of the second layer will be maintained at a certain level when the SWC sum of the shallowest two layers falls below a certain level.These certain levels are reflected using two empirical coefficients related to management levels.

The soil water balance module includes three outputs: daily evapotranspiration (ET),daily water stress factor,and daily multi-layer soil water contents.These three factors are mutually restricted in the water balance module,and the daily water stress factor (fw) is used to constrain the photosynthesis rate in the photosynthesis module.The computation of daily ET and multi-layer SWC is presented in Appendix A of the current study and Baietal.(2018).fwis computed using the following equations (Baietal.2018; Zhangetal.2018a).

where superscriptldenotes soil layer index;l_wstdenotes the index of the wettest soil layer.w,SWSF,andθwrepresent the weight,soil water stress factor,and wilting point of the soil layer,respectively.SWC(l)denotes the SWC in soil layerl,in mm.θmdenotes the SWC in the case of maximum stomatal conductance.Yrepresents the proportion of root biomass within a soil layer (0–1).d(l) denotes the depth of the lower boundary of soil layerl,in cm.βdenotes the root extinction coefficient,set at 0.961 for crops (Baietal.2018).

Approach to approximating irrigation in PRYM– WheatThe approach to approximating irrigation proposed by Baietal.(2018) was based on multiple crops over the globe.This approach is not entirely suitable for our research study as we focused only on winter wheat over NCP.Thus,in this study,we developed a new approach to approximating irrigation (Fig.3) using the information on the growing stages of winter wheat.In this study,the whole life of winter wheat is categorized into ten growing stages: sowing,emergency,overwintering,returning green,jointing,booting,heading,flowering,filling,and maturity.However,we only focused on the growing stages from returning green to maturity to simplify the computation process,as these stages are critical for photosynthetic organic matter forming (Renetal.2006; Wangetal.2019).The new approach to approximating irrigation used growing stages from returning green to filling (Fig.3).

Fig.3 Proposed new approach to approximating irrigations in PRYM–Wheat.IPPs,irrigation pattern parameters; RS,remote sensing.SOSRS,EOSRS,and LOSRS,RS-based returning green date,maturity date,and length from returning green to maturity,respectively; DOY,the day of the year; length,the length of a growing stage,and it is set as ten days in this study; SWCsum,the sum of SWC in three soil layers; FC,filed capacity.Observed phenology includes six observed growing stages,i.e.,Si, obs (i=4,5,6,…,9),representing returning green,jointing,booting,heading,flowering,and filling dates,respectively.Si, RS (i=4,5,6,…,9) denotes the corresponding growing stage estimated based on RS phenology and observed phenology.j values are equal to i values according to IPPs; for example,if the irrigation frequency and irrigation timing in IPPs are twice at returning green and booting stages,then j=i=4 and j=i=6.

The Si,RSin the inequality in Fig.3 is calculated as follows:

where Si,RSand Si,obsrepresent theith growing stage of winter wheat estimated using RS phenology and observed phenology and phenology,respectively;i=4,5,6,…,9,which denotes returning green,jointing,booting,heading,flowering,and filling stages of winter wheat.SOSRS,EOSRS,and LOSRSdenote the start of the growing season,the end of the growing season,and the length of the growing season estimated using RS data,referring to the Day of Year (DOY) of returning green date,the maturity date of winter wheat and days from returning green to maturity,respectively.SOSobs,EOSobs,and LOSobsare observed returning green date,maturity date,and length from winter wheat returning green to maturity.kiis defined as the ratio of the difference between Si,obsand SOSobsto LOSobs.The steps of determination ofkiare listed as follows:

(1) Collect the original argo-meteorological station data and extract the observed phenology information (in date format) from the original data for each site.

(2) Convert the phenology information from the date format to DOY.

(3) Calculate LOSobsusing returning green DOY (as SOSobs) and maturity DOY (as EOSobs) as in eq.(7).

Calculateki(i=4,5,6,…,9) using eq.(6).

Thekiat the returning stage is set as an extreme minimum value of 0.01 to avoid null or irregular results during the calculation ofki.The calculatedkiparameters are listed in Table 1.

Table 1 Determined ki for each stage in this study1)

Modeling cases with different lPPsIn this study,IPPs indicate the irrigation pattern defining the frequency and timing.It should be noted that the defined IPPs are not the real irrigation frequency and timing of winter wheat.In this study,we design multiple modeling cases running with different IPPs to determine IPPs for the whole study area.Each set of IPPs represents an irrigation scenario,which is designed according to the following principles:

(1) I-0 with no irrigation is used as a control experiment.

(2) Irrigation occurs only at certain growing stages,such as returning green,jointing,or other growing stages before maturity.

(3) Irrigation will not occur in the two neighboring growing stages.

Based on these principles,21 irrigation scenarios were designed (Table 2) to drive PRYM–Wheat.These modeling cases were categorized into four groups based on the different irrigation scenarios with different IPPs,I-0,I-1,I-2,and I-3,indicating zero,once,twice,and three times irrigations,respectively.In all these cases,irrigation is implemented using the approach mentioned above.

Table 2 Modeling cases under different irrigation scenarios with different irrigation pattern parameters1)

PRYM–Wheat was run to simulate winter wheat yield in reference years (2010–2015) under each modeling case.The optimal IPPs were determined by comparing the correlation coefficient (R),and the root mean standard error (RMSE) between simulated and official statistical yields from statistical yearbooks.

2.3.Data

RS data(1) MODIS NDVI/EVI.MODIS NDVI/EVI from MOD13A2 (2001–2019) and MYD13A2 (2003–2019) products were used in this study.Both were downloaded from https://search.earthdata.nasa.gov/search and have a 1-km spatial resolution and a 16-day temporal resolution.We obtained NDVI data with an 8-day temporal resolution using the combination of the two NDVI products with a 16-day temporal resolution.

(2) LAI data.The MODIS LAI products were not used in the study because they were mostly underestimated for crops (Yangetal.2006,2015; Fuetal.2016).Instead,we built an empirical LAI estimation (eq.(9)) based on the site-measured LAI data of two years and RS EVI.These measured LAI data,observed at the Yucheng Experiment Station under the China Ecosystem Research Network during 2000–2001,were obtained by digitizing the first figure in the research of Wangetal.(2007).The RS EVI values were extracted from MODIS EVI images.These data were linearly interpolated to a temporal scale consistent with the observed LAI.

where LAI denotes the simulated LAI,and EVI denotes the corresponding EVI.This study used this equation to calculate the LAI of winter wheat.The determination coefficient (R2) between the measured and estimated LAI is 0.75.

(3) Land use and land cover data.In this study,MODIS land use and land cover data,MCD12Q1 product (500 m,yearly) were also used to drive PRYM–Wheat.The density of impervious pixels over the study area was estimated using a fine-resolution (30 m) land-cover map of China (Zhangetal.2019),downloaded from http://data.casearth.cn.

(4) Winter wheat distribution data.The method proposed in the study by Zhangetal.(2018b) was used to extract winter wheat planted areas during 2000–2019.The extracted winter wheat planted areas are validated using statistical cultivated areas at city and county levels,as shown in Appendix E.

(5) RS phenology data.The crop phenology information (returning green up and maturity dates) was derived using MOD13A2 and MYD13A2 NDVI products and the TIMESAT v.3.3 package for MATLAB.The key parameters required for applying TIMESAT are shown as follows: fitting method,double-logistic; start of season method,seasonal amplitude; season start/end value,0.1/0.4.

The observed values at agro-meteorological sites were used to validate the estimated returning green and maturity dates.The validation result is shown in Appendix F.

Climate dataClimate data,such as daily air temperature (Ta),precipitation (Prec),surface net radiation (Rn),global solar radiation (Rsg),vapor pressure deficit (VPD),and wind speed (WS),were used to drive PRYM–Wheat.The Ta,Rn,and WS data were from the ERA-interim reanalysis dataset,available at https://apps.ecmwf.int/datasets/data/interim-full-daily/levtype=sfc/.The VPD was calculated using the dewpoint temperature in the ERA-Interim reanalysis dataset.Variables from ERA were in 12-h temporal and 0.125-arc-degree spatial resolutions.The average of two 12-h values in one day was regarded as the daily value of one variable.The precipitation data were interpolated from meteorological site data,which was available at http://www.nmic.cn/.Daily Rsg was calculated from the interpolated siteobserved daily temporal range (TR) and daily sun hours (HrS).Daily Rsg was not directly retrieved from the site scale observations because the sites observing Rsg were too sparse to be interpolated.Daily Rsg was calculated referring to an equation proposed by Chenetal.(2004):

wherea,b,c,andddenote empirical coefficients.The averaged values for the four coefficients over China were used (Chenetal.2004),i.e.,a=0.04,b=0.48,c=0.83,andd=0.11.

Soil property dataThe dataset ‘Global Gridded Surface of Selected Soil Characteristics’ of the International Geosphere-Biosphere Program–Data and Information System (IGBP–DIS),available at https://daac.ornl.gov/,was used to retrieve soil property data,such as field capacity,wilting point,and saturated water content.In this study,the soil property data have a spatial resolution of 5×5 arcminutes and are measured within the depth interval of 0–100 cm.

Agro-meteorological station dataIn this study,we used 92 agro-meteorological station data (Fig.1),including their latitudes,longitudes,crop species,sowing dates,and maturity dates.These data were obtained from the AGME_AB2_CHN_TEN dataset,which was available at http://data.cma.cn/.The sowing and maturity dates of winter wheat were converted to DOY and then used to calculatekiin eq.(6) and validate the estimated returning green and maturity dates of winter wheat.

Observed winter wheat yieldsWe considered projects previously conducted by our team to collect 24 site-years observed winter wheat yields at 12 sites (Fig.1 and Table 3) in the study area.The observed yield data were used to validate the winter wheat yields estimated by the developed PRYM–Wheat at the site scale in this study.

Table 3 Sites information for observed winter wheat yields

Statistics of winter wheat planted area and yieldIn this study,we used the official statistics of winter wheat cultivated area and yield data during 2000–2019 at city and county levels (Table 4).The yearbooks of each province (https://www.yearbookchina.com/index.aspx) were used to collect the official data that were then used to validate the estimated winter wheat areas (results are shown in Appendix E) and yields.The estimated winter wheat-cultivated areas and yields in these cities or counties were averagely aggregated and compared with the corresponding statistical winter wheat yields.

Table 4 Collected statistical winter wheat cultivated areas and yield data in this study

3.Results

3.1.Determination of IPPs

As shown in Fig.4,theRand RMSE of statistical yieldsvs.estimated yields during reference years (2010–2015) under 21 tests (T0–T20) over the whole study area were calculated.Compared to the T0 case,Rand RMSE from other modeling cases exhibit an increase and decrease of 0.05 (in T6)–0.15 (in T18) and 0.33 t ha–1(in T6)–0.91 t ha–1(in T18),respectively.Rincreased from 0.41 in T0 to 0.56 in T18,while RMSE decreased from 2.20 t ha–1in T0 to 1.30 t ha–1in T18.Both the maximumRand minimum RMSE appeared in T18.Fig.4 also shows that theR(RMSE) values of T1–T6 with once irrigation (I-1),T7–T16 with twice irrigations (I-2),and T17–T20 with a three-time irrigation strategy (I-3) are higher (lower) than that of T0 (I-0).This finding demonstrated that the proposed irrigation strategy could accurately estimate regional winter wheat yield.Fig.4 displays that T17–T20 simulations with three-time irrigation obtain higher (lower)R(RMSE) values than T1–T6 with one-time irrigation and T7–T16 with two-time irrigations.This observation indicated that one- and two-time irrigations are not the best strategy for representing the pattern of winter wheat irrigations and for yield estimation over the NCP.In addition,theRvalues of T7,T8,and T9 with twice irrigations are similar to that of T18 (R=0.56) with threetime irrigations,while the RMSE of T18 is the lowest at 1.30 t ha–1.Thus,the IPPs used in T18 (three-time irrigations at returning green,booting,and filling stages,respectively) were considered as the IPPs for winter wheat yield simulation over the NCP.

Fig.4 The correlation of Yield_Est vs.Yield_Sta data pairs (R; A) and root mean standard error (RMSE; B) of estimated yields under T0–T20 vs.statistical yields for the whole study area.I-0,I-1,I-2,and I-3 represent test groups with zero,once,twice,and three times irrigation,respectively.The total number of Yield_Est and Yield_Sta data pairs is 1 211.Yield_Est,the yields estimated using PRYM–Wheat in this study; Yield_Sta,the statistical yields collected from the yearbooks.

3.2.Validating PRYM–Wheat with optimized IPPs over the NCP

The PRYM–Wheat with optimized IPPs was used to simulate winter wheat yields over the NCP across all validation years (2000–2009 and 2016–2019).In brief,only a one-year simulation is shown (Fig.5).The highyielding areas are mainly distributed in the southeast of Henan Province and the western part of Shandong Province.This result is consistent with previous findings (Wangetal.2011; Yangetal.2017).

Fig.5 Spatial distribution of estimated winter wheat yields in 2019.

The estimated yields of winter wheat were validated using site-observed yields and statistical yields at city and county levels (Fig.6).As the winter wheat planted areas around the site may be smaller than 1 km×1 km,the pixel where the site is located was recognized as a non-winter wheat pixel in the extracted winter wheat areas with a spatial resolution of 1 km.Accordingly,we averaged the pixels within a circle having a radius of 5 km around the site to represent the estimated site-level yield.TheR(RMSE and mean absolute error,MAE) of observedvs.estimated winter wheat yields are 0.80 (0.62 and 0.47 t ha–1),as shown in Fig.6-A.We conducted the validation for each province as the statistical data in different provinces are on different spatial scales.As shown in Fig.6-B–D,theR(RMSE and MAE) of statisticalvs.estimated winter wheat yields of Henan,Hebei,and Shandong provinces are 0.73 (0.97 and 0.79 t ha–1) on a county level,0.61 (0.91 and 0.79 t ha–1) on a county level,and 0.55 (0.75 and 0.62 t ha–1) on a city level,respectively.This observation demonstrates that the new approach to approximating irrigations used in PRYM–Wheat performs relatively satisfactorily for estimating multi-year winter wheat yields.

Fig.6 Scatterplots of observed vs.estimated yields on a site level (A),a county level for Hebei Province (B),a county level for Henan Province (C),and a city level for Shandong Province (D),respectively.Yield_Obs,observed yields; Yield_Sta,the statistical yields collected from the yearbooks; Yield_Est,the yields estimated using PRYM–Wheat in this study,and they were aggregated on a city level for Shandong Province and a county level for Hebei Province and Henan Province.R,the correlation of Yield_Est vs.Yield_Sta data pairs; RMSE,root mean standard error; MAE,mean absolute error.

We showed the model performance statistics (R,RMSE,MAE,slope (a),intercept (b),and sample size (n)) from regression (y=ax+b) of estimated yields (asy) on statistical yields (asx) for each year and each province in Appendix G.Generally,the PRYM–Wheat with IPPs could get a higher (lower)R(RMSE and MAE) value in Henan than in Hebei and Shandong.

4.Discussion

4.1.Further understanding of IPPs

In this study,the proposed IPPs refer to a set of parameters,including irrigation frequency and timing,that define the irrigation pattern implemented during winter wheat cultivation.Notably,the IPPs are pattern parameters,not the actual irrigation frequency and timing.These parameters could be easily coupled into a processbased and RS-driven model to present approximate irrigation of winter wheat over the study area.

As shown in Fig.4,the highestR(0.56) and the lowest RMSE (1.30 t ha–1) for the whole study area (Hebei,Henan,and Shandong provinces) appear in T18.In particular,theR(RMSE and MAE) under T18 for Hebei Province is 0.66 (1.50 and 1.31 t ha–1) on a county level,that for Henan Province is 0.63 (0.98 and 0.79 t ha–1) on a county level,and that for Shandong Province is 0.46 (0.99 and 0.76 t ha–1) on a city level.

We determined the IPPs for each province using the data of reference years (2010–2015) to further understand the role of IPPs in PRYM–Wheat,and the obtained results are shown in Table 5.Based on our findings,the performances of province-specific IPPs were slightly different from that of the IPPs for the whole NCP with respect to the estimation of winter wheat yield.The maximumRappeared in T6 and T15–T20,containing T18,for Hebei Province.For Henan Province,the maximumRappeared in T18,as interpreted in Fig.4 for the whole study area.The maximumRappeared in T17 for Shandong Province,different from that interpreted in Fig.4.Our findings showed insignificant differences in the maximumRvalues for Shandong Province under T17 (0.47) and T18 (0.46).

Table 5 Modeling case with maximum R,minimum RMSE (t ha–1),and minimum MAE (t ha–1) among 21 tests for each province1)

At the city level,we attained the modeling cases with maximumR,minimum RMSE,and minimum MAE among 21 tests for each city using the same data (2010–2015 statistical and estimated yields) (Table 6).As shown in Fig.7,the frequency numbers of regions with maximumR,minimum RMSE,and minimum MAE on each test were counted.Our findings showed that modeling cases with maximumR,minimum RMSE,and minimum MAE for these city-scale regions show more variations,particularly in the case of modeling cases with maximumR.In the whole study area,10 of 45 cities obtained maximumRamong 21 tests on T0,while 31/28 of 45 cities attained the minimum RMSE/MAE among all tests on T18.These imply that PRYM–Wheat using IPPs set in T18 could get lower RMSE and MAE than other tests when used for winter wheat yield simulation over the NCP.Fifteen regions among the total 45 city-scale regions over the whole study area exhibited anRdifference (the maximum absoluteR-value among 21 tests–Rattained on T18) larger than 0.10 (Table 7).This finding suggests that using IPPs for the whole study could result in a precision loss but with no significant effect on the estimation accuracy.This also indicates that the estimation accuracy can be improved by determining IPPs for each subregion.In the future,we intend to optimize IPPs in a smaller region to obtain more accurate simulations.

Table 6 Modeling case with maximum R,minimum RMSE,and minimum MAE among 21 tests for each city1)

Table 7 R difference for each city-scale region1)

Fig.7 Frequency number of regions with maximum R,minimum RMSE,and minimum MAE among 21 tests.R,the correlation of estimated yields vs.statistical yields; RMSE,root mean standard error; MAE,mean absolute error.

4.2.Uncertainty analysis

Uncertainties in yield estimationA few uncertainties were observed in RS input data and biological parameters.The current LAI products are not directly used as model inputs for crop carbon,ET,and yield estimation because they generally are underestimated for crops (Yangetal.2006,2015; Fuetal.2016; Fangetal.2019; Zhangetal.2020).Most of the previous studies focused on regional crop yield simulation generally used LAI yielded by the authors rather than using the current LAI products as one of the model inputs (Lietal.2017; Fangetal.2018).In the present study,we computed LAI using an empirical equation established by field-observed LAI and RS EVI to drive the developed PRYM–Wheat model.The simulated LAI in this study was also underestimated from the booting to flowering stages.However,the use of the estimated LAI enabled it to alleviate the underestimation of MODIS LAI to some extent (Fig.8).A previous study showed that an acceptable error of less than 1 t ha–1in yield estimation was obtained using LAI estimated based on different empirical equations (Zhang 2019).Our findings showed that relatively high RMSE values are yielded,although theRvalues are comparable.Therefore,providing crops LAI with high accuracy is indicated for yield estimation on a regional scale.

Fig.8 Estimated leaf area index (LAI) in this study and observed LAI at a measured site in Raoyang County,Hengshui City,and Hebei Province,China (115.69°E,38.18°N) in 2009.MODIS,moderate resolution imaging spectroradiometer; DOY,the Day of Year.

Other biological parameters are a key photosynthesis parameter (Vm,25),the maximum stomatal conductance to CO2(gsmax),and the restrictive function of nitrogen (N) supply (fN).In this study,Vm,25,gsmax,andfNwere fixed at 180 μmol m–2s–1,0.005 m s–1,and 0.7,respectively (Appendix C).The previously published studies were referred to obtain the values ofVm,25andgsmax(Liuetal.1997; Zhangetal.2018a; Guoetal.2015; Huangetal.2019).The study by Zhang (2019) was referred to obtain the value offNusing a trial-and-error method.These values were set as constant for winter wheat in the whole study area and were effective in estimating the winter wheat yield over the NCP (Zhang 2019).A higher accuracy can be obtained by considering the spatial variations of these biological parameters in yield estimation.However,there is still uncertainty regarding the depiction of spatial variations of such biological parameters.

Uncertainties in model validationUncertainties can also be introduced by validation data,except for the uncertainties in the model input data.The sources of the estimated yields and the statistical yields used for validation showed a scale mismatch.The former is estimated on a pixel scale,while the latter is obtained using a sample survey method.Wangetal.(2011) reported the challenges in validating the modeled yield using ground measurements because of mixed pixels.They also discussed the mismatch between the estimated pixel-level yield and the ground-measured yield.With respect to regional winter wheat yields estimation,some studies aggregated the estimated pixel-level yields at the city or county level and then validated the yields using official statistical yield data at the corresponding level (Huangetal.2015a; Lietal.2017; Wangetal.2019).Such an approach leads to a lowerR-value in validation at the regional scale than at the site or local scale.Statistical yields are commonly used to validate region crop yield estimates due to their easy availability and continuity.However,the source of statistical yields differs from that of modeled yields.In this study,we excluded some regions (mainly cities,such as Zhangjiakou,Chengde,Tangshan,Qinhuangdao,Beijing,and Tianjin) where the winter wheat planted area are small or scattered distributed to mitigate the effect of mixed pixels on validation and reduce the uncertainty introduced by the mismatch between statistical and estimated yields.This approach is more effective compared to the direct use of all cities or counties to conduct the validation.

4.3.Possible reasons for the underestimation of yield in Hebei

As described in Section 3.2,the winter wheat yields were underestimated by PRYM–Wheat in the middle and south of Hebei Province.The parameterizations of the key photosynthesis parameter (Vm,25) and the harvest index (HI) are unlikely to account for this issue,as their values were fixed on relatively high levels (180 μmol m–2s–1and 0.48,shown in Appendices C and D,respectively) in the simulation.Although the empirically estimated LAI showed more improvements compared to MODIS LAI,it still underestimated the in-situ values (Fig.8).The underestimation in LAI could be the potential reason for the underestimation of winter wheat yields in Hebei Province because LAI is used to upscale the leaflevel evapotranspiration and photosynthesis to a canopy scale in PRYM–Wheat.Hebei Province comprises cities,counties,villages,and roads scattered in the middle and south of the region,enhancing the negative effect of empirically estimated LAI.Thus,we calculated the number of impervious pixels per 50 km2over the NCP and found that the impervious surface pixels were concentrated in Beijing Municipality,Tianjin Municipality,and the mid-south part of Hebei Province,as shown in Fig.9.Impervious was also concentrated in the midsouth and mid-east parts of Shandong Province; however,these regions were not the main production areas of winter wheat (Fig.1).The above analysis implies that the dense impervious pixels negatively affect the accuracy of simulating winter wheat yields to some extent.Large numbers of 1-km resolution mixed pixels were formed due to the concentrated distribution of impervious pixels (30 m×30 m spatial resolution) in Beijing Municipality,Tianjin Municipality,and the mid-south part of Hebei Province.These mixed pixels partly account for the underestimation of winter wheat yield in Hebei Province.In addition,the concentrated impervious surface influences the reflectance of the neighboring crop pixels,leading to the deviation of the retrieved spectral vegetation index and leaf area index.

Fig.9 Distribution of impervious pixels (A) and the sum of impervious pixels’ number in per 50 km2 (B) over the North China Plain.The impervious data are with 30-m spatial resolution,downloaded from http://data.casearth.cn.The impervious pixels’ number per 50 km2 is calculated using Focal Statistics Tool in ArcGIS Software.

5.Conclusion

We proposed a new approach to approximating irrigations of winter wheat over the NCP,considering the existence of a relatively stable pattern of irrigations for winter wheat over the NCP.The new approach was incorporated into the PRYM–Wheat model developed in this study to improve the regional estimates of winter wheat yield over the NCP.The main conclusions are as follows:

(1) The optimal IPPs for the NCP’s winter wheat is the three-time irrigation applied at the returning green,booting,and filling stages.

(2) IPPs remarkably contributed to more reasonable estimates of winter wheat yield over the NCP.PRYM–Wheat with the optimal IPPs could accurately simulate winter wheat yield over the NCP,with an increase and decrease inRand RMSE of 0.15 (about 37%) and 0.90 t ha–1(about 41%),respectively.

(3) PRYM–Wheat incorporating the optimal IPPs can effectively estimate winter wheat yields over the NCP and multiple years with reasonable performances.TheR(RMSE) values were 0.80 (0.62 t ha–1) on a site level,0.61 (0.91 t ha–1) for Hebei Province on a county level,0.73 (0.97 t ha–1) for Henan Province on a county level,and 0.55 (0.75 t ha–1) for Shandong Province on a city level.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (42101382 and 41901342),the Shandong Provincial Natural Science Foundation (ZR2020QD016),the National Key Research and Development Program of China (2016YFD0300101).

Declaration of competing interest

The authors declare that they have no conflict of interest.

Appendicesassociated with this paper are available on https://doi.org/10.1016/j.jia.2023.02.036