Estimating evapotranspiration at high spatio-temporal resolution based on the Bayesian model averaging method in Baoding Plain, China

Abstract

Evapotranspiration (ET) is a core parameter of the hydrology and carbon cycles, and its accurate estimation is crucial for water resource management. Satellite-based ET products provide an effective means for large-scale monitoring. However, due to limitations in the spatial and temporal resolution, the use of these products at regional and field scales is limited. In this study, daily ET in the Baoding Plain—a key groundwater resource recharge area—was estimated at a 500 m spatial resolution using the Bayesian Model Averaging (BMA) method. The model was driven by a synthesis of remote sensing datasets, reanalysis products, and interpolated data from meteorological stations. Validation results from in-situ observations indicated that the BMA ET had better performance (R=0.83, RMSE=1.25 mm/d) than each model in the BMA scheme. The spatiotemporal analysis revealed that the average annual BMA ET in the Baoding Plain was 683 mm/year from 2000 to 2019. Seasonal and monthly variations in the BMA ET captured the irrigation and water consumption patterns of the local crop rotation systems. A significant increasing trend of BMA ET (2.40 mm/year²) was observed in the Baoding Plain over the study period. At the regional scale, ET over more than 50% of the plain exhibited a significant positive trend. Further analysis identified water availability, solar radiation, and temperature as the primary drivers of ET variation. The BMA ET product generated in this study is characterized by high spatiotemporal resolution and accuracy. This reliable, high-resolution dataset offers valuable support for precision agricultural water management and hydrological studies, including groundwater investigations, in this predominantly agricultural region.

Graphical abstract

Keywords

Evapotranspiration / Satellite-based model / Random Forest model / Climate change / Water resources

Introduction

Evapotranspiration (ET) refers to the transfer of water in vapor form from liquid and solid states at the land surface back to the atmosphere. ET is a critical parameter for assessing the global water cycle process and water resources evaluation, and it can be used to track the interaction between the land surface and the atmosphere. Meanwhile, ET is not only crucial to the water cycle but also plays a significant role in the carbon and energy cycles (Zheng et al. 2025). Accurate assessment of ET can provide the basis for studies of climate change, terrestrial ecosystem dynamics, weather prediction, and hydrological modeling, as well as agricultural irrigation and planning (Elzain et al. 2024).

Over the last few years, with the development of Earth observation science and technology, a variety of ET models based on remote sensing have been applied at broader regional scales. As a result, a number of well-established data products have been developed. A wide variety of ET models are currently available, which can be broadly categorized into empirical methods, process methods based on the principle of energy balance, and other models such as machine learning and trigonometric algorithms, etc. However, the uncertainty in different ET models is close to 50% (Mu et al. 2011). Furthermore, despite advances in understanding the temporal dynamics of ET through modeling, the capacity to represent fine-scale spatial heterogeneity is often hampered by the limited resolution of current datasets. For example, GLEAM ET has a high temporal resolution (daily) but suffers from a coarse spatial resolution (0.1°), whereas MODIS has a fine spatial resolution (500m) but a low temporal resolution (8-day). Consequently, high-precision and high spatiotemporal resolution ET estimation in a specific area continues to be a formidable challenge.

The Bayesian Model Averaging (BMA) method combines multi-model probabilistic forecasts and post-processing to estimate uncertainty based on model selection, combination, and prediction. BMA has been proven that it could produce reliable and accurate simulations, and it can outperform individual models. The results from many studies using the BMA method for integrated ET assessment show that the estimations of this method are generally better than those of every single model (Yang et al. 2021). In addition to reducing the root mean square error and improving the correlation coefficient of individual ET datasets, BMA has shown robust performance in integrating diverse models for water level prediction and ecological footprint analysis. This makes it an effective technique for mitigating model uncertainty in ET estimation.

The Baoding Plain, a typical semi-arid agricultural area, located in the middle of the Hebei Plain, plays a multifaceted role as a key development zone for the Beijing-Tianjin-Hebei region, and a crucial recharge area for groundwater resources (Cao et al. 2020). The influence of the Taihang Mountains in the west of the Baoding Plain leads to the uneven spatial distribution of precipitation in this region (Bai et al. 2023). ET here exhibits distinct spatiotemporal variations, driven by diverse meteorological factors and human activities. Recent years have witnessed significant alterations in the local water cycle and ET fluxes, attributable to the operation of the South-to-North Water Diversion (SNWD) Project, measures to control groundwater over-exploitation, sustained ecological river recharge, and consistently abundant precipitation (Cao et al. 2020). Therefore, this study selects the Baoding Plain as a representative area of the North China Plain (NCP), to develop accurate ET evaluation methods at a fine-resolution scale, thereby providing a scientific reference for regional water cycle and water resources assessment.

This study aims to generate high spatiotemporal resolution daily ET data for the Baoding Plain through multi-model integration. The prime objectives are: (1) to estimate daily ET using the BMA method; (2) to validate the BMA-integrated results; and (3) to analyze the spatiotemporal patterns and the driving factors of ET derived from the BMA approach in the Baoding Plain.

1 Materials and Methods

1.1 Study area

