Mapping and Assessment of Evapotranspiration over Different Land-Use/Land-Cover Types in Arid Ecosystem

Abstract

Evapotranspiration (ET) is an essential process for defining the mass and energy relationship between soil, crop and atmosphere. This study was conducted in the Eastern Region of Saudi Arabia, to estimate the actual daily, monthly and annual evapotranspiration (ETa) for different land-use systems using Landsat-8 satellite data during the year 2017/2018. Initially, six land-use and land-cover (LULC) types were identified, namely: date palm, cropland, bare land, urban land, aquatic vegetation, and open water bodies. The Surface Energy Balance Algorithm for Land (SEBAL) supported by climate data was used to compute the ETa. The SEBAL model outputs were validated using the FAO Penman-Monteith (FAO P-M) method coupled with field observation. The results showed that the annual ETa values varied between 800 and 1400 mm.year−1 for date palm, 2000 mm.year−1 for open water and 800 mm.year−1 for croplands. The validation measure showed a significant agreement level between the SEBAL model and the FAO P-M method with RMSE of 0.84, 0.98 and 1.38 mm.day−1 for date palm, open water and cropland respectively. The study concludes that the ETa produced from the satellite data and the SEBAL model is useful for water resource management under arid ecosystem of the study area.

1. Introduction

Evapotranspiration (ET) is an essential component of the water cycle that connects hydrologic and biological processes. It is directly affected by water and land management, land-use change and climate variability [1, 2]. Therefore, estimation of ET in vast areas using efficient tools is important for optimum water resources management over different land-use systems [3].

The knowledge of ETa and its spatial distribution can have a great potential to develop new cost-effective indicators of irrigation performance and increase water use efficiency [4]. The ET can be estimated using many methods and techniques such as lysimeters, sap flow, eddy covariance, Bowen ratio, and scintillometer, which were accurate and efficient at the field scale [5]. However, these techniques cannot be used for large regional scale ET mapping due to prohibitive cost and logistical limitations [6]. Accordingly, remote sensing and biophysical modelling are adequate techniques for evaluating the ET patterns over large scale regions.

SEBAL was selected in this study to estimate the ETa because it can measure the ET without requiring quantifying the other complicated hydrological processes [7]. Also, SEBAL can identify the satellite image’s dry and wet pixels by showing a linear relationship between the surface temperature and the temperature gradient difference [8]. The FAO P-M method was adopted to validate the ETa obtained from the SEBAL model, since it used the actual climate data of the study area for quantifying the ETa process.

1.1 Application of LULC systems in hydrology

The study of land use and land cover (LULC) systems directly affects hydrology at different scales. Many studies showed that the LULC changes have clear impacts on the soil surface runoff, evapotranspiration, groundwater recharge, streamflow, and water balance over agricultural areas.

The impacts of LULC changes on hydrology were investigated in the Loess Plateau of China by Lu et al. [9] using hydrological modelling during the period 1995–2010. There was the transformation of farmland into forests, grassland, and built-up land. They showed slight increases in average annual potential evapotranspiration, actual evapotranspiration, and water yield at the basin scale, but soil water decreased between the two intervals. However, in sub-basins, obvious LULC changes did not have clear impacts on hydrology, and the impacts may be affected by precipitation conditions. The streamflow was also affected by the LULC changes in the Dinder and Rahad Rivers basins located in Ethiopia and Sudan [10]. The LULC of the catchment indicated a significant decrease and increase of the woodland and croplands were observed between 1972 and 2011. The effect of LULC change on streamflow was significant during 1986 and 2011, which could be attributed to the severe drought during the mid-1980s and the recent large expansions of cropland.

The hydrological process of a periurban catchment was quantified using high-resolution satellite images (0.50–2.50 m) in the Yzeron district of France [11]. The produced land covers maps of the district categorised into sub-catchments dominated by vegetation and imperviousness areas. According to the image processing and images characteristics, the calculated imperviousness rates were different, and lead to significant differences in the hydrological response. Wolfe et al. [12] compiled climate, geological, topographical, and land-cover data from the Prairie in Canada, and conducted a classification of watersheds using hierarchical clustering of principal components. Their analysis resulted in 7 classes based on the clustering of watersheds. The important defining variables identifying the watersheds clustering were climate, elevation, surficial geology, wetland distribution, and land cover. The authors indicated that developing management strategies is essential to prevent watersheds from future change.

Estimating the actual evapotranspiration and crop coefficients of an almond and pistachio orchard was explored in Central Valley, California, USA during an entire growing season by combining a simple crop evapotranspiration model with remote sensing data [13]. The authors used vegetation index NDVI derived from Landsat-8 to estimate the basal crop coefficient (Kcb), or potential crop water use. Their results showed that the model indicated a difference of 97 mm in transpiration over the season between both crops. However, the soil evaporation accounted for an average of 16% and 13% of the total actual evapotranspiration for almonds and pistachios, respectively. They concluded that the combination of crop evapotranspiration models with remotely-sensed data helps up-scaling irrigation information from plant to field scale and thus may be used by farmers for making day-to-day irrigation management decisions.

