# Estimation of Potential Evapotranspiration across Sri Lanka Using a Distributed Dual-Source Evapotranspiration Model under Data Scarcity

Authors: Udara Senatilleke; Himasha Abeysiriwardana; Randika K. Makubura; Faisal Anwar; Upaka Rathnayake
Journal: Advances in Meteorology
 (2022)

---

## Abstract

Evapotranspiration estimations are not common in developing countries though most of them have water scarcities for agricultural purposes. Therefore, it is essential to estimate the rates of evapotranspiration based on the available climatic parameters. Proper estimations of evapotranspiration are unavailable to Sri Lanka, even though the country has a significant agricultural contribution to its economy. Therefore, the Shuttleworth–Wallace (S-W) model, a process-based two-source potential evapotranspiration (PET) model, is implemented to simulate the spatiotemporal distribution of PET, evaporation from soil (ETs), and transpiration from vegetation canopy (ETc) across the total landmass of Sri Lanka. The country was divided into a grid with6km×6km cells. The meteorological data, including rainfall, temperature, relative humidity, wind speed, net solar radiation, and pan evaporation, for 14 meteorological stations were used in this analysis. They were interpolated using Inverse Distance Weighting (IDW), Universal kriging, and Thiessen polygon methods as appropriate so that the generated thematic layers were fairly closer to reality. Normalized Difference Vegetation Index (NDVI) and soil moisture data were retrieved from publicly available online domains, while the threshold values of vegetation parameters were taken from the literature. Notwithstanding many approximations and uncertainties associated with the input data, the implemented model displayed an adequate ability to capture the spatiotemporal distribution of PET and its components. A comparison between predicted PET and recorded pan evaporations resulted in a root mean square error (RMSE) of 0.75 mm/day. The model showed high sensitivity to Leaf Area Index (LAI). The model revealed that both spatial and temporal distribution of PET is highly correlated with the incoming solar radiation fluxes and affected by the rainfall seasons and cultivation patterns. The model predicted PET values accounted for 80–90% and 40–60% loss of annual mean rainfall, respectively, in the drier and wetter parts of the country. The model predicted a 0.65 ratio of annual transpiration to annual evapotranspiration.


---

## Body

## 1. Introduction

Evapotranspiration (ET) is an integral part of the global hydrological cycle and regional water budget, as ET represents the loss of the surface and soil water to the atmosphere as water vapor by the combined actions of the two processes: evaporation from the surface water bodies; bare soil and other surfaces that intercept rainwater, and transpiration from plants [1]. The primary driver of ET is solar radiation, which provides the latent heat requirement for water vaporization [2]. In addition, several other factors such as wind speed, humidity, air temperature, soil type, crop type, land use type [2, 3], and rainfall [4] affect the rate of ET. Because of the large number of influencing factors and their heterogeneity over a watershed, accurate estimation of ET is challenging, specifically under data scarcities. Thus, estimation of potential evapotranspiration (PET)—the potential amount of water that could evaporate and transpire from a vegetated landscape with unlimited water supply to the surface [5], has been mostly utilized instead of actual ET estimation. [1, 6].A few dozen empirical and semiempirical PET models, such as Penman [2], Thornthwaite [7], Priestley-Taylor [8], Monteith [9], and Shuttleworth and Wallace [10], have been developed over the last few decades. These models are based on different sets of assumptions and inherit various limitations; hence, the results generated are highly inconsistent [6]. Some of the PET models, such as Penman [2] and Priestley-Taylor [8], estimate potential evaporation over water surfaces but do not account for transpiration from vegetation cover. Monteith [9] developed the Penman–Monteith (P-M) model, one of the widely used models to estimate PET, assuming the vegetation canopy as a single uniform cover or ‘big-leaf,’ thus accounting for the transpiration process. One of the drawbacks of the P-M model is that it neglects sparse vegetation. Several researchers, i.e., Shuttleworth and Wallace (S-W) (1985), Choudhury and Monteith [11], Mo et al. [12], have developed PET models by extending the P-M model by incorporating sparse canopy. These extended models assume two-source (the crop and the substrate soil) schemes and balance the energy exchanged between soil, canopy, airspace between soil and canopy, and the atmosphere above the canopy. It has been found that two-source models better predict PET than “big leaf” models [13, 14]. Shuttleworth and Wallace [10] and Choudhury and Monteith (1988) have used a resistance network and estimated PET as the summation of transpiration from vegetation and evaporation from substrate soil. Mo et al. [12] have modified these two models by incorporating evaporation from intercepted storage.In addition, the usage of satellite data and remote sensing techniques were used in the estimation of potential evapotranspiration [15–18]. These techniques were highly useful for remote areas with meteorological data scarcities. As many of the models were originally developed for specific regions, the empirical relationships between evaporation and influencing factors may not be necessarily the same for other regions [19]. Further, not all meteorological data required by models have often been measured at all meteorological stations [1, 19]. Under these conditions, it is of considerable interest to evaluate PET models for their applicability and reliability in different regions and climatic zones.Although many such attempts have been taken worldwide at different regional scales [1, 6, 12, 20], comprehensive studies aiming to evaluate or develop PET models are still scant in Sri Lanka. Sri Lanka is characterized by high temperatures, high humidity, and unevenly distributed (both temporally and spatially) rainfall (Imbulana et al., [21]) and thus, considered mild to semiarid. Earlier studies have shown ET has accounted for 30% to more than 60% loss of total rainfall in semiarid regions [22]. Hence, ET can be considered the major water loss pathway in Sri Lanka’s water budget. This emphasizes the need for accurate quantification of ET, especially in the dry zone, the semiarid region of Sri Lanka, for sustainable water management and efficient irrigation. In addition, the meteorological data scarcity (other than rainfall) in the country is high and the available data are expensive. Therefore, a high necessity can be identified to develop accurate models to estimate ET.Both P-M and S-W models have been widely used [20]. However, considering the limitations in the P-M model, this study employs the S-W model, which has not been applied in the Sri Lankan context. While understanding the research gap, this paper aims to (1) apply the S-W model for simulation of evapotranspiration over the entire country and (2) illustrate the temporal and spatial variations of evapotranspiration over the entire country.


## 2. Study Area