The Baoding Plain is situated in central Hebei Province, northern China (114°30'–116°20'E, 38°20'–39°60'N) (Fig. 1a). This is a gently sloping piedmont terrain, with a gradient of less than 2‰. The overall terrain is higher in the west and lower in the east, transitioning from the Taihang Mountains in the west to a flat alluvial plain. While most of the plain lies below 30 meters in elevation, a small western section reaches 80-200 meters (Fig. 1b). The region has a northern temperate, semi-arid monsoon climate (Bai et al. 2023). Dominated by cropland (Fig. 1c), its primary crops are wheat and corn, with some rice and soybean cultivation. The drainage system originates from the eastern foot of the Taihang Mountains and belongs to the Daqing River system in the Haihe River Basin. The river network forms a fan-shaped pattern, with three main channels converging into Baiyangdian Lake in the south, center, and north (Fig. 1d).

1.2 Data

1.2.1 Site Data

Six flux sites around the Baoding Plain were selected to examine the ET product performance (Fig. 1a). Due to the lack of available flux sites in the study area, we selected these six sites surrounding the study area with similar surface characteristics and climate conditions. The six sites are CN-Bed, DX2, GT, HL, MY, and YC. Five of them were covered by cropland. Details of all sites are shown below (Table 1):

The core advantage of the BMA method lies in its ability to learn the performance weights of each model under different climatic conditions from the limited synchronous observation data, and it has spatiotemporal transferability (Shen et al. 2025). The model weights derived from these sites can be spatiotemporally extrapolated. Although the observation time range of the flux sites cannot fully cover the time range of the simulation (Table 1), the verification results and simulation results of the BMA model are reliable.

