Introduction

Evapotranspiration (ET) is a physical aerodynamics process in which water moves from liquid to gaseous stage, whereas bringing from the soil to the atmospheric surface1. It denotes to both evaporation from vegetation and soil fields and transpiration from plants. Two distinct processes (evaporation and transpiration) happen concurrently, and there is no alternative way of differentiating one from the other. ET is one of the basic elements of the water cycle, and its estimation is necessary to drought mitigation and management as well as other fields, including agro-meteorology, hydrology, climatology, and environmental studies2,3,4,5. For this reason, many drought indices such as Reconnaissance Drought Index (RDI)6, Standardized Precipitation Evapotranspiration Index (SPEI)7, and the Water Surplus Variability Index (WSVI)8 are based on the ET. Besides, two closely associated terms and concepts are potential evapotranspiration (ETp) and reference evapotranspiration (ETo) that estimate the atmospheric evaporation demand. ETp is defined as the rate of water transpired in a specific time by a crop, fully shading the ground, of constant height with adequate soil water setting in the outline9,10. On the other hand, ETo is expressed as the ET rate from a reference crop surface, where the reference crop surface is a theoretical grass or alfalfa with accurate and recognized characteristics1,11. However, the definition of ETo is more precise and specific than the ETp. The application of the terms ETp and ETo have been puzzled for several decades. One ideal example is that Hargreaves and Samani12 used the term “ETp” whereas again Hargreaves and Samani13 applied the term “ETo”. Under the well-known background of global warming in recent decades, the term ETo has been broadly applied in hydrology14,15,16, agronomy field17,18,19, irrigation engineering4,20,21 and meteorological field22,23. The application of ETo is also used in the studies of crop water demand. Therefore, knowledge of ETo is of great importance in agricultural water management, hydrological field, climate change, and irrigation practice24,25.

FAO recommended Penman–Monteith (FAO56-PM) method is the sole standard method for estimating ETo1,20,26. The main limitation of the FAO56-PM method is the difficulty in obtaining all necessary input data (air temperature, humidity, solar radiation, and wind speed). In such circumstances, simple equations or alternative methods are often used to estimate ETo27. The major advantages of an optimal empirical method are the simplicity, low cost, ease of application, and easy access to a few climatic input data measured in most of meteorological stations28. Thus, an optimal empirical model is vital for ETo estimation in data-limited regions across the globe, including Bangladesh.

The selection of an optimal empirical model for calculating the ETo is significant for agricultural water resources management, hydrological planning and irrigation designing11,24,29. In recent decades, a large number of empirical models have been developed to estimate the ETo, which have been widely reported in the literature (Table S1). These empirical models can be demarcated into five types mainly based on the data requirement: mass transfer-based30,31, temperature-based11, radiation-based13,32,33,34,35,36, combined-based24,33,34,35 and the pan evaporation-based models37,38,39 (Table S1). Earlier studies showed that the performance of various empirical models exhibited spatiotemporal variations, and most of these empirical models might be region-specific that enhanced the uncertainty problem in the identified spatiotemporal patterns. To solve this issue, local adjustment and validation of empirical models are required against the standard FAO56-PM model at various regions with various climatic contexts for accurate estimation of the ETo24. Most of the previous studies revealed that the performance of empirical models for estimating ETo showed significant regional differences40,41. For instance, the calibrated adjusted Hargreaves model performed better than the calibrated Priestley-Taylor model for measuring ETo in Serbia42. Quej et al.43 assessed the performance of the temperature-based ETo models and found that the Hargreaves–Samani model exhibited the best performance in a tropical sub-humid climate. Krishna44 pointed out that the Turc model was an optimal alternative for the estimation of the ETo under a humid subtropical climate, India. Li et al.24 reported that combination-based Valiantzas3 was the best model for estimating ETo in the humid to sub-humid region, China. On the contrary, Pandey and Pandey45 found that the Hargreaves–Samani method had a larger overestimation than the standard FAO56-PM in humid areas of India. These contrasting outcomes can be attributed to variations in regional climate, and geography. However, whether empirical models influence the computation procedure of the ETo against the FAO56-PM model remains uncertain. Therefore, it is crucial to conduct research appraising the performance of the empirical models in Bangladesh to determine an optimal alternative to ETo and their changes shifted overtimes at the regional scale and differ spatially14,15.

Bangladesh, the vast deltaic plain, a low-lying subtropical humid climatic country, is positioned in Southeast Asia, has a total land area of 147,700 square km46 (Supplementary Figure S1). The country has a complex geomorphic setting and complicated hydrologic system which comprises various water bodies, wetlands, floodplain, flood basins, agricultural land, forest, and hilly regions. The elevations of most regions of the country varied from 1 to 60 m above the mean sea level, which forms generally low-lying areas from the east to the west, making a so-called “delta-shaped” landform4. Nevertheless, Bangladesh is not only faced this type of difficulty in obtaining long-term and complete climatic datasets but also this poor country experiences similar phenomena owing to naturals such as a complicated hydro-geographic and climatic setting and humankind e.g., low economic growth, lack of proper knowledge and technological hindering causes. Under this circumstance, for the ETo appraisal of Bangladesh, an alternative empirical model depending on the limited climatic dataset is needed47. Therefore, it is of paramount importance to validate an appropriate alternative model which is easier in calculation procedure with fewer climatic variable requirements and good precision in comparison with the FAO56-PM model in various climatic sub-regions of Bangladesh. To the author's knowledge, so far, a systematic and thorough investigation for choosing an optimal alternative for estimating ETo has not been conducted in Bangladesh till now, particularly at a monthly and regional scale, which in itself is the novelty of this study.

