Providing a comprehensive understanding of missing data imputation processes in evapotranspiration-related research: a systematic literature review

ABSTRACT

This study aimed to review the existing research focalizing on the missing data imputation techniques for the systems enabling actual evapotranspiration calculation (such as eddy covariance, Bowen ratio, and lysimeters) and divergent evapotranspiration related variables, i.e. temperature, wind speed, humidity, and solar radiation. Thus, the Scopus engine was utilized to scan the entire literature and 62 articles were diligently investigated. Results show classical approaches have been widely used by researchers due to their ease of implementation. However, the applicability and validity of these methods heavily rely on assumptions made about the distribution and characteristics of missing data. Hence, advanced imputation techniques produce more accurate outcomes as they handle complex and non-linear problems. Also, current trends embraced by the research community revealed that employing deep learning techniques and incorporating explainable artificial intelligence into imputations have significant potential to make insightful contributions to the body of knowledge.

KEYWORDS:

• missing data imputation
• gap filling
• eddy covariance
• evapotranspiration
• multiple imputation

Editor S. Archfield; Associate Editor (not assigned)

1 Introduction

The quality of measurements is of utmost importance for the majority of research areas, such as hydrology, meteorology, and earth science (Yozgatligil et al. 2013). Despite the notable advancements in technology over the past few decades, the acquisition of highly accurate measurements remains a significant challenge due to the potential for measurement errors. Such errors can occur due to various factors such as inadequate calibration, imprecise measurement readings, and environmental influences (e.g. proximity to man-made or natural structures, air pollution, etc.), which can negatively impact all stages of data analysis (Lin and Tsai 2020). To tackle these challenges, one can either remove missing data or employ imputation methods to assign substitute values (Langkamp et al. 2010). Here, the amount of missing data plays a critical role; missing data accounting for less than 10% or 15% of the entire dataset can be removed (Strike et al. 2001), whereas one should pay considerable attention to addressing missing values exceeding 15%.

The research community has increasingly focused on statistical models with varying degrees of complexity to minimize the efforts required to handle the missing data problem. The aim is to develop methods that require minimal assumptions while producing accurate estimations that closely align with the actual measured data. The pertinent literature recognized the critical role played by the classical methods, e.g. linear regression (LR) analysis, principal component analysis (PCA), and expectation-maximization (EM) (Hron et al. 2010). Despite their practical utilization, these algorithms require some prerequisites that pushed researchers to adopt different methods for accomplishing missing data imputation attempts. In this sense, modern algorithms (mostly concentrated on machine learning applications) are used for a plethora of purposes, such as imputation of river flow data (Ng et al. 2009, Sidibe et al. 2018), imputation of water quality data (Tabari and Talaee 2015, Zhang and Thorburn 2022), imputation of meteorological data (Sattari et al. 2020, Wang and Shi 2021), etc. The pertinent literature indicates that modern modelling techniques tend to produce more accurate estimations compared to classical imputation methods, thus offering promising solutions for addressing the challenges of missing data (Gautam and Ravi 2015, Silva-Ramírez et al. 2015, Zhang et al. 2015). Studies using modern techniques are explained in detail in the following sections.

The current paper focuses on the imputation of missing values observed in various systems to facilitate the computation of actual evapotranspiration (ET_a), along with variables aiding in the calculation of reference evapotranspiration (ET_o). Hence, this research encapsulates a comprehensive investigation of data imputation techniques for the parameters utilized in ET measurements and calculations. In this vein, studies carrying out imputation of actual ET values based on eddy covariance (EC), Bowen ratio, and lysimeter measurements as well as other efforts devoted to reconstructing temperature, humidity, wind speed, and solar radiation to assist in calculating the reference ET were diligently investigated. The present research sheds light on the selection of input parameters for imputation methods, the selection of hyperparameters in algorithms used during the imputation process, and the determination of iteration numbers for multiple imputation techniques. The software tools used in the imputation process were further investigated within the scope of the current study. The outcomes of this study are anticipated to not only enhance researchers’ understanding of the imputation techniques employed in ET-related research for the purpose of completing missing data, and their consequent impact on research outcomes, but also provide valuable insights for future investigations in this domain.

2 Mechanisms and patterns of missing data

The mechanisms and patterns also play a significant role in the resultant research outcomes, which may have a greater impact than the extent of missing data itself. Figure 1 visualizes the missing data mechanisms and patterns.

Figure 1. Mechanisms and patterns of missing data. MAR = Missing at Random; MCAR = Missing Completely at Random; MNAR = Missing Not at Random.

2.1 Missing data mechanisms

In imputing missing data, it is crucial to first identify the structural features of the missing data. In particular, it is necessary to designate whether the missing data in the observations are randomly distributed or follow a certain pattern. In this context, the classification system of missing data that was summarized by Rubin (1976), which treats missing data as a probabilistic process, has been widely used in recent decades. Here, among three different mechanisms, missing completely at random (MCAR) is regarded as the least challenging mechanism for statistical analyses as it grounded on the independent relationship between the missing data and other variables in the dataset. On the other hand, missing at random (MAR) assumes that there is a dependent relationship between the variable having missing values and other variables. In this mechanism, the probability of missing data relies on the observed values of other variables, rather than the other missing values of the target variable. In addition, missing not at random (MNAR) is based on a direct dependent relationship between missing data and the nature of the corresponding variable. Hence, datasets having relatively longer gaps are generally characterized by the MNAR mechanism.

2.2 Missing data patterns

According to Fig. 1, a potential indicator emerges regarding the possibility of data deletion within certain groups in the context of missing data mechanisms. When missing data is deemed negligible, this implies that the likelihood of observing a variable is not contingent upon the value of that variable, while the converse denotes non-negligible missing data (Beaulieu-Jones et al. 2018). On the one hand, in instances where the missing data mechanism is MCAR or MAR, data deletion may result in minimal information loss. On the other hand, in the case of an MNAR mechanism, it is recommended to employ data imputation as the missingness is attributable to the impact of other independent variables.

Furthermore, the utilization of missing data imputation techniques is based on the patterns of missing data as the corresponding patterns indicate the location of missing values in the focalized time series. The existing literature covers many missing data patterns, mainly divided into two groups, i.e. monotone or non-monotone (general) and univariate or multivariate. Let the data be represented by the $n\times p$ matrix. In the presence of missing data, $Y$ is partially observed. Notation $Y_j$ is the $j^{th}$ column in $Y$, and $Y_{-j}$ indicates the complement of $Y_j$ – that is, all columns in $Y$ except $Y_j$. The missing data pattern of $Y$ is the $n\times p$ binary response matrix $R.$ A missing data pattern can be regarded as univariate if there is only one variable with missing data, whereas, once a set of series are all observed or missing on the same set of cases, the corresponding pattern is called multivariate. In addition, a missing data pattern is regarded as monotone if the variables $Y_j$ can be ordered such that if $Y_j$ is missing then all variables $Y_k$ with $k>j$ are also missing. However, if the pattern is not monotone, it is called non-monotone or general.