Furthermore, the up-scaling of the daily and seasonally ET using multisource remote sensing images was explored by Cha et al. [14] in the agricultural lands of the Kai-Kong River Basin, Xinjiang, China. They proposed a trapezoidal and a sinusoidal method to upscale daily ET values to seasonal ET. Moreover, the actual ET over LULC types in India’s Malaprabha River Basin was estimated using Landsat 8 data and surface energy balance models [15]. Their results demonstrate the challenge in actual ET estimation at a fine spatial resolution and highlight the importance of choosing a suitable algorithm.

1.2 Using of SEBAL model for ET estimation

The algorithms that use remote sensing products to estimate the energy balance and ET have become increasingly common [16]. Examples of these models are the two-source energy balance (TSEB), developed by Norman et al. [17]; the surface energy balance algorithm for land (SEBAL), formulated by Bastiaanssen (1995); and the mapping of evapotranspiration at high resolution with internalised calibration (METRIC), developed by Allen et al. [18].

SEBAL model is the most widely used algorithm for estimating ET throughout the world that involves applications for agricultural and water resources management, urbanisation impacts, aquifer recharge, and water balances [16].

The SEBAL model applied to determine the distribution of the ETa for analysing water use patterns over a large basin in Kenya [19]. In Indonesia’s Java Province, a method was developed using SEBAL to detect and describe the spatial variability of lowland rice ET [20]. The water used to assess irrigation system performance and management was pointed in the Indus Basin, Pakistan, using the ET estimated by SEBAL and MODIS data [21]. Remote sensing and biophysical modelling were used in recent studies to estimate the ET in different Saudi Arabia regions. Mahmoud and Alazba [22] estimated the ETa in the western and southern regions of Saudi Arabia during 1992–2014 using the SEBAL model, MODIS data and field observations. However, the application of SEBAL in Mara Basin of East Africa indicates that the ETa is measurable over different land-use types in data-scarce regions [23]. Daily ET monitoring using SEBAL model also found to be possible for improving water resources decision support over an oasis in the desert ecosystem [24]. The efficiency of the SEBAL model in estimating ET of pistachio crop in the Semnan province of Iran was investigated. The model shows good efficiency for estimating the actual ET of the pistachio product [25]. Moreover, SEBAL was used to calculate ET during the cultivation and harvesting of wheat crops in the Ilam province, Iran [26]. Thus the evaluation using SEBAL and the FAO-Penman–Monteith method showed that SEBAL has sufficient accuracy for estimating ET.

The Kingdom of Saudi Arabia (KSA) suffers a continuing water scarcity and almost around 90% of the agricultural sector’s water budget [27]. Groundwater is the primary source of water in the KSA, considering the limited precipitation and high agricultural demands. Moreover, the increased population of Saudi Arabia resulted in a significant increase in water use [28]. Also, the diversity of the LULC in the arid region is critical to water consumption. Accordingly, for effective water resources management in these regions, the impacts of LULC on hydrology need to be assessed. Hence, precise information of the ETa is crucial for policymakers and water planners to develop and formulate strategies for agricultural water utilisation. Therefore, the objective of this study is to assess the potential of Landsat-8 data and SEBAL model for estimating the daily, monthly and annual ETa under different LULC systems in Al-Ahsa region of Saudi Arabia.

2. Materials and methods

2.1 Study area

The study area is located in the eastern province of Saudi Arabia and lies about 300 km from Riyadh and 70 km west of the Arabian Gulf (Figure 1). The study area consists of Al-Ahsa Oasis, one of the leading agricultural centres in Saudi Arabia and dominated by date palm plantation [29]. The climate in the study area is hot and dry during the summer, with a relatively high air temperature that might exceed 45 °C. However, the winter is wet, with a minimum air temperature that might reach 5 °C. The annual rainfall is less than 250 mm and occurs mainly during the winter [30]. The dominated soils of the study area are sandy to sandy loam soil. Groundwater is the primary water source in the oasis and used mainly for irrigation, domestic and industrial purposes [31]. The land-use system in the study area is predominated by date palm plantation as the main agricultural activity. However, the cropping of rice and vegetables is also practiced in the area.

Figure 1.

Location of the study area.

2.2 Data

The data used in this study consist of remote sensing, climate information and field observations. They collected during April 2017–March 2018 to include the summer and the winter season of the study area. The summer season is considered for Aril–September, while the winter is for October–March.

A total number of 22 Landsat–8 images (path/row is 164/042) acquired from the United States Geological Survey (USGS-https://earthexplorer.usgs.gov/) to cover the entire study area, and the characteristics of these data are shown in (Table 1). The obtained Landsat-8 images have a cloud cover of less than 10%, and they have been geometrically and radiometrically corrected. All images bands were resampled into a pixel size of 30 m × 30 m using the nearest neighbour method. A global digital elevation model (DEM) is generated from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) known as ASTER GDEM obtained from the USGS website. It is a 30 m grid size DEM produced by the National Aeronautics and Space Administration (NASA) and the Ministry of Economy, Trade, and Industry of Japan (METI).

Sensor	Bands Type	Wavelength (μm)	Spatial Resolution (m); Resolution
Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS)	Band 1, Coastal aerosol	0.43–0.45	30
	Band 2, Blue	0.45–0.51	30
	Band 3, Green	0.53–0.59	30
	Band 4, Red	0.64–0.67	30
	Band 5, Near Infrared (NIR)	0.85–0.88	30
	Band 6, Short-wave Infrared 1 (SWIR-1)	1.57–1.65	30
	Band 7, Short-wave Infrared 2 (SWIR-2)	2.11–2.29	30
	Band 8, Panchromatic	0.50–0.68	15
	Band 9, Cirrus	1.36–1.38	30
	Band 10, Thermal Infrared 1 (TIRS-1)	10.60–11.19	100
	Band 11, Thermal Infrared 2 (TIRS-2)	11.50–12.51	100