Sri Lanka (refer to Figure1) is a tropical island in the Indian Ocean, located between 5°N and 10°N latitudes and 79°E to 82°E longitudes, with a total geographical area of 65,610 km2, comprising 62,705 km2 area of land and 2,905 km2 area of water. These water bodies comprise 103 distinct natural river basins and an extensive network of tanks and reservoirs (about 13,000). Approximately two-thirds of the country's landmass is low lands with elevations less than 100 m above the mean sea level. Highlands, elevations varying from 100 m to 2500 m approximately (highest mountain peak 2525 m), lie in the country's central part (Imbulana et al., [21].Figure 1
Map of Sri Lanka indicating three major climatic zones and topographic zonation (created by the authors).The only precipitation method, rainfall, has an unequal spatial and temporal distribution with a mean annual rainfall of 1861 mm, while the rainfall distribution is governed by the two major monsoon seasons: southwest monsoon (SWM) from March to September and northeast monsoon (NEM) from December to February. The country is divided into three major climatic zones, i.e., the wet zone, intermediate zone, and dry zone, based on the rainfall received and distribution. The wet zone is separated by the 2000 m annual average rainfall isohyet. In the Wet zone, rainfall ranges from 2000 mm to over 5000 mm, with an annual average rainfall of about 2,400 mm. But in the dry zone, the annual average rainfall is about 1450 mm with a minimum of lower than 1000 mm. In addition to two major monsoons and intermonsoon rains, tropical depressions that originate in the Bay of Bengal frequently enter Sri Lanka resulting in extreme rainfall events which sometimes may exceed 500 mm/day. Rainfall by all the methods counts for a mean rainfall of 1861 mm over the country per annum. Nearly 35%–45% of annual rainfall contributes to annual surface runoff. However, in most dry zone river basins, the runoff percentage is less than 35%, with the rest of the rainfall lost as evaporation and groundwater recharge.Mean annual temperature in lowlands and highlands varies between 26.5 and 28.5oC and 14.7–17.1 oC, respectively. Pan evaporation values show considerable temporal and spatial variations—varying between 1900 and 795 mm/year, with higher values recorded in the hotter dry zone. The climate of the country is characterized by high relative humidity, generally ranging between 75% and 95%. (Imbulana et al., [21]).As it was stated in the introduction, the measured meteorological data are expensive in Sri Lanka and also there is a scarcity of measured data. One of the major limitations was obtaining solar radiation data, as solar radiation is not measured at all 22 main meteorological stations that the Department of Meteorology of Sri Lanka is maintaining across the country. Therefore, only 14 stations were selected (refer to Figure1) as they mostly cover all climatic (wet, intermediate, and dry) and topographic (hilly areas and lowlands) zones of the country.


## 3. Methodology

### 3.1. Evapotranspiration Model

#### 3.1.1. Basic Governing Equations to the Evapotranspiration Model

In S-W model, total PET is computed as the summation of two major evapotranspiration components: soil evapotranspiration and transpiration from the dry canopy [10]:(1)Et=Es+Ec,where Et is the total PET (mm), Es is the evaporation from soil (mm), and Ec is the transpiration from the dry canopy (mm). Ec and Es can be expressed as follows [10]:(2)Ec=1λ⋅ΔRnc+ρCpD0/racΔ+γ1+rc/rac,(3)Ec=1λ⋅ΔRns−G+ρCpD0/rasΔ+γ1+rs/ras,where λ is the latent heat of vaporization (MJ kg −1), Δ is the slope of saturation vapour pressure curve (kPa oC −1), Rnc is the net radiation absorbed by the canopy (MJ m −2), ρ is the air density (kg m−3), Cp is the air specific heat at constant pressure (=1.013 × 10−3 MJ kg −1oC −1), D0 is the water vapour deficit at the canopy height (kPa), rac is the bulk boundary-layer resistance of the canopy (s m−1), γ is the psychrometric constant (kPa oC −1), rc is the canopy resistance (s m−1), Rns is the net radiation at the substrate surface (MJ m −2), G is the soil heat flux (MJ m −2), ras is the aerodynamic resistance between the soil surface and canopy air space (s m−1), and rs is the soil resistance (s m−1). D0 can be expressed as follows [10]:(4)D0=D+ΔRn−G−Δ+γλEraρCp,where D is the water vapor deficit at the reference height (kPa), Rn is the net incoming radiation (kPa), and ra is the aerodynamic resistance between canopy source and reference level (s m−1). By substituting equations (2)–(4) in equation (1) and manipulating, Shuttleworth and Wallace [10] have derived the following equation for the PET model.(5)λEt=CcPMc+CsPMs.More explanations of these parameters are given in Appendix equations (A.1)–(A.7). Chow et al. [23] have also expressed some of the parameters, and these are given in Appendix equations (A.8)–(A.12).


#### 3.1.2. Estimation of Net Radiation (Rn)

Estimation of radiation is important. The net radiationRn received at the Earth surface can be subdivided (refer to equation (6)) into the net radiation absorbed by the canopy Rnc and the net radiation absorbed by the soil Rns [12].(6)Rn=Rnc+Rns.Shuttleworth and Wallace [10] have expressed the relationship of Rn and Rns as in equation (7), where Cr is the extinction coefficient of the vegetation for net radiation. All radiation terms are in MJ m−2. Cr was taken as 0.5 [12, 20]. Cr was taken as 0.7 [10].(7)Rns=Rnexp−CrLAI.LAI is the Leaf Area Index, a dimensionless parameter that characterises vegetation cover. There are various definitions of LAI; hence, different sets of equations can be found in the literature to estimate LAI. In our study, we utilized the method followed by Zhou et al. [20], which is given in Appendix equations (A.13)–(A.15). In addition, the maximum LAI values for different vegetation types are given in Table 1. Numbers 1–3 are for the tall vegetation, while the others are for the shorter vegetation. The complete table is given in the Appendix as Table 2.Table 1
Maximum LAI values for different vegetation types [20].
CodeLand use typeClassification as per the literatureLAImax1CoconutEvergreen needle leaf forests5.52RubberEvergreen broadleaf forests73Forest, unclassifiedMixed forests5.74Homesteads/gardenOpen shrub lands35Shrublands36Tea37GrasslandsGrasslands1.88Marshy landsPermanent wetlands69ChenaCroplands710Other cultivations711Paddy712Urban and built-upUrban and built-up013Barren landBarren or sparsely vegetated0.314Water bodiesWater bodies0Table 2
Comparison of averaged annual PET (simulated), pan evaporation, and rainfall at the fourteen selected stations.
StationCell numberClimate zonePET (mm/yr)Pan evaporation (mm/yr)Rainfall (mm/yr)PET/rainfall (%)PanEvap./rainfall (%)Jaffna15Dry zone1430.161414.451297.77110.2109.0Vavuniya2381320.011252.801555.8584.880.5Anuradhapura4611506.761250.361624.7692.777.0Puttalam6021516.561503.501224.44123.9122.8Polonnaruwa6901671.671488.141784.1593.783.4Hamabantota17861585.071524.471075.55147.4141.7Kurunegala954Intermediate zone1552.391416.712075.2774.868.3Badulla1303955.34866.581896.4650.445.7Bandarawela14121095.88982.441650.8066.459.5Katugastota1069Wet zone1233.441143.911934.2263.859.1Nuwara-eliya1298968.67889.841846.5852.548.2Colombo13191041.391319.972493.0441.852.9Rathnapura14731025.53921.474072.7025.222.6Galle18091058.68963.872450.2643.239.3


#### 3.1.3. Estimation of Water Vapour Deficit (D)

Estimation of water vapor deficit at the reference height (D) can be found in the following equation [20]:(8)D=es−ea,where es and ea are the saturation vapor pressure (kPa) and the ambient vapor pressure (kPa), respectively, and the expressions for es and ea are given in Appendix equations (A.16) and (A.17) [23].


#### 3.1.4. Estimation of Aerodynamic Resistance between Canopy Source and Reference Level (ra)

The aerodynamic resistance between the canopy source and reference level was calculated using the following equation (9), which was adapted by Shuttleworth and Gurney [24]:(9)ra=1ku∗lnza−d0hc−d0+hcηKhexpη1−Z0+dphc−1,where k is the von Karman’s constant (k = 0.41), u∗ is the friction velocity (m s−1), za is the reference height (m), d0 is the zero-plane displacement of the canopy (m), hc is the canopy height (m) (refer to Table 2), η is the eddy diffusivity decay constant of the vegetation, Kh is the eddy diffusion coefficient at the top of the canopy (m2 s−1), Z0 is the ‘‘preferred’’ roughness length (m), and dp is the ‘‘preferred’’ zero plane displacement. The equations used to compute each term are given in Appendix equations (A.18)–(A.27).


#### 3.1.5. Estimation of Bulk Boundary-Layer Resistance of the Canopy (rac)

Shuttleworth and Gurney [24] have estimated bulk boundary layer resistance of the canopy by assuming that energy transfer only occurs by molecular diffusion through a laminar layer around leaves and using the following equation, which was used in our model:(10)rac=100η⋅luh⋅1−exp−η2−1⋅12LAI,where l is the canopy characteristic leaf width (m), and uh is the wind speed at the top of the canopy (m s−1). uh was computed using equation (11) [23], where u∗ is the shear velocity (m s−1), which was calculated assigning wind speed values recorded at the reference height and reference height, respectively, for uh and hc. Then uh was calculated, assigning the respective u∗ values. l was calculated using equation (12), where lmax is the maximum leaf width (m) (refer to Table 2).(11)uhu∗=1klnhcz0,(12)l=lmax,for perennial vegetation,lmax1−exp−0.6LAI,for annual vegetation.


#### 3.1.6. Estimation of Aerodynamic Resistance between the Soil Surface and Canopy Air Space (ras)

Estimation of the aerodynamic resistance between the soil surface and canopy air space requires complex formulations asras is affected by many factors. ras can be calculated using (notations were previously defined)(13)ras=hcexpηηKhexp−ηzoghc−exp−ηZ0+dphc.


#### 3.1.7. Estimation of Canopy Resistance (rc)

The canopy resistance was estimated using the Jarvis canopy resistance model [20, 25]. Jarvis model is given by the following equation:(14)rc=rstminLAIefRnfTkfDfθ,where rstmin is the minimum stomatal resistance (s m−1) (refer to Table 2) and LAIe is the effective LAI. Expressions for LAIe, fRn, fD, fTk, and fθ are given in Appendix equations (A.28)–(A.32).


#### 3.1.8. Estimation of Soil Resistance (rs)

Due to the complexities in acquiring data required for accurate assessments of the soil resistance, it was setrs=0sm−1 at saturation point and rs=200sm−1 at wilting point, and this was suggested by Shuttleworth and Wallace [10]. Then interpolate between two extremes to estimate rs at average soil moisture content.


#### 3.1.9. Estimation of Soil Heat Flux (G)

Different methods are available in the literature to estimateG. Many researchers consider G is 30% of the Rn [20]; however, Mo et al. [12] have suggested the following equation, which was used in this study:(15)G=0.183Rnexp−0.299LAI.


### 3.2. Data Collection and Preprocessing

#### 3.2.1. Topographic Data

Elevation data are required to calculate the atmospheric pressure (equation (A.12)), which is then used to compute ρ (air density-equation (A.10)). Elevation values were extracted from a digital elevation model (DEM) with a spatial resolution of 30m×30m, which was obtained from the Department of Survey, Sri Lanka. The DEM resolution was adjusted in ArcGIS to match the model grid size (6km×6km) of Figure 2(a).Figure 2
Model input parameters (created by the authors). (a) For elevation. (b) For land-use type. (c) For NDVI. (d) For soil type. (e) For soil moisture. (f) For mean temperature. (g) For wind speed. (h) For average daily net solar radiation. (i) For average relative humidity. (j) For monthly average rainfall.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)