3 Methods for measuring and calculating ET

3.1 The eddy covariance method

The EC method is widely used to measure fluxes in ecosystem studies across the globe. Flux refers to the amount of matter passing through a closed surface per unit area in a unit of time. The EC method describes the airflow in the atmosphere as consisting of vertical movements that transport materials through turbulence. The difference between the number of gas molecules transported upwards and those transported downwards at any given moment results in net convection. In the EC method, the concentration of the gas being measured and the covariance of the vertical wind speed are the fundamental principles used for estimating gas exchange between the atmosphere and the surface (Burba et al. 2013). Equation (1)(1) $LE={\rho }_a\overline{w^{\prime }{q}^{\prime }}$(1) shows how the evaporation latent heat flux (LE) is calculated:

(1) $LE={\rho }_a\overline{w^{\prime }{q}^{\prime }}$(1)

Here, LE represents the latent heat of evaporation (W m⁻²), ρ_a is the air density (kg m⁻³), w′ represents the deviation of the vertical wind speed from the mean (m s⁻¹), q′ represents the deviation of the water vapor concentration from the mean (kg m⁻³), and $\overline{w^{\prime }{q}^{\prime }}$ represents the covariance of these deviations (Jiang et al. 2014).

3.2 Bowen ratio energy balance method

The Bowen ratio energy balance (BREB) method is used to determine the actual ET based on the ratio of distribution of the convective flux between the evaporative latent heat flux and the sensible heat flux (Oke 2002). The BREB method involves measuring the temperature and humidity gradient, as well as net radiation and soil heat flux. The BREB method is based on the energy balance equation (Bowen 1926):

(2) $R_n-G=LE+H$(2)

Here, R_n represents net radiation (W m⁻²), G represents the soil heat flux (W m⁻²), LE represents evaporative latent heat flux (W m⁻²), and H represents sensible heat flux (W m⁻²). The method is based on the Bowen ratio calculation introduced by Bowen, which shows the ratio of distribution of the receiving energy between LE and H (Bowen 1926).

(3) $\beta =\gamma \left[\frac{T_1-T_2}{e_{a_1-}{e}_{a_2}}\right]=\frac{H}{LE}$(3)

where β represents the Bowen ratio, γ indicates the psychrometric coefficient (kPa °C⁻¹), T₁ represents the air temperature (°C) measured at the lower level, T₂ represents the air temperature (°C) measured at the upper level, e_a1 represents the actual vapor pressure (kPa) of the air at the lower level, and e_a2 represents the actual vapor pressure (kPa) of the air at the upper level. With the determination of the Bowen ratio, LE (W m⁻²) can be calculated using the expression given in Equation (4)(4) $LE=\frac{R_n-G}{1+\beta }$(4) .

(4) $LE=\frac{R_n-G}{1+\beta }$(4)

3.3 Weighing lysimeter method

The water budget equation was used for ET measurements in lysimeters. For a typical lysimeter tank insulated from the environment, the water budget equation can be written as in Equation (5)(5) $\left(P+I+C_p+R_{fc}\right)-\left(ET+D_p+R_{fo}\right)=\mp \mathrm{\Delta }\mathrm{W}$(5) .

(5) $\left(P+I+C_p+R_{fc}\right)-\left(ET+D_p+R_{fo}\right)=\mp \mathrm{\Delta }\mathrm{W}$(5)

where P is precipitation and I is irrigation water. Also, C_p and R_fc denote capillary rise and runoff, respectively. In Equation (5), ET, D_p, and R_fo indicate evapotranspiration, deep infiltration or drainage water, and runoff, respectively. ΔW shows the changes in soil water content. It is worth mentioning that all units are in millimetres. In general, C_p,R_fc, and R_fo values are zero in lysimeters. Thus, the equation for ET ends up in the form of Equation (6)(6) $ET=P+I-D_p\mp \mathrm{\Delta }\mathrm{W}$(6) :

(6) $ET=P+I-D_p\mp \mathrm{\Delta }\mathrm{W}$(6)

3.4 FAO-56 Penman-Monteith method

The FAO Penman-Monteith method is an equation derived from the original Penman-Monteith equation along with the aerodynamic and surface resistance equations (Allen et al. 1998).

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

where ET_o is the reference evapotranspiration [mm day⁻¹], R_n is the net radiation at the crop surface [MJ m⁻² day⁻¹], G is the soil heat flux density [MJ m⁻² day⁻¹], T is the mean daily air temperature at 2 m height [°C], U₂ is the wind speed at 2 m height [m s⁻¹], e_s is the saturation vapor pressure [kPa], e_a is the observed vapor pressure [kPa], e_s – e_a refers to the saturation vapor pressure deficit [kPa], Δ represents the slope of vapor pressure curve [kPa °C⁻¹], and γ denotes the psychrometric constant [kPa °C⁻¹].

4 Literature review rationale

In this study, missing data imputation studies of the parameters used to obtain actual and reference ET values were examined. The Scopus search engine was used to conduct the comprehensive literature survey. To ensure relevance, the search was limited to articles published in English between 2002 and 2022, reflecting the most recent research trends in the field. A total of 282 articles were subjected to the abstract-level investigation, after which the number of studies was reduced to 110. Finally, the articles deemed directly relevant to the missing data imputation of the variables taken into account were further analysed through full-text investigation, which resulted in a total of 62 articles. The entire process of the literature review is illustrated in Fig. 2.

Figure 2. The steps of the literature review conducted via Scopus.

5 Results

Three methods, i.e. EC, Bowen ratio, and lysimeter systems, that are used to measure actual ET values were examined in the present research. For reference ET calculation, the parameters (i.e. temperature, wind speed, pressure, insolation, radiation, and humidity) of the FAO-56 Penman-Monteith method, which is one of the most commonly adopted empirical approaches in the literature, were considered. It is also worth noting that the review of the pertinent literature revealed that there is no study on the imputation of missing data during the measurement of pressure and insolation times.

Along with the content-based investigation, the present work also aimed to highlight the scientometrics of the relevant studies. Therefore, a total of 62 articles extracted through the Scopus engine were subjected to the VOSviewer software and a keyword-based analysis was performed to explore the knowledge distribution over years. In Fig. 3, elements’ size is directly related to the occurrence frequency of the keywords, while the colour represents distinct time periods associated with the keywords, and the links between the nodes symbolize the co-occurrence relationship between the words. The corresponding analysis revealed that the early studies (represented by deep blue) mainly concentrated on the ET calculations were linked with climatic data, data obtained from EC and Bowen-ratio systems. However, most studies primarily focused on filling gaps to perform further processes rather than directly embracing the missing data imputation objective. In addition, recent attempts tend to discover the contribution of air temperature, land surface temperature (LST), and wind speed, while the research community also focalized on the integration of classical techniques in gap filling. Still, few studies have acknowledged the role of advanced modelling techniques such as machine learning approaches.