Table 1.

Characteristics of Landsat-8 data used in this study [32].

The climate data was collected from two meteorological stations located within the study area. These data include air temperature, relative humidity, wind speed, net radiation, precipitation and vapour pressure, and all data collection on an hourly and daily basis. However, the field observations include the identification of the main land-use system in the study area. Besides field notes, site descriptions, and terrestrial photographs were taken to relate the site location to scene features.

2.3 Analysis of LULC

A field survey was conducted throughout the study area to identify the LULC classes during the study period. Global Positioning System (GPS) instrument was used to obtain accurate location point data for each LULC class included in the classification process. A total number of 115 ground control points (GCPs) were collected. A supervised maximum likelihood classification (MLC) method was employed to classify images. Based on the study objectives, the supervised classification applied in this study does not compare different classifiers. Therefore, the MLC was adopted to be the only classification method for this study. MLC widely used in remote sensing for image classification [33, 34, 35]. The accuracy assessment of the classified images was performed using 30% of the collected GCPs. Also, visual interpretation of the unclassified satellite images supported with the field observations was used to validate the LULC maps. However, to reduce bias, the stratified random sampling method was adopted to classify images [36]. Finally, the overall accuracy, user’s and producer’s accuracies, and the Kappa statistic were derived from the error matrices [37].

2.4 SEBAL method

The SEBAL model developed by Bastiaanssen et al. [38] was used to calculate the actual evapotranspiration from Landsat-8 satellite images. The model’s key inputs were the satellite measurements of surface albedo, normalised difference vegetation index (NDVI) and surface temperature (Ts). Also, the DEM and land-use map were used as additional input data. The DEM was applied for topographic and atmospheric correction [39]. However, the land-use map was used mainly to differentiate between LULC types exist in the study area. In addition to the satellite data, the SEBAL model requires minimum inputs of routine weather data (see the Data section). Figure 2 shows a flowchart that describes the SEBAL model process. The SEBAL model scripts were formed using the Spatial Modeller Tool of the ERDAS IMAGINE 9.2 software, and the ArcGIS 10.2 software was used for data mapping and visualisation.

Figure 2.

A flowchart explains the SEBAL model process.

The SEBAL algorithm computes the latent heat flux as the residue of the energy balance Equation [38, 40, 41]:

$$\mathit{\lambda ET}={\mathrm{R}}_{\mathrm{n}}-\mathrm{G}-\mathrm{H}$$E1

where ${\mathrm{R}}_{\mathrm{n}}$ is the net radiation over the surface (W/m²), $\mathrm{G}$ is the soil heat flux (W/m²), $\mathrm{H}$ is the sensible heat flux (W/m²), $\mathrm{\lambda ET}$ is the latent heat flux (W/m²) and $\mathrm{\lambda }$ is the latent heat of vaporisation (J/Kg).

2.4.1 Estimation of the Main surface parameters

The main surface parameters of the SEBAL include: the surface albedo, land surface emissivity (ε) and land surface temperature (Ts).

Surface albedo is defined as the ratio of reflected radiation from the surface to the incident shortwave radiation [40]. Thus the albedo is a single value that represents the integrated reflectance across the entire shortwave spectrum as represented by bands 2–7 of Landsat-8 data. Accordingly, the surface albedo (α) was estimated using the following formula [42]:

$$\alpha =0.2453{\mathrm{\alpha }}_2+0.0508{\mathrm{\alpha }}_3+0.1804{\mathrm{\alpha }}_4+0.3081{\mathrm{\alpha }}_5+0.1332{\mathrm{\alpha }}_6+0.05221{\mathrm{\alpha }}_7+0.0011$$E2

where α2, α3, α4, α5, α6, α7 represent the albedo of band2, band3, band4, band5, band6, and band7, respectively.

Surface emissivity is the thermal energy ratio radiated by the surface to the thermal energy radiated by a blackbody at the same temperature [42]. The NDVI is vital for calculating the land surface emissivity since it is used to estimate the vegetation coverage [43]. The calculation of the NDVI and vegetation coverage was as follows [8, 44]:

$$\mathit{\text{NDVI}}=\frac{\left(\mathit{NIR}-\mathit{RED}\right)}{\left(\mathit{NIR}+\mathit{RED}\right)}$$E3

$$P_v=\frac{\left(\mathit{\text{NDVI}}-{\mathit{\text{NDVI}}}_{\mathit{min}}\right)}{\left({\mathit{\text{NDVI}}}_{\mathit{max}}+{\mathit{\text{NDVI}}}_{\mathit{min}}\right)}$$E4

where $\mathrm{NIR}$ is the reflectance in the near-infrared band, which corresponds to band 5 in lansat-8 data, while the RED reflectance corresponds to band 4. ${\mathrm{P}}_{\mathrm{v}}$ is the vegetation coverage.

