Applicability of Different Assimilation Algorithms in Crop Growth Model Simulation of Evapotranspiration

Abstract

Remote sensing spatiotemporal fusion technology can provide abundant data source information for assimilating crop growth model data, enhancing crop growth monitoring, and providing theoretical support for crop irrigation management. This study focused on the winter wheat planting area in the southeastern part of the Loess Plateau, a typical semi-arid region, specifically the Linfen Basin. The SEBAL and ESTARFM were used to obtain 8 d, 30 m evapotranspiration (ET) for the growth period of winter wheat. Then, based on the ‘localization’ of the CERES-Wheat model, the fused results were incorporated into the data assimilation process to further determine the optimal assimilation method. The results indicate that (1) ESTARFM ET can accurately capture the spatial details of SEBAL ET (R > 0.9, p < 0.01). (2) ESTARFM ET can accurately capture the spatial details of SEBAL ET (R > 0.9, p < 0.01). The calibrated CERES-Wheat ET characteristic curve effectively reflects the ET variation throughout the winter wheat growth period while being consistent with the trend and magnitude of ESTARFM ET variation. (3) The correlation between Ensemble Kalman filter (EnKF) ET and ESTARFM ET (R² = 0.7119, p < 0.01) was significantly higher than that of Four-Dimensional Variational data assimilation (4DVar) ET (R² = 0.5142, p < 0.01) and particle filter (PF) ET (R² = 0.5596, p < 0.01). The results of the study provide theoretical guidance to improve the yield and water use efficiency of winter wheat in the region, which will help promote sustainable agricultural development.

1. Introduction

The crop growth model converts the processes of photosynthesis, respiration, soil moisture (SM), and energy balance of the plant–soil–atmosphere system into mathematical expressions and assembles them into a program that can simulate crop growth on a computer [1]. CERES-Wheat, widely utilized in agricultural production worldwide, was developed by the Agricultural Research Service of the United States Department of Agriculture specifically for wheat crops [2,3,4]. However, in regions with complex terrain and significant spatial variation, adaptability issues may occur, leading to delays in detecting and responding to changes in crop growth over time [5]. With the advancement of remote sensing technology and the increase in data availability [6], coupled with the growing demand for in-depth research and widespread agricultural applications [7,8], the integration of remote sensing with crop growth models through data assimilation methods has become one of the key focuses and hot topics in agricultural research. Despite advances, current ET simulation methods struggle to capture the variability in semi-arid areas with complex terrain, which affects irrigation planning and water resource allocation. The core concept of assimilating remote sensing into crop growth model shares common ground with land surface process assimilation. It is defined as the method of coupling multi-source remote sensing observations with model simulation values during the simulation process of the crop growth model to enhance simulation accuracy [3,9,10]. EnKF is able to handle high-dimensional, nonlinear, and complex systems without explicit linearization of the model. PF does not rely on the assumption of a Gaussian distribution for the system and therefore works well with highly nonlinear, non-Gaussian systems. Through global optimization (usually a gradient-based optimization method), 4DVar is able to approximate the optimal solution very well and is particularly suitable for time-dependent observations. State variables, such as SM and leaf area index (LAI), are critical for crop models as they capture real-time environmental conditions. Ji et al. employed EnKF to assimilate time-series Sentinel-2 into the CASA-WOFOST coupled model for wheat yield estimation. Field validation results demonstrated that this method exhibited good accuracy and significant potential for field yield estimation. Additionally, the EnKF method was proved more effective in reducing computation time [11]. Zhang et al. utilized cross-wavelet transform to determine the relationship between assimilated variables PF SM (0–20 cm) and PF LAI at multiple time scales, enhancing the extraction of assimilation parameter characteristics and improving the accuracy of yield estimation [12]. Liu et al. used multi-source Sentinel data to assimilate LAI and SM using 4DVar and EnKF and conducted winter wheat yield estimation. The findings indicated that 4DVar-LAI can more accurately identify the phenological period of winter wheat, which aligns with the actual growth and development period [13]. Accurate ET simulation methods can improve water conservation practices for winter wheat farming in semi-arid areas, addressing regional water scarcity challenges.

At present, state variables derived from satellite remote sensing data, such as SM, LAI, and ET, are typically chosen as optimal state variables in crop growth models. LAI can reflect stress and irrigation information to some extent, making it a crucial comprehensive parameter for monitoring crop growth [14]. SM, as a vital state variable, is associated with soil, atmosphere, and vegetation, playing a significant role in hydrological, ecological, and biological processes [15]. Remote sensing LAI is typically acquired through vegetation indices, which exhibit a time-lagged response to drought and are unsuitable for crop yield estimation under drought conditions [16]. The accuracy of remote sensing SM retrieval largely relies on in situ measurements, leading to errors and making it difficult to extrapolate to large regional scales [17]. The ET in farmland refers to the sum of crop transpiration and soil evaporation, constituting the total flux of water vapor transferred from vegetation and the ground to the atmosphere [18]. Currently, remote sensing single-source estimation SEBAL possesses advantages such as easy data acquisition, clear physical concepts, high accuracy, and broad adaptability. It is suitable for calculating surface ET over large-scale regions and long timeframes. Moreover, ET computed based on the SEBAL can fully consider spatial heterogeneity in arid and semi-arid areas, leading to its widespread application both domestically and internationally [19,20,21]. Vazifedoust et al. employed surface energy balance algorithms to integrate MODIS LAI and relative MODIS ET in the Borkhar region of Isfahan, Iran. They implemented an agricultural hydrological model in a distributed manner with a grid size of 250 m. Utilizing a constant-gain Kalman filtering data assimilation algorithm, they corrected internal variables of the distributed model, significantly enhancing the accuracy of winter wheat yield estimation in a small production area of Iran [22]. Zhuo et al. utilized EnKF to assimilate MOD16 ET PET⁻¹ into the WOFOST model for winter wheat yield estimation in the North China Plain from 2008 to 2018. The results showed that assimilating MOD16 ET PET⁻¹ into the WOFOST model reduced model simulation uncertainty and improved yield estimation accuracy [23].