Figure 3. Keyword networks generated by VOSviewer based on the Scopus database.

5.1 Missing data imputation techniques

The existing literature recognized the critical role played by divergent approaches (e.g. statistical, data-based, machine learning techniques, etc.) in missing data imputation to acquire precise estimations. Hence, in the present research we diligently reviewed the pertinent literature and extracted the commonly utilized methodologies to provide a better overview regarding the definitions, pros, and cons of the corresponding techniques (Table 1).

Table 1. Commonly applied missing data imputation methods.

Method	Basis	Rationale	Assumption	Complexity	HP	CC	Explanation
Deletion	N/A	N/A	N/A	Low	N/A	Low	It removes observations containing missing values from the dataset. Advantage: Easy to implement. Disadvantage: Information loss.
Simple imputation	Statistical	N/A	N/A	Low	N/A	Low	It replaces missing values with the estimations derived from the quantitative attribute of the available values. Advantage: Easy to implement. Disadvantage: Prone to change the original distribution of the variable.
Mean imputation	Statistical	N/A	N/A	Low	N/A	Low	It replaces the missing value with the arithmetic mean of the remaining dataset. Advantage: Easy to implement. Disadvantage: Potentially distorts imputations due to unequal data distribution and the influence of extreme values on the mean.
Regression imputation	Statistical	N/A	Normal distribution	Low	N/A	Low	It estimates the relationship between input and output variables, or between a data point and its corresponding variable. Advantage: Easy to implement. Disadvantage: Requires many assumptions due to its parametric nature.
Hot-deck imputation	Data based	N/A	N/A	Mid	N/A	High	It replaces missing data with values from complete data cases. Advantage: Being robust to misidentification and interaction effects as the imputations are generated using an implicit model. Disadvantage: Underestimates standard errors due to reduced variability that results from the selection of completed data from already observed data. In addition, when the dataset is not MCAR, data imputations can result in unpredictable bias.
Expectation–maximization	Statistical	N/A	Distribution	Mid	N/A	High	It is an iterative approach utilized to compute the maximum likelihood estimation of variables within statistical models. Advantage: Better than mean imputations as it preserves the relationship with other variables. Disadvantage: It does not provide standard error estimations.
Multiple imputation	Statistical	N/A	N/A	Mid	N/A	High	The regression assignment is re-operated more than once and each dataset is stored for later use, while the obtained errors are calculated. Then, the real dataset is attained by taking the average of the corresponding data. Advantage: It produces statistically valid estimates under the assumption of MAR and yields asymptotically consistent, efficient, and normally distributed estimates. Disadvantage: Due to its inherent variability, it often results in different estimates each time, leading to different outcomes among researchers working with the same dataset.
Linear interpolation	Statistical	Distance-based	Linearity	Mid	N/A	Low	It performs interpolation between a data point immediately preceding the gap and another data point shortly following the gap. Advantage: When properly utilized, it provides improved estimates for statistical inference. Disadvantage: It provides limited accuracy where continuous long missing data exists.
Kriging	Statistical	Distance-based	Stationarity and isotropy	Mid	N/A	Mid	It is a spatial estimation method that utilizes nearby data points to predict the values of unsampled points and uses semi-variogram structural features to make unbiased estimations of positional changes at unsampled locations. Advantage: Measuring the uncertainty of predictions and ensuring uncertainties to the user, enabling more reliable and informed decision making processes. Disadvantage: Explaining spatial relationships using a linear model, which may result in an inability to capture complex or non-linear relationships present in the dataset.
Inverse distance weighting	Statistical	Distance-based	N/A	Mid	N/A	Low	It is grounded on the estimation of the values of unsampled points by utilizing the values of sampled points. The estimated values are a function of both distance and the size of neighbouring points, and the importance and impact on the cell being estimated decrease as distance increases. Advantage: When used properly, it enhances the accuracy of statistical inference. Disadvantage: If the data has low density or is unevenly distributed, it can lead to inaccurate estimation of missing values in certain regions.
K-nearest neighbours	Machine learning	Distance-based	N/A	Low	Yes	Low	It works by identifying the closest k (referring to the number of neighbours) number of data points to the data point being evaluated, and assigning it to the same class or estimating its value based on the values of those k data points. Advantage: Easy to implement. Disadvantage: It does not work well with large datasets as calculating the distance between data instances is costly. Similarly, it is not compatible with high-dimensional data as the calculation process becomes more complex for each dimension. Additionally, the KNN algorithm is sensitive to noisy and missing data.
Decision trees	Machine learning	Tree-based	N/A	Mid	Yes	High	It uses a sequential approach to perform the classification or regression process, aiming to minimize the standard error to obtain an optimal decision tree. Advantage: Easy to understand, interpret, and visualize and does not require normalization or scaling of data. Disadvantage: It is prone to overfitting.
Artificial neural networks	Machine learning	Black box	N/A	High	Yes	High	Inspired by biological neural networks, ANNs are a set of mathematical operations aiming to solve inherent non-linear problems. Advantage: Non-linear and flexible, able to work with incomplete information. Disadvantage: High computational costs, requiring huge amounts of data, and limited explanations regarding its internal processes.
Adaptive neuro-fuzzy inference system	Machine learning	Black box	N/A	High	Yes	Very high	It defines the relationship between inputs and outputs in fuzzy space, while neural networks try to find the most suitable parameters of the model. Advantage: Less complex compared to other decision making systems. Disadvantage: Poses significant challenges in estimations having large number of inputs, including the curse of dimensionality, interpretability of rules, and parameter training.
Support vector machine	Machine learning	Distance-based and rule-based	N/A	High	Yes	High	It performs optimal classification with a single differential decision surface and is based on the principle of structural risk minimization. Advantage: It is more effective when there is a clear margin of separation between classes. Disadvantage: It has limited suitability for high-dimensional spaces and requires choosing an optimal kernel, which is a challenging task.

HP: hyperparameter; CC: computational cost.

5.2 Missing data imputation attempts concentrated on Bowen ratio data

The pertinent literature illustrated that studies investigated to impute the missing values in Bowen ratio measurements account for 5% of the total number of studies. The corresponding attempts mostly concentrated on the utilization of classical approaches, such as regression and linear interpolation techniques. For instance, Guo et al. (2007) employed two distinct imputation methods for completing missing LE and H values. Initially, the Priestly-Taylor method was utilized to estimate the missing LE values, followed by the completion of missing H values using the energy balance equation. This sequential imputation approach was effective in simulating both daytime evaporation and night-time condensation. The imputed fluxes closely matched the quality-controlled data for each flow scenario in terms of their temporal distribution. Similarly, in other studies using the Bowen ratio method, linear interpolation was the most widely adopted method, as can be seen from Table 2 (Balogun et al. 2009, Sexstone et al. 2016, Burakowski et al. 2018).

Table 2. Details regarding the missing value imputation studies utilizing the Bowen-Ratio technique.