Based on the aforementioned research gaps, in this research, 13 widely employed empirical models selected for performance evaluation, including four temperature-based model (Hargreaves–Samani, Hargreaves M1, Hargreaves M2, and Berti), one mass transfer-based models (WMO) and eight radiation-based models (Makkink, Priestly–Taylor, Jensen–Haise, Abtew, Irmak1, Irmak 2, Turc, and Tabari), based on the extensive literature review of meteorological variables, climatic regional differences, and their universal applicability. Subsequently, this study seeks three hypotheses: first, the various empirical models will generate considerably various outcomes for estimating the ETo at a monthly and regional scale; second, identifying the importance degree of empirical models that can indicate which model is outperformed against the FAO56-PM model at the regional scale and third, the simple linear regression can efficiently validate the 13 models against the FAO56-PM model in different sub-regions and whole Bangladesh. The specific objectives of this study are: (1) to analyze the spatiotemporal changes and the trends of ETo in Bangladesh for the period of 1980–2017 at a monthly scale, (2) to compare the performances of 13 empirical models against the FAO56-PM model for ETo estimation in climatic sub-regions of Bangladesh, (3) to choose an optimal alternative of the FAO56-PM ETo model, which will be easier in ETo quantification and apply few meteorological variables, (4) to identify the most outperformed empirical model against FAO56-PM model using the heuristic random forest method, and (5) to validate 13 empirical models using linear regression to opt an alternative empirical model against FAO56-PM model. The novelty of this research lies in employing 13 empirical models with a heuristic random forest model for the first time in Bangladesh that enables us to find an optimal alternative used for important environmental implication from the most reliable and outperformed equations for ETo computation in different climatic sub-regions of Bangladesh.

Data and methods

Study area description and data sources

Bangladesh, situated in Southeast Asia, geographically it encompasses between 20° 30′ N and 26° 45′ N latitudes and 88° 0′ E to 92° 45′ E longitudes (Fig. 1). Banglapedia48 divided Bangladesh into seven climatic sub-regions based on climatology and geography as shown in the Supplementary Material of Figure S2. The seven sub-regions are (1) south-eastern zone; (2) north-eastern zone; (3) northern part of the northern zone; (4) north-western zone; (5) western zone; (6) south-western zone and (7) south-central zone. Bangladesh experiences a sub-tropical humid monsoon climate with seasonal differences4.

Figure 1
figure 1

Map showing the geographical location of study area, prepared by ArcGis 10.5 (www.esri.com).

Full size image

Western Bangladesh has usually become drier compared to other regions in Bangladesh49. Here, climatic variability is a regular scenario. Long-term daily average relative humidity, minimum temperature, maximum temperature, wind speed (at 2 m height), net radiation, evapotranspiration across the country are, respectively, 80%, 21.39 °C, 29.94 °C 1.32 ms−1, 10.44 MJm−2day−1, and 3.72 mm day−1.

Bangladesh Meteorological Department (BMD) runs only 43 meteorological stations across the country. The meteorological stations are unevenly distributed all over the country and most of the stations are located in the south-eastern regions. These meteorological stations are available for the climatic dataset, although some of the stations are newly established after the 1990s in Bangladesh and they do not have long-term data records (www.bmd.gov.bd). When more climatic variables are required, the dataset from a smaller number of stations was available. Due to these drawbacks, 20 stations were chosen for ETo estimation over the 38 years from 1980 to 2017. These selected 20 stations embody the seven climatic sub-regions of the country. Daily minimum (Tmin) and maximum temperature (Tmax) (°C), mean relative humidity (Hr) (%), wind speed (Uz) (Knots) and sunshine hour (h day−1) datasets of 20 stations were sourced from the BMD. Net radiation (Rn) and wind speed at 2 m height (U2) cannot directly be measured by weather stations. Daily Rn and U2 were estimated using the procedures recommended by Allen et al.1 with the available meteorological datasets. A brief geographical and meteorological description of the selected stations is found in the Supplementary Material of Table S2. However, missing data in almost all the 20 stations was found. After the initial screening test, missing data of the 20 stations were less than 5% for the period of 1980–2017. Missing data for each station were filled by the existing records for the respective days from the adjacent neighbor stations. It is worthy to note that sunshine hour dataset in this study is continuous with no missing data. More details about the fill-up of missing meteorological datasets is given in the Supplementary Material (Table S3). The BMD follows the guideline of World Meteorological Organization (WMO) for weather dataset collection and record archiving. However, quality control of the dataset was primarily undertaken thoroughly by checking namely, positive values of parameters, for example, Tmin is lower than Tmax, and humidity is less than 100%. The homogeneity tests of the dataset were conducted to exhibit any anomaly in the dataset50. All of the datasets were passed through the quality control by the staff of the BMD.

FAO56 Penman–Monteith model (FAO56-PM model)

The FAO56-PM equation is used for estimating daily ETo of this study. This model is well-known as the standard model for estimating ETo across the whole world, which was proposed by Allen et al.1 The original form of FAO56-PM model is expressed by the following Eq. (1):

$${\text{ET}}_{{\text{o}}} = \frac{{0.408\Delta \left( {{\text{R}}_{{\text{n}}} - {\text{ G}}} \right) + {\upgamma }\frac{900}{{{\text{T}} + 273}}{\text{U}}_{2} \left( {{\text{e}}_{{\text{s}}} - {\text{e}}_{{\text{a}}} } \right)}}{{\Delta + {\upgamma }\left( {1{ } + { }0.34{\text{U}}_{2} } \right)}}$$
(1)

where, ETo is the reference evapotranspiration (mm day−1),\({\mathrm{R}}_{\mathrm{n}}\) isthe net radiation atcrop surface (MJm−2 day−1), Gis the soil heat flux density (MJ m−2 day−1), T is the average daily air temperature at 2 m height (°C), \({\mathrm{U}}_{2}\) is the wind speed at 2 m height (ms−1), \({\mathrm{e}}_{\mathrm{s}}\) is the saturation vapour pressure (kPa),\({\mathrm{e}}_{\mathrm{a}}\) is the actual vapour pressure (kPa), \({\mathrm{e}}_{\mathrm{s}}-{\mathrm{e}}_{\mathrm{a}}\) is the saturation vapour pressure deficit (kPa), \(\Delta\) is the slope of vapour pressure curve (kPa °C−1), γ is the psychrometric constant (kPa °C−1). Allen et al.1 recommended G = 0. The detailed procedures of ETo estimation is found in FAO 56 paper1.