The surface emissivity was conditionally determined based on the NDVI values using Eq. (5) [45]:

$${\epsilon }_0={\mathrm{\epsilon }}_{\mathrm{v}}{\mathrm{P}}_v+{\mathrm{\epsilon }}_s\left(1-P_v\right)+{\mathrm{C}}_{\epsilon }$$E5

where ${\mathrm{\epsilon }}_0$ is the land surface emissivity, ${\mathrm{\epsilon }}_{\mathrm{v}}$ and ${\mathrm{\epsilon }}_{\mathrm{s}}$ are the vegetation and soil emissivity, respectively, and ${\mathrm{C}}_{\mathrm{\epsilon }}$ represents the surface roughness (${\mathrm{C}}_{\mathrm{\epsilon }}$ = 0 for homogenous and flat surfaces) taken as a constant value of 0.005 [46].

The ${\mathrm{P}}_{\mathrm{v}}$ values are conditioned with the NDVI ones. The land cover is classified as water when the NDVI ˂0. For the NDVI values range between 0 and 0.2, the land is considered covered with soil. The NDVI values of 0.2–0.5 the land cover are considered mixtures of soil and vegetation. However, when the NDVI >0.5, the land is considered covered with vegetation [43].

The land surface temperature was derived from the thermal bands. This step needs the spectral radiance to be converted into a sufficient brightness temperature. That means it has the black body temperature, assuming that the Earth’s surface is a black body. Consequently, brightness temperature was determined using the formula [47]:

$${\mathrm{T}}_{\mathrm{b}}=\frac{{\mathrm{K}}_2}{ln\left(\frac{{\mathrm{K}}_1}{{\mathrm{L}}_{\mathrm{\lambda }}}+1\right)}-273.15$$E6

where, ${\mathrm{T}}_{\mathrm{b}}$ is the satellite brightness temperature (C°), Lλ spectral radiance at top of the atmosphere, K1 and K2 are satellite calibration constants from the image metadata. The absolute zero value of −273.15C° was added to obtain results in Celsius [48].

The land surface temperature (${\mathrm{T}}_s$) was computed using Eq. (7) [49]:

$${\mathrm{T}}_s=\frac{{\mathrm{T}}_{\mathrm{b}}}{1+\left(\mathrm{\lambda }.\frac{{\mathrm{T}}_{\mathrm{b}}}{\mathrm{\rho }}\right).ln{\epsilon }_0}$$E7

where, λ is the wavelength of emitted radiance (11.5 μm for Landsat 5 &7; 10.9 μm for Landsat 8), ${\epsilon }_0$ is the land surface emissivity, $\mathrm{\rho }=\mathrm{h}.\frac{\mathrm{c}}{\mathrm{\sigma }}=1.438\times {10}^{-2}\mathrm{mK},\mathrm{\sigma }$ is the Boltzmann constant $\left(1.38\times {10}^{-23}\mathrm{J}/\mathrm{K}\right)$, $\mathrm{h}$ is the Planck’s constant $\left(6.626\times {10}^{-34}\mathrm{J}\mathrm{s}\right)$, and $\mathrm{c}$ is the velocity of light $(2.998\times {10}^8\mathrm{m}/\mathrm{s}$).

2.4.2 Determination of heat fluxes

The net radiation (${\mathrm{R}}_{\mathrm{n}}$) was calculated using surface reflectance and surface temperature (Ts):

$${\mathrm{R}}_{\mathrm{n}}={\mathrm{R}}_{\mathrm{s}}\downarrow -\propto {\mathrm{R}}_{\mathrm{s}}\downarrow +{\mathrm{R}}_{\mathrm{L}}\downarrow -{\mathrm{R}}_{\mathrm{L}}\uparrow -\left(1-{\mathrm{\epsilon }}_0\right){\mathrm{R}}_{\mathrm{L}}\downarrow$$E8

where ${\mathrm{R}}_{\mathrm{s}}\downarrow$ is the incoming short wave radiation (W/m²) (solar radiation), $\propto$ surface albedo (dimensionless), ${\mathrm{R}}_{\mathrm{L}}\downarrow$ is the incoming long wave radiation (W/m²), $\mathrm{RL}\uparrow$ is the outgoing long wave radiation (W/m²), and ε0 is the surface thermal emissivity (dimensionless). The calculation of these radiations was performed in Eqs. (9)–(11) [50]:

$$\mathit{RS}\downarrow =G_{\mathit{sc}}.cos\theta .r.{\tau }_{\mathit{sw}}$$E9

$$\mathit{RL}\uparrow ={\epsilon }_0.\sigma .T_s^4$$E10

$$\mathit{RL}\downarrow ={\epsilon }_{\propto }.\sigma .T_{\propto }^4$$E11

where ${\mathrm{G}}_{\mathrm{sc}}$ is the solar constant, 1367 W/m² and $cos\mathrm{\theta }$ is the cosine of the solar incidence angle, $\mathrm{r}$is the Earth-Sun distance (dimensionless), and ${\mathrm{\tau }}_{\mathrm{sw}}$ is atmospheric transmissivity. $\mathrm{RSI}\downarrow$ values range from 200–1000 W/m², depending on the image’s time and location and the local weather conditions [42]. $\mathrm{\sigma }$ is the Stefan-Boltzmann constant (5.67 × 10⁻⁸ W.m⁻².K⁻⁴), ${\mathrm{T}}_{\mathrm{s}}$ is the surface temperature (K), ${\mathrm{\epsilon }}_{\mathrm{\alpha }}$ is the atmospheric emissivity, and ${\mathrm{T}}_{\propto }$ is the atmospheric temperature (K). The empirical Eq. (6) was used for calculating the ${\mathrm{\epsilon }}_{\propto }$ [42].