ID	Reference	Imputation methods	Method/data	Gap size	Performance metric
1.	Burakowski et al. (2018)	LI	Bowen Ratio/RSWin	<0.01%	N/A
2.	Sexstone et al. (2016)	LI	Bowen Ratio/Rn	N/A	N/A
3.	Balogun et al. (2009)	LI	Bowen Ratio/Rn	15%	N/A
4.	Guo et al. (2007)	PT, LR	Bowen Ratio/LE	1-7%	R^2

LI: Linear Interpolation; PT: Priestley-Taylor; LR: Linear Regression; RSWin: Incoming Shortwave Radiation; LE: Latent Heat Flux; R2: Coefficient of Determination.

5.3 Missing data imputation attempts concentrated on EC data

EC measurements are significant for determining ET as they provide direct and continuous measurements of the turbulent fluxes of water vapor, heat, and momentum between the Earth’s surface and the atmosphere. Table 3 presents the studies devoted to attempts regarding the imputation of missing values encountered in EC measurements. For instance, Tang et al. (2019) performed the classical LR method to impute missing data in LE and H values that were measured using the EC method over four years in a tropical peat forest in Sarawak, Malaysia. Concerning the software that can be directly used in missing data imputation attempts, Deb Burman et al. (2019) implemented the REddyProc software package compiled with the R program to impute missing data in H, LE, Rn, short wave incoming solar radiation (R_SWin), G, and photosynthetic photon flux densities (PPFD) provided by the EC method in a tropical forest located in Northeast India. Similar attempts further illustrated the effectiveness of the REddyProc package regarding the imputation of missing data in LE values (Zitouna-Chebbi et al. 2018, Foltýnová et al. 2020).

Table 3. Details regarding the missing value imputation studies utilizing the eddy-covariance technique.

ID	Reference	Imputation methods	Method/data	Gap size	Performance metric
5.	Goodrich et al. (2021)	PLS, NN, KNN	Eddy covariance/N2O	15%	R^2
6.	Tang et al. (2019)	LR	Eddy covariance/H, LE	N/A	N/A
7.	Deb Burman et al. (2019)	REddyProc package	Eddy covariance/H, LE, Rn, RSWin, G, PPFD	N/A	N/A
8.	Abudu et al. (2010)	ANN	Eddy covariance/LE	5, 10, 20, 30, 40%	R^2, RMSE, MAE, MRE
9.	Chen et al. (2012)	MLR, KNN	Eddy covariance/LE	5, 10, 20, 40%	RMSE
10.	Hui et al. (2004)	MI	Eddy covariance/NE, LE, H	6–51%	R^2
11.	Alavi et al. (2006)	KF	Eddy covariance/LE	4–17%	R^2, RMSE, Rbias, R
12.	Kunwor et al. (2017)	Traditional, static window method, moving window method, parameter prediction method	Eddy covariance/NEE	35%	RMSE
13.	Schmidt et al. (2008)	NN	Eddy covariance/ CO2	13%	RMSE
14.	Mahabbati et al. (2021)	ANN, LR, RF, SVM, XGB, FFNN	Eddy covariance/NE, LE, H	1–32%	R^2, RMSE, MBE, VR, AI
15.	Li et al. (2007)	LI	Eddy covariance/H, LE	N/A	N/A
16.	Eamus et al. (2013)	Self-organizing linear output, ANN	Eddy covariance/E, C	N/A	RMSE
17.	Foltýnová et al. (2020)	REddyProc package	Eddy covariance/LE	1–100 half-hour	RMSE, MBE
18.	Boudhina et al. (2018)	REddyProc package, LR, MLR, evaporative fraction	Eddy covariance/LE	N/A	R^2, RMSE
19.	Zitouna-Chebbi et al. (2018)	REddyProc package	Eddy covariance/LE	10–70%	B, RB, RRMSE, R^2

PLS: Partial Least Square; NN: Neural Network; KNN: K-Nearest Neighbors; H: Sensible Heat Flux; G: Soil Heat Flux; PPFD: Photosynthetic Photon Flux Densities; ANN: Artificial Neural Networks; RMSE: Root Mean Square Error; MAE: Mean Absolute Error; MRE: Mean Relative Error; MLR: Multiple Linear Regression; MI: Mean Imputation; KF: Kalman Filter; Rbias: Relative Bias; RF: Random Forest; SVM: Support Vector Machine; XGB: Extra Gradient Boosting; FFNN: Feed Forward Neural Network; MBE: Mean Bias Error; VR: Ratio of Variance; AI: Index of Agreement; E: Expectation Stage; Rrmse: Relative Root Mean Square Error; NEE: Net Ecosystem Exchange.

In addition, there are many studies in which missing LE data were imputed with machine learning algorithms, such as artificial neural networks (ANN) (Schmidt et al. 2008, Abudu et al. 2010, Eamus et al. 2013, Mahabbati et al. 2021), k-nearest neighbours (KNN) (Chen et al. 2012) and Kalman filter (KF) (Alavi et al. 2006). Kunwor et al. (2017) used three non-linear regression approaches to fill artificially created gaps of different sizes via light and temperature response curves: (1) “traditional” fixed monthly window, (2) moving window and (3) moving window with parameter estimation using physiological drivers. The results of the incomplete data-imputation simulations showed that the net ecosystem exchange (NEE) estimates made with the moving window and parameter estimation methods were closer to those of observed NEE, while the traditional method had yielded lower overall accuracy.

5.4 Missing data imputation attempts concentrated on lysimeter data

Lysimeter measurements are a valuable tool for understanding the water balance of a specific area and can provide important insights into the complex processes that drive ET. They ensure direct measurement of the water balance of a certain region, and can also be used to validate and improve remote sensing-based estimates of ET. Despite their importance, limited studies have been conducted due to the scarcity of the corresponding systems across the globe. Still, Huang et al. (2020) (which is the 62^nd study investigated within the scope of the present research) artificially generated random gaps for the three time series attained by the lysimeter systems from Basel (Switzerland), Rheindahlen (Germany), and Rietholzbach (Switzerland) to explore the utility of four techniques: potential ET, ratio, FAO-based water balance and HYDRUS modelling methods. Interestingly, the results suggested that the ratio method was well suited for imputing missing values of lysimeter data for Basel (Switzerland), while the HYDRUS method outperformed its counterparts concerning Rheindahlen (Germany). For Rietholzbach (Switzerland), all methods yielded very similar results but greater root mean square error (RMSE) values were observed with the FAO method.

5.5 Missing data imputation attempts concentrated on humidity data