Before the flux data were applied, quality control was conducted to ensure the accuracy of the observations. The vegetation index, including Leaf Area Index (LAI) and Normalized Difference Vegetation Index (NDVI) at the site scale, was obtained from the relevant website ( https://daac.ornl.gov/MODIS/). These datasets were screened using Quality Assurance (QA), and continuous NDVI and LAI were generated based on the linear interpolation method.

1.2.2 Regional data

This study employed a combination of regional meteorological and remote sensing datasets. Meteorological data were obtained from The National Meteorological Information Center ( https://data.cma.cn/en). Net radiation (R_n) was acquired from the GLDAS 2.1 reanalysis product. To achieve the target spatial resolution, the meteorological variables were interpolated using the ANUSPLIN software.

Remote sensing data included NDVI, LAI, land cover type, Digital Elevation Model (DEM), and fraction of Photosynthetically Active Radiation (fPAR). With the exception of the DEM (from SRTMv4.1, https://srtm.csi.cgiar.org/), all these datasets were sourced from MODIS products accessed via the Level-1 and Atmosphere Archive & Distribution System Distributed Active Archive Center (LAADS DAAC). Specifically, NDVI was derived from MOD13Q1, LAI and fPAR from MOD15A2, and land cover information from MCD12Q1 ( Table 2).

1.3 Methods

1.3.1 ET model

Three processes-based ET models were chosen in this study, including the PM-MOD16 model (Mu et al. 2011), the PT-JPL model (Fisher et al. 2008) and the improved RS-PM model (RRS-PM model) (Yuan et al. 2010). Chen et al. (2014) compared eight ET models in China. The results showed that these three models had better model performance. The models with good performance could increase the accuracy of BMA estimates (Chen et al. 2015). Therefore, these three models were used in this study within the BMA framework. To extend the models involved in the BMA process, a machine learning approach, i.e. the Random Forest model (Breiman, 2001), was used in this study for the integration process.

PM-MOD16 is the most commonly used remote sensing model, which was modified from the Penman-Monteith (P-M) equation. The PM-MOD16 model is a typical algorithm to couple the P-M equation with remote sensing data. The parameterization of effective resistance of land surface evaporation and canopy transpiration is the core of the PM-MOD16 model algorithm. In this model, relative humidity (R_h) and water Vapor Pressure Deficit (VPD) are used to describe the moisture content in the soil, and specific conductance terms are used to describe the water transport process between the land and the atmosphere. The specific calculation process consists of three main parts, as shown in the following equations:

(1)$ ET={ET}_{i}+{ET}_{t}+{ET}_{s} $

Where: ET_i represents the evaporation of rainwater trapped by the canopy; ET_t stands for vegetation transpiration; ET_s represents the soil evaporation. The total ET is composed of these three parts.

The PT-JPL model is a simple remote-sensing ET model developed by Fisher et al. (2008), based on the Priestley-Taylor (PT) method. Due to its high simulation accuracy, the model is already widely applied in the estimation of site, regional, and global ET. It uses meteorological and remote sensing data to drive the model and estimate actual ET from potential ET. In this model, the total ET consists of three parts: Canopy evaporation, soil evaporation, and canopy interception, as shown below:

(2)$ ET={ET}_{c}+{ET}_{s}+{ET}_{i} $

Where: ET_c represents canopy transpiration; ET_s represents soil evaporation; ET_i represents canopy interception evaporation.

The RRS-PM model is revised from the method proposed by Cleugh et al. (2007). The model has been widely used in regional ET estimation because of its simple parameterization scheme. The RS-PM model improved by Yuan et al. (2010), also known as RRS-PM, added the soil evaporation component, and upgraded the canopy conductance by using LAI. At the same time, the RRS-PM method distributed R_n between the soil surface and the canopy using the Beer-Lambert law and used the temperature and humidity limits of stomatal conductance (Yuan et al. 2010). The temperature limiting equation used to control stomatal conductance in this model is:

(3)$ {R}_{nc}={R}_{n}-{R}_{ns} $

(4)$ {R}_{ns}={R}_{n}exp (-{k}_{Rn}LAI) $

Where: R_nc represents the part of the net radiation at the canopy, R_ns is the part of the net radiation that reaches the soil, and k_Rn is the extinction coefficient.

The equation controlling the temperature limit of stomatal conductance is:

(5)$ {m}_{c}=exp\left(-{\left(\frac{{T}_{m}-{T}_{f}}{{T}_{f}}\right)}^{2}\right) $

Where: T_m represents air temperature, and T_f is the optimum temperature.

The Random Forest model is a tree-based non-parametric machine learning algorithm for classification and regression developed by Breiman (2001). It is composed of multiple prediction trees. As an ensemble method, it is less sensitive to outliers, effectively reduces the risk of overfitting, and has high prediction accuracy and stability (Mohtaram et al. 2025). In this study, the Random Forest model was constructed by using site data, where the dataset was divided into an 80% training set and a 20% test set to rigorously evaluate the performance of each model.

1.3.2 BMA model

The BMA method obtains a weighted average of all model results by assembling the results of other candidate models to explain the uncertainty of the model and improve the simulation (Yang et al. 2021). This is a statistical treatment method for probability prediction with multi-model integration. Based on four ET models, the BMA method is used to improve the accuracy of the estimated daily ET. In this method, the dependent variable y is used to represent the predicted value based on the BMA, that is, the estimated ET; y_T represents the training data with training length T; X{x₁, x₂, ···, x_k} is the set of predictions of k different models. According to the conditional probability sum formula, the posterior probability distribution of the model can be obtained as follows:

(6)$ p\left(y|{x}_{1},{x}_{2},{x}_{3},\cdots ,{x}_{k}\right)=\sum \limits_{k=1}^{K}p(y|{x}_{k})\cdot p({x}_{k}|{y}_{T}) $

Where: p(y|x_k) refers to the prediction Probability Density Function (PDF) obtained through the simulation of the x_k model; p(x_k|y_T) is the posterior probability of x_k obtained by the model simulation; p(x_k|y_T) is a statistical weight w_k, reflecting the matching degree between x_k and y_T. If p(x_k|y_T) is represented by w_k, it can be expressed as follows:

(7)$ p\left({x}_{k}|{y}_{T}\right)=\sum \limits_{k=1}^{K}{w}_{k}=1 $

Then, by substituting equation (7) into equation (6), we get:

(8)$ p\left(y|{x}_{1},{x}_{2},{x}_{3},\cdots ,{x}_{k}\right)=\sum \limits_{k=1}^{K}p(y|{x}_{k})\cdot {w}_{k} $

Assuming that p(y|x_k) obeys a Gaussian distribution, whose mean is x_k and variance is w_k², the Gaussian distribution and the parameter vector can be expressed as:

(9)$ p\left(y|{x}_{k}\right)=g(y|{\theta }_{k}) $

(10)$ {\theta }_{k}=\left\{{x}_{k}\right.,\left.{{{w}_{k}}}^{2}\right\} $

Where: g represents a Gaussian distribution.

By combining formula (8) and (9), we can get:

(11)$ p\left(y|{x}_{1},{x}_{2},{x}_{3},\cdots ,{x}_{k}\right)=\sum \limits_{k=1}^{K}g(y|{\theta }_{k})\cdot {w}_{k} $

The calculation process based on the log-likelihood function is simpler than the likelihood function, so converting the above formula into the log-likelihood function can be approximated as:

(12)$ l({\theta }_{1},{\theta }_{2}{,\cdots ,\theta }_{k})=\sum \limits_{\left(s,t\right)}log \left[\sum \limits_{k=1}^{K}g\left({y}_{s,t}|{\theta }_{k}\right)\cdot {w}_{k}\right] $

Where: y_s,t is the target data of time and place, ∑_(s,t) represents the sum of ET observed values at each point. Due to the need to accurately estimate the weight and variance of each model in the set when using this method, Raftery et al. (2005) estimated the weight and variance of BMA by using the expectation-maximization method. It is also important to note that the application of the BMA should be consistent with the data used by the Gaussian distribution, so this paper applies the Box-Cox method to preprocess the data.

1.3.3 Statistical analysis

The statistical indicators in this study include the correlation coefficient (R), centralized Root Mean Square Deviation (RMSD), Root-Mean-Square Error (RMSE), bias error, and Nash-Sutcliffe efficiency coefficient (NSE). NSE is a comprehensive index that considers the R and RMSE of the model. The calculation formulae of all the above evaluation indicators are expressed as:

(13)$ {\sigma }_{0}=\sqrt{\sum \limits_{t=1}^{n}{\left({{O}_{t}}-\overline{O}\right)}^{2}/n} $

(14)$ {\sigma }_{M}=\sqrt{\sum \limits_{t=1}^{n}{\left({{M}_{t}}-\overline{M}\right)}^{2}/n} $

(15)$ RMSD=\sqrt{\frac{1}{n}\sum \limits_{t=1}^{n}{\left[\left({O}_{t}-\overline{O}\right)\left({M}_{t}-\overline{M}\right)\right]}^{2}} $

(16)$ R=\frac{\displaystyle\sum \limits_{t=1}^{n}\left({O}_{t}-\overline{O}\right)\left({M}_{t}-\overline{M}\right)}{n{\sigma }_{0}{\sigma }_{M}} $

(17)$ RMSE=\sqrt{\frac{1}{n}\sum \limits_{t=1}^{n}{[\left({O}_{t}-{M}_{t}\right)}^{2}} $

(18)$ bias=\overline{O}-\overline{M} $

(19)$ NSE=1-\frac{\displaystyle\sum \limits_{t=1}^{n}{({{O}_{t}}-{{M}_{t}})}^{2}}{\displaystyle\sum \limits_{t=1}^{n}\left({{M}_{t}}-\overline{M}\right)^{2}} $

Where: n is the total number of observations; O_t is the observed value of ET; $ \overline{O} $ represents the mean of the observed values; M_t represents the predicted value of ET; $ \overline{M} $ is the average of the predicted ET; ${\sigma }_{0} $ and ${\sigma }_{{\mathrm{M}}} $ represent the deviations between the observed and predicted ET values, respectively; bias is a deviation error.

In addition, NSE represents the variance of the residual, and the closer its value is to 1, the better the simulation performance. If the value of NSE is less than 0, the performance of the model is not acceptable. Bias is used to describe the degree to which the evaluated value deviates from the observed value, and the optimal value is 0. The closer the value to 0, the better the performance of the model.

1.3.4 Uncertainty analysis

In this study, we employed a site-grouping approach and conducted Monte Carlo simulation analyses for each site individually to assess the uncertainty of the BMA ET relative to the observational ET. The Monte Carlo sampling method, grounded in statistical sampling theory, utilizes computer-based large-scale random sampling to derive statistical characteristics of relevant random variables (Yan et al. 2018). For each site, this method was applied to compute key statistical metrics—including the 90% confidence interval, mean, and standard deviation—of the BMA ET predictions, thereby quantifying the associated uncertainty.

To evaluate the confidence intervals, two statistical indicators were used: The Containing Ratio (CR) and the average Relative Deviation amplitude (RD). The CR, ranging from 0 to 1, represents the proportion of observed values falling within the confidence interval. A higher CR indicates better performance, with CR=1 denoting a perfect simulation. The RD is a dimensionless metric quantifying the deviation between estimated and observed values, ranging from 0 to +∞. An RD value of 0 signifies a perfect simulation (Huang et al. 2019).

2 Results

2.1 Model and climate data validation

To assess the reliability and quantify the uncertainty of the regionally interpolated datasets, we validated them against the flux tower measurements (Table 3). The interpolation results showed strong agreement with the observations, with R values of 0.87 for R_n, 0.99 for T_m, 0.93 for R_h and 0.82 for atmospheric pressure (Ps) (all statistically significant at P<0.05). Further supporting their accuracy, the bias and RMSE values were consistently low for all variables. The NSE of these climate variables were all close to 1 (NSE_Rn = 0.60, NSE_Tm = 0.98, NSE_Rh = 0.84, NSE_Ps = 0.65). These comprehensive metrics collectively demonstrate that the interpolated datasets are scientifically robust and reliable for subsequent analysis.

The performance of the BMA model was benchmarked against four individual models by comparing their estimates with site-scale observations. As shown in Fig. 2, the scatter points of the BMA method were more densely clustered along the 1:1 line than those of the four single models. Quantitatively, the BMA model achieved the highest NSE value, which was closest to 1, alongside the strongest correlation (R=0.83) and the lowest root mean square error (RMSE=1.25 mm/d) and bias. Although the four individual models participating in the integration were themselves stable and reliable (all with NSE ≥ 0.5), the BMA integration consistently outperformed each of them across all metrics.

The BMA-based ET estimates were validated against observations from 6 Eddy Covariance (EC) sites and demonstrated superior overall performance compared to the four individual models (Fig. 3). The box plots showed that the medians of the BMA and PM-MOD16 simulations were closer to the observed median (Fig. 3a). Both the BMA simulations and the in-situ observations exhibited positively skewed distributions. While the Random Forest model achieved the smallest RMSE and a bias closer to 0 at some single sites (Fig. 3b, d), the BMA method performed most consistently across sites. Furthermore, the BMA ensemble ET showed significantly higher R values at most sites compared to the single models (Fig. 3c). The NSE values also confirmed the good performance of the BMA method (Fig. 3e). In summary, when comprehensively considering the RMSE, R value, bias, and NSE, the BMA approach yielded the most robust and accurate results.

2.2 Spatial and temporal distribution characteristics of ET

Fig. 4 illustrates the spatial distribution of the annual average ET across the Baoding Plain, as estimated by the BMA method. The regional mean annual ET was 683 mm. Lower ET values were predominantly observed in the central, northeastern, and southwestern parts of the plain, which correlated with grassland areas. In contrast, the highest ET values were concentrated in the northwestern and southeastern regions, characterized by extensive croplands. Overall, ET exhibited a distinct mosaic pattern of high and low values throughout the region (Fig. 4).

Fig. 5 depicts the spatial distributions of seasonal ET across the Baoding Plain. As expected, the lowest ET values occurred in winter (December-February), with nearly the entire plain exhibiting values below 50 mm (Fig. 5d). In contrast, summer (June-August) showed the peak ET, with values exceeding 300 mm across the region (Fig. 5b). Spring (March-May) and autumn (September-November) displayed intermediate values, although ET in autumn was generally lower than in spring (Fig. 5a, 5c). Notably, the spatial heterogeneity of ET was most pronounced in spring and summer. Meanwhile, the ET values showed an increasing trend from spring to summer. Consequently, the spatial pattern of the annual total ET was predominantly determined by the distributions during these two seasons.

The spatial and temporal variations of monthly ET across the Baoding Plain are influenced by the winter wheat-summer corn rotation system and meteorological conditions (Fig. 6). Overall, ET exhibited a distinct seasonal pattern, with low values in winter months, a gradual increase in spring months, and a decline in autumn months. During January and February, ET remained very low (< 20 mm) throughout the region (Fig. 6a, 6b) due to low temperatures. From March to May, ET increased steadily (Fig. 6c, 6d, 6e), reflecting the growth of winter wheat and rising temperatures. Compared to May, ET values in June showed a decrease across most of the region, notably in the central and southwestern parts (Fig. 6f). This is caused by the harvest of winter wheat. The highest ET occurred in July and August (Fig. 6g, 6h), driven by the growth of summer corn, strong solar radiation, higher precipitation, and higher atmospheric demand. After September, ET gradually decreased, returning to low levels in late autumn and winter (October-December) (Fig. 6i-6l).

2.3 Spatial-temporal Long-term change trend of ET

The interannual and seasonal variability of the BMA ET over the Baoding Plain from 2000 to 2019 is shown in Fig. 7. The annual ET series exhibited a statistically significant increasing trend (P<0.05) at a rate of 2.40 mm/year² (Fig. 7a), with values ranging from 643.06 mm/year to 723.99 mm/year across the study period. Notably, sharp increases in ET occurred in 2008 and 2011, while pronounced decreases were observed in 2009 and 2014. The seasonal cycle followed a distinct double-peak pattern, with peaks occurring in May and August (Fig. 7b), which is consistent with the spatial patterns shown in Fig. 6. The ET values increased steadily from January to May and decreased steadily from August to December.

The spatial distributions of the BMA ET trend are shown in Fig. 8. From 2000 to 2019, 54.25% of the total area experienced a significant increase in ET. The largest positive trends (>5 mm/year²) were concentrated in the eastern and southern parts. In contrast, only 3.46% of the total area showed a decreasing trend (ranging from −5 mm/year ² to 0 mm/year²). The remaining 42.29% of the area showed no statistically significant trend.

2.4 Contributions of climate and vegetation variables to changes in ET

The spatial distributions of the key climate variables, including R_n, T_m, air water conditions (VPD and R_h), surface moisture (Pre) and wind speed (W_s), are presented in Fig. 9. Relatively lower T_m and higher R_n were observed in the northern areas, whereas the southern areas exhibited the opposite pattern (Fig. 9a, 9b). The spatial patterns of VPD and R_h were relatively uniform across most of the plain (Fig. 9c, 9d). Areas proximal to the mountains received higher Pre and experienced lower W_s; conversely, regions farther from the mountains showed lower Pre and higher W_s (Fig. 9e, 9f). The highest NDVI and LAI were concentrated in the southwestern and central parts of the region (Fig. 9g, 9h).

Spatial trends of key climatic and vegetation variables from 2000 to 2019 are mapped in Fig. 10. Pronounced declining trends in R_n were concentrated in the central and southwestern Baoding Plain (Fig. 10a). T_m and VPD exhibited significant increasing trends in the southern regions (Fig. 10b, 10c), whereas R_h showed an opposite (i.e., decreasing) trend to VPD in these areas (Fig. 10d). Pre trends were largely non-significant across the study area (Fig. 10e). W_s declined significantly in the northern and eastern zones (Fig. 10f). In contrast, the vegetation indices (NDVI and LAI) showed significant and extensive increasing trends across nearly the entire plain (Fig. 10g, 10h).

The interannual trends of key variables from 2000 to 2019 are shown in Fig. 11. The trend in R_n (−0.21 MJ/m²/day/year) was marginally significant (P = 0.051), while T_m and Pre showed no statistically significant trend over the study period (Fig. 11a, 11b, 11e). VPD, NDVI, and LAI increased significantly at rates of 0.045 hPa/year, 0.003/year, and 0.038 m²/m²/year, respectively (Fig. 11c, 11g, 11h). In contrast, R_h and W_s exhibited significant decreasing trends (Fig. 11d, 11f).

The spatial patterns of Pearson's correlation coefficients between environmental variables and ET are presented in Fig. 12. Analysis revealed that ET variations were primarily controlled by radiation, temperature, and moisture conditions. Specifically, in the southwest part of the plain, ET was significantly positively correlated with R_n and significantly negatively correlated with T_m (Fig. 12a, 12b). A strong negative correlation between ET and VPD was widespread across most of the region (Fig. 12c), while a positive correlation with R_h was observed in the south (Fig. 12d). Pre was positively correlated with ET over a large portion of the study area (Fig. 12e). In contrast, W_s, NDVI, and LAI showed no strong correlation patterns with ET (Fig. 12f, 12g, 12h). In summary, water availability, radiation, and temperature were the dominant factors controlling ET dynamics in the Baoding Plain.

2.5 Uncertainty analysis

The uncertainty analysis (Table 4) indicated that the DX2 station recorded the smallest RD value and the largest CR value, reflecting the BMA model's excellent performance and high reliability at this location. Conversely, the GT site exhibited the smallest mean error and standard deviation, as well as the narrowest uncertainty range, indicating stable model performance with consistently low uncertainty. Overall, metrics including standard deviation, confidence interval width, and RD collectively capture the inherent uncertainty of the model. Notably, the BMA model showed distinctlydifferent behavior at the CN-Bed site compared to others, a pattern potentially influenced by site-specific factors including underlying surface characteristics and local meteorological variability.

3 Discussion

3.1 Comparison of ET distribution and variation with other studies

Research focusing specifically on ET in the Baoding Plain remains limited. To contextualize our findings, we therefore compare them with studies conducted across the broader NCP. For instance, Zhang et al. (2021) reported an increasing trend in ET across the NCP from 2008 to 2017, based on a daily 1 km resolution dataset. Similarly, Yang and Lei (2022) documented increasing ET trends for major cropping systems from 2005 to 2020, respectively. Although these studies encompass a larger area, their findings are highly relevant for our work. The Baoding Plain is a constituent part of the NCP, and these regional studies primarily focus on croplands that share identical crop types and agricultural practices with our study area. Consequently, the reported regional trends provide a valuable benchmark for validating the temporal patterns observed in our study. It is noteworthy, however, that a direct comparison of spatial distributions is challenging, as the larger study extent of regional analyses inherently leads to greater value ranges and diminishes the visual clarity of local heterogeneity in published maps. This confirms that our integrated dataset not only faithfully represents the spatial heterogeneity but also delineates more granular spatial features, thanks to its enhanced spatiotemporal resolution.

3.2 The response of ET to the control factors

Our analysis demonstrates that ET in the Baoding Plain is predominantly controlled by a suite of climatic variables: Pre, R_n, T_m, R_h, and VPD. This conclusion is strongly supported by regional-scale studies conducted across the NCP and Northern China, which have consistently identified the same factors as key drivers of ET spatiotemporal variability (Su et al. 2022). The tight coupling between ET and climate is clearly reflected in the interannual record, where sharp changes in ET in 2008, 2009, 2011, and 2014 (Fig. 7a) correspond to notable climate anomalies. For instance, the pronounced ET increases in the La Niña years of 2008 and 2011 ( https://psl.noaa.gov/,National Oceanic and Atmospheric Administration Physical Sciences Laboratory, NOAA PSL) provide direct evidence of this linkage. In contrast, in the years 2009 and 2014, which were El Niño years (NOAA PSL), the BMA ET decreased markedly.

These results not only show the reliability of the BMA ET products but also indicate the crucial role of ET in the climate change processes. ET is both a major influence on near-surface climate through land-atmosphere interactions and a responsive component of the terrestrial water cycle to climatic changes. Unraveling the complex, scale-dependent feedbacks between ET and climate change remains a pivotal challenge for future research.

3.3 The effects of human activities

The Baoding Plain is a representative agricultural region in China, dominated by farmland and built-up area (Fig. 13). Key anthropogenic drivers of environmental change in the region encompass agricultural practices—notably irrigation and fertilization—as well as land use and land cover (LULC) change (Ren et al. 2020). Studying ET in cropland systems is critical not only globally but also at local scales. The Baoding Plain thus serves as an ideal and typical area for exploring these interlinked processes.

LULC change affects ET by altering surface energy and moisture budgets (Chilukoti and Xue, 2020). While modeling studies generally agree on the biophysical impacts of LULC change on ET at high and low latitudes, mid-latitude responses are more complex and location-specific (Cherubini et al. 2018). The Baoding Plain, a mid-latitude region, is characterized by three ecosystem types: Croplands, grasslands, and cropland/natural vegetation mosaics (Fig. 1c). Although statistically significant, the observed trends in land cover proportions were minimal: Grassland area increased at a rate of only 0.012%/year, while cropland area decreased at a commensurate rate of 0.012%/year (Fig. 14a, 14b). The area of cropland/natural vegetation mosaics remained unchanged (Fig. 14c). Consequently, LULC change exerted a negligible influence on regional ET dynamics. This finding is consistent with broader regional studies of the NCP and Haihe River Basin, which also concluded that land use change had only a limited impact on ET variations (Sun et al. 2024).

Spatially, our results show that lower ET values in the central, northeastern, and southwestern Baoding Plain contrast with higher values in the northwestern and southeastern regions. This pattern is attributable to greater precipitation (Fig. 9e) and more intensive irrigation in the piedmont area of the Taihang Mountains. The research of Peng et al. (2020) also indicated that the precipitation in the eastern foothills of the Taihang Mountains was higher, and the Baoding Plain belongs to this region. However, natural precipitation still falls short of meeting the water needs for crop growth, and substantial irrigation is still required. These practices directly influence ET, which in turn interacts with irrigation scheduling through soil water dynamics. The additional water and enhanced vegetation lead to greater ET, driven primarily by increased evaporation from wetter soil conditions and, to a lesser extent, by transpiration from more vigorous vegetation (Fan et al. 2025).

Monthly ET follows a bimodal temporal pattern, with peaks occurring in May and August (Fig. 7b). This pattern corresponds to the irrigation water demand of the winter wheat-summer corn rotation system in the Baoding Plain. The initial ET increase in March is linked to the start of irrigation, which boosts soil moisture and promotes ET (Koch et al. 2020). The second rapid rise in June aligns with the winter wheat harvest (Zhang et al. 2021), after which the vigorous growth of summer corn (June-October) maintains high ET, thereby forming the observed double-peak pattern. Intriguingly, groundwater use from self-contained wells and for agriculture exhibit biannual peaks, coinciding with the irrigation periods of winter wheat (March-May) and summer maize (July-August) (Long et al. 2020). Groundwater serves as a vital resource in semi-arid regions, constituting approximately 40% of irrigation water consumption (Siebert et al. 2010). As an agricultural water‐dominated region (Feng and Sun, 2024), the ET of the Baoding Plain is actually influenced by human activities involving water resource allocation, such as the SNWD Project. As the world's largest water transfer diversion project, the SNWD Project has been shown to enhance surface water supply and the recovery of groundwater storage in North China (Zheng et al. 2025). Fig. 7a shows that ET increased from a lower to a higher annual mean between 2000 and 2014, after which the values plateaued at a high level. This corresponds to the operation of the central SNWD Project in December 2014. However, due to the complex interactions between ET and its controlling factors, this study cannot provide a quantitative assessment of the impact of the SNWD Project on ET. Consequently, further research that integrates the BMA ET data from this study with ancillary datasets, such as soil moisture, irrigation, and groundwater, is warranted.

3.4 Future studies and uncertainty

This study quantified ET in the Baoding Plain by applying the BMA method to integrate multiple individual models. The BMA-derived ET demonstrated higher accuracy than any single model within the ensemble. Nonetheless, several limitations inherent to the datasets and methods introduce uncertainties. First, the scarcity of flux towers within the study area constrained model validation. Although we utilized sites from both the Baoding Plain and its immediate vicinity to better represent spatial heterogeneity, the limited number of sites remains a source of uncertainty, underscoring the need for expanded observational networks in the future. Second, while the 500 m resolution of the final BMA ET product captured broad spatial patterns, higher-resolution data at the sub-grid scale would depict regional heterogeneity more accurately. Future work should therefore incorporate higher spatiotemporal resolution satellite data to refine ET estimates. Third, although the high-resolution ET dataset produced here provides valuable support for research in climate change, groundwater dynamics, and water resource management, accurately disentangling the respective impacts of meteorological factors and human activities remains a challenge, warranting further investigation.

4 Conclusion

This study has successfully generated a high-accuracy, 500-meter resolution daily ET dataset for the Baoding Plain using the BMA method. The main conclusions are as follows: (1) the BMA framework proved superior to individual models, enhancing ET estimation accuracy and robustness; (2) the dataset uncovered critical spatiotemporal patterns: A widespread increasing ET trend in over 50% of the region, with hotspots in croplands, which underscores the imperative for enhanced water conservation; (3) ET dynamics in the Baoding Plain are predominantly controlled by water availability, solar radiation, and temperature. This finding provides a scientific basis for targeted water resource management. Consequently, the high-precision BMA ET dataset generated in this study serves as a valuable resource for subsequent research in water cycle dynamics (including groundwater), climate change impacts, and agricultural water conservation.

References

[1] Bai H, Yang HF, Meng RF, et al. Chemical characteristics and evolution of groundwater in Baoding Plain. Geological Review, 2023, 69(06): 2216-2228.

[2] Breiman L. Random forests. Machine learning, 2001, 45: 5-32.

[3] Cao WG, Yang HF, Gao YY, et al. Prediction of groundwater quality evolution in the Baoding Plain of the SNWDP benefited regions. Journal of Hydraulic Engineering, 2020, 51(08): 924-935.

[4] Chen Y, Xia JZ, Liang SL, et al. Comparison of satellite-based evapotranspiration models over terrestrial ecosystems in China. Remote Sensing of Environment, 2014, 140: 279-293.

[5] Chen Y, Yuan WP, Xia JZ, et al. Using bayesian model averaging to estimate terrestrial evapotranspiration in China. Journal of Hydrology, 2015, 528: 537-549.

[6] Cherubini F, Huang B, Hu XP, et al. Quantifying the climate response to extreme land cover changes in europe with a regional model. Environmental Research Letters, 2018, 13(7): 074002.

[7] Chilukoti N, Xue YK. An assessment of potential climate impact during 1948–2010 using historical land use land cover change maps. International Journal of Climatology, 2020, 41(1): 295-315.

[8] Cleugh HA, Leuning R, Mu Q, et al. Regional evaporation estimates from flux tower and modis satellite data. Remote Sensing of Environment, 2007, 106(3): 285-304.

[9] Elzain HE, Abdalla OA, Abdallah M, et al. Innovative approach for predicting daily reference evapotranspiration using improved shallow and deep learning models in a coastal region: A comparative study. Journal of Environmental Management, 2024, 354: 120246.

[10] Fan Y, Yang Z, Lo MH, et al. Deciphering the capricious precipitation response: Irrigation impact in the north China plain. Npj Climate and Atmospheric Science, 2025, 8(1): 211.

[11] Feng Z, Sun L. Water conservation implications based on tempo-spatial characteristics of water footprint in the water-receiving areas of the South-to-North Water Diversion Project, China. Sustainability, 2024, 16(3): 1270.

[12] Fisher JB, Tu KP, Baldocchi DD. Global estimates of the land–atmosphere water flux based on monthly avhrr and islscp-ii data, validated at 16 fluxnet sites. Remote Sensing of Environment, 2008, 112(3): 901-919.

[13] Huang HP, Liang ZM, Li BQ, et al. Combination of multiple data-driven models for long-term monthly runoff predictions based on bayesian model averaging. Water Resources Management, 2019, 33(9): 3321-3338.

[14] Koch J, Zhang WM, Martinsen G, et al. Estimating net irrigation across the north China plain through dual modeling of evapotranspiration. Water Resources Research, 2020, 56(12): e2020WR027413.

[15] Long D, Yang W, Scanlon BR, et al. South-to-north water diversion stabilizing beijing's groundwater levels. Nature Communications, 2020, 11(1): 3665.

[16] Mohtaram A, Shafizadeh-Moghadam H, Ketabchi H. A flexible multi-scale approach for downscaling grace-derived groundwater storage anomaly using lightgbm and random forest in the tashk-bakhtegan basin, Iran. Journal of Hydrology: Regional Studies, 2025, 57: 102086.

[17] Mu QZ, Zhao MS, Running SW. Improvements to a modis global terrestrial evapotranspiration algorithm. Remote Sensing of Environment, 2011, 115(8): 1781-1800.

[18] Peng H, Jia Y, Zhan C, et al. Topographic controls on ecosystem evapotranspiration and net primary productivity under climate warming in the Taihang mountains, China. Journal of Hydrology, 2020, 581: 124394.

[19] Raftery AE, Gneiting T, Balabdaoui F, et al. Using bayesian model averaging to calibrate forecast ensembles. American Meteorological Society, 2005, 133(5): 1155-1174.

[20] Ren W, Banger K, Tao B, et al. Global pattern and change of cropland soil organic carbon during 1901-2010: Roles of climate, atmospheric chemistry, land use and management. Geography and Sustainability, 2020, 1(1): 59-69.

[21] Shen H, Li J, Wu G, et al. Can cmip6 models accurately reproduce terrestrial evapotranspiration across China. International Journal of Climatology, 2025, 45(6): e8794.

[22] Siebert S, Burke J, Faures JM, et al. Groundwater use for irrigation – a global inventory. Hydrology and Earth System Sciences, 2010, 14(10): 1863-1880.

[23] Su T, Sun S, Wang S, et al. Spatiotemporal variation of actual evapotranspiration and its relationship with precipitation in northern china under global warming. Remote Sensing, 2022, 14(18): 4554.

[24] Sun SB, Chen BZ, Yan JW, et al. Potential impacts of land use and land cover change (LUCC) and climate change on evapotranspiration and gross primary productivity in the haihe river basin, China. Journal of Cleaner Production, 2024, 476: 143729.

[25] Yan LL, Cheng YB, Chen XH, et al. Verification of Type A Evaluation Methods of Uncertainty Based on Monte Carlo Method. Journal of Henan University of Science and Technology (Natural Science), 2018, 39(5): 17-20.

[26] Yang C, Lei HM. Climate and management impacts on crop growth and evapotranspiration in the north china plain based on long-term eddy covariance observation. Agricultural and Forest Meteorology, 2022, 325: 109147.

[27] Yang Y, Sun HW, Xue J, et al. Estimating evapotranspiration by coupling bayesian model averaging methods with machine learning algorithms. Environmental Monitoring and Assessment, 2021, 193(3): 156-156.

[28] Yuan WP, Liu SG, Yu GR, et al. Global estimates of evapotranspiration and gross primary production based on modis and global meteorology data. Remote Sensing of Environment, 2010, 114(7): 1416-1431.

[29] Zhang C, Liu JG, Shang JL, et al. Improving winter wheat biomass and evapotranspiration simulation by assimilating leaf area index from spectral information into a crop growth model. Agricultural Water Management, 2021, 255: 107057.

[30] Zheng M, Zha Y, Bian J, et al. Analysis of temporal and spatial changes in regional evapotranspiration in china and the factors affecting them, 2003–2020. Theoretical and Applied Climatology, 2025, 156(11): 601.

RIGHTS & PERMISSIONS

Journal of Groundwater Science and Engineering Editorial Office