#### 3.2.2. Land Use and Vegetation Data

Land cover data were obtained from the Department of Survey, Sri Lanka. The raw data obtained were in vector format. The raw data were then converted to raster format with a resolution of6km×6km. There are altogether 46 land use classes in the original vector file. They were recategorized into ten classes, as shown in Table 1 and Figure 2(b). Further, the vegetation cover was subcategorised as short-tall and annual-perennial (Table 1). Collecting the land cover threshold parameters required by the model, i.e., LAImax, lmax, hc,Fcl, rstmin, NDVI98%, and zog was the greatest challenge, as estimation of said parameters was extremely difficult. Therefore, aforesaid parameters were gathered from Zhou et al. [20] and are given in Table 1.


#### 3.2.3. NDVI

Remotely sensed Landsat images from 2009 to 2019, with a spatial resolution of30m×30m, were acquired from USGS [26]. A total of 9 Landsat tiles were required to cover entire Sri Lanka. However, due to the high dense cloud cover (>5%), Landsat images from a previous or a subsequent month were utilized for some months. This may result in some errors in the actual condition, thus considered a potential limitation of the study. The NDVI values were computed as the difference between near-infrared (NIR) and red (RED) reflectance divided by their sum (equation (16)) in ArcGIS. Finally, the calculated NDVI layers’ resolution was converted to 6km×6km (Figure 2(c)).(16)NDVI=NIR−RedNIR+Red.


#### 3.2.4. Soil Type and Soil Moisture

Sri Lanka soil type map was obtained from the Department of Survey, Sri Lanka. The soil map of Sri Lanka was then recategorized into six soil textural classes based on the classifications of Moormann and Panabokke [27] (Figure 2(d)). Much of the dry zone consists of clay loam, and most of the wet zone contains sandy loam. Soil saturation and wilting point values were assessed based on the soil textural classes. The reclassified soil vector layer was converted to a raster layer with a 6 km × 6 km grid.Due to the limitation imposed by the unavailability of field measured soil moisture data in the country, monthly root zone soil moisture values were retrieved from NASA [28] from 2009 to 2019 (Figure 2(e)). When data for a month was missing, data from the previous or subsequent month (depending on the rainfall seasons) were utilized. The study may face limitations incurred by this alternative method. Retrieved data has a spatial resolution of 30m×30m; thus, it was transformed to the model’s spatial resolution. Incur.


#### 3.2.5. Meteorological Data

Monthly time series data of maximum and minimum daily temperatures, relative humidity (day and night), daily precipitation, solar net radiation, wind speed, and evaporation from 2009 to 2019 were obtained from the Department of Meteorology of Sri Lanka. The model requires only the temperature, solar radiation, humidity, and wind speed data (Figures2(f)–2(j)). Average values of maximum and minimum temperatures and day and night relative humidity data were fed into the model. Precipitation and evaporation data were used for trend analysis and model validation. The data was preprocessed in ArcGIS. The spatial distributions of monthly averaged temperature, relative humidity, and evaporation were generated by Inverse Distance Weighted (IDW) interpolation. The Universal Kriging interpolation method was used to obtain the spatial distribution of solar net radiation and wind speed. The spatial distributions of monthly precipitation were generated by applying the Thiessen polygon method.


### 3.3. Overall Methodology

Processed data showcased in Figures2(a)–2(j) were used to simulate the potential evapotranspiration for the whole of Sri Lanka. All the equations were modeled in a Microsoft Excel office package and then extracted to ArcGIS to develop the graphical presentation of potential evapotranspiration over the country. The results were validated by a comparison analysis using recorded pan evapotranspiration and the predicted potential evapotranspiration. RMSE was calculated between the predictions and recorded values.


## 3.1. Evapotranspiration Model

### 3.1.1. Basic Governing Equations to the Evapotranspiration Model

In S-W model, total PET is computed as the summation of two major evapotranspiration components: soil evapotranspiration and transpiration from the dry canopy [10]:(1)Et=Es+Ec,where Et is the total PET (mm), Es is the evaporation from soil (mm), and Ec is the transpiration from the dry canopy (mm). Ec and Es can be expressed as follows [10]:(2)Ec=1λ⋅ΔRnc+ρCpD0/racΔ+γ1+rc/rac,(3)Ec=1λ⋅ΔRns−G+ρCpD0/rasΔ+γ1+rs/ras,where λ is the latent heat of vaporization (MJ kg −1), Δ is the slope of saturation vapour pressure curve (kPa oC −1), Rnc is the net radiation absorbed by the canopy (MJ m −2), ρ is the air density (kg m−3), Cp is the air specific heat at constant pressure (=1.013 × 10−3 MJ kg −1oC −1), D0 is the water vapour deficit at the canopy height (kPa), rac is the bulk boundary-layer resistance of the canopy (s m−1), γ is the psychrometric constant (kPa oC −1), rc is the canopy resistance (s m−1), Rns is the net radiation at the substrate surface (MJ m −2), G is the soil heat flux (MJ m −2), ras is the aerodynamic resistance between the soil surface and canopy air space (s m−1), and rs is the soil resistance (s m−1). D0 can be expressed as follows [10]:(4)D0=D+ΔRn−G−Δ+γλEraρCp,where D is the water vapor deficit at the reference height (kPa), Rn is the net incoming radiation (kPa), and ra is the aerodynamic resistance between canopy source and reference level (s m−1). By substituting equations (2)–(4) in equation (1) and manipulating, Shuttleworth and Wallace [10] have derived the following equation for the PET model.(5)λEt=CcPMc+CsPMs.More explanations of these parameters are given in Appendix equations (A.1)–(A.7). Chow et al. [23] have also expressed some of the parameters, and these are given in Appendix equations (A.8)–(A.12).


### 3.1.2. Estimation of Net Radiation (Rn)

Estimation of radiation is important. The net radiationRn received at the Earth surface can be subdivided (refer to equation (6)) into the net radiation absorbed by the canopy Rnc and the net radiation absorbed by the soil Rns [12].(6)Rn=Rnc+Rns.Shuttleworth and Wallace [10] have expressed the relationship of Rn and Rns as in equation (7), where Cr is the extinction coefficient of the vegetation for net radiation. All radiation terms are in MJ m−2. Cr was taken as 0.5 [12, 20]. Cr was taken as 0.7 [10].(7)Rns=Rnexp−CrLAI.LAI is the Leaf Area Index, a dimensionless parameter that characterises vegetation cover. There are various definitions of LAI; hence, different sets of equations can be found in the literature to estimate LAI. In our study, we utilized the method followed by Zhou et al. [20], which is given in Appendix equations (A.13)–(A.15). In addition, the maximum LAI values for different vegetation types are given in Table 1. Numbers 1–3 are for the tall vegetation, while the others are for the shorter vegetation. The complete table is given in the Appendix as Table 2.Table 1
Maximum LAI values for different vegetation types [20].
CodeLand use typeClassification as per the literatureLAImax1CoconutEvergreen needle leaf forests5.52RubberEvergreen broadleaf forests73Forest, unclassifiedMixed forests5.74Homesteads/gardenOpen shrub lands35Shrublands36Tea37GrasslandsGrasslands1.88Marshy landsPermanent wetlands69ChenaCroplands710Other cultivations711Paddy712Urban and built-upUrban and built-up013Barren landBarren or sparsely vegetated0.314Water bodiesWater bodies0Table 2
Comparison of averaged annual PET (simulated), pan evaporation, and rainfall at the fourteen selected stations.
StationCell numberClimate zonePET (mm/yr)Pan evaporation (mm/yr)Rainfall (mm/yr)PET/rainfall (%)PanEvap./rainfall (%)Jaffna15Dry zone1430.161414.451297.77110.2109.0Vavuniya2381320.011252.801555.8584.880.5Anuradhapura4611506.761250.361624.7692.777.0Puttalam6021516.561503.501224.44123.9122.8Polonnaruwa6901671.671488.141784.1593.783.4Hamabantota17861585.071524.471075.55147.4141.7Kurunegala954Intermediate zone1552.391416.712075.2774.868.3Badulla1303955.34866.581896.4650.445.7Bandarawela14121095.88982.441650.8066.459.5Katugastota1069Wet zone1233.441143.911934.2263.859.1Nuwara-eliya1298968.67889.841846.5852.548.2Colombo13191041.391319.972493.0441.852.9Rathnapura14731025.53921.474072.7025.222.6Galle18091058.68963.872450.2643.239.3