Measuring the humidity for a region plays a crucial role in calculating the vapor pressure deficit (VPD), which is a key factor in the rate of ET and is important in understanding the overall water balance of a region. For instance, Costa et al. (2021) used the multiple imputations by chained equations (MICE) algorithm to impute missing humidity data. Despite the usage of classical approaches, related studies mostly rely on the application of machine learning techniques, such that along with the adoption of Kalman smoothing and autoregressive integrated moving average (ARIMA) in their study, Afrifa-Yamoah et al. (2020) further benchmarked them with multiple linear regression (MLR) methods for missing humidity data estimation. In a similar vein, Hasanpour Kashani and Dinpashoh (2012) used ANN, adaptive neuro fuzzy inference system (ANFIS), and genetic programming (GP) methods, while Körner et al. (2018) and Lara-Estrada et al. (2018) discovered the performance of gradient boosting and Bayesian network methods, respectively (Table 4).

Table 4. Details regarding the studies devoted to imputing missing humidity data.

ID	Reference	Imputation methods	Method/data	Gap size	Performance metric
20.	Costa et al. (2021)	MICE	Humidity	17–54%	RMSE
21.	Hasanpour Kashani and Dinpashoh (2012)	AA, ID, NR, single best estimator, MLR, UK traditional method, CSM, ANN, ANFIS, GP,	Humidity	N/A	MAE, NSE, SS
22.	Afrifa-Yamoah et al. (2020)	Kalman smoothing, ARIMA, MLR	Humidity	10%	MAE, RMSE, SMAPE
23.	Körner et al. (2018)	GTB	Humidity	50%	RMSE, MAE, R^2
24.	Anjomshoaa and Salmanzadeh (2019)	LI, CI	Humidity	1–374 h	VR, MAE, RMSE
25.	Lara-Estrada et al. (2018)	Bayesian network	Humidity	20%	RMSE, Bias
26.	Lee et al. (2010)	LI	Humidity	0.1%	N/A

MICE: Multiple Imputations by Chained Equations; NSE: Nash–Sutcliffe Model Efficiency Coefficient; SS: Skill Score; AA: Simple Arithmetic Averaging; ID: Inverse Distance Interpolation; NR: Normal Ratio; CSM: Closest Station Method; ANFIS: Adaptive Neuro-Fuzzy Inference System; GP: Genetic Programming; ARIMA: Autoregressive Integrated Moving Average; SMAPE: Symmetric Mean Absolute Percentage Error; GTB: Gradient Tree Boosting; CI = Cubic Splines Interpolation.

5.6 Missing data imputation attempts concentrated on solar radiation data

Measurements pertaining to the solar radiation variable help researchers understand the ET mechanisms as they provide insights into the energy available for photosynthesis, transpiration, and soil evaporation, and can be used to estimate the net radiation balance. Among the missing data imputation studies in solar radiation values (Table 5), Zainudin et al. (2015) assessed the utility of one of the commonly adopted machine learning algorithms, KNN, by comparing its performance with three classical techniques (i.e. linear interpolation, spline interpolation, and cubic interpolation). Likewise, Kim et al. (2019) employed not only different conventional techniques (including linear interpolation (LI), mode imputation, and MICE) but also the KNN to complete missing solar radiation data. The authors reported that LI produced optimal results when the proportion of missing data was below 15%. However, as the proportion of missing data increased, the KNN method, producing lower error rates, outperformed its counterparts. Still, relevant research has generally focused on classical imputation methods, such as the LR method (Dimas et al. 2011, Turrado et al. 2014, Žliobaitė et al. 2014). Demirhan and Renwick (2018) created a comprehensive comparison scheme encompassing 36 models with different variants of interpolation, Kalman filters, persistence, weighted moving average (WMA), and random sample techniques. From a different perspective, Ernst and Gooday (2019) utilized an equation consisting of direct normal irradiance (DNI) and diffuse horizontal irradiance (DHI) values instead of interpolation methods to complete missing solar radiation data.

Table 5. Details regarding the studies devoted to imputing missing solar radiation data.

ID	Reference	Imputation methods	Method/data	Gap size	Performance metric
27.	Zainudin et al. (2015)	LI, SI, CI, KNN	Solar radiation	10, 20, 30, 40, 50%	RMSE, R^2
28.	Ogunsola and Song (2014)	SSA, SASR, TBA	Solar radiation	5, 10, 15, 20, 25 days	ME, MAE, RMSE
29.	Žliobaitė et al. (2014)	LR	Solar radiation	1–25%	MSE
30.	Kim et al. (2019)	LI, MI, MICE, KNN	Solar radiation	10, 15, 20%	RMSE, MRE, RMSD, MRD
31.	Turrado et al. (2014)	MICE, IDW, MLR	Solar radiation	0.002–1%	RMSE, MAE
32.	Demirhan and Renwick (2018)	36 imputation methods	Solar radiation	5, 10, 25, 50%	rMAE, rRMSE, MASE
33.	Schwandt et al. (2014)	Basic gap filling, satellite-based gap filling	Solar radiation	0.001–1.6%	NA
34.	Ernst and Gooday (2019)	DNI and DHI	Solar radiation	17%	NA

SI: Simple Imputation; SSA: Singular Spectrum Analysis; SASR: Statistically Adjusted Solar Radiation; TBA: Temperature-Based Approach; MSE: Mean Square Error; IDW: Inverse Distance Weighting; DNI: Direct Normal Irradiance; DHI: Horizontal Irradiance; MASE: Mean Absolute Standard Error.

5.7 Missing data imputation attempts concentrated on wind speed data

Wind speeds impact the rate of air movement over the surface, which in turn affects the rate of water vapour transfer from the surface to the atmosphere. As shown in Table 6, the use of machine learning methods predominates in studies related to completing wind speed missing data. For instance, the ANFIS approach was adopted by Hocaoglu et al. (2009) and Zhiling Yang et al. (2011) for Turkey and China, respectively. To illustrate the effectiveness of classical and modern techniques, Liu et al. (2018) combined the multiple imputation (MPI) and Gaussian process regression (GPR) approaches to perform predictions regarding various datasets containing missing values. In this regard, a plethora of research further validated the superiority of the machine learning rationale for similar purposes, such as back propagation neural network (BPNN) (Basheer Shukur and Hisyam Lee 2015, Mao and Jian 2016), generalized regression neural network (GRNN) (Du et al. 2019), and recurrent neural network (RNN) (Li et al. 2019), in the pertinent literature. Similar to the NN-based methods, other machine learning techniques have been proven to be practical tools; for example, Chen et al. (2016) integrated the least-squares support vector machine (LSSVM) model with comprehensive learning particle swarm optimization (CLPSO) to reveal that the data tested with the optimized LSSVM model has better imputation sensitivity and the interpolation results of the missing wind speed data are suitable for the variation of the wind speed data at the adjacent level. de Oliveira et al. (2021) developed machine learning models to predict tropospheric ozone concentrations and completed missing values in the input parameters using the MissForest algorithm. Wind speed, which is one of the parameters in the input set, had varying degrees of missing data in different seasons, with the smallest and largest missing rates being 1.69% and 43.75%, respectively. Prior to imputing missing data, the statistical characterization of the dataset was ensured with PCA, and following the imputation process, further controls showed no significant change in the dataset characteristics, leading the researchers to conclude that the imputation process was successful.