Due to constraints such as satellite revisit cycles [24] and adverse weather conditions (clouds, rain, and snow) [25,26], the acquisition of remote sensing data is often limited [27]. Optical remote sensing data may contain significant gaps in data records, leading to considerable discontinuity in SEBAL’s estimation of ET, thereby preventing continuous monitoring of crop growth during the growing season [28]. The spatiotemporal fusion model combines high temporal resolution data for capturing temporal phase information with high spatial resolution data for acquiring spatial texture details. This makes it a key tool for obtaining continuous and accurate remote sensing observations of crop growth throughout the growing season [29]. Bai et al. utilized ESTARFM to acquire high spatiotemporal resolution spectral parameters and, based on this, performed the inversion of regional high spatiotemporal resolution ET and analyzed the interannual variations in crop ET proportions since the implementation of large-scale irrigation water-saving transformations [30]. He et al. utilized a spatiotemporal fusion model to acquire continuous temporal scale Landsat ET and validated the adaptability of the spatiotemporal fusion model in this region through measured ET [31]. Therefore, integrating product data with complex remote sensing parameters within the same framework has become a development direction for addressing this issue. In this study, for the winter wheat planting area in Linfen Basin, Shanxi Province, ET was used as the state variable for assimilation, and high spatial and temporal resolution ET data (8 d, 30 m) were obtained by ESTARFM spatial and temporal fusion modeling of SEBAL ET and MOD16A2 ET, which provided a rich resource for remote sensing data assimilation for crop growth modeling. Selecting the optimal solution from three assimilation algorithms, EnKF, PF and 4DVar, aims to improve the yield and water use efficiency of winter wheat in the region by improving the science and precision of irrigation management, which helps to promote the sustainable development of agriculture and enhance ecological environment protection.

2. Materials and Methods

2.1. Study Area

The Linfen Basin is situated in the southeastern part of the Loess Plateau in the semiarid zone of the North Temperate Zone. It borders the Taiyue Mountains to the north, with the Taihang Mountains and the Lvliang Mountains surrounding its southeast and northwest sides, respectively (Figure 1). The Fenhe Basin runs through the entire area, converging with the Yellow River in the southwestern part of the basin [31]. The climate type is a semiarid temperate continental monsoon climate, characterized by simultaneous heat and rainfall. The annual average precipitation is approximately 500–600 mm, with rainfall mainly concentrated from June to September. The annual average temperature is around 11.7–13.2 °C, with higher temperatures in summer and lower temperatures in winter. The Linfen Basin is an important grain production base in North China [32,33], suitable for the cultivation of crops such as wheat and corn, with winter wheat accounting for 58.36% of the cultivated land area. Due to its typical semi-arid conditions characteristic of the Loess Plateau, spring and midsummer droughts are frequent in this area. High rates of soil evaporation, coupled with fragmented terrain and water scarcity, are the primary limiting factors for agricultural development in the region [33].

2.2. Data Source