### 3.1.3. Estimation of Water Vapour Deficit (D)

Estimation of water vapor deficit at the reference height (D) can be found in the following equation [20]:(8)D=es−ea,where es and ea are the saturation vapor pressure (kPa) and the ambient vapor pressure (kPa), respectively, and the expressions for es and ea are given in Appendix equations (A.16) and (A.17) [23].


### 3.1.4. Estimation of Aerodynamic Resistance between Canopy Source and Reference Level (ra)

The aerodynamic resistance between the canopy source and reference level was calculated using the following equation (9), which was adapted by Shuttleworth and Gurney [24]:(9)ra=1ku∗lnza−d0hc−d0+hcηKhexpη1−Z0+dphc−1,where k is the von Karman’s constant (k = 0.41), u∗ is the friction velocity (m s−1), za is the reference height (m), d0 is the zero-plane displacement of the canopy (m), hc is the canopy height (m) (refer to Table 2), η is the eddy diffusivity decay constant of the vegetation, Kh is the eddy diffusion coefficient at the top of the canopy (m2 s−1), Z0 is the ‘‘preferred’’ roughness length (m), and dp is the ‘‘preferred’’ zero plane displacement. The equations used to compute each term are given in Appendix equations (A.18)–(A.27).


### 3.1.5. Estimation of Bulk Boundary-Layer Resistance of the Canopy (rac)

Shuttleworth and Gurney [24] have estimated bulk boundary layer resistance of the canopy by assuming that energy transfer only occurs by molecular diffusion through a laminar layer around leaves and using the following equation, which was used in our model:(10)rac=100η⋅luh⋅1−exp−η2−1⋅12LAI,where l is the canopy characteristic leaf width (m), and uh is the wind speed at the top of the canopy (m s−1). uh was computed using equation (11) [23], where u∗ is the shear velocity (m s−1), which was calculated assigning wind speed values recorded at the reference height and reference height, respectively, for uh and hc. Then uh was calculated, assigning the respective u∗ values. l was calculated using equation (12), where lmax is the maximum leaf width (m) (refer to Table 2).(11)uhu∗=1klnhcz0,(12)l=lmax,for perennial vegetation,lmax1−exp−0.6LAI,for annual vegetation.


### 3.1.6. Estimation of Aerodynamic Resistance between the Soil Surface and Canopy Air Space (ras)

Estimation of the aerodynamic resistance between the soil surface and canopy air space requires complex formulations asras is affected by many factors. ras can be calculated using (notations were previously defined)(13)ras=hcexpηηKhexp−ηzoghc−exp−ηZ0+dphc.


### 3.1.7. Estimation of Canopy Resistance (rc)

The canopy resistance was estimated using the Jarvis canopy resistance model [20, 25]. Jarvis model is given by the following equation:(14)rc=rstminLAIefRnfTkfDfθ,where rstmin is the minimum stomatal resistance (s m−1) (refer to Table 2) and LAIe is the effective LAI. Expressions for LAIe, fRn, fD, fTk, and fθ are given in Appendix equations (A.28)–(A.32).


### 3.1.8. Estimation of Soil Resistance (rs)

Due to the complexities in acquiring data required for accurate assessments of the soil resistance, it was setrs=0sm−1 at saturation point and rs=200sm−1 at wilting point, and this was suggested by Shuttleworth and Wallace [10]. Then interpolate between two extremes to estimate rs at average soil moisture content.


### 3.1.9. Estimation of Soil Heat Flux (G)

Different methods are available in the literature to estimateG. Many researchers consider G is 30% of the Rn [20]; however, Mo et al. [12] have suggested the following equation, which was used in this study:(15)G=0.183Rnexp−0.299LAI.


## 3.1.1. Basic Governing Equations to the Evapotranspiration Model

In S-W model, total PET is computed as the summation of two major evapotranspiration components: soil evapotranspiration and transpiration from the dry canopy [10]:(1)Et=Es+Ec,where Et is the total PET (mm), Es is the evaporation from soil (mm), and Ec is the transpiration from the dry canopy (mm). Ec and Es can be expressed as follows [10]:(2)Ec=1λ⋅ΔRnc+ρCpD0/racΔ+γ1+rc/rac,(3)Ec=1λ⋅ΔRns−G+ρCpD0/rasΔ+γ1+rs/ras,where λ is the latent heat of vaporization (MJ kg −1), Δ is the slope of saturation vapour pressure curve (kPa oC −1), Rnc is the net radiation absorbed by the canopy (MJ m −2), ρ is the air density (kg m−3), Cp is the air specific heat at constant pressure (=1.013 × 10−3 MJ kg −1oC −1), D0 is the water vapour deficit at the canopy height (kPa), rac is the bulk boundary-layer resistance of the canopy (s m−1), γ is the psychrometric constant (kPa oC −1), rc is the canopy resistance (s m−1), Rns is the net radiation at the substrate surface (MJ m −2), G is the soil heat flux (MJ m −2), ras is the aerodynamic resistance between the soil surface and canopy air space (s m−1), and rs is the soil resistance (s m−1). D0 can be expressed as follows [10]:(4)D0=D+ΔRn−G−Δ+γλEraρCp,where D is the water vapor deficit at the reference height (kPa), Rn is the net incoming radiation (kPa), and ra is the aerodynamic resistance between canopy source and reference level (s m−1). By substituting equations (2)–(4) in equation (1) and manipulating, Shuttleworth and Wallace [10] have derived the following equation for the PET model.(5)λEt=CcPMc+CsPMs.More explanations of these parameters are given in Appendix equations (A.1)–(A.7). Chow et al. [23] have also expressed some of the parameters, and these are given in Appendix equations (A.8)–(A.12).


## 3.1.2. Estimation of Net Radiation (Rn)

Estimation of radiation is important. The net radiationRn received at the Earth surface can be subdivided (refer to equation (6)) into the net radiation absorbed by the canopy Rnc and the net radiation absorbed by the soil Rns [12].(6)Rn=Rnc+Rns.Shuttleworth and Wallace [10] have expressed the relationship of Rn and Rns as in equation (7), where Cr is the extinction coefficient of the vegetation for net radiation. All radiation terms are in MJ m−2. Cr was taken as 0.5 [12, 20]. Cr was taken as 0.7 [10].(7)Rns=Rnexp−CrLAI.LAI is the Leaf Area Index, a dimensionless parameter that characterises vegetation cover. There are various definitions of LAI; hence, different sets of equations can be found in the literature to estimate LAI. In our study, we utilized the method followed by Zhou et al. [20], which is given in Appendix equations (A.13)–(A.15). In addition, the maximum LAI values for different vegetation types are given in Table 1. Numbers 1–3 are for the tall vegetation, while the others are for the shorter vegetation. The complete table is given in the Appendix as Table 2.Table 1
Maximum LAI values for different vegetation types [20].
CodeLand use typeClassification as per the literatureLAImax1CoconutEvergreen needle leaf forests5.52RubberEvergreen broadleaf forests73Forest, unclassifiedMixed forests5.74Homesteads/gardenOpen shrub lands35Shrublands36Tea37GrasslandsGrasslands1.88Marshy landsPermanent wetlands69ChenaCroplands710Other cultivations711Paddy712Urban and built-upUrban and built-up013Barren landBarren or sparsely vegetated0.314Water bodiesWater bodies0Table 2
Comparison of averaged annual PET (simulated), pan evaporation, and rainfall at the fourteen selected stations.
StationCell numberClimate zonePET (mm/yr)Pan evaporation (mm/yr)Rainfall (mm/yr)PET/rainfall (%)PanEvap./rainfall (%)Jaffna15Dry zone1430.161414.451297.77110.2109.0Vavuniya2381320.011252.801555.8584.880.5Anuradhapura4611506.761250.361624.7692.777.0Puttalam6021516.561503.501224.44123.9122.8Polonnaruwa6901671.671488.141784.1593.783.4Hamabantota17861585.071524.471075.55147.4141.7Kurunegala954Intermediate zone1552.391416.712075.2774.868.3Badulla1303955.34866.581896.4650.445.7Bandarawela14121095.88982.441650.8066.459.5Katugastota1069Wet zone1233.441143.911934.2263.859.1Nuwara-eliya1298968.67889.841846.5852.548.2Colombo13191041.391319.972493.0441.852.9Rathnapura14731025.53921.474072.7025.222.6Galle18091058.68963.872450.2643.239.3