$${\epsilon }_{\propto }=0.85\times {\left(-\mathit{ln}{\tau }_{\mathit{sw}}\right)}^{0.09}$$E12

The soil heat flux (G) is the rate of the heat flux stored or released into the soil and vegetation due to conduction. The ratio $\mathrm{G}/{\mathrm{R}}_{\mathrm{n}}$ was computed using Eq. (13) developed by [51]:

$$\mathrm{G}/{\mathrm{R}}_{\mathrm{n}}=\frac{T_s}{\propto }\left(0.0038+0.0074{\propto }^2\right)\left(1-0.98{\mathit{\text{NDVI}}}^4\right)$$E13

where ${\mathrm{T}}_{\mathrm{s}}$ is the surface temperature (C°), $\propto$ is the surface albedo, and NDVI is the Normalised Difference Vegetation Index (ranged between −1 and + 1). NDVI values between 0 and 0.2 correspond to bare soil or very sparse vegetation, the values greater than 0.2 represent vegetated areas. The typical estimates of $\mathrm{G}/{\mathrm{R}}_{\mathrm{n}}$ assumed to be 0.5 for water, 0.05–0.4 for agriculture and 0.2–0.4 for bare soil [42].

The sensible heat flux (H) is the rate of heat loss to the air by convection and conduction, due to a temperature difference. H was determined using the aerodynamic based heat transfer equation as follows [51]:

$$H=\frac{\left(\rho \times {\mathrm{C}}_{\mathrm{p}}\times \mathit{dT}\right)}{{\mathrm{r}}_{\mathrm{ah}}}$$E14

where $\mathrm{\rho }$ is air density (kg/m³), ${\mathrm{C}}_{\mathrm{p}}$ is air specific heat (1004 J/kg/K), $\mathrm{dT}$ (K) is the temperature difference (T1 – T2) between two heights (z1 and z2), and ${\mathrm{r}}_{\mathrm{ah}}$ is the aerodynamic resistance to heat transport (s/m). The ${\mathrm{r}}_{\mathrm{ah}}$ is computed for neutral atmospheric stability conditions as:

$${\mathrm{r}}_{\mathrm{ah}}=\frac{\mathit{ln}\left[\frac{Z_2}{Z_1}\right]}{u\ast \times k}$$E15

where ${\mathrm{Z}}_1$ and ${\mathrm{Z}}_2$ are heights in meters above the zero plane displacement of the vegetation, u* is the friction velocity (m/s) which quantifies the turbulent velocity fluctuations in the air, and k is von Karman’s constant (0.41).

2.4.3 Estimating the evapotranspiration

The instantaneous value of ET in equivalent evaporation depth was computed as:

$${\mathrm{ET}}_{\text{inst}}=3600\frac{\mathrm{\lambda ET}}{\mathrm{\lambda }}$$E16

where ${\mathrm{ET}}_{\text{inst}}$ is the instantaneous ET (mm/hr), 3600 is the conversion from seconds to hours, $\mathrm{\lambda ET}$ is the latent heat flux (W/m) consumed by ET, $\mathrm{\rho w}$ is the density of water (1000 kg/m³), and $\mathrm{\lambda }$ is the latent heat of vaporisation (J/kg) and was computed as:

$$\mathrm{\lambda }=\left[2.501-0.00236\left({\mathrm{T}}_{\mathrm{s}}-273.15\right)\times {10}^6\right]$$E17

The reference ET fraction (${\mathrm{ET}}_0\mathrm{F}$) or crop coefficient (kc) was calculated based on ${\mathrm{ET}}_{\text{inst}}$ for each pixel and ${\mathrm{ET}}_0$ was obtained from local ground weather stations.

$${\mathrm{ET}}_0\mathrm{F}={\mathrm{ET}}_{\text{inst}}/{\mathrm{ET}}_0$$E18

The daily values of ET (ET24) (mm/day) for each pixel was calculated as follows:

$${\mathrm{ET}}_{\mathrm{a}}={\mathrm{ET}}_0\mathrm{F}\times {\mathrm{ET}}_{0}24$$E19

where ${\mathrm{ET}}_0\mathrm{F}$ is the reference ET fraction, ${\mathrm{ET}}_{0}24$ is the cumulative alfalfa reference for the day (mm/day), and ${\mathrm{ET}}_{\mathrm{a}}$ is the actual evapotranspiration for the entire 24-hour period (mm/day).

The actual monthly and annual ET was calculated using daily ET data as follows [42]:

$${\mathrm{ET}}_{\mathrm{a},\text{period}}=\sum _{i=m}^n{\mathrm{ET}}_0\mathrm{F}\times {\mathrm{ET}}_{0}24$$E20

$${\mathrm{ET}}_{\mathrm{a},\text{annual}}=\sum {\mathit{ET}}_{a,\text{period}}$$E21

