Phenology-Based Remote Sensing Assessment of Crop Water Productivity

Abstract

The assessment of crop water productivity (CWP) is of practical significance for improving regional agricultural water use efficiency and water conservation levels. The remote sensing method is a common method for estimating large scale CWP, and the assessment errors in CWP by remote sensing originate mainly from remote sensing inversion errors in crop yield and evapotranspiration (ET). The phenological period is the important factor in crop ET and yield estimation. The crop coefficient (Kc) and harvest index (HI), which are closely related to different phenological periods, are considered during the processes of crop ET and yield estimation. The crop phenological period is detected from enhanced vegetation index (EVI) curves using Moderate Resolution Imaging Spectroradiometer (MODIS) data and Sentinel-2 data. The crop ET is estimated using the surface–energy balance algorithm for land (SEBAL) model and Penman‒Monteith (P-M) equation, and the crop yield is estimated using the dry matter mass–harvest index method. The CWP is calculated as the ratio of the crop yield to ET during the growing season. The results show that the daily ET and crop yield estimated from remote sensing images are consistent with the measured values. It is found from the variation in daily ET that the peaks appear at the heading period of wheat and maize, which are in good agreement with the rainfall and growth characteristics of the crop. The relationship between crop yield and ET shows a negative parabolic correlation, and that between CWP and crop yield shows a linear correlation. The average CWPs of wheat and maize are 1.60 kg/m³ and 1.39 kg/m³, respectively. The results indicate that the phenology-based remote sensing inversion method has a good effect on the assessment of CWP in Lixin County.

1. Introduction

Scarcity of water resources is one of the main constraints to achieving more food production [1]. CWP is an important indicator for assessing agricultural water resources, reflecting whether agricultural production levels and irrigation methods are scientific and reasonable. It also provides a reliable basis for water resource management [2]. CWP is the ratio of crop yield to its evapotranspiration during the growing season [3,4]. There are different water requirements and water contributions to yield at different stages of crop growth and development. Therefore, the crop phenological stage is an important factor in the remote sensing retrieval process of ET and yield, which is closely related to the estimation of CWP.

The remote sensing method of retrieving crop ET and yield from satellite images is a common method for estimating large-scale CWP [5,6]. Cai et al. [7] estimated the actual ET by the Surface Energy Balance System (SSEB) remote sensing inversion model and crop yield from NDVI products of MODIS and then evaluated the CWP of rice in the Indo-Gangetic River basin. Zwart et al. [8] developed a water productivity model (WATPRO) on the basis of remote-sensing-derived input datasets, which was a combination of Monteith’s theoretical framework for crop yield and an energy balance model to assess ET. Nana Yan et al. [9] estimated crop ET using ETWatch and yield using the Carnegie–Ames–Stanford Approach (CASA) model, with the purpose of analyzing crop water productivity in the Haihe basin. Bastiaanssen and Steduto [1] assessed the global water productivity of wheat, rice, and maize by using the SEBAL model for ET and the dry matter mass–harvest index method. In terms of crop ET estimation, the SEBAL model as a typical single-source remote sensing model simplifies input parameters compared with dual-source models, and is more flexible for processing remote sensing images [10]. In terms of yield estimation, the dry biomass–harvest index method is a widely used crop yield estimation model by remote sensing [11,12,13]. Although remote sensing inversion models of ET and yield are widely applied in CWP assessments on a large scale, the remote sensing estimation accuracy still needs to be improved.

The assessment errors in CWP by remote sensing originate mainly from remote sensing inversion errors in crop yield and ET, which may propagate to errors in CWP as high as 108.2% [14]. Yuan et al. [15] found that the EC-LUE model has good performance in simulating gross primary production (GPP). However, HI is a fixed value without consideration of different phenology characteristics, which leads to errors in crop yield simulation. Kalma et al. [16] considered that the temporal scale expansion has an impact on evapotranspiration estimation. Errors exist especially in the temporal scale expansion from daily to the whole growth stage, and the main factor leading to this situation may be that different crops in different places have different phenological characteristics. Therefore, the accuracy of CWP may be improved by considering phenological characteristics in crop ET and yield estimation.

Due to crop phenology variability, applying fixed dates of the crop calendar to crop ET and yield estimation may affect the assessment accuracy of CWP. Kc and HI are parameters that are closely related to the phenological period. Kc is an important parameter in the temporal scale expansion of ET, and HI is an important parameter for crop yield estimation. However, Kc and HI are often defined as constants, which may lead to the low accuracy of CWP estimation [1,17]. Zhao Nana et al. [18] found that Kc varies greatly at different phenological stages of the summer maize growing season. M. Moriondo calculated HI by combining NDVI and crop phenological characteristics and applied it to agricultural water use assessment [19]. The accuracy of CWP estimation may be improved by acquiring the variability in Kc and HI on the basis of crop phenology.