## 3.1.3. Estimation of Water Vapour Deficit (D)

Estimation of water vapor deficit at the reference height (D) can be found in the following equation [20]:(8)D=es−ea,where es and ea are the saturation vapor pressure (kPa) and the ambient vapor pressure (kPa), respectively, and the expressions for es and ea are given in Appendix equations (A.16) and (A.17) [23].


## 3.1.4. Estimation of Aerodynamic Resistance between Canopy Source and Reference Level (ra)

The aerodynamic resistance between the canopy source and reference level was calculated using the following equation (9), which was adapted by Shuttleworth and Gurney [24]:(9)ra=1ku∗lnza−d0hc−d0+hcηKhexpη1−Z0+dphc−1,where k is the von Karman’s constant (k = 0.41), u∗ is the friction velocity (m s−1), za is the reference height (m), d0 is the zero-plane displacement of the canopy (m), hc is the canopy height (m) (refer to Table 2), η is the eddy diffusivity decay constant of the vegetation, Kh is the eddy diffusion coefficient at the top of the canopy (m2 s−1), Z0 is the ‘‘preferred’’ roughness length (m), and dp is the ‘‘preferred’’ zero plane displacement. The equations used to compute each term are given in Appendix equations (A.18)–(A.27).


## 3.1.5. Estimation of Bulk Boundary-Layer Resistance of the Canopy (rac)

Shuttleworth and Gurney [24] have estimated bulk boundary layer resistance of the canopy by assuming that energy transfer only occurs by molecular diffusion through a laminar layer around leaves and using the following equation, which was used in our model:(10)rac=100η⋅luh⋅1−exp−η2−1⋅12LAI,where l is the canopy characteristic leaf width (m), and uh is the wind speed at the top of the canopy (m s−1). uh was computed using equation (11) [23], where u∗ is the shear velocity (m s−1), which was calculated assigning wind speed values recorded at the reference height and reference height, respectively, for uh and hc. Then uh was calculated, assigning the respective u∗ values. l was calculated using equation (12), where lmax is the maximum leaf width (m) (refer to Table 2).(11)uhu∗=1klnhcz0,(12)l=lmax,for perennial vegetation,lmax1−exp−0.6LAI,for annual vegetation.


## 3.1.6. Estimation of Aerodynamic Resistance between the Soil Surface and Canopy Air Space (ras)

Estimation of the aerodynamic resistance between the soil surface and canopy air space requires complex formulations asras is affected by many factors. ras can be calculated using (notations were previously defined)(13)ras=hcexpηηKhexp−ηzoghc−exp−ηZ0+dphc.


## 3.1.7. Estimation of Canopy Resistance (rc)

The canopy resistance was estimated using the Jarvis canopy resistance model [20, 25]. Jarvis model is given by the following equation:(14)rc=rstminLAIefRnfTkfDfθ,where rstmin is the minimum stomatal resistance (s m−1) (refer to Table 2) and LAIe is the effective LAI. Expressions for LAIe, fRn, fD, fTk, and fθ are given in Appendix equations (A.28)–(A.32).


## 3.1.8. Estimation of Soil Resistance (rs)

Due to the complexities in acquiring data required for accurate assessments of the soil resistance, it was setrs=0sm−1 at saturation point and rs=200sm−1 at wilting point, and this was suggested by Shuttleworth and Wallace [10]. Then interpolate between two extremes to estimate rs at average soil moisture content.


## 3.1.9. Estimation of Soil Heat Flux (G)

Different methods are available in the literature to estimateG. Many researchers consider G is 30% of the Rn [20]; however, Mo et al. [12] have suggested the following equation, which was used in this study:(15)G=0.183Rnexp−0.299LAI.


## 3.2. Data Collection and Preprocessing

### 3.2.1. Topographic Data

Elevation data are required to calculate the atmospheric pressure (equation (A.12)), which is then used to compute ρ (air density-equation (A.10)). Elevation values were extracted from a digital elevation model (DEM) with a spatial resolution of 30m×30m, which was obtained from the Department of Survey, Sri Lanka. The DEM resolution was adjusted in ArcGIS to match the model grid size (6km×6km) of Figure 2(a).Figure 2
Model input parameters (created by the authors). (a) For elevation. (b) For land-use type. (c) For NDVI. (d) For soil type. (e) For soil moisture. (f) For mean temperature. (g) For wind speed. (h) For average daily net solar radiation. (i) For average relative humidity. (j) For monthly average rainfall.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)


### 3.2.2. Land Use and Vegetation Data

Land cover data were obtained from the Department of Survey, Sri Lanka. The raw data obtained were in vector format. The raw data were then converted to raster format with a resolution of6km×6km. There are altogether 46 land use classes in the original vector file. They were recategorized into ten classes, as shown in Table 1 and Figure 2(b). Further, the vegetation cover was subcategorised as short-tall and annual-perennial (Table 1). Collecting the land cover threshold parameters required by the model, i.e., LAImax, lmax, hc,Fcl, rstmin, NDVI98%, and zog was the greatest challenge, as estimation of said parameters was extremely difficult. Therefore, aforesaid parameters were gathered from Zhou et al. [20] and are given in Table 1.


### 3.2.3. NDVI

Remotely sensed Landsat images from 2009 to 2019, with a spatial resolution of30m×30m, were acquired from USGS [26]. A total of 9 Landsat tiles were required to cover entire Sri Lanka. However, due to the high dense cloud cover (>5%), Landsat images from a previous or a subsequent month were utilized for some months. This may result in some errors in the actual condition, thus considered a potential limitation of the study. The NDVI values were computed as the difference between near-infrared (NIR) and red (RED) reflectance divided by their sum (equation (16)) in ArcGIS. Finally, the calculated NDVI layers’ resolution was converted to 6km×6km (Figure 2(c)).(16)NDVI=NIR−RedNIR+Red.


### 3.2.4. Soil Type and Soil Moisture

Sri Lanka soil type map was obtained from the Department of Survey, Sri Lanka. The soil map of Sri Lanka was then recategorized into six soil textural classes based on the classifications of Moormann and Panabokke [27] (Figure 2(d)). Much of the dry zone consists of clay loam, and most of the wet zone contains sandy loam. Soil saturation and wilting point values were assessed based on the soil textural classes. The reclassified soil vector layer was converted to a raster layer with a 6 km × 6 km grid.Due to the limitation imposed by the unavailability of field measured soil moisture data in the country, monthly root zone soil moisture values were retrieved from NASA [28] from 2009 to 2019 (Figure 2(e)). When data for a month was missing, data from the previous or subsequent month (depending on the rainfall seasons) were utilized. The study may face limitations incurred by this alternative method. Retrieved data has a spatial resolution of 30m×30m; thus, it was transformed to the model’s spatial resolution. Incur.


### 3.2.5. Meteorological Data

Monthly time series data of maximum and minimum daily temperatures, relative humidity (day and night), daily precipitation, solar net radiation, wind speed, and evaporation from 2009 to 2019 were obtained from the Department of Meteorology of Sri Lanka. The model requires only the temperature, solar radiation, humidity, and wind speed data (Figures2(f)–2(j)). Average values of maximum and minimum temperatures and day and night relative humidity data were fed into the model. Precipitation and evaporation data were used for trend analysis and model validation. The data was preprocessed in ArcGIS. The spatial distributions of monthly averaged temperature, relative humidity, and evaporation were generated by Inverse Distance Weighted (IDW) interpolation. The Universal Kriging interpolation method was used to obtain the spatial distribution of solar net radiation and wind speed. The spatial distributions of monthly precipitation were generated by applying the Thiessen polygon method.