Table 6. Details regarding the studies devoted to imputing missing wind speed data.

ID	Reference	Imputation methods	Method/data	Gap size	Performance metric
35.	Yang et al. (2011)	ANFIS	Wind speed	N/A	RMSE, R^2, MAE
36.	Liu et al. (2018)	GPR, MI	Wind speed	5, 10, 15, 20, 25, 30%	RMSE, MAE
37.	Mao and Jian (2016)	NN	Wind power	N/A	N/A
38.	Du et al. (2019)	NN	Wind speed	N/A	MSE, MAPE, MAE
39.	Li et al. (2019)	RNN	Wind speed	N/A	N/A
40.	He et al. (2013)	Mahalanobis distance	Wind speed	N/A	N/A
41.	Hocaoglu et al. (2009)	ANFIS	Wind speed	N/A	RMSE, R^2
42.	Xie and Sun (2018)	Low-rank matrix approximation, KNN, MPI, MI	Wind speed	N/A	RMSE
43.	Chen et al. (2016)	LSSVM	Wind speed	10%	MAPE, RMSE,
44.	Wesonga (2015)	MICE	Wind speed	28%	N/A
45.	Lin and Wang (2014)	CEM	Wind speed	1–12 h	MBE, MAPE, RMSE, r
46.	de Oliveira et al. (2021)	MissForest	Wind speed	1.69–43.75%	Mean_RMSE

GPR: Gaussian Process Regression; RNN: Recurrent Neural Network; MPI: Multiple Imputation; LSSVM: Least-Squares Support Vector Machine; CEM: Correlated Echelon Model; MAPE: Mean Absolute Percentage Error; r: Correlation Coefficient

5.8 Missing data imputation attempts concentrated on temperature data

Table 7 demonstrates the attempts focused on the imputation of missing temperature data in the relevant literature. For instance, Xu et al. (2013) introduced a neoteric approach, called Biased Sentinel Hospitals-based Area Disease Estimation (P-BSHADE) point estimation model, and benchmarked it with the kriging, inverse distance weighting (IDW), and spatial regression testing (SRT). The effectiveness of the technique was empirically evaluated in terms of the annual temperature dataset from 1950 to 2000 in China and the scholars concluded that the P-BSHADE is applicable to the annual Chinese temperature datasets.

Table 7. Details regarding the studies devoted to imputing missing temperature data.

#	Reference	Imputation methods	Method/data	Gap size	Performance metric
47.	Xu et al. (2013)	P-BSHADE, IDW, SRT	Temperature	N/A	RMSE, MAE
48.	Coulibaly and Evora (2007)	MLP, RBF, RNN, TLFN, TDRNN, CFNN	Temperature	2.1–65.2%	MAE, r
49.	Park et al. (2019)	LSTM	Temperature	6, 12, 24 h	RMSE, MBE
50.	López et al. (2021)	NN, KNN	Temperature	N/A	R^2
51.	Thevakaran and Sonnadara (2018)	Between-station technique	Temperature	0.55–1.22%	N/A
52.	Altan and Ustundag (2012)	Wavelet	Temperature	50%	MSE, MAE
53.	Katipoğlu (2022)	SVM, ANFIS, DT	Temperature	20%	RMSE, MAE, R^2
54.	Shtiliyanova et al. (2017)	Kriging	Temperature	16, 17.69%	NSE, MAPE, MAE, BIAS
55.	Cheval et al. (2020)	DINEOF	Land surface temperature	1.1–58.4%	R^2
56.	Kadow et al. (2020)	DL	Temperature	12.5%	P, RMSE
57.	Dumitrescu et al. (2020)	MLR, GAM	Land surface temperature	N/A	MAE, RMSE, r
58.	Zhou et al. (2020)	MLR	GIS/land surface temperature	N/A	RMSE
59.	Williams et al. (2018)	EMD, Bayesian compressive sensing	GIS/temperature	1, 5, 10, 20%	RMSE
60.	Lompar et al. (2019)	Debiasing procedure	GIS/temperature	5 days	RMSE
61.	Koksal et al. (2017)	ANN	Air temperature, relative humidity	30%	R^2, RMSE, MAE, MAPE

P-BSHADE: Biased Sentinel Hospitals-Based Area Disease Estimation; SRT: Spatial Regression Testing; MLP: Multi-Layer Perceptron; RBF: Radial Basis Function; TLFN: Time-Lagged Feedforward Network; TDRNN: Time Delay Recurrent Neural Network; CFNN: Counter Propagation Fuzzy-Neural Network; LSTM: Long Short-Term Memory; DT: Decision Trees; DINEOF: Data Interpolating Empirical Orthogonal Functions; DL: Deep Learning; GAM: Generalized Additive Model; EMD: Empirical Mode Decomposition; P: Precipitation.

Among other attempts devoted to imputing missing temperature values, several studies have been carried out using machine learning approaches for different regions across the globe, such as multi-layer perceptron (MLP), radial basis function (RBF), RNN, time-lagged feedforward network (TLFN), time delay recurrent neural network (TDRNN), and counter propagation fuzzy-neural network (CFNN) methods for Canada (Coulibaly and Evora 2007), the LSTM method in South Korea (Park et al. 2019), and NN and KNN methods in Chile. In a study using the temperature data of four neighbouring stations with complete temperature data records pertaining to the same climatological area, Thevakaran and Sonnadara (2018) estimated the missing temperature data at the target station and performed accurate predictions due to the inter-stational correlation. Regarding Europe, Shtiliyanova et al. (2017) applied a kriging-based interpolation in the spatial dimension to fill the data gaps in air temperatures and a data-driven machine learning-based approach for the temporal dimension.

In addition, Geographic Information System (GIS)-related data plays a pivotal role in estimating ET by providing information on the spatial heterogeneity of environmental factors, including topography, land cover, soil type, and meteorological variables. Moreover, GIS data facilitate the assessment of catchment surface area and characteristics that are fundamental in ET calculations. The GIS data can also be utilized to evaluate energy balance components of the land surface, which are necessary for determining the energy available for ET. Hence, several efforts have been devoted to acquiring complete records regarding GIS-related datasets (Table 7). For instance, to accomplish the imputation of raw images containing between 1.1% and 58.4% incomplete LST data, Cheval et al. (2020) utilized the data interpolating empirical orthogonal functions (DINEOF) algorithm via the R software. A void-filled LST dataset with high spatial resolution was used to investigate LST climatology on a monthly, seasonal, and annual scale over an urban area in Bucharest, Romania. In addition, to reduce uncertainties in climate records, Kadow et al. (2020) proposed a new deep learning approach, while Dumitrescu et al. (2020) used two class regression models (i.e. MLR and generalized additive model [GAM]) characterized by differing complexity and flexibility to impute the missing LST data.