The Sentinel-2 L2A level data used in this study were obtained from the European Space Agency (ESA) (https://dataspace.copernicus.eu/ (accessed on 13 November 2023)). The object-based image analysis (OBIA) within eCognition 9.0 software was used to identify and extract winter wheat areas in the Linfen Basin [34]. The Landsat 8 data were sourced from the Geographic Spatial Data Cloud (http://www.gscloud.cn/ (accessed on 13 November 2023)) and the United States Geological Survey (USGS) data center (https://earthexplorer.usgs.gov/ (accessed on 13 November 2023)). Based on the growth period of winter wheat in the Linfen Basin, Landsat 8 covering the entire growth period of winter wheat with good quality (cloud cover < 10%) were downloaded. The data were as follows: 12 November 2020 (tillering stage); 14 December 2020 (overwintering stage); 15 March 2021 (re-greening stage); 31 March 2021 (jointing stage); 16 April 2021 (tasseling-irrigation stage); 7 May 2021 (tasseling-irrigation stage); 3 June 2021 (maturity stage), totaling 7 images. The MOD16A2 ET were obtained from the NASA (https://ladsweb.modaps.eosdis.nasa.gov/ (accessed on 13 November 2023)). The MOD16A2 ET product is derived from the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor on the Terra satellite. This product has been widely used to assess the spatial and temporal dynamics of ET at different spatial scales. The MOD16 ET algorithm employs a multi-source approach using the Penman–Monteith equation to calculate vegetation transpiration and soil evapotranspiration. Currently, the sixth version of the MOD16A2 ET product has modified its ET algorithm by increasing the spatial resolution (from 1000 m to 500 m) and using an updated land cover product [35,36]. In addition, MOD16A2 ET has been validated in different ecological regions around the globe, such as the Savannah pastoral region in South Africa [37], the shrublands of northwestern Mexico [38], the irrigated areas of southern Italy [39], the humid regions of China [37], and the arid regions of the United States [40]. In summary, it can be seen that MOD16A2 ET data have good accuracy on a global scale and are widely used.

The Digital Elevation Model (DEM) data for the study area utilized the ASTER GDEM 30 m product, which was sourced from the Geographic Spatial Data Cloud Platform (http://www.gscloud.cn/ (accessed on 13 November 2023)). Temperature and wind speed observations were obtained from the China Meteorological Data Network (http://data.cma.cn/ (accessed on 13 November 2023)). Twenty-three meteorological stations within and around the Linfen Basin were selected for the years 2020–2021. The ANUSPLIN model is a tool for interpolating multivariate data using ordinary thin-disk and local thin-disk spline functions, which are widely used for time-series meteorological data [41]. Utilizing the meteorological data interpolation tool ANUSPLIN [41], spatial interpolation of the meteorological data was conducted using latitude and longitude as independent variables and elevation as a covariate. The pixel size for interpolation was set at 30 m.

2.3. ET Calculation

2.3.1. Localization of the CERES-Wheat

The CERES-Wheat model requires four sets of input parameters: meteorological data, soil data, management practices, and genetic parameters. The meteorological data include maximum temperature, daily minimum temperature, daily precipitation, and daily solar radiation, all provided by agricultural meteorological departments. The soil data encompasses sampling depth, soil bulk density, organic carbon, total nitrogen, total phosphorus, total potassium, cation exchange capacity, and pH, which are determined through field and laboratory experiments. Management parameters include irrigation and fertilization amounts, obtained through field surveys. The genetic parameters for wheat varieties include vernalization characteristics, photoperiod characteristics, grain-filling characteristics, grain characteristics, potential grain-filling rate, potential single-stem spike weight at flowering, and leaf emergence interval parameters. These are adjusted to local conditions in the Linfen Basin by inputting the meteorological, soil, and management data. Sample points were chosen based on soil composition, elevation, and proximity to primary wheat-growing regions to ensure they represent typical conditions in the Linfen Basin. Taking Daoxiaoli Village in Houma City as an example, the specific data are as follows: geographic coordinates of the Xiaoli Village site are 111°35′ E, 35°65′ N (the center of the sample area was recorded with a hand-held GPS); date of sowing: 20,280; date of maturity: 21,160; date of irrigation: 20,352 (winter irrigation), 21,079 (spring irrigation); irrigation amount: 60 mm (winter irrigation), 60 mm (spring irrigation); fertilizer application date: 20,280; fertilizer amount: 175 kg/hm² (N), 135 kg/hm² (P), where 20,280, 20,352, 21,079, 21,160 represent 6 October 2020, 17 December 2020, 20 March 2021, 9 June 2021, respectively.

The specific steps are as follows: First, verify whether the simulated maturity date aligns with the observed maturity date. Second, check if the simulated total growth duration of winter wheat matches the actual duration. Third, compare the simulated LAI curve with the actual LAI curve, and determine if the date of maximum LAI occurrence is close. Fourth, assess whether the simulated LAI at a specific time and the final yield are accurate. For this, a canopy analyzer is used to measure the LAI of winter wheat in five 60 cm × 60 cm sample plots, taking the average as the observed LAI for that location. Additionally, on the date of winter wheat harvest, 60 random stalks are selected from the five plots, and the grains are collected, weighed, and recorded. The samples are then dried in an oven at a constant temperature (60–70 °C) until they reach a constant weight, cooled, and weighed again. Finally, the weight of the grain bags is subtracted, and the actual yield for that sample point is calculated.

2.3.2. ESTARFM

The ESTARFM model assumes that remote sensing data from different sources are correlated over a period of time, and the systematic bias between them is stable. Based on at least two pairs of high and low-resolution images before and after the prediction date, as well as a low-resolution image on the prediction date itself, calculate the weights and transformation coefficients between corresponding pixels to simulate high spatiotemporal resolution data on the prediction date [42]. The calculation formula is:

$${ET}_k\left(x_{\omega /2},y_{\omega /2},t_p\right)={ET}_o\left(x_{\omega /2},y_{\omega /2},t_k\right)+\sum _{i=1}^N{W}_i{V}_i\left[M\left(x_i,y_i,t_p\right)-M\left(x_i,y_i,t_k\right)\right]$$

(1)

where ET_k is a high resolution predicted ET at time k; ET₀ is the SEBAL ET at time k; M is the MODISET for the corresponding period; ω is the sliding window’s size; (x_ω/2, y_ω/2) is the sliding window center pixel location; N is the number of similarity pixels in the sliding window; (x_i, y_i) is the location of similarity pixel i; V_i is the conversion coefficient, which is the regression coefficient of the two groups of SEBAL and MODIS ET on the corresponding date; and W_i is the weights of similarity pixels i. The forecast ET at m and n times were then utilized to compute the final forecast ET. The calculation formula is:

$$ET\left(x_{\omega /2},y_{\omega /2},t_p\right)=T_m\times {ET}_m\left(x_{\omega /2},y_{\omega /2},t_p\right)+T_n\times {ET}_n\left(x_{\omega /2},y_{\omega /2},t_p\right)$$

(2)

where ET is the ET of the final prediction date; ET_m is the predicted ET at time m; ET_n is the predicted ET at time n; T_m and T_n are the weight factors at times m and n.

$$T_k={\displaystyle \frac{1/\left|{\sum }_{j=1}^{\omega }{\sum }_{i=1}^{\omega }M\left(x_i,y_i,t_k\right)-{\sum }_{j=1}^{\omega }{\sum }_{i=1}^{\omega }M\left(x_i,y_i,t_p\right)\right|}{\left(1/\left|{\sum }_{j=1}^{\omega }{\sum }_{i=1}^{\omega }M\left(x_i,y_i,t_k\right)-{\sum }_{j=1}^{\omega }{\sum }_{i=1}^{\omega }M\left(x_i,y_i,t_p\right)\right|\right)}}$$

(3)

2.4. Assimilation Methods

2.4.1. 4DVar

The 4DVar algorithm constructs a cost function to describe the discrepancy between simulated state variables and observed values. Within a time assimilation window T, it iteratively adjusts the simulated state variables of the model using all available external observational data until the simulated trajectory of the model closely matches the trajectory of the observational data within the assimilation window, thereby achieving optimization of the simulated state variables of the model [43]. The objective function of 4DVar is:

$$J\left(X_0\right)={\displaystyle \frac{1}{2}}{\left[X_0-X_0^b\right]}^T{B}^{-1}\left[X_0-X_0^b\right]+{\displaystyle \frac{1}{2}}{\sum }_{i=0}^n{\left[X_i-X_i^{obs}\right]}^T{R}_i^{-1}\left[X_i-X_i^{obs}\right]$$

(4)

where X₀ and X₀^b are the initial field and background field of model state variables, respectively, and where X₀^b is the a priori probability estimate within X₀ model state variables; B and R are the covariance matrices of background field error and observation field error, respectively; X_i and X_i^obs are the forecast value and remote sensing retrieval values of the model at time i, respectively.

2.4.2. EnKF

The core of the EnKF is to weight external observations and model-simulated data within the framework of the mechanistic model, computing the optimal estimate of the current state variables and substituting them for the model-simulated values. It then proceeds to the next time step and repeats this process until all external observations have been assimilated into the model, completing the optimization of the model simulated trajectory [44,45]. The objective function of EnKF is:

$$x_k^f=Q\left(x_{k-1}^a\right)$$

(5)

$$x_k^a=x_k^f+P_k{H}^T{\left(H{P}_k{H}^T+R_k\right)}^{-1}\left(D_k-H{x}_k^f\right)$$

(6)

where x_k^f and x_k^a are the forecast and analysis matrix of the ensemble of state variables at time k, respectively; P_k and R_k are the covariance matrices of the forecast and observation matrix, respectively; H is the observation operator; D_k is the observation matrix; Q is the state transition equation.

2.4.3. PF

The basic idea of PF is as follows: during the prediction phase, the simulated state variables (x_k) at time k are perturbed with random noise from any distribution to generate a particle swarm xⁱ_k (i = 1, 2, …, 200). Then, these 200 particles are substituted into the non-linear model and run to time k + 1 to obtain the predicted state of the model xⁱ_k+1. During the update phase, the xⁱ_k+1 at time k + 1, the remote sensing observations (y_k+1) and the initial importance probability density (q(xⁱ_k+1|xⁱ_k, y_k+1)) are used to compute the importance weights (w^ˆi_k+1) for each particle. In this study, it is assumed that q(xⁱ_k+1|xⁱ_k, y_k+1) = p(xⁱ_k+1|xⁱ_k), where p(xⁱ_k+1|xⁱ_k) represents the prior probability distribution function. Then, assuming that the deviation of the predicted particles from the observed values follows a normal distribution, w^ˆi_k+1 is computed using the probability density function of the normal distribution as follows:

$$w_{k+1}^{\widehat{i}}=1/\left(200\sqrt{2\pi O_{k+1}}\right)exp\left(-{\left(x_{k+1}^i-y_{k+1}\right)}^2/\left(2O_{k+1}\right)\right)$$

(7)

where O_k+1 is the standard deviation of the observations (Figure 2).

3. Results

3.1. Accuracy Evaluation of ESTARFM ET

Using Landsat 8 ET and MODIS ET as data sources, the ESTARFM was employed to obtain ET data for the winter wheat planting area in the Linfen Basin with a spatial resolution of 30 m and a temporal resolution of 8 d. Comparing ESTARFM ET, MODIS ET, and the corresponding SEBAL model-retrieved ET at each time point (Figure 3 and Figure S1), it can be observed that from 12 November 2020 to 3 June 2021, there was a clear overall increasing trend in ET across the winter wheat planting area of the Linfen Basin. Spatially, ET gradually decreases from the Fenhe Basin towards both sides, with high ET values mainly distributed on or near the Fenhe Basin sides. The hue of MODIS ET and SEBAL ET images showed significant differences on 16 April 2021 and 7 May 2021, while the remaining time periods exhibited relatively consistent trends. Due to the influence of mixed pixels, MODIS ET exhibited insufficient spatial detail response to land cover types in the Linfen Basin, resulting in a smaller range of ET values compared to the corresponding SEBAL ET values at each time point. MODIS ET was particularly affected by cloud and rainfall during the tasseling-irrigation stage of winter wheat, which exacerbated the simulation errors of ET, leading to significant simulation errors on 16 April 2021 and 7 May 2021.

Further calculations of the correlation coefficient (R), bias, and RMSE between ESTARFM ET and SEBAL ET were carried out and the results are shown in Figure 4. The maximum and minimum values of R between ESTARFM ET and SEBAL ET were 0.98 and 0.281, respectively, with bias ranging from −0.061 to 0.055 and RMSE ranging from 0.015 to 0.145. In addition, by examining the difference histogram, it was observed that ESTARFM ET had higher simulation errors prior to the re-greening stage, with a scatter plot showing R as 0.281 and a difference distribution range of [−0.3, 0.3]. The overall evaluation indicated that the R between ESTARFM ET and SEBAL ET exceeded 0.9 (p < 0.01), with an RMSE of 0.010. At the same time, the range of the difference distribution between ESTARFM ET and Landsat 8 ET is [−0.1, 0.1]. Based on the above quantitative statistical analysis, it was found that ESTARFM ET has a good applicability in the winter wheat production area of the Linfen Basin.

Selected regions within the Fenhe Basin in the inner part of the Linfen Basin, namely the upstream, midstream, and downstream winter wheat planting areas, are depicted in the overview map of the study area. Comparative analysis of MODIS ET, ESTARFM ET, and SEBAL ET (Figure 5) revealed that all ET in the upstream, midstream, and downstream winter wheat planting areas exhibited an increasing trend over time, peaking in May to June, with relatively consistent temporal patterns. Compared to SEBAL ET, MODIS ET had lower spatial resolution and was not able to capture the spatial details of the winter wheat planting areas. Prior to the winter-wheat tasseling stage, MODIS ET was influenced by mixed pixels, resulting in a narrower range of value changes compared to SEBAL ET and a less timely response to detail changes. During periods of high water demand, such as tillering, MODIS ET exhibited inaccuracies in simulating lower thresholds, where the lower threshold was lower than the winter MODIS ET lower threshold. Therefore, MODIS ET data were less suitable for the rugged terrain of the Linfen Basin. Compared to MODIS ET, ESTARFM ET data provided better spatial detail information, effectively expressing spatial differences between smaller features, which mitigated inaccuracies in MODIS ET values and reduced discrepancies between MODIS ET and SEBAL ET. Additionally, the hue of the ESTARFM ET remained highly consistent with SEBAL ET at different times, while its range of values generally aligned with SEBAL ET.

3.2. Localization of the CERES-Wheat ET

The relative error (RE) between calculated CERES-Wheat LAI and CERES-Wheat yield compared to measured LAI and yield at different sample points is presented in

Table 1. CERES-Wheat accuracy validation.

Location	LAI RE (%)	Yield RE (%)
Xiqiaozhuang Village in Hongtong County	2.41	12.57
Tundou Village in Yaodu District	0.77	8.30
Cui Village in Xiangfen County	1.30	11.16
Xinanjie Village in Quwo County	7.18	1.97
Daoxiaoli Village in Houma City	11.17	1.92
Duanjiazhuang Village in Xinjiang County	10.03	4.35
Shangbai Village in Jishan County	9.05	5.71
Cangtou Village in Hejin City	8.01	14.71

. The RE between CERES-Wheat LAI and measured LAI for all sample points was 6.24%. Among them, Tundou Village in Yaodu District had the lowest RE of 0.77%, while Daoxiaoli Village in Houma City had the highest RE of 11.17%. The RE for the remaining sample points, from lowest to highest, were as follows: Cui Village in Xiangfen County (1.30%), Xiqiaozhuang Village in Hongtong County (2.41%), Xinanjie Village in Quwo County (7.18%), Cangtou Village in Hejin City (8.01%), Shangbai Village in Jishan County (9.05%), and Duanjiazhuang Village in Xinjiang County (10.03%). The relative error (RE) between the calibrated CERES-Wheat yield and the measured yield was 7.59%, which was within the acceptable range of error. The RE between CERES wheat yield and measured yield for Xiqiaozhuang Village in Hongtong County, Tundou Village in Yaodu District, Cui Village in Xiangfen County, Xinanjie Village in Quwo County, Daoxiaoli Village in Houma City, Duanjiazhuang Village in Xinjiang County, Shangbai Village in Jishan County, and Cangtou Village in Hejin City were 12.58%, 8.30%, 11.17%, 1.97%, −1.92%, 4.35%, 5.71%, and −14.71%, respectively. Therefore, the calibrated CERES-Wheat model demonstrated high simulation accuracy, making it suitable for winter wheat production areas in the Linfen Basin.

The calibrated CERES-Wheat was used to simulate the variation of ET during the winter-wheat growing period at eight sample points in the Linfen Basin for the years 2020–2021. Here, a-h represent Xiqiaozhuang Village in Hongtong County, Tundou Village in Yaodu District, Cui Village in Xiangfen County, Xinanjie in Quwo County, Daoxiaoli Village in Houma City, Duanjiazhuang Village in Xinjiang County, Shangbai Village in Jishan County, and Cangtou Village in Hejin City, respectively. The results are shown in Figure 6. During the major growth stages of winter wheat, including re-greening, jointing, and tasseling irrigation, there was a significant increasing trend with ET peaking at DOY200 (tasseling-irrigation stage). Due to abrupt changes in meteorological factors at the daily scale, the CERES-Wheat showed strong variability in daily scale ET at the winter wheat sample points, with noticeable noise and large fluctuations at the daily scale, but without significant overall trend changes. Therefore, S-G filtering was applied in this study to smooth and remove the noise. After S-G filtering, ET values at different sample points consistently exhibited significant increases during the major growth stages (re-greening, jointing, tasseling irrigation, and maturity), which agrees well with the actual growth conditions of winter wheat. Due to pre-planting irrigation practices aimed at improving soil moisture conditions in the Linfen Basin, the filtered ET curve of winter wheat exhibited a small wave peak at the early planting stage, specifically at DOY = 25. During the major growth stages of winter wheat, characterized by rapid stem and leaf growth, water use increases significantly and peaks during the tasseling-irrigation stage, represented by a large wave peak at DOY = 225. Dryland sample points (Xiqiaozhuang Village in Hongtong County, Tundou Village in Yaodu District and Cangtou Village in Hejin City) and irrigated sample points showed a consistent increasing trend during the major growth stages, starting from the jointing stage. The highest growth rate of ET occurred between DOY 98 and DOY 106, reaching its maximum value at the end of the tillering stage and the beginning of the tasseling-irrigation stage, followed by a decreasing trend. The average ET values for winter wheat in the irrigated area during the re-greening, jointing, tasseling-irrigation, and maturity stages were 0.901/mm, 2.217/mm, 3.507/mm and 2.058/mm, respectively. In the rainfed area, the ET values at different growth stages of winter wheat are all lower than those in the irrigated area (0.718/mm, 1.944/mm, 3.056/mm, 2.01/mm), and the increase in ET from the jointing stage to the tasseling-irrigation stage was significantly higher in the irrigated sample points than in the dryland sample points. Except for the re-greening stage, the mean ET values for winter wheat in the irrigated area were higher than the median, indicating a greater proportion of high values during the jointing, tasseling-irrigation, and maturity stages. The differences in median and mean ET values at different growth stages of winter wheat in the rainfed area were not significant. Affected by mixed pixels, the 8 d MODIS ET was significantly lower than the winter wheat ET curve throughout the growing period. It also exhibited a delayed response to changes in vegetation growth information, characterized by relatively small fluctuations between adjacent time points. This phenomenon was particularly pronounced during the main stages of winter wheat sowing (late September to early to mid-October), emergence, tillering (early November), and overwintering (late November to early December). Taking Xiqiaozhuang Village in Hongtong County as an example, during the 0–120 d after sowing, the 8 d MODIS ET showed no significant change, with ET values remaining relatively constant. While the irrigated sample points showed a more pronounced trend in MODIS ET value changes compared to the dryland sample points, there was still a phenomenon of delayed response to vegetation growth information due to pre-planting irrigation. During the main growth stages of winter wheat (from re-greening to maturity), there was significant improvement in the trend of MODIS ET changes compared to earlier stages. However, there still remained a considerable discrepancy with the characteristic ET curve of winter wheat. ESTARFM ET, building upon the temporal information provided by MODIS ET, utilized the spatial information from the higher resolution SEBAL ET to mitigate the impact of mixed pixels on crop growth information. Its trend and magnitude of change were more consistent with the existing ET characteristic curve, thus better capturing the characteristics of ET variation throughout the winter wheat growing period.

3.3. Assimilation ET

The 4DVar, EnKF, and PF methods were used to calculate the assimilated ET of the winter wheat growth period simulated by CERES-Wheat on a daily scale and the ESTARFM ET, respectively (Figure 7). Throughout the entire growth period of winter wheat, all three assimilated ET methods exhibited consistent trends with simulated ET. That is, the minimum ET values occurred during the sowing and re-greening stages, with the highest rate of increase observed during the jointing stage, and the maximum values occurred during the wheat tasseling-irrigation stage. Comparing 4DVar ET, EnKF ET, and PF ET revealed that PF ET exhibited inferior assimilation performance at rainfed sites (Xiqiaozhuang Village in Hongtong County, Tundou Village in Yaodu District, and Cangtou Village in Hejin City). Specifically, this was manifested by an overestimation of the maximum ET values during the tasseling-irrigation stage. Additionally, in Hongtong County, PF ET showed a trend discrepancy between days 100 and 175 after sowing compared to the growth characteristics of winter wheat during this period. At the sampling points in Xiangfen County and Houma City, PF ET showed a delayed response to sudden ET changes before day 50 after sowing and between days 150 and 210. Furthermore, in Quwo County, there was a significant overestimation of ET after 100 days post-sowing. 4DVar ET performed good fitting at most of the sampling points, with only a reverse trend observed between days 175 and 210 after sowing at the Houma point, attributed to the poor quality of the ESTARFM ET during this period. At the remaining points, 4DVar ET demonstrated satisfactory simulation results and accurately captured the trend of winter wheat ET throughout the growth period. Except for Xiqiaozhuang Village in Hongtong County and Cangtou Village in Hejin City, EnKF ET showed satisfactory performance in simulating the variation trend of winter wheat ET and responding to sudden moisture information changes at other sampling points. However, at Xiqiaozhuang Village in Hongtong County and Cangtou Village in Hejin City, ESTARFM ET significantly underestimates between days 100 and 175 after sowing, leading to considerable simulation errors for EnKF ET during this period.

Based on the above research, the correlation between ESTARFM ET and assimilated ET was calculated using 4DVar, EnKF, and PF methods (Figure 8). It is evident that the observed ESTARFM ET data correlate well with assimilated ET, with R² values of 0.5142, 0.7119, and 0.5596, respectively, all passing significance tests at p < 0.01. This indicated that assimilated ET from different assimilation methods could effectively incorporate the remote sensing information parameters expressed by ESTARFM ET. Among them, the correlation between ESTARFM ET and EnKF ET was the highest with the smallest systematic error (intercept = 0.5024).

4. Discussion

ET-flux data typically require expensive instruments and equipment, with corresponding capital investment for maintenance and calibration, and measurements from specific locations may not accurately reflect evapotranspiration over a larger area, leading to spatial inconsistencies. ESTARFM ET effectively preserves the temporal variations of MODIS while reflecting the spatial details of higher resolution data. The ESTARFM is based on reconstruction and assumes that the reflectance of pixels in low spatial resolution data can be regarded as a linear combination of reflectance in high spatial resolution pixels. The algorithm sets the lower limit of ET to 0 without constraining the upper limit, resulting in the fused ET having a lower limit consistent with SEBAL ET for the corresponding time but an upper limit smaller than the corresponding SEBAL ET. ESTARFM ET maintains consistency in tone and trend with SEBAL ET for the corresponding time, indicating that higher quality data sources and sufficient data quantity can reduce the uncertainty of fusion models, especially when underlying surface heterogeneity is strong. However, due to the need for two sets of corresponding high and low spatial resolution data, the ESTARFM fusion method is susceptible to large biases in weighting calculation when there are significant changes in land cover abundance between the two input time periods. This can lead to substantial errors in fusion results, particularly during long periods of missing high-resolution data (such as poor data quality of Landsat 8 during winter). Although MOD16A2 data are processed through cloud detection and filtering, remote sensing data may be incomplete or inaccurate under severe weather conditions (e.g., frequent cloud cover), which can affect the reliability of the data. Therefore, we use the ESTARFM fusion model to improve the data completeness and accuracy [46]. MODIS has a spatial resolution of 250–1000 m, which may not adequately capture detailed changes in small-scale study areas and may result in some small-scale changes being overlooked. Future research could explore the use of high-resolution remotely sensed data (e.g., Sentinel-2) to improve ET estimates, especially at local scales and in areas of complex topography.

Due to limitations imposed by satellite revisit cycles and adverse weather conditions (clouds, rain, and snow), Landsat 8 only provided seven observational datasets during the winter wheat growing period. The considerable dispersion of these datasets hinders continuous monitoring of crop growth throughout the growing period. Additionally, remote sensing data can capture the instantaneous physical condition of the vegetation canopy [31]. By incorporating ESTARFM ET, timely adjustments can be made to simulated ET, enhancing responsiveness to changes in instantaneous meteorological factors. However, the PF assimilation method demands high requirements for particle initialization and resampling strategies. Improper selection may lead to rapid convergence or divergence of results, resulting in significant assimilation biases. Moreover, particle degeneration issues may restrict weight sampling, causing the loss of some instantaneous information and delaying responses to abrupt changes [47]. The 4DVar assimilation method involves using remote sensing information to construct a cost function that is then optimized through adjustment of model initial parameters using optimization algorithms to minimize the cost function, thereby achieving optimization objectives. However, this assumption can lead to significant simulation errors when remote sensing information is inaccurate. Additionally, this algorithm requires the setting of an assimilation time window during the assimilation process, where the cost function is minimized within this window to optimize the initial ET. Consequently, when setting the assimilation window to 10 d, there may be a lack of assimilated ET values for the subsequent 10 d. The Kalman filter assumes that noise follows a Gaussian distribution and is based on the minimum variance estimation criterion, performing well in linear system simulations [44,45]. The sequential filtering algorithms, represented by EnKF, integrate continuous canopy observations into model simulations to reduce errors in assimilated state variables, thereby enhancing the accuracy of model simulations. Compared to other assimilation algorithms, the EnKF formula is concise, computationally efficient, and more adept at capturing the nonlinear, high-dimensional characteristics of crop growth models. Previous studies have demonstrated that under conditions of high-quality observational data for model calibration, EnKF can achieve higher simulation accuracy [10]. In the Linfen Basin, after the winter wheat enters the tasseling-irrigation stage, there is a noticeable increase in rainfall, leading to significant errors in SEBAL ET inversion due to cloud cover. This results in systematic errors. Additionally, during this period, ET values are relatively high. When ESTARFM ET exceeds 4, EnKF ET shows a divergent trend, leading to a significant decrease in simulation accuracy. The correlation between ESTARFM ET and PF ET is second only to that of ESTARFM ET and EnKF ET. When ESTARFM ET exceeds 2, the divergence trend of PF ET increases, indicating that excessive reliance on initial state estimates can cause rapid convergence or divergence of assimilation results. EnKF is an ensemble-based Kalman filtering method, which is very effective in dealing with linear or weakly nonlinear systems, but may lead to filter divergence under highly nonlinear or complex dynamic systems, especially in the presence of extreme conditions such as high evaporative transpiration. Local Particle Filtering (L-PF), on the other hand, is better able to cope with nonlinear systems by introducing random particles to approximate the a posteriori probability distribution, thus providing more flexible and accurate state estimation. By combining the two, the fast computational advantage of EnKF and the powerful nonlinear processing capability of particle filtering complement each other, which can effectively improve the stability and accuracy of the model under extreme conditions. Future research focuses on exploring how to further optimize the computational efficiency of hybrid assimilation methods and incorporate high-performance computing techniques for wider application in complex environmental models. In addition, combining more advanced filtering algorithms, such as variational methods and deep learning, is considered to further improve the accuracy and computational efficiency of the model, especially for applications under complex and extreme climate conditions.

After an initial increase during the sowing period, winter wheat ET diminishes during the tillering and overwintering stages, primarily due to soil evaporation, with temperature serving as the primary limiting factor for transpiration during this period. Reduced temperatures slow the growth of winter wheat organs such as leaves, stems, and roots, leading to less biomass accumulation and consequently lower water demand, resulting in a decreasing trend in ET [31]. Due to the compensatory growth effect of crops, water stress occurs before and after the re-greening stage, followed by normal irrigation during the reproductive stage. The impact of water deficit on plant height, final tiller number, leaf area, and dry matter accumulation is minimal during the re-greening stage, indicating low sensitivity to water stress. During the rapid growth stage to the mid-reproductive stage, ample sunlight ensures simultaneous growth of winter wheat organs such as leaves, stems, and roots. The length and volume of leaves and stems increase exponentially or tens of times, while the root system rapidly extends, leading to more apparent absorption and utilization of soil moisture in deeper layers. Plant transpiration increases, while soil evaporation tends to increase due to rising temperatures. ET reaches its maximum during the tasseling-irrigation stage. Insufficient soil moisture supply during this growth stage can affect the efficiency and speed of photosynthesis and reduce the synthesis of starch, protein, and organic matter, significantly decreasing winter wheat grain weight and ultimately impacting its growth and yield [48]. It is important to ensure that sufficient water is available to meet the crop’s growth requirements at this stage, when winter wheat is most sensitive to water deficit [49]. From the heading to the tasseling-irrigation stage to the maturity stage, there is a marked decrease in ET and winter wheat growth comes to a virtual standstill. Leaves become dry and yellow, indicating the end of the process of transporting nutrients to the grain. Leaves begin to turn yellow, photosynthetic capacity decreases, plants begin to wither, roots begin to die, and plant transpiration decreases, leading to a reduction in ET. In addition, late irrigation during the maturity stage can cause new tillers to emerge from the base of old stems, consuming nutrients and leading to a reduction in thousand kernel weight. Therefore, reducing field irrigation during this growth stage can help to increase yield.

5. Conclusions

Incorporating remote sensing fusion techniques into the assimilation of remote sensing data with crop growth models, various assimilation algorithms (4DVar ET, EnKF ET, and PF ET) were coupled to determine the optimal assimilation algorithm suitable for the winter wheat planting area in the Linfen Basin. By improving the accuracy of monitoring changes in ET over the full life span of winter wheat, policy makers can better predict changes in water demand as the crop grows.

(1)

Utilizing the SEBAL and ESTARFM, an 8 d and 30 m time series of ET was obtained for the winter wheat growing season. This approach effectively preserved the temporal variation information of the MODIS while spatially reflecting the detailed information of the Landsat 8.

(2)

Using the calibrated CERES-Wheat, the variation curve of the entire growth period ET of winter wheat was obtained. The trend and magnitude of the ESTARFM ET variation closely matched the characteristic curve of CERES-Wheat ET after S-G filtering, effectively capturing the variations in winter wheat ET throughout the growing period.

(3)

The correlation between EnKF ET and ESTARFM ET (R² = 0.7119, p < 0.01) was significantly higher than that of 4DVar ET (R² = 0.5142, p < 0.01) and PF ET (R² = 0.5596, p < 0.01). EnKF ET demonstrated a more ideal response to changing trends and abrupt moisture information. However, due to the significant underestimation of ESTARFM ET from 100 to 175 days after planting, there was notable simulation error in EnKF ET at sample points in Hongtong County and Hejin City during this period.