## 3.2.1. Topographic Data

Elevation data are required to calculate the atmospheric pressure (equation (A.12)), which is then used to compute ρ (air density-equation (A.10)). Elevation values were extracted from a digital elevation model (DEM) with a spatial resolution of 30m×30m, which was obtained from the Department of Survey, Sri Lanka. The DEM resolution was adjusted in ArcGIS to match the model grid size (6km×6km) of Figure 2(a).Figure 2
Model input parameters (created by the authors). (a) For elevation. (b) For land-use type. (c) For NDVI. (d) For soil type. (e) For soil moisture. (f) For mean temperature. (g) For wind speed. (h) For average daily net solar radiation. (i) For average relative humidity. (j) For monthly average rainfall.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)


## 3.2.2. Land Use and Vegetation Data

Land cover data were obtained from the Department of Survey, Sri Lanka. The raw data obtained were in vector format. The raw data were then converted to raster format with a resolution of6km×6km. There are altogether 46 land use classes in the original vector file. They were recategorized into ten classes, as shown in Table 1 and Figure 2(b). Further, the vegetation cover was subcategorised as short-tall and annual-perennial (Table 1). Collecting the land cover threshold parameters required by the model, i.e., LAImax, lmax, hc,Fcl, rstmin, NDVI98%, and zog was the greatest challenge, as estimation of said parameters was extremely difficult. Therefore, aforesaid parameters were gathered from Zhou et al. [20] and are given in Table 1.


## 3.2.3. NDVI

Remotely sensed Landsat images from 2009 to 2019, with a spatial resolution of30m×30m, were acquired from USGS [26]. A total of 9 Landsat tiles were required to cover entire Sri Lanka. However, due to the high dense cloud cover (>5%), Landsat images from a previous or a subsequent month were utilized for some months. This may result in some errors in the actual condition, thus considered a potential limitation of the study. The NDVI values were computed as the difference between near-infrared (NIR) and red (RED) reflectance divided by their sum (equation (16)) in ArcGIS. Finally, the calculated NDVI layers’ resolution was converted to 6km×6km (Figure 2(c)).(16)NDVI=NIR−RedNIR+Red.


## 3.2.4. Soil Type and Soil Moisture

Sri Lanka soil type map was obtained from the Department of Survey, Sri Lanka. The soil map of Sri Lanka was then recategorized into six soil textural classes based on the classifications of Moormann and Panabokke [27] (Figure 2(d)). Much of the dry zone consists of clay loam, and most of the wet zone contains sandy loam. Soil saturation and wilting point values were assessed based on the soil textural classes. The reclassified soil vector layer was converted to a raster layer with a 6 km × 6 km grid.Due to the limitation imposed by the unavailability of field measured soil moisture data in the country, monthly root zone soil moisture values were retrieved from NASA [28] from 2009 to 2019 (Figure 2(e)). When data for a month was missing, data from the previous or subsequent month (depending on the rainfall seasons) were utilized. The study may face limitations incurred by this alternative method. Retrieved data has a spatial resolution of 30m×30m; thus, it was transformed to the model’s spatial resolution. Incur.


## 3.2.5. Meteorological Data

Monthly time series data of maximum and minimum daily temperatures, relative humidity (day and night), daily precipitation, solar net radiation, wind speed, and evaporation from 2009 to 2019 were obtained from the Department of Meteorology of Sri Lanka. The model requires only the temperature, solar radiation, humidity, and wind speed data (Figures2(f)–2(j)). Average values of maximum and minimum temperatures and day and night relative humidity data were fed into the model. Precipitation and evaporation data were used for trend analysis and model validation. The data was preprocessed in ArcGIS. The spatial distributions of monthly averaged temperature, relative humidity, and evaporation were generated by Inverse Distance Weighted (IDW) interpolation. The Universal Kriging interpolation method was used to obtain the spatial distribution of solar net radiation and wind speed. The spatial distributions of monthly precipitation were generated by applying the Thiessen polygon method.


## 3.3. Overall Methodology

Processed data showcased in Figures2(a)–2(j) were used to simulate the potential evapotranspiration for the whole of Sri Lanka. All the equations were modeled in a Microsoft Excel office package and then extracted to ArcGIS to develop the graphical presentation of potential evapotranspiration over the country. The results were validated by a comparison analysis using recorded pan evapotranspiration and the predicted potential evapotranspiration. RMSE was calculated between the predictions and recorded values.


## 4. Results and Discussion

### 4.1. Impact of Data Scarcity on Analysis

The model results of fourteen grid cells (refer to Table3), in which the selected meteorological stations lie, compared with the recorded pan evaporation at the respective meteorological stations. They could lessen the error resulting from interpolation. The model was not calibrated in this study due to a lack of data and a shorter data period. Coupled with the limitations exerted by NDVI and soil moisture data, further challenges to model calibration were exerted by the uncertain canopy and soil parameters obtained from the literature. However, the root means square error (RMSE) of model prediction to observation, 0.75 mm/day, suggests that the model can predict PET and its component over Sri Lanka with moderate accuracy.Table 3
Predicted annual averaged PET, ETs, and ETc from different land-use types.
PET (mm/year)ETs (mm/year)ETc (mm/year)Paddy1391.65505.03886.63Rubber1292.72168.201124.52Coconut1276.18206.751069.42Chena (arid cultivation)1255.75331.32924.43Tea1199.56295.68903.87Grass lands663.07429.49233.58Shrub lands1173.70285.13888.57Home gardens872.00300.53571.47Urban and built-up areas405.46405.46Historical solar radiation measurements are not commonly available at many meteorological stations in Sri Lanka. This reduced the number of meteorological stations we could use for this study. The density of selected meteorological stations is notably lesser in the dry zone. Since there is no station in the South-East (SE) part of the country, observed values at Badulla and Bandarawela seem to influence the interpolation outputs in the SE region, irrespective of the interpolation method. As per the actual situation, in the SE low land, some meteorological parameters such as rainfall, temperature, and humidity show a drastic difference from those at Badulla and Bandarawela. Therefore, the accuracy and precision of the model predictions in the SE region are questionable. However, the model can be applied to the rest of the dry zone and intermediate zone with moderate confidence and the wet zone with high confidence.


### 4.2. Comparison of Predicted PET, Pan Evaporation, and Precipitation

The model generated PET values are greater than the observed, except at the Colombo station, at which the simulated results are significantly smaller than recorded values by an average of 32%. The disparity between simulated and observed values at Colombo can be attributed to several factors. The model is highly sensitive to LAI, and there is a positive correlation between LAI and simulated results, which has been confirmed by Zhou et al. [20]. Colombo is the capital city of Sri Lanka and is densely populated; thus, the vegetation cover is significantly lower than the rest of the wet zone. Hence, the predicted transpiration from vegetation cover decreases. As a result of the large extent of built-up areas, the evaporation from soil is also limited. Consequently, it is reasonable to assume that pan evaporation values are higher than the actual evapotranspiration in Colombo and the simulated results closely represent the actual evapotranspiration than pan evaporation. These discussion points are showcased in Table 3, which explicitly averaged annual pan evaporation, annual simulated PET, and annual rainfall at the 14 stations which were analyzed.At Puttalam, Polonnaruwa and Hambantota averaged annual simulated figures and observed results differ only by 4–6 mm. However, at all three locations, predicted values are higher than those observed during the rainy season and lower than those observed during dry months. The highest variations between the average annual figures of simulated and pan evaporation can be seen in the wet zone, which could be due to good vegetation coverage throughout the year. In Jaffna, Puttalam, and Hambantota, generally the hottest areas in Sri Lanka, annual evapotranspiration is higher than the annual precipitation. Deficit water must be coming from the groundwater, as the aquifers in these areas are being abstracted excessively (Imbulana et al., [21]). In the rest of the dry zone, although the evapotranspiration is less than the rainfall, 80–93% of precipitation loses to the atmosphere as evapotranspiration, as per the predictions, whereas observed values indicate 75–85% evaporation loss. The reason for this massive evaporation loss might be due to the extensive network of small to large scale man-made irrigation tanks, which provides a constant supply of water over a large land area. In addition, paddy, the major cultivation in these areas, relies on flood irrigation, keeping the soil saturated and aiding both the evaporation and transpiration processes. In the wet zone, the ratio of PET to rainfall ranges between 40 and 65%, and the higher percentages are in the highlands. The predicted annual PET over the entire country is about 2 280 674 mm, with an average of 1243 mm/year. Approximately 65% of PET is from canopy transpiration.