The study conducted by Koksal et al. (2017) focused on estimating incomplete hourly climate data in the semi-humid Bafra plain of northern Turkey, which plays a significant role in terms of agricultural facilities. They utilized the ANN technique and found that it not only accurately mapped the missing ET data but also estimated the data for areas where the ET was not recorded. Along with the adoption of machine learning algorithms, signal processing approaches have also been employed in the imputation of GIS-related data. For instance, Williams et al. (2018) incorporated two signal processing methods into the missing temperature values attained from remote sensing products. They only used time series with missing data and disregarded the information of nearby stations. The results showed that the BCS-on-Signal method outperformed the bayesian compressive sensing (BCS)-on-intrinsic mode functions (IMF) method, highlighting the effectiveness of incorporating signal processing methods to estimate missing temperature data from remote sensing products.

6 Discussion

Missing data is a common problem in data analysis, and a plethora of algorithms have been introduced to overcome the corresponding challenge. In addition, the reliability and consistency of the imputation attempts should be diligently examined in various aspects. Hence, this research assessed the 62 existing studies based on three major phases, called pre-imputation (in terms of assumptions, mechanism and pattern control, and standardization), during-imputation (in terms of software used and the number of iterations), and post-imputation (model run times, statistical distribution checks, time series variance, and validation), that must be taken into account to pinpoint the crucial issues in implementing data imputation strategies and ensuring the accuracy and validity of the imputed datasets. The following subsections briefly discuss the corresponding elements.

6.1 Assessment of the pre-imputation stage

The pre-imputation stage is a crucial step in missing data imputation, as the foundation for the entire imputation process is set in this stage. Hence, scholars are advised to evaluate the corresponding step concerning different aspects including the pattern and extent of missing data, and assess the causes of missing data, statistical assumptions and data standardization. Failing to adequately address these issues during the pre-imputation stage can lead to biased and inaccurate imputation results, which can have serious consequences for subsequent analysis and the interpretation of the adopted models. This research assessed the pertinent literature based on whether the existing works concerned the statistical assumptions, checked mechanisms and patterns, and performed the standardization of datasets.

Concerning the statistical assumptions, it was found that the vast majority of studies utilizing the classical linear regression-based imputation techniques lack a proper reporting of statistical assumptions. Specifically, Kunwor et al. (2017) rigorously tested the normality of the data and homogeneity of the variances in imputing their missing EC measurements. Furthermore, Afrifa-Yamoah et al. (2020) suggested the use of non-parametric regression models for the imputation of humidity data. This aspect was further emphasized by Rezvan et al. (2015), who reported that such imputations lack reliability and the direct use of data following imputation may lead to spurious results. Therefore, it is imperative that researchers acknowledge the importance of adhering to appropriate statistical assumptions prior to the missing data imputation in order to ensure the validity and accuracy of subsequent analyses.

The present study also examined the existing studies in terms of the mechanisms of the missing data. It is of utmost importance that researchers possess a prior understanding of the mechanism underlying missing data before imputation procedures. Such knowledge is critical for the development of input sets utilized in the imputations (Ali et al. 2023, Başakın et al. 2023). Here, incorporating parameters believed to be associated with missingness as auxiliary variables into the model can positively influence imputation results. Various methods exist that enable the statistical testing of the corresponding mechanisms, with some R programming language packages offering ready-made solutions (Jamshidian et al. 2014). However, only six of the reviewed studies (Hui et al. 2004, Alavi et al. 2006, Lo Presti et al. 2010, Demirhan and Renwick 2018, Liu et al. 2018, Kim et al. 2019, Afrifa-Yamoah et al. 2020) provided information regarding the mechanism underlying missing data. Similarly, the investigation of missing data patterns is of vital importance in missing data imputation studies. However, only three of the reviewed studies (Lo Presti et al. 2010, Lin and Wang 2014, Goodrich et al. 2021) provided information regarding the pattern of missing data. It is therefore crucial for researchers to consider the mechanism and pattern underlying missing data before imputation and toincorporate such information in their analyses to ensure accurate and valid results.

The last facet inspected under the umbrella of the pre-imputation stage is the standardization of the parameter values incorporated into the imputations. To define the parameters involved in the estimations, it is essential to standardize the input variables (Sovilj et al. 2016, Hasan et al. 2021). Multiple imputation procedures generally focus on solving a regression problem, and reducing the dataset to a dimensionless and evenly distributed structure (such as [0,1]) is regarded as a routine procedure in almost every regression model in the literature. Despite the vitality of standardization, only limited attempts (Lee and Chou 2004, Coulibaly and Evora 2007, Schmidt et al. 2008, Kim and Ahn 2009, Chen et al. 2016, Mao and Jian 2016) reported information regarding the standardization of parameters among the investigated articles.

6.2 Assessment of the during-imputation stage

The existing literature further recognized the critical role played by the essential points during the imputation of missing data. This paper examined the literature based on two aspects, i.e. the software used and the number of iterations in multiple imputations. Concerning the software utilized, the investigated studies reported the use of various programs to perform their imputation attempts. Among them, it was concluded that the most widely used tool is the R programming language (Wesonga 2015, Demirhan and Renwick 2018, Körner et al. 2018, Williams et al. 2018, Zitouna-Chebbi et al. 2018, Meher and Das 2019, Cheval et al. 2020, Costa et al. 2021). In addition, Schmidt et al. (2008) employed both STATISTICA and MATLAB. Likewise, Miró et al. (2017) utilized MATLAB, whereas many other software programs, such as Neuto-Solution (Abudu et al. 2010), SAS (Hui et al. 2004, Alavi et al. 2006), GeneXPro (Hasanpour Kashani and Dinpashoh 2012), C# (Altan and Ustundag 2012), and SPSS (Zainudin et al. 2015), were also adopted to perform the imputation of missing values that exist in the records of the variables related to ET.

Furthermore, in multiple imputation algorithms, the imputation process is repeated for different numbers of iterations. The appropriate number of iterations can vary depending on the size of the dataset, the number of independent variables, and the pattern of missingness. For instance, the default value for the number of iterations in the MICE package is 5. When examining studies that utilized multiple imputations, Hui et al. (2004), Turrado et al. (2014), Xie and Sun (2018), Kim et al. (2019), and Costa et al. (2021) selected the iteration value as “m = 5.” However, Miró et al. (2017) chose a value of 20. When considering datasets with varying sizes, it may be beneficial to explore different combinations of iteration values.

6.3 Assessment of the post-imputation stage