Crop phenological characteristics and key time points of the crop calendar can be extracted from satellite images on the basis of remote sensing. The normalized difference vegetation index (NDVI) and EVI can reflect the change in vegetation properties in the phenological period [20]. In comparison with NDVI, EVI, which integrates information from NIR (near-infrared band) and red and blue bands, provides a better optimization for vegetation phenological signals by reducing the influence of atmospheric attenuation and canopy background uncertainty; thus, crop phenological characteristics are more accurately extracted from EVI by remote sensing than NDVI [21,22]. Two extraction methods of crop phenology from vegetation indices are the threshold method and the slope method. The threshold method needs to be manually adjusted according to crop types and area range, while the slope method occurs through the extraction of crop phenology information directly from the slope of the EVI curve [23,24,25].

In this paper, crop phenological characteristics are extracted from EVI curves according to the slope method, and the MODIS data and Sentinel-2 data are fused to improve the spatial resolution. After retrieving crop phenology by remote sensing, the cumulative ET during the growing season was obtained by using the temporal scale expansion of Kc based on crop phenology. The HI for crop yield estimation is calculated according to key phenological time points. Thus, Kc and HI based on crop phenological characteristics are introduced into the remote sensing assessment of CWP, which may improve the assessment accuracy.

Lixin County, the study area, is located in the Huaibei Plain, which is one of the important commercial grain production bases in China. The objectives of this study are to: (1) detect the crop phenological periods; (2) estimate ET, yield, and CWP of wheat and maize; (3) analyze the spatial distribution of ET, yield, and CWP; (4) analyze the results of phenological periods and ET; and (5) analyze the relationships among ET, yield, and CWP. These results will provide a reference for regional farmland water consumption and a decision-making basis for the rational and efficient allocation of regional water resources.

2. Study Area and Data Acquisition

2.1. Study Area

Lixin County (N 32°51′–33°27′, E 115°54′–116°31′) is in the northwest of Anhui Province, which is shown in Figure 1. The total area is approximately 2005 km², with a population of 1.236 million. The terrain slopes gently, with an average elevation of 28 m. Arable land resources are abundant, with 1.44 ha of arable land per capita. Lixin County is located in a semi-humid monsoon climate zone, with an average annual temperature of 14.8 °C and an average annual precipitation of 824 mm. The winter is dry and less rainy, and the summer is hot and rainy. The summer harvesting crop is mainly wheat, and the autumn harvesting crop is mainly maize. The classification maps of wheat and maize in Figure 1 are extracted from Landsat images using the decision tree classification method.

2.2. Data Acquisition

2.2.1. Satellite Data

The time of the satellite data is from 2020 to 2021. The information regarding the data is shown in

Table 1. Information of satellite data.

Datasets	Source	Path, Row/Tile	Spatial Resolution	Time Resolution
MODIS(MOD09A1)	https://ladsweb.modaps.eosdis.nasa.gov/ (accessed on 9 December 2021)	h27v05	500 m and 250 m	8 d
Landsat7 and 8	https://earthexplorer.usgs.gov/ (accessed on 15 December 2021)	122,037	30 m	16 d
Sentinel-2	https://scihub.copernicus.eu/ (accessed on 17 December 2021)	T50SLB and T50SMB	10 m	5 d

and

Table 2. Acquisition times of Landsat7 and 8 and Sentinel-2, represented by day of year (DOY).

Datasets	2020	2021
Landsat7	249,297,313	123,155,251
Landsat8	177,241	019,051,083,099,211,243,275,307
Sentinel-2	315,320,335,355,365	009,014,019,029,049,084,094,099,109,124,129,184,254,264,274,279,304,314

.

MOD09A1 product is the surface reflectance dataset of MODIS. The spatial resolution of the red and NIR bands is 250 m, and others are 500 m.

Sentinel-2 includes 13 bands, the highest spatial resolution of which is 10 m and the lowest spatial resolution of which is 60 m. In this paper, the blue, red, and NIR bands were chosen to calculate the EVI, with a spatial resolution of 10 m.

2.2.2. Meteorological Data