Rn is calculated by the Eqs. (2–11):

$${\text{R}}_{{\text{n}}} = {\text{R}}_{{{\text{ns}}}} - {\text{R}}_{{{\text{nl}}}}$$
(2)
$${\text{R}}_{{{\text{ns}}}} = ({1} - \alpha ){\text{R}}_{{\text{s}}}$$
(3)
$${\text{R}}_{{\text{s}}} = \left[ {{\text{a}}_{{\text{s}}} + {\text{b}}_{{\text{s}}} \frac{{\text{n}}}{{\text{N}}}} \right]{\text{ R}}_{{\text{a}}}$$
(4)
$${\text{Ra}} = \frac{{24{ }\left( {60} \right)}}{\pi } {\text{G}}_{{{\text{sc}}}} {\text{d}}_{{\text{r}}} [\omega_{{\text{s}}} {\text{sin}}(\varphi ){\text{sin}}(\delta ) + {\text{cos}}(\varphi ){\text{cos}}(\delta ){\text{sin}}(\omega_{{\text{s}}} )]$$
(5)
$${\text{d}}_{{\text{r}}} = {1 }0.0{33}\;{\text{cos}}\left( {\frac{2\pi }{{365}}{\text{J}}} \right)$$
(6)
$$\delta = 0.{4}0{\text{9 sin}}\left( {\frac{2\pi }{{365}}{\text{J}} - 1.39} \right)$$
(7)
$$\omega_{{\text{s}}} = {\text{arccos }}[ - {\text{tan }}(\varphi ){\text{ tan }}(\delta )]$$
(8)
$${\text{Radians}} = \pi /{18}0 \, \left( {\text{decimal degrees}} \right)$$
(9)
$${\text{R}}_{{{\text{nl}}}} = \sigma \left[ {\frac{{{\text{T}}_{{{\text{max}}}} {\text{k}}^{4} + {\text{T}}_{{{\text{min}}}} {\text{K}}^{4} }}{2}} \right]\left( {0.{34}{-}0.{14}\surd {\text{e}}_{{\text{a}}} } \right)\left[ {{1}.{35}\frac{{{\text{R}}_{{\text{s}}} }}{{{\text{R}}_{{{\text{so}}}} }} - 0.{35}} \right]$$
(10)
$${\text{R}}_{{{\text{so}}}} = (0.{75} + {2} \times { }10^{ - 5} {\text{Z}}){\text{ R}}_{{\text{a}}}$$
(11)

U2 is calculated from the following Eq. (12) recommended by Allen et al.1,

$${\text{U}}_{{2}} = {\text{U}}_{{\text{z}}} \frac{4.87}{{{\text{In}}\left( {67.8{\text{z}} - 5.42} \right)}}$$
(12)

where, Rns is the net solar or shortwave radiation (MJ m−2 day−1), Rnl is the net outgoing longwave radiation (MJ m−2 day−1), Rs is the global solar or shortwave radiation (MJ m−2 day−1), N and n are, respectively, the maximum and actual possible sunshine duration, Ra is the extraterrestrial radiation (MJ m−2 d−1), Gsc is the solar constant (0.0820 MJ m−2 min−1), dr is the inverse relative distance Earth-Sun, ωs is the sunset hour angle (rad), φ is latitude (rad), δ is solar declination (rad), J is the number of the day in the year between 1 (1 January) and 365 or 366 (31 December), σ is the Stefan-Boltzmann constant (4.903 × 10−9 MJ K−4 m−2 day−1), α is albedo (α = 0.23), Tmaxk and Tmink are, respectively, the maximum and minimum absolute temperatures during 24-h, and Rso is the clear sky solar radiation (MJ m−2 day−1). Allen et al.1 recommended 0.25 for as and 0.50 for bs. Uz is measured wind speed at Zm above ground surface (ms−2) and z is respective station elevation above sea level (m).

According to Allen et al.1, Saturation Vapour Pressure (es), Actual Vapour Pressure (ea), Slope Vapour Pressure Curve (∆) and Psychrometric Constant (γ) are calculated by the following Eqs. (13–19), respectively:

$${e}_{s}=\frac{{e}^{0}\left({T}_{max}\right) +{e}^{0}({T}_{min})}{2}$$
(13)
$${e}^{0}\left({T}_{max}\right)=0.6108\;exp\;\left[\frac{17.27{T}_{max}}{{T}_{max}+237.3}\right]$$
(14)
$${e}^{0}\left({T}_{min}\right)=0.6108\;exp\;\left[\frac{17.27{T}_{min}}{{T}_{min}+237.3}\right]$$
(15)
$${e}_{a}=\frac{Hr(mean)}{100}\left[\frac{{e}^{0}\left({T}_{max}\right) +{e}^{0}({T}_{min})}{2}\right]$$
(16)
$$\Delta =\frac{4098\left[0.6108 \;exp\;\left(\frac{17.27 T}{T+237.3}\right)\right]}{{\left(T+237.3\right)}^{2}}$$
(17)
$$\gamma = \frac{CpP}{{\varepsilon \lambda }} = 0.{665} \times {1}0^{{ - {3}}} {\text{P}}$$
(18)
$${\text{P}} = 101.3\left( {\frac{293 - 0.0065Z}{{293}}} \right)^{5.26}$$
(19)