Allen et al. [18] showed that one cloud-free satellite image per month is sufficient to develop ${\mathrm{ET}}_0\mathrm{F}$ curves for seasonal ${\mathrm{ET}}_{\mathrm{a}}$ estimations.

2.4.4 Validation of the SEBAL model evapotranspiration

The produced ETa from Landsat-8 images and SEBAL model was validated using the FAO P-M method [52]. The FAO P-M was used to calculate the reference crop evaporation (${\mathrm{ET}}_0$) from the actual climate data in the study area based on Eq. (22):

$${\mathrm{ET}}_0=\frac{0.408∆\left(\mathrm{Rn}-\mathrm{G}\right)+\mathrm{\gamma }\frac{900}{\mathrm{T}+273}\mathrm{U}2\left(\mathrm{es}-\mathrm{ea}\right)}{∆+\gamma \left(1+0.34\mathrm{U}2\right)}$$E22

where ${\mathrm{ET}}_0$ is reference evapotranspiration (mm/day), $∆$ is slope vapour pressure curve (kPa/°C), $\mathrm{\gamma }$ is psychrometric constant (kPa/°C-1), T is mean daily air temperature at 2 m height (°C), $\mathrm{U}2$ is wind speed at 2 m height (m/s), es is saturation vapour pressure (kPa), ea. is actual vapour pressure (kPa), $\left(\mathrm{es}-\mathrm{ea}\right)$ represents the saturation vapour pressure deficit (kPa).

The crop coefficient (Kc) for the different croplands and the open water determined based on Allen et al. [52]. The ${\mathrm{ET}}_0$ obtained from the FAO P-M method and the kc were used to calculate the ${\mathrm{ET}}_{\mathrm{a}}$ depending on actual weather data as follows:

$${\mathrm{ET}}_{\mathrm{a}}={\mathrm{ET}}_0\times \mathrm{kc}$$E23

The ${\mathrm{ET}}_{\mathrm{a}}$ resulted from the FAO P-M method was used to validate the ${\mathrm{ET}}_{\mathrm{a}}$ obtained from SEBAL model.

A linear correlation and the root mean square error (RMSE) between the measured (FAO P-M) and the SEBAL daily ETa was computed [13]. The RMSE was calculated as follows [26]:

$$\text{RMSE}=\sqrt{\frac{{\sum }_{i=1}^N{\left(O_i-P_i\right)}^2}{N}}$$E24

where ${\mathrm{O}}_{\mathrm{i}}$ represents the observed values of the FAO P-M method as the standard model; ${\mathrm{P}}_{\mathrm{i}}$ represents the estimated values from the SEBAL algorithm; and $\overline{{\mathrm{O}}_{\mathrm{i}}}$ and $\overline{{\mathrm{P}}_{\mathrm{i}}}$ are the mean values from the FAO P-M method and SEBAL model, respectively.

3. Results and discussions

3.1 LULC mapping

The LULC map of the study area showed that the main identified classes were the date palm, cropland, bare land, urban land, aquatic vegetation, and water (Figure 3). The area occupied by each LULC type within the oasis boundaries is shown in Table 2. The date palm covers about 131 km² of Al-Ahsa Oasis, and it is the most important land-use class for the local and national economy. Croplands used only 144 km² of the oasis area, and it is dominated by rice and vegetables. The bare land class occupies most of the oasis area with 4759 km². Bare lands dominated by desert and rock outcrops also occurred in this class. Most of the urban land occurs on the oasis periphery, as most oasis land is under agricultural use. The aquatic vegetation and water classes occupy together an area of about 17 km². Al-Dakheel [53] reported that date palm covered about 92% within the oasis boundary.

Figure 3.

LULC map of the study area.

LULC	Area		Classification Accuracy (%)
	Km^2	%	User’s	Producer’s
Date Palm	131	2.5	95	94
Cropland	144	2.7	83	83
Bare Land	4759	89.9	95	94
Urban Land	241	4.5	82	86
Aquatic Vegetation	7	0.1	86	92
Water	10	0.2	100	100
Overall Accuracy	90%
Kappa Statistic	88%

Table 2.

Areas and accuracy assessment of the LULC classes.

The overall classification accuracy of the LULC map was 90%, with a kappa index of 88%, while the user’s and producer’s accuracies differed along with LULC types (Table 2). This accuracy level indicates that the classification method adopted in this study effectively produced a compatible LULC map over the study area.

3.2 Analysis of land surface parameters

The statistical means values of land surface albedo show that it was raining between 0.46 and 0.51 during Apr. 2017-Mar. 2018 (Figure 4a). However, the spatial distribution of land surface albedo indicates higher in the bare lands areas than the other LULC types (Figure 5a). Nevertheless, the seasonal variation shows that the surface albedo was higher during the summer (April, July) compared to the winter (November, February). The surface albedos levels were found lower in the vegetation-covered areas than the exposed soil [8].

Figure 4.

Mean values of land surface parameters: (a) albedo; (b) emissivity and (c) surface temperature. Bars denoted standard error.

Figure 5.

Spatial distributions of land surface parameters. (a) Albedo; (b) emissivity and (c) surface temperature.

The land surface emissivity means values clearly show that it does not vary along the study period (Figure 4b). However, the land surface emissivity’s spatial patterns indicated higher in the open water and lowered in the date palm and croplands classes (Figure 5b). The lowest values of the land surface emissivity were observed in bare lands. These results are inconsistent with that reported by Kong et al. [8].