### 4.3. Seasonal Changes of PET

Five stations, i.e., Vavuniya, Kurunegala, Bandarawela, Katugastota, and Galle, which give the smallest RMSEs, were selected from the dry zone, intermediate zone lowlands, intermediate zone high lands, wet zone high lands, and wet zone lowlands, respectively, to graphically illustrate the temporal variations of average monthly (1) PET, (2) ETs, (3) ETc, (4) rainfall, and (5) average daily Net radiation (refer to Figure3).Figure 3
Seasonal patterns of monthly PET, ETs, ETc, rainfall, and averaged daily net solar radiation (all parameters are averaged for the period of 2009–2019). (a) For Vavuniya. (b) For Kurunegala. (c) For Bandarawela. (d) For Katugastota. (e) For Galle.
(a)(b)(c)(d)(e)Figure3 indicates that the total monthly evapotranspiration follows a pattern similar to the seasonal fluctuations of net radiation. Nevertheless, a few minor deviations correlated with rainfall and cultivation seasons can be observed. Within a year, monthly PET values fluctuate between 125 and 70 mm in the Vavuniya (84%), Katugastota (69%), and Bandarawela (50%) and Kurunegala, whilst Galle (46.3%) and Kurunegala (43%) are seeing much more minor variations around 125 mm and 100 mm, respectively. Vavuniya experiences the highest PET, ETc values in the period from June to August, while peaks arise either in July or August. The reason is the high incoming solar radiation and the irrigated cultivation, which facilitates the continuous supply of water despite the dry spell the dry zone experiences from June to August. The two stations from the intermediate zone show two distinct peaks of PET, one in March or April and the other one in July or August. Similarly, two peaks are seen in the wet zone, either in March or April and October. A decrease in incoming radiation fluxes has been observed all over the country from October to December. Relating to this observation, PET shows a decreasing trend from October to December in all three climatic zones, with a mean rate of 10 mm/month. ETs, as predicted, is associated with the rainy seasons. The highest amount of rainfall has been recorded from October to November, and this period is a rainy season for the entire island. Thus, a peak in ETs can be observed everywhere from October to December, and it decreases from January to mid-year until ETs record their lowest in June or August. When annual sums are considered in all the climate zones, ETc is superior to ETs, unless at bare lands, build-up lands, and water bodies. In Vavunia and Bandarawela, ETc is always greater than ETc, indicating that the dry substrate has limited evaporation from the soil.


### 4.4. Spatial Distribution of PET

The model outputs, PET, ETs, and ETc, were averaged from 2009 to 2019 and these were exported to the ArcGIS environment to visualize the spatial distribution of PER, ETs, and ET covering the entire country. These distributions are presented in Figure4.Figure 4
Spatial distribution of averaged annual ETc, ETs, and PET in 6km×6km grid (created by the authors). (a) For ETc. (b) For ETs. (c) For PET.
(a)(b)(c)The spatial distribution of these parameters shows significant heterogeneity over the small landmass of Sri Lanka. Low to medium PET occurs in the high lands. ETs also account for the evaporation from water bodies. It is shown that higher ETs occur in the North-central part, the Anuradhapura region, most evidently due to a large number of irrigation tanks in the region. The wet zone also has shown high ETs, as expected, due to high precipitation keeping the soil moist throughout the year.Table4 showcases the predicted values of annual PET, ETs ETc for different land uses. As discussed earlier, nine land uses, including paddy, rubber, coconut, Chena, tea, grasslands, shrublands, home gardens, and urban and built-up areas, were considered in generating the predicted results.Table 4
Maximum LAI values for different vegetation types [20].
CodeLand use typeClassification as per the literatureLAImaxhc (m)lmax (m)Fclrstmin (sm−1)NDVI98%zog (m)1CoconutEvergreen needle leaf forests5.5170.00111500.6890.022RubberEvergreen broadleaf forests7300.0501500.6110.023Forest, unclassifiedMixed forests5.7200.040.51500.7210.024Homesteads/gardenOpen shrub lands310.0111000.6740.025Shrublands310.0111000.6740.026Tea310.0111000.6740.027GrasslandsGrasslands1.80.80.0101150.6740.018Marshy landsPermanent wetlands610.010650.6740.019ChenaCroplands70.60.010900.6740.0510Other cultivations70.60.010900.6740.0511Paddy70.60.010900.6740.0512Urban and built-upUrban and built-up000000.6740.0213Barren landBarren or sparsely vegetated0.30.050.0111200.6740.0114Water bodiesWater bodies000000.6740.001As shown in Table4, PET from paddy lands, which has been predicted to have an annual average of 1392 mm and an annual maximum of 2035 mm, is the highest among the PET from cultivated lands. This is related to flood irrigation. It is also shown that PET from paddy fields is highest in the dry zone and lowest in the wet zone. This spatial distribution mainly correlates with the net solar radiation flux, which is higher in the dry zone. Rubber, which is planted only in the wet zone on gentle slopes, transpires water at an average rate of 1293 mm/year and seconds only to paddy. The evergreen broadleaf plant, growing on year-long moist soil, explains the second-highest transpiration rate of rubber. The lowest PET values occur in the urban built-up areas, grasslands, and home gardens with an average of 405 mm/year, 665 mm/year, and 872 mm/year, respectively.


## 4.1. Impact of Data Scarcity on Analysis

The model results of fourteen grid cells (refer to Table3), in which the selected meteorological stations lie, compared with the recorded pan evaporation at the respective meteorological stations. They could lessen the error resulting from interpolation. The model was not calibrated in this study due to a lack of data and a shorter data period. Coupled with the limitations exerted by NDVI and soil moisture data, further challenges to model calibration were exerted by the uncertain canopy and soil parameters obtained from the literature. However, the root means square error (RMSE) of model prediction to observation, 0.75 mm/day, suggests that the model can predict PET and its component over Sri Lanka with moderate accuracy.Table 3
Predicted annual averaged PET, ETs, and ETc from different land-use types.
PET (mm/year)ETs (mm/year)ETc (mm/year)Paddy1391.65505.03886.63Rubber1292.72168.201124.52Coconut1276.18206.751069.42Chena (arid cultivation)1255.75331.32924.43Tea1199.56295.68903.87Grass lands663.07429.49233.58Shrub lands1173.70285.13888.57Home gardens872.00300.53571.47Urban and built-up areas405.46405.46Historical solar radiation measurements are not commonly available at many meteorological stations in Sri Lanka. This reduced the number of meteorological stations we could use for this study. The density of selected meteorological stations is notably lesser in the dry zone. Since there is no station in the South-East (SE) part of the country, observed values at Badulla and Bandarawela seem to influence the interpolation outputs in the SE region, irrespective of the interpolation method. As per the actual situation, in the SE low land, some meteorological parameters such as rainfall, temperature, and humidity show a drastic difference from those at Badulla and Bandarawela. Therefore, the accuracy and precision of the model predictions in the SE region are questionable. However, the model can be applied to the rest of the dry zone and intermediate zone with moderate confidence and the wet zone with high confidence.


## 4.2. Comparison of Predicted PET, Pan Evaporation, and Precipitation

The model generated PET values are greater than the observed, except at the Colombo station, at which the simulated results are significantly smaller than recorded values by an average of 32%. The disparity between simulated and observed values at Colombo can be attributed to several factors. The model is highly sensitive to LAI, and there is a positive correlation between LAI and simulated results, which has been confirmed by Zhou et al. [20]. Colombo is the capital city of Sri Lanka and is densely populated; thus, the vegetation cover is significantly lower than the rest of the wet zone. Hence, the predicted transpiration from vegetation cover decreases. As a result of the large extent of built-up areas, the evaporation from soil is also limited. Consequently, it is reasonable to assume that pan evaporation values are higher than the actual evapotranspiration in Colombo and the simulated results closely represent the actual evapotranspiration than pan evaporation. These discussion points are showcased in Table 3, which explicitly averaged annual pan evaporation, annual simulated PET, and annual rainfall at the 14 stations which were analyzed.At Puttalam, Polonnaruwa and Hambantota averaged annual simulated figures and observed results differ only by 4–6 mm. However, at all three locations, predicted values are higher than those observed during the rainy season and lower than those observed during dry months. The highest variations between the average annual figures of simulated and pan evaporation can be seen in the wet zone, which could be due to good vegetation coverage throughout the year. In Jaffna, Puttalam, and Hambantota, generally the hottest areas in Sri Lanka, annual evapotranspiration is higher than the annual precipitation. Deficit water must be coming from the groundwater, as the aquifers in these areas are being abstracted excessively (Imbulana et al., [21]). In the rest of the dry zone, although the evapotranspiration is less than the rainfall, 80–93% of precipitation loses to the atmosphere as evapotranspiration, as per the predictions, whereas observed values indicate 75–85% evaporation loss. The reason for this massive evaporation loss might be due to the extensive network of small to large scale man-made irrigation tanks, which provides a constant supply of water over a large land area. In addition, paddy, the major cultivation in these areas, relies on flood irrigation, keeping the soil saturated and aiding both the evaporation and transpiration processes. In the wet zone, the ratio of PET to rainfall ranges between 40 and 65%, and the higher percentages are in the highlands. The predicted annual PET over the entire country is about 2 280 674 mm, with an average of 1243 mm/year. Approximately 65% of PET is from canopy transpiration.