where, es is the mean saturation vapour pressure (kPa), \({e}^{0}\left({T}_{max}\right) \mathrm{and}{ e}^{0}({T}_{min})\) are the saturation vapor pressure at maximum and minimum temperature, respectively. ea is the actual vapour pressure function (kPa) and Hr is the mean relative humidity. Tave, Tmax and Tmin are the mean, maximum and minimum air temperature, respectively, in °C and exp [·] is 2.7183 (base of natural logarithm) raised to the power [··]. P is the atmospheric pressure (kPa), λ is the latent heat of vaporization (2.45 MJ kg−1), Cp is the specific heat at constant pressure (1.013 × 10–3 MJ kg−1 °C−1), ε is the ratio molecular weight of water vapour/dry air (0.622).

Empirical models

A primary survey of literature clearly showed that the 13 ETo empirical model performed usually well in various sub-regions worldwide. Abtew51, Jensen and Haise52, Irmak53, Makkink54, Priestley-Taylor55, Hargreaves-Samani13, Berti56, WMO30, Tabari40, and Turc57 models were chosen to compare to the FAO56-PM model. The 13 empirical models were chosen based on the available input meteorological variables, universal acceptance and their applicability worldwide (Table S1). The Hargreaves–Samani (HS), Hargreaves M1 (HM1), Hargreaves M2 (HM2) and Berti models used in this study, as the HS, HM1 and HM2 models require only the temperature and extraterrestrial radiation datasets and Berti model requires only temperature data, making these models less complex. Therefore, the 13 empirical models used in the present study can be classified into the three classes: four temperature-based model (Hargreaves–Samani, Hargreaves M1, Hargreaves M2, and Berti), one mass transfer-based models (WMO), eight radiation-based models (Makkink, Priestly–Taylor, Jensen–Haise, Abtew, Irmak1, Irmak2, Tabari, and Turc). The performances and application of these models had never been validated in Bangladesh so far. The studied models, input parameters, computed equations with references are outlined in Table 1.

Table 1 The original form of the 13 empirical models associated with the input parameters.
Full size table

Performance evaluation of 13 empirical models

Performance evaluation of 13 empirical models, based on the accuracy of each model for estimating ETo, was undertaken by six statistical criteria. The six statistical criteria were the mean bias error (MBE)58; mean absolute error (MAE), correlation of determination (R2), root mean square error (RMSE)43; relative error (RE), Nash–Sutcliffe efficacy coefficient (NSE)59 expressed by the following Eqs. (20–25):

$$MAE=\frac{1}{n}\sum_{i=1}^{n}\left|{ET}_{o,s}-{ET}_{o,i}\right|$$
(20)
$$MBE=\frac{1}{n}\sum_{i=1}^{n}{ET}_{o,i}-{ET}_{o,s}$$
(21)
$$NSE=1-\frac{\sum_{i=1}^{n}{({ET}_{o,s}-{ET}_{o,i})}^{2}}{\sum_{i=1}^{n}{({ET}_{o,s}-\stackrel{-}{{ET}_{o,s})}}^{2}}$$
(22)
$${R}^{2}=\frac{\sum_{i=1}^{n}{({ET}_{o,s}-{ET}_{o,i})}^{2}}{\sum_{i=1}^{n}{({ET}_{o,s}-{\stackrel{-}{ET}}_{o,i})}^{2}}$$
(23)
$$RE=\frac{{ET}_{o,i}-{ET}_{o,s}}{{ET}_{o,s}}$$
(24)
$$RMSE=\sqrt{\frac{1}{n}\sum_{i=1}^{n}{({ET}_{o,s}-{ET}_{o,i})}^{2}}$$
(25)

where, \({\mathrm{ET}}_{\mathrm{o},\mathrm{s}}\), \({\mathrm{ET}}_{\mathrm{o},\mathrm{i}}\) and n are the observed ETo (estimated by FAO56-PM), estimated ETo (estimated by empirical models) and total observations, respectively.

Modified Mann–Kendall test

Modified Mann–Kendall (MMK) test is a non-parametric test which was applied for detecting the increasing and decreasing trend of ETo,s in Bangladesh during 1980–201760. To carry out this MKK test, it is imperative to confirm the serial autocorrelation of the time series dataset. Hence, the serial autocorrelation should be excluded before employing the MMK test. To exclude the serial autocorrelation, the trend free pre-whitening approach proposed by Yue and Wang61 has been utilized. The original form of MK test62,63 statistics (S) is as followed:

$$S=\sum_{i=1}^{n-1}\sum_{j=i+1}^{n}sgn({X}_{j}-{X}_{i})$$
(26)
$$sgn\left( \theta \right) = \left\{ {\begin{array}{*{20}l} 1 \hfill & {if\;\theta > 0} \hfill \\ 0 \hfill & {if\;\theta = 0} \hfill \\ { - 1} \hfill & {if\;\theta = 0} \hfill \\ \end{array} } \right.$$
(27)

Direction of increasing or decreasing trend is indicated by S. The variance64 of S is followed by the Eq. (29):

$$V\left(S\right)=\frac{n\left(n-1\right)\left(2n+5\right)-\sum_{i=1}^{n}{t}_{i}i\left(i-1\right)(2{t}_{i}+5)}{18}$$
(28)

V*(S), is the modified variance61 given by following Eq. (30):

$${V}^{*}\left(S\right)=V\left(S\right)\cdot \frac{n}{{n}^{*}}$$
(29)

n/n* is termed as correction factor and denoted65 by following Eq. (31):

$$\frac{n}{{n^{*} }} = 1 + 2 \cdot \mathop \sum \limits_{k = 1}^{n - 1} \left( {1 - \frac{k}{n}} \right) \cdot \rho_{k}$$
(30)

Test statistic Z is calculated by following Eq. (32):

$$Z = \left\{ {\begin{array}{*{20}l} {\frac{{S - 1}}{{\sqrt {V^{*} (S)} }}} &\quad {S > 0} \\ 0 &\quad {S = 0} \\ {\frac{{S + 1}}{{\sqrt {V^{*} (S)} }}} &\quad {S < 0} \\ \end{array} } \right.$$
(31)

Positive Z statistic indicates increasing trend of ETo,s and negative Z statistic indicates decreasing trend of ETo,s in Bangladesh.

Sen’s slope of estimator

Sen’s slope of estimator66 was applied for calculating the change of ETo,s in Bangladesh per decade and the statistics is as followed:

$${\text{Q}} = xj-xkj$$
(32)

Q is denoted as the slope between xj and xk.

The spatial distributions of the monthly ETo and its trends are mapped.

Spatial distributions of the monthly meteorological variables; ETo,s, ETo,i, trends and the other examined variables are mapped by the inverse distance weighted interpolation model in ArcGIS 10.5 software.

Random forest (RF) model

The RF is a heuristic decision tree-based supervised machine learning model67 that is appropriate for addressing the existence of the over-fitting problem to the decision trees, and other machine learning algorithm68. The RF is most robust, can handle numerous heterogeneous covariates, and has been effectively employed into the hydrological field69, genetic engineering field70 and hydro-meteorological field71. The RF model has been benefited from the two more powerful algorithms e.g., bagging and random binary trees, which are called the powerhouse of this model. For developing the RF model, the number of trees and features in each split is essential. RF is a classifier which comprises of an assortment of classifier trees \({f}_{m}(x)\) for m = 1, …, M which relies on the parameters and every single tree casts a unit vote for input x71. Each tree generates an individual class which then combined and the majority vote predicts the final results. Present study optimized its accuracy with 100 trees, 1 execution slot, 5 seeds and with maximum depth 1. As a tree-based ensemble learning model, this model has extensively used to evaluate the importance degree of any climatic dataset in various regions25,72. To the best of author’s knowledge, the RF model has not yet been employed to explore the importance degree of 13 empirical models against the FAO56-PM model in Bangladesh. The RF model is used to know which model is most reliable and dominant for estimating ETo. More detailed about the RF model can be found elsewhere46,71.

Results

Spatial distribution of meteorological variables

Figure 2 represents the distribution of multi-year mean meteorological variables of Tave, Tmin, Tmax, Rn, U2, and Hr from 1980 to 2017. Distribution of Tave (Fig. 2a), Tmin (Fig. 2b) and Tmax (Fig. 2c) showed almost similar results. Sub-region VI showed the higher values of historical temperature, while sub-regions II and III entirely showed the lowest temperature. Sub-region V showed the lowest temperature for the distribution of Tave and Tmin and moderate temperature for Tmax. Sub-regions I, IV and VII showed the moderate values. The higher rate of net radiation was observed in the sub-region I and lower rate of net radiation was found in the sub-regions II and III (Fig. 2d). Average values of Rn were seen in the sub-regions IV, V, VI, and VII. The highest wind speed was found in the sub-region of VI and the lowest in the sub-regions II, III and V (Fig. 2e). Sub-region I experienced a comparatively higher rate of Hr as it located near the Bay of Bengal from where this region took available moisture and sub-region V experienced the lower rate of Hr (Fig. 2f).

Figure 2
figure 2

Spatial distribution of multi-year mean meteorological variables in Bangladesh, prepared by ArcGis 10.5 (www.esri.com).

Full size image

Spatial and temporal patterns of ETo,s and ETo,i

Figure 3 represents the long-term multi-year mean monthly trends of ETo,s (FAO56-PM). Most parts of sub-regions I and VII showed the higher values of MMK-Z statistic and the sub-regions III, IV and V showed the lowest values of MMK-Z statistic (Fig. 3a). Moderate values are shown in the sub-regions II and VI. In general, the rates between increasing and decreasing of ETo,s was from 72.18 to − 72.17 mm per decade (Fig. 3b). Figure 3c shows the nature of the trend whether it was significant or insignificant. The significant increasing trend of ETo,s was detected in Bhola, Cumilla, Feni (\(\alpha = 0.01\)); Rangamati, and Patuakhali (\(\alpha = 0.05\)) in the sub-regions I and VII.

Figure 3
figure 3

Representation of the multi-year mean monthly ETo,s trends of (a) MMK-Z values (mm); (b) Sen’s slope estimation and (c) station wise increasing or decreasing trends of Bangladesh, prepared by ArcGis 10.5 (www.esri.com).

Full size image

All the weather stations of sub-regions II (\(\alpha = 0.1\)); III, IV (\(\alpha = 0.01\)) and V (\(\alpha = 0.05\)) showed a significant decreasing trend of ETo,s. Faridpur, Madaripur, Dhaka and Barishal of sub-regions I and VII showed an insignificant decreasing trend of ETo,s. Cox’s Bazar (\(\alpha = 0.05\)); Teknaf, Sandwip and Chattogram (\(\alpha = 0.01\)) of sub-region I; Jashore (\(\alpha = 0.05\)) of sub-region VI and Mymensingh (\(\alpha = 0.01\)) and Khulna (\(\alpha = 0.1\)) of sub-region VII showed a significant decreasing trend of ETo,s.

Spatial distribution of multi-year mean monthly ETo,s and ETo,i from 1980 to 2017 is presented in Fig. 4. The higher value (4.12 mm) of ETo,s seen in the sub-regions V, VI and some parts of region I. The lower value (3.46 mm) of ETo,s seen in the sub-regions II, III and some part of the sub-regions I, IV and VII. The distribution of ETo,i showed the homogeneous distribution of ETo,s. The high-low values of spatial distribution of ETo,1 (temperature-based); ETo,6; ETo,7; ETo,10 (radiation-based) models were analogous with that of ETo,s. ETo,5 (mass transfer-based) and ETo,13 (radiation-based) models showed the most heterogeneity with ETo,s. Sub-regions V and VI showed the highest value of ETo (computed by most of the empirical models and FAO56-PM) and all the models found that the sub-regions II and III experienced a lower rate of ETo. All the models revealed that moderate ETo was experienced by the sub-regions I, IV and VII. Approximately, similar zonation (based on the highest and lowest values) was observed between ETo,s and ETo,6.

Figure 4
figure 4

Spatial distribution of multi-year mean monthly ETo,s and ETo,i in Bangladesh, prepared by ArcGis 10.5 (www.esri.com).

Full size image

Temporal distribution of long-term monthly ETo,s, and ETo,i in different sub-regions, as well as whole Bangladesh for the period of 1980–2017, is shown in Fig. 5. The highest rate of reference evapotranspiration (ETo,s and ETo,i) occurred in April in all the regions of Bangladesh. The lowest ETo (both ETo,s and ETo,i) found in January and December. Except for sub-regions II and III, the range of the rate of ETo,s, and ETo,i was the same in all the sub-regions. Among the 13 empirical models, ETo,5; ETo13 (elevated the lowest values than ETo,s) and ETo,11 (elevated the highest values than ETo,s) models showed the most deficit values of ETo,i compared to ETo,s identifying least suitable method for estimating ETo. Conversely, the values estimated by ETo,6 showed the closest values to ETo,s demonstrating as the most preferable model for estimating ETo. ETo,1 and ETo,7 also estimated values with the smallest difference with ETo,s in all the regions, confirmed as the preferable method for estimating ETo. Moderate values estimated by ETo,2; ETo,4; ETo,8; ETo,9; ETo,10 and ETo,12 models compared to ETo,s.

Figure 5
figure 5

Temporal distribution of multi-year mean monthly ETo,s and ETo,i in different sub-regions and whole Bangladesh.

Full size image

Long term inter-annual variation of ETo,s, and ETo,i from 1980 to 2017 are represented by Fig. 6. It also demonstrates similar results as shown in Fig. 5. The largest deviation from the ETo,s values occurred by estimating ETo,i (values lower than 1) by both ETo,5 and ETo,13 models. The Fig. 6 shows that the rate of ETo in each sub-region along with whole Bangladesh, estimated by both ETo,s and ETo,i models, declined gradually. Unlike the multi-year monthly distribution, the range (from 0 to 5.5 mm) of ETo values was nearly similar in every sub-regions and Bangladesh. Like the multi-year monthly distribution, ETo,6; ETo,1; ETo,7 and ETo,10 showed the very closest values to ETo,s in each sub-region. Values larger than that of ETo,s found by ETo,11. Based on the spatiotemporal distribution of ETo,i estimated by 13 empirical models in each sub-region and whole Bangladesh, the empirical models can be ranked in ascending order based on the closest values to ETo,s as ETo,6 > ETo,1 > ETo,7 > ETo,10 > ETo,3 > ETo,2 > ETo,12 > ETo,4 > ETo,8 > ETo,9 > ETo,11 > ETo,13 > ETo,5.

Figure 6
figure 6

The inter-annual variations of ETo,s and ETo,i in different sub-regions and whole Bangladesh.

Full size image

Performance appraisal of 13 empirical models for estimating ETo

Figure 7 shows the long term monthly RE of ETo,i by spatial distribution for the period of 1980–2017. RE of the maximum ETo,i covered mutually positive and negative values. The lowest RE was observed in ETo,i calculated by ETo,6. Furthermore, ETo,1; ETo,7; ETo,10 and ETo,12 models also produced lower RE, respectively. The worst performance with the high RE belongs to the ETo,5 and ETo,13 models. The relative error in different sub-regions varied with the variation of different empirical models. The ETo,11, ETo,9, ETo,8, ETo,2, ETo,3 and ETo,4 models also recognized as the worst models, respectively, producing high RE.

Figure 7
figure 7

Spatial distribution of relative error values for multi-year mean monthly ETo,i in Bangladesh, prepared by ArcGis 10.5 (www.esri.com).

Full size image

Table 2 shows the long term mean monthly and annual RMSE for 13 selected ETo,i models for whole Bangladesh. The higher RMSE was produced by the ETo,5; ETo,13 and ETo,11 models, respectively. The ETo,6, ETo,3, ETo,7 and ETo,1 models generated the least RMSE for calculating ETo. The descending order of the other models based on performance was ETo,2 > ETo,9 > ETo,4 > ETo,10 > ETo,8 > ETo,12.

Table 2 Root Mean Square Error (RMSE) values of the 13 empirical models for calculating ETo,i at the monthly and annual timescales for Bangladesh. The bold values indicate the optimal model among other models.
Full size table

Table 3 represents the long term mean monthly and annual MAE values of 13 empirical models used for estimating ETo,i. Lower values of MAE indicates higher accuracy. Respectively, ETo,6 and ETo,1 models produced the smallest MAE values for both the monthly and annual ETo,i values in Bangladesh. Higher MAE values produced by the models (respectively) of ETo,5, ETo,13, ETo,11, ETo,9 and ETo,8 in calculating both the monthly and annual ETo,i.

Table 3 Mean Absolute Error (MAE) values of 13 empirical models for calculating ETo,i at the monthly and annual timescales for Bangladesh. The bold values indicate the optimal model among other models.
Full size table

Table 4 shows the NSE coefficient of empirical models used for calculating ETo,i at the annual and monthly time scale in Bangladesh. Approximately, all the models gave negative NSE value indicating the least correlated method with ETo,s. In both monthly and annual timescales, ETo,6 outperformed as it gave positive NSE value, except for the month of October and November. The performance accuracy of ETo,1 was comparatively higher than that of other models except for ETo,6, ETo,13, ETo,5 and ETo,11 models. In annual timescale, the order of the performance of the models was ETo,6 > ETo,7 > ETo,1 > ETo,10 > ETo,12 > ETo,4 > ETo,2 > ETo,8 > ETo,9 > ETo,3 > ETo,11 > ETo,5 > ETo,13.

Table 4 Nash–Sutcliffe efficiency (NSE) coefficients of the 13 empirical models for calculating ETo,i at the monthly and annual timescales for Bangladesh. The bold values indicate the optimal model among other models.
Full size table

MBE at both long term annual and monthly time scale in calculating ETo,i in Bangladesh is shown in Table 5. ETo,6 and ETo,1 showed higher performance accuracy, producing lower MBE values. Again, ETo,5; ETo,13 and ETo,11 models showed the least performance accuracy for calculating ETo,i. Other models showed moderate performance accuracy in estimating ETo,i. From the above discussions of performance accuracy, ETo,6 revealed as the suitable alternative model for estimating ETo in Bangladesh. Whereas, ETo,5; ETo,13 and ETo,11models explored as the worst performed empirical models inappropriate for estimating ETo in Bangladesh. Figure 8 shows the scatter plots of daily ETo,s vs. ETo,i of Bangladesh during 1980–2017. The ETo,6; ETo,7 and ETo,11 models outperformed among 13 empirical models with higher r-value (0.92). ETo,10 (r-value 0.91) and ETo,12 (r-value 0.90) performed well, respectively, which produced r-values ≥ 0.90.

Table 5 Mean Bias Error (MBE) values of the 13 empirical models for calculating ETo,i at the monthly and annual timescales for Bangladesh. The bold values indicate the optimal model among other models.
Full size table
Figure 8
figure 8

Scatter plots showing the comparisons of daily ETo,i and ETo,s in Bangladesh.

Full size image

All the scatter plots produced r values with significant p-value (< 0.05) indicating a strong correlation between ETo,s and ETo,i. Among 13 empirical models, ETo,5 recognized as the worst model with a lower accuracy and reliability than other models as this model produced the r-value of 0.78 which is less than 0.80. From the comparisons of 13 empirical models by RE, RMSE, MAE, NSE, MBE and scatter plots, ETo,6 model found as the best suitable alternative model against the ETo,s model for calculating ETo for all the sub-regions and whole Bangladesh.

Exploring the most suitable model for ETo computation using the RF algorithm

The importance degree analysis of the empirical models was calculated by the RF method at different sub-regions of Bangladesh. The RF model depicted a suitable significant model followed by other analyses against the FAO56-PM model for all of the sub-regions (Table 6). The best performed model is highlighted by bold face and the least performed model is marked by italic face (Table 6) for estimating ETo. As can be seen from Table 6, the ETo,6 outperformed for ETo estimation in all the sub-regions of Bangladesh as it produced a comparatively higher importance degree. Whereas, ETo,5 found as the worst model producing the lowest importance degree among all the models.

Table 6 Importance degree of 13 empirical models against FAO56-PM model in seven sub-regions in Bangladesh using RF model.
Full size table

Validation of the best alternative model for ETo,s

Validation is very important for determining the most suitable model from several potential models. All the applied empirical models need to validate to ensure whether the pre-analyses gave the accurate results or not and to find the best suitable alternative empirical model against the FAO56-PM. Validation was undertaken by utilizing the linear correlation method. Following Peng et al.11, the linear correlation was calculated by the below Eq. (33):

$${ET}_{o,s}=\frac{{ET}_{o,i}-b}{a}$$
(33)

where, a and b denoted as fitted coefficients.

Table 7 shows the fitted a, b and R2 values of correlation between ETo,s and ETo,i in seven sub-regions and whole Bangladesh. A strong correlation between ETo,s and ETo,i was found in all sub-regions and whole Bangladesh. Values of R2 is greater than 0.8 in every sub-regions and in whole Bangladesh for all the 13 empirical models indicating a strong correlation between ETo,s and ETo,i. The model which is highly correlated with the ETo,s is highlighted with light pink color. Among 13 models ETo,6 performed best, producing greater R2 values than the other models. ETo,6 model is simple among all the models used in this study as this model utilized only Tmax and Rs. Rs is calculated from Tmax and Tmin. Tmax and Tmin are available everywhere in each region and can estimate easily. So, it can be affirmed that the ETo,6 is the best suitable, preferred, accurate, simple and reliable model, highly consisted with the results of pre-analyses (spatial and temporal distribution, performance evaluation and importance degree analysis), for ETo estimation in all the sub-regions and whole Bangladesh.

Table 7 Fitted a, b and R2 values of correlation between ETo,s and ETo,i in seven sub-regions and whole Bangladesh.
Full size table

Discussions

Increasing and decreasing as well as significant and non-significant trends were found across the country from 1980 to 2017. The rates between increasing and decreasing of ETo,s was 72.18 to − 72.17 mm per decade in this study. For example, Bhola, Cumilla, Feni, Rangamati and Patuakhali showed a significant increasing trend of ETo,s and non-significant decreasing trend of ETo,s found in Faridpur, Madaripur, Dhaka and Barishal of sub-regions I and VII (southeastern and south-central regions). The significant decreasing trend of ETo,s showed by the sub-regions of II, III, IV and V. Cox’s Bazar, Teknaf, Sandwip, Chattogram, Jashore, Mymensingh and Khulna showed a significant decreasing trend of ETo,s. From the above results, it is evident that the trend of ETo,s was decreasing gradually in Bangladesh from 1980 to 2017. Rahman et al.15 found that most of the area of Bangladesh showed decreasing trends and some parts of the study area showed an increasing trend of ETo,s which is analogous to the results of this study. Decreasing the trend of ETo,s may be the results of worldwide climate change impacts. Spatial distribution of the long-term mean monthly between ETo,s and ETo,6 (Abtew) in Bangladesh shown a homogenous pattern. By contrary, ETo,5 (WMO) and ETo,13 (Turc) models produced the least close values to ETo,s in terms of spatial distribution. Temporal distribution of long-term monthly ETo,s, and ETo,i in Bangladesh found that the highest and lowest rate of ETo occurred in April and January–December, respectively. Peng et al.11 and Li et al.24 showed that the highest and lowest rate of ETo,s occurred in China in July (dissimilar to this study) and January–December which is similar to this study. Long term inter-annual variation of ETo,s, and ETo,i in Bangladesh revealed that ETo,5 and ETo,13 models produced very lower values, identifying the least performed method for calculating ETo. Previous studies5,11,24,29 along with this study revealed that ETo values obtained by both ETo,s and ETo,i models were very closer to each other in the month of January–December (cold season) and highest discrepancy occurred among them in the hot summer season.

Long term monthly RE (relative error) of ETo,i revealed that the WMO and Turc models were the least suitable with greater RE values and Abtew model was an optimal alternative with lower RE values against the FAO56-PM for calculating ETo in Bangladesh. The values of RMSE, MAE, NSE, and MBE had the strong concurrence with the RE exploring the same results as Abtew (ETo,6) was the best alternative and WMO (ETo,5); Turc (ETo,13) were the least suitable model for estimating ETo in Bangladesh. Gabriela and Irmak73 evaluated the impact of the meteorological variables on the estimates of 13 empirical models in various regions and found that the Doorenbos and Pruitt (ETo,12) ranked top in the three regions of Iran under sub-humid to sub-arid climate conditions which are in disagreement with the results of this study. Identification of the best suitable model against ETo,s varies from region to region and country to country. It might be due to the variation of geographical and meteorological variations from one country to another country and input model combinations. Correlation between long term daily ETo,s and ETo,i in Bangladesh was explored that a very strong correlation aligns with no one line (1:1) existed between ETo,s and ETo,6 (Abtew) with the r values of 0.92. Li et al.24 also found a strong correlation (R2 was 0.972) between daily ETo,s and ETo,13 (Valiantzas 3) in China. Xystrakis and Matzarakis35 found a strong correlation (r-value 0.993) between monthly ETo,s, and ETo,i (Turc) in Greece. There also existed a strong correlation (r-value 0.996) between ETo,s, and ETo,i (Blaney–Criddle) explored by Tabari et al.40 in Iran. Peng et al.11 found the largest correlation between monthly ETo,s, and ETo,653 in China. Present study found the Abtew (ETo,6) model as the most reliable model for estimating ETo, compared to the other empirical models in all the sub-regions of Bangladesh, as this model produced a higher importance degree.

Evaluation of the performance of different empirical models for estimating the ETo,s was finally validated by Eq. (33). This study explored that Abtew (ETo,6) model outperformed other models with R2 values ranged from 0.856 to 0.989 at all the sub-regions of Bangladesh. This is in good agreement with the earlier performance appraisal results in which RMSE, MAE, MBE, and NSE are the lowest in ETo,6 model. The main reason is that this model has high precision, easy, consistent; requiring less climatic datasets and strong association with FAO56-PM. This model provides satisfactory outcomes and generally uses simple computable parameters and has easy model forms. Djaman et al.33 also found a similar result as this study that the Abtew was the best alternative model against the FAO56-PM in New Mexico, USA. Li et al.24 found Valiantzas 3 as the outperformed model among 13 empirical models for estimating ETo with R2 values ranged from 0.882 to 0.993. Peng et al.11 explored ETo,656 as the best-performed model at all the sub-regions and EMC among ten empirical models for calculating ETo with R2 values ranged from 0.87 to 0.99. Shiri5 showed Priestley–Taylor outweighed the other 6 empirical models for estimating ETo with R2 values ranged from 0.636 to 0.792. Mohawesh38 found Penman as the best-performed model in different regions of Jordan for estimating ETo,s with R2 values ranged from 0.66 to 0.7874. Similarly, the mass-transfer-based model was the optimal model in computing ETo compared to the other models in humid regions in Iran40 and forest regions in Greece32. Based on the accuracy, reliability, simplicity and higher correlation with ETo,s, the most suitable method for ETo calculation in Bangladesh is ETo,6 model.

Conclusions

In this study, daily meteorological datasets from 20 weather stations from seven sub-regions in Bangladesh for the period of 1980–2017 were used. A widespread comparison between ETo,s (calculated by FAO56-PM) and ETo,i (calculated by HS, HM1, HM2, BT, WMO, ABT, IR1, IR2, MAK, PT, JH, TAB, TR, respectively) has been carried out. The possible roles of these 13 empirical models against the FAO56-PM were also explored in this study. Out of 20 stations, 5 stations showed an increasing trend of ETo,s; 11 stations showed a decreasing trend and 4 stations showed no trend of ETo,s. Spatiotemporal distribution of ETo,s, and ETo,i revealed that the model proposed by Abtew model showing the closest distribution of ETo,i to ETo,s. RE, RMSE, MAE, MBE, and NSE were employed for evaluating the empirical models which were identified ETo,6 as the outperformed model with the lowest errors for calculating ETo in different sub-regions and whole Bangladesh. By contrast, ETo,5 (WMO) and ETo,13 (Turc) models selected as the poorer alternative models with the higher statistical errors. RF model also confirmed the Abtew as the outperformed model. The linear regression model showed that a strong linear correlation was found between FAO56-PM and Abtew model. Validation by using Eq. (33) explored the similar outcomes that the ETo,6 model outperformed than the other models. This study recommends the model proposed by the Abtew (ETo,6) as the best alternative model with high accuracy, reliability and lowest errors for all the sub-regions and whole Bangladesh for calculating ETo when full climatic datasets for FAO56-PM model are unavailable. Future study should be focused on the evaluation of machine learning ensemble models for estimating daily ETo in Bangladesh. This research is a vital scientific contribution to ETo quantification and influential empirical models in Bangladesh where the large set of meteorological datasets could not be acquired. This study provides an important guidance for agricultural water practices, hydrological processes and irrigation management in Bangladesh, also useful as well as the similar subtropical climate region elsewhere in the world.