The land surface temperature statistics showed that the average minimum value of 299 K was observed in January and February 2018, while the maximum mean value of 332 K occurred on 10 August 2017 (Figure 4c). Figure 5c indicated that the highest land surface temperature values were shown in summer, and they associated with the bare lands. Nevertheless, during the winter, the difference in the land surface temperature between the bare lands and the other LULC types was 10–20 K. The land surface temperature is an essential parameter for quantifying the ET process among the different LULC types of the study area. Accordingly, the land surface temperature estimation is essential for land classification, energy budget estimations, and crop production [54].

3.3 Heat fluxes analysis

The surface heat fluxes estimated over the different LULC types can affect the ET amount measured along the study period at varying scales.

The net radiation flux’s statistical values show that the values in April, May and June are clearly higher than those in October, November, December and January (Figure 6a). As shown in Figure 7a, the open water received the highest amount of net radiation followed by the aquatic vegetation, date palm, and croplands. However, the lowest net radiation found in the bare lands. Moreover, the net radiation flux increased significantly in April and July compared to November and February.

Figure 6.

Mean values of heat fluxes: (a) net radiation; (b) soil heat and (c) sensible heat. Bars denoted standard error.

Figure 7.

Spatial distributions of heat fluxes. (a) Net radiation; (b) soil heat and (c) sensible heat.

The increase of the net radiation in April and May might be due to rising ground-air temperature and the gradual death of the sparse vegetation over bare lands. The change in the surface energy budgets due to irrigation results in increased net radiation over agricultural lands [55]. Also, the significant variation of the net radiation flux during the study period could be due to the region’s heterogeneous LULC types.

The soil heat flux tendency showed higher in the summer recording average value of 91 W/m² on 23 June 2017, while the lowest mean value was 5 W/m² observed on 16 December 2017 (Figure 6b). The spatial distribution of the soil heat flux for the vegetation cover classes was higher than that of the bare lands during the summertime, while the difference was slight in winter (Figure 7b).

The sensible heat fluxes mean values do not vary consistently along with the summer and wintertime. The average highest value of 203 W/m² detected on 04 April 2017, while the lowest was on 10 August 2017 (Figure 6c). According to the spatial distribution in Figure 7c, the bare and urban lands’ sensible heat flux show relatively high values compared to the vegetation lands.

The increasing trend in sensible heat flux during April and May can be attributed to increased net radiation flux under persistent irrigation land used for agriculture and urban areas. However, Amatya et al. [56] indicated that the increase in wind speed and the ground-air temperature difference could increase sensible heat flux.

3.4 Analysis of the actual evapotranspiration

The spatial distribution of the daily ETa over the study area during 20 April 2017–2022 March 2018 is shown in Figure 8. The temporal patterns of the daily ETa showed that the highest values were observed during the peak summertime in July and August. The mean daily ETa values for the different LULC types in the oasis are shown in Figure 9. It is clear that the daily ETa for the water bodies and aquatic vegetation was between 5.6–8.7 mm.day⁻¹ in summer and about 2.3–5.6 mm.day⁻¹ in winter. However, the date palm and croplands showed daily ETa of 3.5–8.0 and 2.0–3.6 mm.day⁻¹ for the winter and summer, respectively.

Figure 8.

The spatial distribution of the daily ETa.

Figure 9.

Temporal variation of mean daily ETa for different LULC systems.

The variability of the daily ETa values for the Date palm and cropland during the summer and winter times mainly attributed to the irrigation water distribution, soil salinity, drainage, and agricultural practices and their impact on moisture and salinity in the root zone. Under Saharan oasis conditions, soil texture, plot size, and farmers’ practices in particular irrigation found to have significant effects on the daily ETa [57]. In Saudi Arabia, the daily ETa of date palm was observed to decrease during winter and increased during summer, depending on the crop’s growth stage [58]. Also, the daily water consumption of major cropping systems in Saudi Arabia varied spatially depending on cropping practises and climatic conditions [22].

The mean daily ETa for the urban land ranged from 1.3 to 4.5 mm.day⁻¹ throughout the study period (Figure 6). The urban land is covered with some lanes and parks that make seasonal variations of the daily ETa within the study area. Nevertheless, the low values of the daily ETa showed in the bare land mainly resulted from the sparse vegetation in this land-use system.

The variation of the daily ETa estimates for the different LULC types under the oasis-arid ecosystem indicates that they change significantly throughout the seasons [24].

Figure 10 shows the spatial pattern of the monthly and annual ETa across the study area. The high rates of the ETa found to be between 80 to 200 mm.month⁻¹ during July, August and September for the water, aquatic vegetation, date palm and cropland. However, at the beginning of the summer in April, May and June, the ETa rates were 60–100 mm.month⁻¹ for the same LULC types. Nevertheless, in the winter during Oct 2017–Mar 2018 the ETa ranged between 40 to 140 mm.month⁻¹ for the water, aquatic vegetation, date palm and cropland land-use systems.

Figure 10.

The spatial patterns of the monthly and annual ETa.