Meteorological data from Bozhou station were obtained from the China Meteorological Data Network (https://data.cma.cn/ (accessed on 15 December 2021)). The location of this station is shown in Figure 1a. The meteorological data in this paper include temperature, wind speed, actual water vapor pressure, saturated water vapor pressure, and so on, from April 2020 to June 2021.

2.2.3. Field Data

The measured ET, as the validation data, were calculated through the Soil Water Atmosphere Plant model (SWAP) combined with meteorological data, soil profile data, and soil water content data. The measured soil moisture data were obtained from the smart soil moisture monitoring instrument (SSMMI), the location of which is shown in Figure 1b.

The measured yields of maize and wheat were measured near the Zhuma Ditch, Chezhe Ditch, and Xihongsi Ditch, and at the same time, the locations of the sampling sites were recorded through the handheld GPS. The locations of the sampling sites are shown in Figure 2. Maize yield was measured in late September 2021. In the field investigation, we selected 50 maize sampling points and 48 wheat sampling points. At every sampling site, 10 maize cobs were randomly selected in an area of 6 m², which were then threshed, dried, and weighed. Wheat yield was measured at the end of May 2021. At every sampling site, some wheat ears were selected in a random area of 1 m² and then were weighed after threshing and drying.

3. Methods

Crop water productivity (CWP) is the ratio of crop yield to ET during the growing season. Crop phenological period is detected from EVI curves based on MODIS and Sentinel-2 data. Crop seasonal ET (ET_season) is calculated by the SEBAL and P-M methods. Crop yield is estimated by the dry biomass–harvest index method. For accurate ET and yield estimation, the dates of crop phenological period are added to the processes of crop ET and yield estimation. The methodology framework is shown in Figure 3.

3.1. Crop Phenological Period Detection Method

Crop phenological periods are detected from the EVI curve by using the slope method. The MODIS and Sentinel-2 data are fused to increase the accuracy. Five phenological periods are obtained, including the sowing period, emergence period, heading period, anthesis, and maturity period.

3.1.1. The EVI Curve Extraction from the Fusion Data of MODIS and Sentinel-2

An EVI similar to NDVI is used to quantify vegetation greenness, which is appropriate for crop phenological period detection [26,27]. Satellite remote sensing image data resources provide spatially and temporally extensive options for estimate EVI values. The revisit period of MODIS is short, and the data quality is stable. However, it does not apply to the small-scale area due to the low spatial resolution [28,29]. The spatial–temporal resolution of Sentinel-2 data is higher than that of MODIS data, but the data quality is unstable [30].

In this study, MODIS and Sentinel-2 data are fused to extract EVI time series with a spatial resolution of 10 m and a time resolution of 8 d. The continuity correction method is used to correct every MODIS image using all the available Sentinel-2 images of the year. The data fusion formula is as follows:

$$V_H\left(t_i\right)=V_M\left(t_i\right)+K·\frac{{{\displaystyle \sum }}_{j=1}^n\left[\omega \left(t_i,t_j\right)\left(V_T\left(t_j\right)-V_M\left(t_j\right)\right)\right]}{{{\displaystyle \sum }}_{j=1}^n\omega \left(t_i,t_j\right)}$$

(1)

$$\omega \left(t_i,t_j\right)=\frac{1}{\left|t_i-t_j\right|+1}$$

(2)

where $V_H\left(t_i\right)$ is the value of the image pixel after fusing. $V_M\left(t_i\right)$ is the value of the MODIS image pixel. $V_T\left(t_j\right)$ is the value of the Sentinel-2 image pixel. n is the number of Sentinel-2 images. $t_i$ is the DOY corresponding to the MODIS data; $t_j$ is the DOY corresponding to the Sentinel-2 data; and $\omega \left(t_i,t_j\right)$ is the weight of the high-resolution image at $t_j$. K is an empirical coefficient and accounts for the influence of missing Sentinel-2 images on data fusion.

The fundamental principle of the above data fusion method was proposed by Liu et al. [31]. In this study, the correction coefficient K is introduced in Equation (2) for correction to consider the effect of the uneven time interval of available Sentinel-2 images.

The data fusion effect is shown in Figure 4. Figure 4a presents the MOD09A1-EVI data image on 5 September 2020 as an example, and the image effect is shown in Figure 4. The spatial resolution of the MODIS image before fusion is 500 m, and the features of the image are vague due to the mixed pixels (Figure 4a). After fusing with all the available Sentinel-2 images, the features become clear (Figure 4b).

The EVI time series curves are extracted from the fusion data of MODIS and Sentinel-2. However, noise exists in the EVI curves, and irregular fluctuations in the EVI curve are not conducive to detecting phonological characteristics [32,33]. To obtain the smooth EVI time series curves, the EVI data are interpolated using the cubic spline method and reconstructed using the S-G filtering method.

3.1.2. Crop Phenological Period Detection

On the basis of crop classification maps, phenological periods of wheat and maize are detected in the crop-planting regions. The crop phenological periods are classified into sowing (D₁), emergence (D₂), heading (D₃), anthesis (D₄), and maturity (D₅) periods. Figure 5 shows an example of the EVI time series curve and the metrics for detecting phenological periods. The sowing, heading, and maturity periods are feature points in EVI time series curves, so they are detected using the slope method from the curves [34,35]. The emergence and anthesis periods were detected according to other studies [36,37]. The approaches for detecting phenological periods can be expressed as:

Sowing period ($D_1$): The cultivated land is generally bare after harvesting the previous crop, so the EVI curve witnesses a downward trend. After sowing, the EVI curve witnesses an upward trend. Therefore, the lowest point of the EVI curve is identified as the sowing period. Specifically, the slope of the point is zero and changes from negative to positive, i.e., $EV{I}_{D_1}=Min\left(EVI\right)$ and $EV{I}_{D_1}^{\prime }=0$.

Emergence period ($D_2$): In this paper, one week after the sowing period is identified as the emergence period [36], i.e., $D_2=D_1+7$.

Heading period ($D_3$): The growth and development of vegetative organs tend to be complete. At the same time, the EVI value of the crop peaks at this time, with the slope of the curve being 0 and the slope changing from positive to negative. Therefore, the period corresponding to the maximum EVI is identified as the heading period, i.e., $EV{I}_{D_3}=Max\left(EVI\right)$ and $EV{I}_{D_3}^{\prime }=0$.

Anthesis ($D_4$): Wheat anthesis starts after the heading period of one week, i.e., $D_4=D_3+7$; maize anthesis starts after the heading period 3 days when the bracts are exposed from the filaments of the female ears [37], i.e., $D_4=D_3+3$.

Maturity period ($D_5$): When crops reach maturity, the leaves begin to gradually senesce and turn yellow because of the decrease in chlorophyll content. At the same time, the EVI decreases at the fastest rate, so the point when the first-order derivative reaches the minimum is judged as the maturity period, i.e., $EV{I}_{D_5}=Min\left(EVI\prime \right)$.

Figure 5 shows an example of the EVI time series curve and the metrics for detecting phenological periods. The X-axis indicates the day of year (DOY). The Y-axis indicates the value of EVI. The length of the crop growing season is the duration between the sowing ($D_1$) and maturity ($D_5$) periods.

3.2. Crop ET

The discrete ET_24-hour images for the study area are calculated by the SEBAL model. ETc at the weather station point is calculated by combining the P-M model and crop phenological periods. Daily ET (ET_daily) images for the study area are estimated according to point ETc and ETRS images. Crop seasonal ET_season is equal to the sum of ET_daily images from the sowing period to the maturity period.

3.2.1. Twenty-four-Hour ET (ET_24-hour)

In the SEBAL model, ET is computed from satellite images and weather data using the surface energy balance. Since the satellite image provides information for the overpass time only, SEBAL computes an instantaneous ET (ET_inst) flux for the image time. The ET flux is calculated for each pixel of the image as a “residual” of the surface energy budget equation [38]:

$$\lambda E{T}_{inst}=R_n-G-H$$

(3)

where Rn is the net radiation flux at the surface, G is the soil heat flux, H is the sensible heat flux, and λET_inst is the latent heat flux $\left(\mathrm{W}·{\mathrm{m}}^{-2}\right)$.

Net radiation flux (Rn) is estimated using albedo, transmittance, and long wave emission. Albedo is calculated by integrating surface reflectivity from bands, and weighting coefficients are applied to each band when estimating albedo. Soil heat flux (G) is predicted using vegetation indices and net radiation. Sensible heat flux (H) is calculated from several factors: surface temperature, wind speed, surface roughness, and surface-to-air temperature differences [38].

The instantaneous evaporation rate ($\Lambda )$ can be expressed as:

$$\Lambda =\frac{\lambda E{T}_{inst}}{R_n-G}$$

(4)

G for 24 h periods is assumed to be nearly zero due to the canceling effect of positive G during daylight and negative G during nighttime. SEBAL computes the ET_24-hour by assuming that the instantaneous evaporation rate (Λ) is the same as the 24 h average. $E{T}_{24\text{-}\mathrm{hour}}$ can be expressed as.

$$E{T}_{24\text{-}\mathrm{hour}}=\frac{86400\Lambda \times \left(R_{n24}-G_{24}\right)}{\lambda }=\frac{86400\Lambda \times R_{n24}}{\lambda }$$

(5)

where λ is the latent heat of vaporization ($\mathrm{MJ}·{\mathrm{kg}}^{-1}$); $R_{n24}$ is the daily net radiation; and $G_{24}$ is the daily soil heat flux.

3.2.2. Daily ET (ET_daily)

For the days without accessible remote sensing images, ET_daily images of the study area can be estimated by applying the conversion coefficient (P) to the field crop ET (ETc). ETc can be obtained using the P-M equation and Kc. P is the ratio of ET_24-hour to ETc [39,40].

• ET_C at the weather station point based on the Penman–Monteith equation.

$ETc$ is calculated by using ET₀ and the Kc. The daily crop reference ET (ET₀) is calculated using the P-M equation recommended by FAO. The Penman–Monteith equation is as follows [41]:

$$E{T}_0=\frac{0.408\Delta \left(R_n-G\right)+r\frac{900}{T+273}{U}_2\left(e_s-e_a\right)}{\Delta +r\left(1+0.34u_2\right)}$$

(6)

$$ETc=K_C\times E{T}_0$$

(7)

where $E{T}_0$ is the reference crop ET ($\mathrm{mm}·{\mathrm{d}}^{-1}$); $ETc$ is the well-watered crop ET ($\mathrm{mm}·{\mathrm{d}}^{-1}$); $R_n$ is the net radiation flux at the surface $\left(\mathrm{MJ}·{\mathrm{m}}^{-2}·{\mathrm{d}}^{-1}\right)$; $G$ is the soil heat flux $\left(\mathrm{MJ}·{\mathrm{m}}^{-2}·{\mathrm{d}}^{-1}\right)$; $T$ is the average daily temperature at an altitude of 2 m (°C); $U_2$ is the wind speed at a height of 2 m $\left(\mathrm{m}/\mathrm{s}\right)$; $e_a$ is the actual water vapor pressure $\left(\mathrm{kPa}\right)$; $e_s$ is the saturated water vapor pressure $\left(\mathrm{kPa}\right)$; $\Delta$ is the slope of the curve of saturated water vapor pressure versus temperature $\left(\mathrm{kPa}{°\mathrm{C}}^{-1}\right)$; $r$ is the hygrometer constant 0.066 $\left(\mathrm{kPa}{°\mathrm{C}}^{-1}\right)$; and $K_C$ is the crop coefficient.

Generally, the $K_C$ in Equation (8) is obtained from field trials. The Kc for this study is taken from Hua et al. and Wang et al. [42,43]. The values for Kc are summarized in

Table 3. Values of Kc for wheat and maize.

Phenological Period	Sowing–Emergence	Emergence–Heading	Heading–Maturity	Growing Season
Kc (wheat)	0.67	0.75	0.92	0.75
Growing Season	0~18%	18~54%	54~86%	86~100%
Kc (maize)	0.1	1.1	1.2	0.6

.

2.

ET_daily images based on ET_24-hour images and ETc at the weather station

ET_daily images are calculated by interpolating between two consecutive remote sensing images. P is calculated as the ratio of ET_24-hour from Equation (5) to ETc from Equation (7). The approach assumes that the ET_24-hour for the study area changes in proportion to the change in ETc.

$$P=\frac{E{T}_{24\text{-}\mathrm{hour}}}{E{T}_C}$$

(8)

The ET_daily is calculated as follows:

$$E{T}_{daily}=E{T}_C\times P$$

(9)

The P between remote sensing image dates is calculated by linearly interpolating the P on the satellite image date. The units for ET_daily will be expressed in mm, whereas ETc is expressed in mm/day.

3.2.3. Seasonal ET (ET_season)

The ET_season during the growing season is calculated as follows:

$$E{T}_{season}={{\displaystyle \sum }}_{D_1}^{D_5}E{T}_{daily}$$

(10)

where $D_1$ is the sowing period and $D_5$ is the maturity period. ET_daily is the daily ET during the growing season.

3.3. Crop Yield

The crop yield is estimated using the dry biomass–harvest index method [44]. The following is the yield estimation formula with remote sensing.

$$Y=Biomass\times HI$$

(11)

where Y is the crop yield (kg/ha) and Biomass is the accumulated above-ground dry biomass during the growing season (kg/ha). HI is the harvest index.

The accumulation of biomass is calculated by the following formula:

$$Biomass={{\displaystyle \sum }}_{D_1}^{D_5}\left(\epsilon \times APAR\right)$$

(12)

where $\epsilon$ is the light use efficiency (g/MJ). $APAR$ is the absorbed photosynthetic active radiation. From Section 3.1.2, $D_1$ and $D_5$ are the sowing period and the maturity period, respectively.

The absorbed photosynthetically active radiation (APAR) is calculated using the ratio of photosynthetically active radiation (PAR) to the fraction of photosynthetically active radiation (FPAR). The FPAR and simple ration vegetation index are linearly related. The light use efficiency (ε) is calculated from three factors: soil moisture, the limit of the low temperature, and the influence of the difference between the environmental temperature and optimum temperature [12,45].

The conversion of biomass into crop yield is achieved through HI. Considering the negative environmental effects, $HI$ is predicted using the average NDVI values from the pre-anthesis to post-anthesis stages [19]. The equations are expressed as:

$$HI=H{I}_{max}-H{I}_{range}\times \left(1-\raisebox{1ex}{$NDV{I}_{post}$}\!\left/ \!\raisebox{-1ex}{$NDV{I}_{pre}$}\right.\right)$$

(13)

$$NDV{I}_{pre}=\frac{{{\displaystyle \sum }}_{D_2}^{D_4}NDVI}{D_4-D_2}$$

(14)

$$NDV{I}_{post}=\frac{{{\displaystyle \sum }}_{D_4}^{D_5}NDVI}{D_5-D_4}$$

(15)

where $H{I}_{max}$ is the optimal value of HI (here chosen to be equal to 0.48); $H{I}_{range}$ is the possible range of HI (here chosen to be equal to 0.18);$NDV{I}_{pre}$ is the mean value of NDVI from emergence to anthesis period; $NDV{I}_{post}$ is the mean value of NDVI from anthesis period to maturity period; and $D_2$, $D_4$, and $D_5$ are the emergence, anthesis, and maturity periods, respectively.

3.4. CWP

Crop water productivity (CWP) is the ratio of crop yield to its evapotranspiration during the growing season, and its unit is (kg/m³). The formula can be expressed as:

$$CWP=\frac{Y}{10\times E{T}_{season}}$$

(16)

where $Y$ is the crop yield (kg/ha) and $E{T}_{season}$ is the accumulated crop daily ET during the growing period (mm).

4. Results

4.1. Crop Phenological Period

The sowing, heading, and maturity periods of winter wheat and summer maize are detected from EVI time series curves. Figure 6 shows the spatial distribution of wheat phenological periods in Lixin County. The dates of the sowing period are mainly from 16 to 31 October 2020 (DOY: 290–305), with an average of 296 (22 October 2020). In the northern region, the sowing period is earlier than that in the middle region. The dates of the heading period are mainly from 3 to 15 April 2021 (DOY: 93–105), with an average of 102 (12 April 2021). The spatial distribution of the heading period is similar to that of the sowing period. The dates of the maturity period are mainly from 3 to 21 May 2021 (DOY: 123–131), with an average of 127 (7 May 2021). The wheat in the eastern region matures later than that in the other regions.

The spatial distribution of maize phenological periods is shown in Figure 7. The dates of the sowing period are mainly from 6 to 12 June 2021 (DOY: 157–163), with an average of 160 (9 June). The sowing period in the northern region is later than that in the south. The dates of the heading period are mainly from 30 July to 16 August 2021 (DOY: 212–228), with an average of 223 (11 August). The heading period of maize in the northwest is later than that in other regions. The dates of the maturity period are mainly from September 8 to 22 (DOY: 251–265), with an average of 257 (14 September).

4.2. Crop ET

The crop ET in the study area is estimated using the SEBAL model and P-M equation. Figure 8 shows the spatial distribution of the crop ET. In Figure 8a, the ET_season of wheat mainly ranges from 350 to 580 mm, with an average of 428.51 mm. The ET_season in the central and northern regions is lower than that in the eastern and southern regions, with some areas less than 200 mm. In Figure 8b, the ET_season of maize mainly ranges from 400 to 570 mm, with an average of 422.85 mm. The ET_season in the eastern region is higher than that in other regions. The ET_season in the central and northern regions is low. This is because the farmland in these areas is close to the construction land, which influences the estimated ET values in the mixed pixels.

To evaluate the prediction accuracy, the coefficient of determination (R²), root mean square error (RMSE), and Nash–Sutcliffe efficiency coefficient (NSE) are chosen to analyze the ET. The ET_24-hour estimated by the SEBAL model and the measured ET obtained on the basis of the smart soil moisture monitoring instrument are shown in Figure 9. The RMSE is 0.56 mm/d, the R2 is 0.96, and the NSE is 0.90. Figure 9 shows that the estimated ET_24-hour results are consistent with the observed values.

4.3. Crop Yield

The crop yield is estimated using the dry biomass–harvest index method. Figure 10 shows the spatial distribution of winter wheat and summer maize in Lixin County, 2021. In Figure 10a, the yield of wheat ranges mainly between 5000 and 10,000 kg/ha, with an average of 7169 kg/ha. The yield in the eastern and southern regions is higher than that in the central region. This is because the wheat planted in the central region is close to construction land, which is not conducive to uniform agricultural management. Meanwhile, the estimation of wheat yield is influenced by the mixed pixels. There are many varieties of wheat in Lixin County, including Huai Mai 28 and Huai Mai 22. Varieties of wheat with different resistances to frost and disease lead to the spatial variability in yields. In Figure 10b, the yield of maize ranges mainly between 4000 and 8000 kg/ha, with an average yield of 6081 kg/ha. The yield in the eastern region is higher than that in the northern region.

The crop yield estimated from remote sensing images and the measured yield obtained at the sampling points in the study area are shown in Figure 11. The RMSE of wheat is 758.92 kg/ha, the R² is 0.90, and the NSE is 0.66. The RMSE of maize is 1359.77 kg/ha, the R² is 0.97, and the NSE is 0.58. Figure 11 shows that the estimated yield results are consistent with the measured values.

4.4. CWP

On the basis of Equation (16), crop ET and yield estimated by remote sensing images were used to calculate the CWP. Figure 12 shows the spatial distribution of wheat and maize. The CWP of wheat ranges mainly from 1.0 to 2.0 kg/m³, with an average of 1.6 kg/m³. The regions with low CWP (less than 1.0 kg/m³) are distributed mainly in the center and south. The regions with high CWP (more than 2.0 kg/m³) are distributed mainly in the north. The CWP of maize ranges mainly from 0.9 to 1.5 kg/m³, with an average of 1.39 kg/m³. The regions with CWP less than 1.0 kg/m³ are distributed in the east and north, and there are fewer regions with CWP greater than 2.0 kg/m³.

5. Discussion

5.1. Analysis of Phenological Period Results

Many scholars have studied crop phenological periods on the Huaibei Plain [46,47]. Table 3 shows the crop phenological periods detected by remote sensing, the measured results from Xinmaqiao Agricultural and Water Comprehensive Experimental Station, and the multiyear average of wheat obtained from Yang et al. [46] and maize obtained from Li et al. [47]. Xinmaqiao Agricultural and Water Comprehensive Experimental Station is located in Xinmaqiao Town, Bengbu City, Anhui Province, and this station is affiliated with the Anhui and Huaihe River Institute of Hydraulic Research. The crop phenological periods at Xinmaqiao Station were obtained from 2015–2016. The ranges of crop phenological periods in this paper, presented in

Table 4. Statistics of crop phenological periods.

	Other Studies		Xinmaqiao Station		Predicted Value
	Wheat	Maize	Wheat	Maize	Wheat	Maize
Sowing period	10.13	6.01	10.12–11.3	6.6–7.4	10.10–11.20	5.31–6.31
Heading period	4.21	8.10	4.8–4.24	7.31–8.15	4.2–5.2	7.26–8.16
Maturity period	5.31	9.30	4.25–5.25	8.16–9.21	5.2–5.31	8.31–10.5

, are from minimum to maximum values. The crop phenological periods predicted in this paper are in good agreement with other results.

Figure 13 and Figure 14 show the growing season frequency distribution and spatial distribution of wheat and maize. The duration of the wheat growing season ranged mainly from 187 to 206 days, with an average of 197 days. In the western region, the wheat growing season is shorter than that in the northern and southern regions. The spatial distribution of the wheat growing season is opposite of that of the wheat sowing period and is similar to that of the wheat maturity period. The inverse relation indicates that early-sown wheat matured later than others. The duration of the maize growing season ranged mainly from 89 to 106 days, with an average of 97 days. In the northern region, the maize growing season is shorter than that in the southern region. Similar to the relationship among the growing season, sowing period, and maturity period of wheat, earlier-sown maize matured later.

5.2. Analysis of ET Results Based on Phenology Results

Figure 15 and Figure 16 show the variation in daily ET as well as the comparison between the rainfall and temperature in 2021 and the average of the last five years. Winter wheat and summer maize rotation is the main cropping pattern in Lixin County. The variation in daily ET exhibits an M shape, and the peaks appear at the heading period of wheat and maize. From November 2020 to February 2021, the temperature and daily ET were low. After the period of seedling establishment, the daily ET increases reached the peak value in the heading period of wheat. In the maturity period of wheat and the sowing period of maize, the daily ET was low. In the heading period of maize, the temperature was high, the rainfall was copious, and the daily ET reached the other value peak. In the maturity period of maize, the temperature began to fall, the daily ET also gradually decreased.

In this paper, the ET_season of wheat and maize were similar, which was caused by abnormal precipitation. As Figure 16 shows, the rainfall in the heading period of maize was above the average of the last five years. In the heading period, maize was sensitive to moisture, and its growth was fast [48]. Therefore, the ET_seasn of maize is above the multiyear average.

Figure 17 shows the variation in the ratio of actual to potential ET for different phenological periods. The ratios in the heading period were higher than those in the other periods. In the heading and maturity periods of wheat and maize, rainfall was abundant, and the actual ET was high. Especially in the heading period of maize, rainfall was sufficient to meet the needs of crops.

5.3. Relationships among CWP, Yield, and ET

It is very important in water-saving agriculture to understand the relationship among CWP, crop yield, and ET. By comparing the spatial distribution of yield and CWP, it was found that some areas have high yield but low CWP, and some areas have low yield but high CWP. For agriculture, saving water increases CWP but does not decrease crop yield.

Figure 18 illustrates diagrams between crop yield and ET required to achieve that yield. The linear fitting curves show that the relationship between yield and ET is positive. Compared with the linear fitting method, the polynomial fitting method is better. The polynomial fitting curves show that there is a negative parabolic relationship between the crop yield and ET. The crop yield increases before reaching the peak and then drops as the ET increases. This means that more irrigation does not necessarily result in a higher yield [6]. Therefore, high crop yield, supported by high water supply, will not lead to high CWP.

In Figure 19, there is a close linear relationship between CWP and crop yield. The CWP generally increases as crop yield increases. This conclusion is similar to the study results of Li et al. [6]. Therefore, an important direction of CWP development is to improve the yield in Lixin County.

Improving crop yield is an important research direction to improve CWP in Lixin County. According to a local field investigation, there is little irrigation of wheat in Lixin County, rainfall is concentrated mostly in summer, and adequate irrigation is necessary in the heading period. The variety of wheat is an important factor affecting the final crop yield. According to the investigation, the wheat yield of Woyang 101 is higher than that of other varieties, but many farmers are not familiar with the varieties of wheat. Hence, the crop yield and CWP can be increased by promoting the fine varieties of wheat. Due to the abundant rainfall during the maize growing season, preventing waterlogging is an important factor increasing maize yield. According to statistics, the yield of maize affected by waterlogging is half that of the others.

6. Conclusions

In this paper, the crop phenological periods were detected from EVI curves based on MODIS and Sentinel-2 data. The SEBAL model and P-M equation were applied to estimate the daily ET and accumulated ET during the growing season of winter wheat and summer maize in Lixin County. The crop yields of wheat and maize were estimated using the dry biomass–harvest index method. The CWP of wheat and maize was calculated as the ratio of crop yield to the accumulated ET during the growing season. The main conclusions are as follows:

• The phenological periods detected from remote sensing images are in good agreement with other studies. The average duration of the growing season for wheat and maize were 197 days and 97 days, respectively.
• The daily ET estimated by the SEBAL model and P-M equation is consistent with the measured ET, with R², NSE, and RMSE values of 0.96, 0.90, and 0.56 mm/day, respectively. The average ET_season of wheat and maize from 2020 to 2021 were 428.51 mm and 422.85 mm, respectively.
• The average yields of wheat and maize were 7169 kg/ha and 6081 kg/ha, respectively, and the yield estimated by the dry mass–harvest index was consistent with the measured yield. Based on the analysis of yield and ET_season, there is a negative parabolic relationship between crop yield and ET_season.
• The average CWPs of wheat and maize were 1.6 kg/m³ and 1.39 kg/m³, respectively, and there was a close linear relationship between CWP and crop yield.

In this paper, the detection of crop phenological periods has been effectively applied in the estimation of CWP in Lixin County. However, the conditions of crop growth are complex. Research on Kc and HI needs further study, which will improve the accuracy of the CWP estimated by remote sensing methods.