The evaluation of results obtained regarding imputation processes should essentially be carried out based on several aspects. Various methods can be used to evaluate the imputed data, including comparison with the observed data and comparison of the imputed data with the data from a complete case analysis. Furthermore, visual inspection of the imputed data can be helpful in detecting any systematic patterns or outliers that may indicate problems with the imputation model. Along with these essential aspects, this review further underpinned the role of run time evaluations, statistical distribution checks, assessing time series variance, and validation of the obtained results. Concerning the first aspect, one can conclude that the run time of imputation algorithms may vary depending on the size of the dataset as well as the utilized algorithms and computer capabilities. However, it is well known that the classic methods typically have shorter run times compared to modern techniques such as machine learning algorithms. Here, machine learning methods may require longer run times due to the optimization of hyperparameters. This can pose a challenge for researchers who need to make quick decisions in time-sensitive situations. Therefore, sharing information regarding the run times of different imputation techniques is crucial to provide a better overview for interested researchers in similar domains. Among the papers investigated within the scope of the current research, only three studies (Lin and Wang 2014, Turrado et al. 2014, Körner et al. 2018) have comparatively reported the computational costs of the utilized algorithms.

This study also assessed the previous work based on statistical distribution checks. Scholars who perform imputation attempts typically aim to minimize error regardless of the chosen method. Error levels can be measured through various approaches, providing important quantitative results regarding the imputation success. However, these quantitative metrics do not necessarily account for the distributional properties of the imputed values, which ideally should align with the distribution of the original dataset (Zhao et al. 2022). Here, the existing literature highlighted that testing the similarity between distributions can be achieved through statistical tests (Özger et al. 2020). Therefore, if the imputed values do not correspond to the dataset distribution, the imputation procedure should be reviewed and/or re-operated in order to achieve accurate imputation outcomes. We found only eight studies (Hui et al. 2004, Lo Presti et al. 2010, He et al. 2013, Wesonga 2015, Miró et al. 2017, Demirhan and Renwick 2018, Körner et al. 2018, Liu et al. 2018) that have investigated the similarity between imputed values and dataset distributions, highlighting the need for more research in this area.

Furthermore, when performing imputations for time series data, only the values of the parameter under investigation can be used for the imputation process. Therefore, it is crucial to ensure that the variance of the time series remains the same after imputations (Gao et al. 2018). Among 62 investigated papers, only Anjomshoaa and Salmanzadeh (2019) and Kunwor et al. (2017) checked the variance of the time series. The first used a quantitative metric to measure the variance values of imputed and gapped time series, while the latter provided a detailed analysis of the variance values after the imputation process. In order to reliably impute missing data, another important step that should diligently be taken into account is to conduct a training and validation process using the complete parts of the dataset used (Aghili et al. 2022). This process serves to demonstrate the accuracy of the imputation model and allows for the conclusion that the completed values are close to real values.

Using the entire dataset to obtain the imputation model without subjecting the models to any validation process, i.e. using the entire dataset for training, can lead to overfitting and result in completed values that are underestimated and/or overestimated. To tackle the overfitting phenomenon, cross-validation is one of the most commonly used methods in the pertinent literature (Yang et al. 2021, Eum and Yoo 2022). Of the studies examined within the scope of current research, 25% of the focalized research considered embracing a validation rationale, while others did not report the differentiation of training and validation stages.

The assessment of whether missing data imputation algorithms achieve optimal performance can be determined by utilizing various statistical evaluation metrics. These metrics are used to assess the degree of agreement between the predicted values and the actual values. Some of these metrics can be used to express the adequacy of a single model/algorithm quantitatively, while others are utilized to compare the results of multiple models (Zouzou and Citakoglu 2023). In this sense, metrics such as mean square error (MSE), root mean square error (RMSE), mean absolute error (MAE), and mean bias (MB) are used to measure which of two or more models makes predictions with lower errors, while some other metrics (such as coefficient of determination [R²], Nash-Sutcliffe efficiency [NSE], and Kling-Gupta efficiency [KGE]) measure the performance of a single model. Metrics that measure individual model performance typically range between [−1,1] or [0,1] and are able to measure the goodness of fit independently of the unit used (Bayram and Çıtakoğlu 2023). Hence, in evaluating the performance of these algorithms, it is crucial to diligently consider the choice of metric and its interpretability, as well as the potential implications of the missing data imputation process for subsequent statistical analyses and modelling outcomes.

It is worth mentioning that nearly 20% of the studies examined do not provide any information regarding the metric used for determining imputation accuracy, while 25% of them paid regard to only one metric. On the other hand, nearly half of the pertinent literature considered two or more metrics in assessing the accuracy of the imputation attempts. Diligent investigation demonstrated that the most commonly used metric is RMSE, with 34 occurrences, followed by MAE and R² with 17 and 16 occurrences, respectively. One can conclude that using different performance evaluation indicators helps increase the reliability of the obtained imputation results due to the various data structures representing divergent data-gathering processes.

7 Conclusions

Missing data is a common problem in meteorological data for many reasons, including instrument and equipment malfunction, errors in data collection, natural events, and lack of resources. In addition, the existence of missing data poses serious challenges for many modelling approaches in the hydrology domain, especially considering the importance of having complete datasets that allow for accurate and reliable analysis of water-related phenomena. Hence, the present research carried out an investigation of the existing literature with reference to calculating ET, which is one of the significant components of the hydrologic cycle. Among nearly 20,000 studies directly or indirectly related to the ET, a total of 62 studies were extracted to draw a comprehensive representation regarding the efforts devoted in the last two decades. As a result, the following conclusions and implications are drawn:

• The applicability and validity of the classical/parametric imputation methods rely heavily on various assumptions; hence, it is essential to approach these methods with caution and rigorously assess their underlying assumptions to avoid inaccurate conclusions or biased results.

• Assessment of the distribution of data is essential prior to commencing imputation procedures, while optimizing the run time of the imputation model and fine-tuning the model parameters allow efficient performance to be achieved in the during-imputation stage. Furthermore, within the post-imputation step, it is important to verify the stability of the imputed values; hence, assessing the constancy of variance and similarity of the distributions can help identify potential issues and provide insight into the robustness of the imputation model.

• Although various conventional techniques have been employed to handle missing data due to their ease of implementation, advanced imputation techniques that are particularly leveraged by the introduction of machine learning algorithms have emerged as more reliable alternatives to address missing data.

• Although attempts regarding the incorporation of deep learning techniques are limited thus far, these techniques show great potential for enhancing the performance of imputation models, particularly in datasets with a large number of observations.

• Interpretability of machine learning algorithms is an emerging concept in recent years; however, there was no study exploring the comprehensibility of missing data imputation procedures. Hence, integration of explainable artificial intelligence (e.g. SHapley Additive exPlanations [SHAP], Local Interpretable Model-Agnostic Explanations [LIME], etc.) can enhance the interpretability and transparency of the imputation results.

Furthermore, it is especially worth mentioning that, despite numerous studies examined within the scope of the current research, Scopus was the only database used to extract the relevant studies that exist in the literature. However, subsequent attempts could include different research databases to extend the scope of missing data imputation studies by the inclusion of conference papers, reports, and dissertations. Overall, a thorough investigation of the ET-related literature and the conclusions made in the current study are expected to assist the pertinent research community in understanding the dynamics of missing data imputation strategies and provide valuable insights regarding decision making processes for downstream analyses.

Disclosure statement

No potential conflict of interest was reported by the authors.