The mean annual ETa produced by SEBAL model for the different LULC types in the study area is shown in Figure 11. The ETa rates of date palm trees ranged from 800 to 1400 mm.year⁻¹ during Apr. 2017 to Mar. 2018. The annual water consumption for date palm was highly variable. This might be attributed to the type of irrigation system and the age variations of date palm trees along the oasis.

Figure 11.

Mean annual ETa produced by SEBAL for the different LULC types in the study area during April 2017–March 2018. Bars denoted standard error.

The open water evaporation lost was around 2000 mm.year⁻¹, while an average of 1600 mm.year⁻¹ was evaporated from aquatic vegetation. Nevertheless, croplands showed lower annual ETa of 800 mm.year⁻¹ compared to the date palm. The main crops like rice and vegetables can be cultivated during a particular time of the year in the oasis; therefore, croplands showed relatively low annual ETa compared to date palm areas.

The annual ETa of urban lands was 400 mm.year⁻¹. Urban lands are affected by the irrigation of trees, lanes and parklands, which resulted in consumption of a large amount of the oasis groundwater. The annual evaporation from bare lands was very low (200 mm.year⁻¹) and less than the long-term average rainfall. Bare lands equipped most of the study region areas, and they covered with sand dunes (Figure 11). Moreover, the bare lands characterised by low vegetation coverage levels and low water contents in the soil surface [41]. Also, very low rates of ETa from bare soil observed in the western and southern parts of Saudi Arabia [22].

3.5 Validation, limitation and uncertainty of SEBAL model

The FAO P-M was used as a standard method to validate the SEBAL model [8, 59]. The validation measurement for ETa between the SEBAL model and FAO P-M method for the different land-use system is shown in Figure 12. A significantly high level of agreement can be observed between the two methods for the selected LULC types. The RMSE for the most validated LULC system in the study area found to be less than 1.0 mm.day−1. However, the RMSE was slightly higher for cropland areas (Figure 12d), mainly due to the method used for calculating the crop coefficient (kc) for the different crop types within the croplands system. The SEBAL model does not require kc information because the model biophysical properties estimated kc as part of the SEBAL process. However, the FAO P-M computed kc based on the characteristics and climatic regions for the different crops. Accordingly, the kc of vegetables (tomato and cucumber) ranged between 0.5–1.15, and for the rice, it was 1.0–1.35 [52]. Nevertheless, the date palm’s kc was in the range of 0.9–0.95, while it was 1.05 for the open water [52].

Figure 12.

Linear correlation between the FAO P-M and SEBAL model for the different LULC types. (a) Open water; (b) aquatic vegetation; (c) date palm; (d) cropland.

Moreover, it seems that SEBAL underestimates the ETa for croplands since they were diverse in terms of crop type and growing season. Also, the kc used to calibrate the ETo was made only for a few experimental plots. Therefore, the variability in the ETa values predicted by SEBAL and measured by the FAO P-M method was slightly higher. However, the low values of ETa in winter do not affect the outputs of seasonal ETa over croplands, as farmers use a little water for irrigation. The kc of crops can vary during the growing season, depending on their growth stage [60]. Rahimi et al. [3] reported that the application of SEBAL for estimating the ETa over agricultural land in the Tajan catchment of Iran resulted in RMSE of 1.49 mm.day-1. However, the daily estimated ETa for both cropland and grassland in the Midwestern United States using SEBAL contributed to RMSE ranging between 1.74 and 2.46 mm.day⁻¹ [6].

SEBAL model utilises satellite imagery and a small set of surface weather information, including wind speed and air temperature at a reference height of 2 m to solve the energy balance [38]. However, the steps used for the SEBAL process are complicated. The model requires an experienced modeller and a substantial number of working hours for image processing. Also, the large model number of mathematical formulas increases the possibility of human errors [42]. Furthermore, the procedure used by SEBAL for selecting the “hot” and “cold” pixels endmembers involves not only an analysis of the surface temperature but also an understanding of other biophysical properties, such as vegetation indices, surface albedo, and LULC type [16, 61]. Therefore, identifying these endmembers requires modeller intervention and knowledge of the biophysical parameter values and ranges.

4. Conclusions

This study demonstrates the power of remote sensing data and the biophysical modelling for quantifying the ETa process over an arid ecosystem in Saudi Arabia. The estimated mean annual ETa was 2000 mm.year⁻¹ for open water and varied between 800 and 1400 mm.year⁻¹ for date palm. However, it was 1600 mm.year⁻¹ for the aquatic vegetation, while an average of 800 mm.year⁻¹ was observed in croplands. The validation measure showed a significant agreement level between the SEBAL model and the FAO P-M method with RMSE of 0.62, 0.84, 0.98 and 1.38 mm.day⁻¹ for aquatic vegetation, date palm, open water and cropland respectively.

The obtained ETa information will help Saudi Arabia formulate strategies to reduce the gap between the water supply and demand in the irrigated areas. Furthermore, the ETa patterns mapped over the diverse LULC systems can be used as a baseline framework for sustainable water resources management and agrometeorological services in the different regions of Saudi Arabia. However, conducting long-term ETa studies using remote sensing data coupled with the implementation of different models and field tools may improve the assessment of the ETa dynamic process in arid regions.

Acknowledgments

The author’s acknowledge the Deanship of Scientific Research at King Faisal University, for the financial support under the annual research project (Grant No. 186002).