## 4.3. Seasonal Changes of PET

Five stations, i.e., Vavuniya, Kurunegala, Bandarawela, Katugastota, and Galle, which give the smallest RMSEs, were selected from the dry zone, intermediate zone lowlands, intermediate zone high lands, wet zone high lands, and wet zone lowlands, respectively, to graphically illustrate the temporal variations of average monthly (1) PET, (2) ETs, (3) ETc, (4) rainfall, and (5) average daily Net radiation (refer to Figure3).Figure 3
Seasonal patterns of monthly PET, ETs, ETc, rainfall, and averaged daily net solar radiation (all parameters are averaged for the period of 2009–2019). (a) For Vavuniya. (b) For Kurunegala. (c) For Bandarawela. (d) For Katugastota. (e) For Galle.
(a)(b)(c)(d)(e)Figure3 indicates that the total monthly evapotranspiration follows a pattern similar to the seasonal fluctuations of net radiation. Nevertheless, a few minor deviations correlated with rainfall and cultivation seasons can be observed. Within a year, monthly PET values fluctuate between 125 and 70 mm in the Vavuniya (84%), Katugastota (69%), and Bandarawela (50%) and Kurunegala, whilst Galle (46.3%) and Kurunegala (43%) are seeing much more minor variations around 125 mm and 100 mm, respectively. Vavuniya experiences the highest PET, ETc values in the period from June to August, while peaks arise either in July or August. The reason is the high incoming solar radiation and the irrigated cultivation, which facilitates the continuous supply of water despite the dry spell the dry zone experiences from June to August. The two stations from the intermediate zone show two distinct peaks of PET, one in March or April and the other one in July or August. Similarly, two peaks are seen in the wet zone, either in March or April and October. A decrease in incoming radiation fluxes has been observed all over the country from October to December. Relating to this observation, PET shows a decreasing trend from October to December in all three climatic zones, with a mean rate of 10 mm/month. ETs, as predicted, is associated with the rainy seasons. The highest amount of rainfall has been recorded from October to November, and this period is a rainy season for the entire island. Thus, a peak in ETs can be observed everywhere from October to December, and it decreases from January to mid-year until ETs record their lowest in June or August. When annual sums are considered in all the climate zones, ETc is superior to ETs, unless at bare lands, build-up lands, and water bodies. In Vavunia and Bandarawela, ETc is always greater than ETc, indicating that the dry substrate has limited evaporation from the soil.


## 4.4. Spatial Distribution of PET

The model outputs, PET, ETs, and ETc, were averaged from 2009 to 2019 and these were exported to the ArcGIS environment to visualize the spatial distribution of PER, ETs, and ET covering the entire country. These distributions are presented in Figure4.Figure 4
Spatial distribution of averaged annual ETc, ETs, and PET in 6km×6km grid (created by the authors). (a) For ETc. (b) For ETs. (c) For PET.
(a)(b)(c)The spatial distribution of these parameters shows significant heterogeneity over the small landmass of Sri Lanka. Low to medium PET occurs in the high lands. ETs also account for the evaporation from water bodies. It is shown that higher ETs occur in the North-central part, the Anuradhapura region, most evidently due to a large number of irrigation tanks in the region. The wet zone also has shown high ETs, as expected, due to high precipitation keeping the soil moist throughout the year.Table4 showcases the predicted values of annual PET, ETs ETc for different land uses. As discussed earlier, nine land uses, including paddy, rubber, coconut, Chena, tea, grasslands, shrublands, home gardens, and urban and built-up areas, were considered in generating the predicted results.Table 4
Maximum LAI values for different vegetation types [20].
CodeLand use typeClassification as per the literatureLAImaxhc (m)lmax (m)Fclrstmin (sm−1)NDVI98%zog (m)1CoconutEvergreen needle leaf forests5.5170.00111500.6890.022RubberEvergreen broadleaf forests7300.0501500.6110.023Forest, unclassifiedMixed forests5.7200.040.51500.7210.024Homesteads/gardenOpen shrub lands310.0111000.6740.025Shrublands310.0111000.6740.026Tea310.0111000.6740.027GrasslandsGrasslands1.80.80.0101150.6740.018Marshy landsPermanent wetlands610.010650.6740.019ChenaCroplands70.60.010900.6740.0510Other cultivations70.60.010900.6740.0511Paddy70.60.010900.6740.0512Urban and built-upUrban and built-up000000.6740.0213Barren landBarren or sparsely vegetated0.30.050.0111200.6740.0114Water bodiesWater bodies000000.6740.001As shown in Table4, PET from paddy lands, which has been predicted to have an annual average of 1392 mm and an annual maximum of 2035 mm, is the highest among the PET from cultivated lands. This is related to flood irrigation. It is also shown that PET from paddy fields is highest in the dry zone and lowest in the wet zone. This spatial distribution mainly correlates with the net solar radiation flux, which is higher in the dry zone. Rubber, which is planted only in the wet zone on gentle slopes, transpires water at an average rate of 1293 mm/year and seconds only to paddy. The evergreen broadleaf plant, growing on year-long moist soil, explains the second-highest transpiration rate of rubber. The lowest PET values occur in the urban built-up areas, grasslands, and home gardens with an average of 405 mm/year, 665 mm/year, and 872 mm/year, respectively.


## 5. Implications on Water Resources Management

As the S-W model accounts for leaf area index, with proper calibration, the model can be used in agricultural and irrigation water management, which was evident from the observations of this work: (1) peaks in plant growing seasons (Figure3) and (2) significantly higher PETs in extensively cultivated areas (Figure 4). However, given the complexity of the model, the model may use for long-term irrigation planning instead of seasonal planning. Since the model shows an excellent correlation with pan evaporation, which has been used in water management in Sri Lanka, the S-W model could be utilized for future predictions of evapotranspiration for climate change resilience attempts.


## 6. Summary and Conclusions

This study aims to develop the S-W model and evaluate the model’s efficiency in predicting PET over entire Sri Lanka. The model demands many meteorological data, soil, and vegetation parameters. Since some required meteorological data are not readily available, especially in data-scarce countries like Sri Lanka, the model holds many limitations. In addition, most of the soil and vegetation parameters are impossible or hard to measure. These facts restrict the extent to which the model can be calibrated and validated and may incur significant yet inevitable uncertainties in the model outputs.However, even under these limitations, the model developed gives an RMSE of 0.75 mm/day when compared with historic pan evaporation values. It can be concluded that a moderate to high confidence can be placed on model predictions, except in the SE quadrant of the country. The lower density of weather stations in the SE region caused this variation. The annual ETc, ETs, and PET were presented in 6km×6km grid to showcase the spatial distribution of these important parameters. Therefore, the country now has its own spatial distribution maps for the PET.The model shows that the PET over Sri Lanka is primarily governed by the incoming radiation flux, followed by rainfall and cultivation seasons. The model is highly sensitive to LAI, thus the vegetation cover. In the dry zone, 80–90% of rainwater escapes to the atmosphere through evaporation, while the percentage loss in the wet zone varies within 40–65%. Evaluation of annual figures indicates that approximately 65% of the total evapotranspiration amount is coming from vegetation canopy through transpiration. The model can be used as a guide to develop new agricultural policies in the changing climate of Sri Lanka. In addition, this study can be used as an example for tropical countries like Sri Lanka.



---

*Source: 6819539-2022-04-23.xml*