Long-Time Variation and Mechanism of Surface Energy Budget Over Diverse Geographical Regions in Pakistan

—Earth’s energy budget is a major force that drives global climate. The long-term pattern of land surface energy bud-get with pronounced biophysical effects on climate was normally ignored at a regional scale, particularly in Pakistan. In this article, the land surface energy budget from 2001 to 2018 was estimated and analyzed from Moderate Resolution Imaging Spectroradiometer (MODIS) and Cloud and the Earth’s Radiant Energy System (CERES) observations over three geographical regions (northern highlands, Indus plains, and Bolochistan plateau) in Pakistan. Biophysical and energy budget parameters, such as Land Surface Temperature (LST), albedo, emissivity, and Normalized Difference Vegetation Index (NDVI), were obtained from MODIS, whereas the downward shortwave solar and longwave thermal radiation were obtained from CERES satellite data. Spatiotemporal trends of three energy budget parameters, i.e., net radiation, latent heat ﬂux, and sensible heat ﬂux, and three biophysical parameters, i.e., albedo, NDVI, and LST, were investigated from 2001 to 2018. The latent heat ﬂux showed a signiﬁcant increase with a trend of 0.24, whereas a decrease in sensible heat ﬂux with a trend of − 0.21 was observedoverPakistan.Netradiationshowedanignorableincreasewithatrendof0.054overthewholeofPakistan.Asigniﬁcant negativerelationwasfoundbetweennetradiationandsensibleheatﬂuxwithalbedo,whereasasigniﬁcantpositiverelationwasfound betweenlatentheatﬂuxandNDVI.Biophysicalparameters,suchas


I. INTRODUCTION
T HE surface and atmospheric energy budget is the principal driving force of the earth's climate.The geographical distribution of the energy budget is more important at the surface and decides the energy distribution between different fluxes.The surface energy budget drives the hydrological and heat cycle in a particular location [1].It is important to monitor the spatiotemporal patterns of the energy fluxes to understand the local climate better.Changing biophysical parameters, such as clearing forests or transforming one land cover into another, alters the energy budget patterns [2], [3].Also, variations in energy flux over different latitudes are of prime importance.Different responses were observed for different latitudinal levels [4], [5].The energy budget at the surface is written as where R N is net radiation and LE, SH, and G are latent heat flux, sensible heat flux, and ground heat flux, respectively.R N is described as a balance between the radiations which the earth's surface receives and the radiations it reflects as well emits into the atmosphere.It is a basic quantitative measure to analyze any perturbation in climate [6], [7].Anthropogenic activities such as land use and land cover change, emission of greenhouse gases, and aerosol optical depth influence R N [8]- [12].R N comprises albedo-driven net shortwave radiation and Land Surface Temperature (LST)-driven net longwave radiation.Both albedo and LST are influenced by natural and anthropogenic changes [13].
It is evident from recent literature that clearing of forests into cropland or any other land cover may vary surface albedo and LST, which ultimately resulted in significant variations in R N [4], [7], [14].For the impact assessment of land cover changes on the surface energy budget, LST and Normalized Difference Vegetation Index (NDVI) are often used as investigative measures [5], [7], [15]- [17].
The right-hand side of (1) presents the nonradiative energy fluxes in which G is the available energy while LE and SH are This work is licensed under a Creative Commons Attribution 4.0 License.For more information, see https://creativecommons.org/licenses/by/4.0/turbulent fluxes [18], [19].LE describes an isothermal process that involves the vaporization of water (change of phase) at the expense of available energy.LE is directly related to NDVI as the higher value of NDVI denotes more vegetation cover, which evaporates more water and hence causes a net cooling effect [7], [20], [21].SH is a process in which no change of phase occurs, and solar radiations are directly converted into measurable heat.Transferred energy from the surface to the aloft atmosphere depends upon temperature difference [1], [22].The impact of land cover change is equally important for radiative as well as nonradiative processes [17], [23], [24].Two atmospheric parameters, i.e., cloud cover and aerosols, largely influence the surface radiation budget.The influence of clouds on surface net radiation is generally referred to as cloud radiative effects, the difference of energy between the actual energy and the energy if there would be clouds [25].Aerosols influences surface net radiation by either scattering or absorbing radiations.Moreover, aerosols have indirect effects associated with the interactions of clouds.A recent study reported the significant impact of aerosols on global net radiation [12].
Despite the above-mentioned significant atmospheric parameters, this study focuses only on surface parameters to investigate the mechanism of energy budget variations.Concerning the surface factors, the variations in energy budget are the response of land cover change, anthropogenic activities, and latitudinal variations.These factors exist at their peak in Pakistan, a highly heterogeneous geographic area including high mountains in the north and the coastal region in the south [26]- [28].
Bastiaanssen et al. [29] studied energy fluxes over the whole Indus basin, including Pakistan, India, China, and Afghanistan.They proposed a new ETLook model to compute LE using microwave remote sensing data.Liaqat et al. [30] analyzed the spatiotemporal distribution of actual evapotranspiration (ET) using a surface energy balance system over the Indus Basin region that belongs to Pakistan.Liaqut et al. [31], Waseem and Khayyam [32], and Hussain and Karuppannan [33] studied the impact of land use and land cover changes on LST over different regions of Pakistan.However, long-term climate change analysis using energy budget parameters is yet widely ignored in Pakistan.According to the authors' knowledge, this is the first-ever comprehensive analysis of all energy budget parameters in Pakistan.The novelty of this study is the spatiotemporal analysis of long-term energy budget patterns over Pakistan using satellite data.
Despite their accuracy, ground-based measurements limit our knowledge of certain time and space events, whereas satellite observations overcome these limitations and prove their worth for climate change studies [34]- [36].In this study, satellite data from Moderate Resolution Imaging Spectroradiometer (MODIS) and Cloud and the Earth's Radiant Energy System (CERES) sensors were used to analyze three energy budget parameters, i.e., R N , LE, and SH, and three biophysical parameters.i.e., albedo, LST, and NDVI, over three regions of Pakistan, i.e., northern highlands (NH), Indus plains (IP), and Bolochistan plateau (BP).The energy transformation from one medium to the other (atmosphere to the land surface and then backward) is explained based on cause and their produced effect.Here, LST, albedo, and NDVI were causes or driving factors that alter the R N , SH, and LE.Spatiotemporal, temporal, and seasonal patterns of all energy fluxes are computed and find the appropriate justifications for these patterns based on the aforementioned causes.Overall, the following three objectives are served in this study: 1) long-term spatiotemporal trend analysis of energy budget parameters over Pakistan; 2) comparison of energy budget parameters over various geographical regions of Pakistan; 3) explaining the impact of biophysical parameters on energy budget parameters.

II. STUDY AREA
Pakistan is situated in south Asia between 24°-38°N latitude and 61°-78°E longitude.The total area of Pakistan is 881 913 km 2 (340 509 square miles), with a total population of 214.2 million, which makes it the sixth-largest country in the world.Pakistan has four seasons, i.e., dry and cold winters (December to February), dry spring (March to May), wet summer with the eastern monsoon (June to August), and partially wet autumn (September to November).About 65% of the total annual rainfall occurs in the monsoon [30].Pakistan has two crop cycles per year, including wet Kharif (summer to autumn) followed by dry Rabi (winter to spring) [29].A major part of Pakistan has an arid and semiarid climate, but few regions are classified as subhumid to humid [28].
Geographically Pakistan is divided into three broader regions: NH, IP, and BP.NH is the higher topographical region, including many mountain ranges such as the Hindu Kush, Karakorum, and the Himalayas.Snow cover peaks, vast grasslands, and conifer forests are the major land covers in this region.The annual average temperature is approximately 15-20 °C, and the annual precipitation is more than 200 mm.IP is a diversified region that is abundant with cultivated lands, however, the Cholistan, Thar, and Thal deserts contribute to the versatility of this region.A rich fertile land of IP makes it the main source of vegetation production in Pakistan and serves as the prime producer of wheat, rice, cotton, and sugar cane.The Indus River is the main feature that crosses vertically throughout the IP.Topographically, IP varies from Himalayan Piedmont in the north to the sea level in the south.On average, IP has a moderate climate while monsoon is the major contributor to shaping its climate.Gradually, a decrease in precipitation while an increase in temperature is observed from the north to the south of IP.BP is also spread vertically from the Sulaiman range in the northeast to the Gawader coastal area in the south.The region is predominantly classified as a dry arid climate with barren lands as the major land cover.Few regions in BP have extreme summer with a maximum temperature of 55 °C [26].The geographical regions of Pakistan and land cover and elevation map are shown in Fig. 1.The land cover map is prepared from MODIS land cover product MCD12Q1.SRTM 30-m Digital Elevation Model is used to show the topography of Pakistan.With highly climatic and topographical variability, it is important to study the surface energy budget of Pakistan.In recent times, this area has remained ignored for such a study.

A. Cloud and the Earth's Radiant Energy System (CERES) Data
CERES sensors are mounted on Terra, Aqua, Suomi National Polar-Orbiting Partnership and National Oceanic and Atmospheric Administration-20 satellites [37].Energy Balance and Filled (EBAF) surface product is a reliable remote sensing product for global radiation data.This study used the CERES EBAF_Surface_Ed4.1 product for downward shortwave and longwave flux under clear sky conditions.Each CERES-EBAF measures filtered radiances in shortwaves from 0.3 to 5.0 μm, a window from 8 to 12 μm, and a total from 0.3 to 200 μm wavelengths [38].Radiance to flux conversion uses empirical distribution models [39].Uncertainties in shortwave and longwave downward clear-sky flux are 4 W•m −2 and 6 W•m −2 , respectively [40].These data were obtained from the NASA Langley Research Center CERES ordering tool at. 1

B. Moderate Resolution Imaging Spectroradiometer (MODIS) Data
MODIS provides multiple reliable products related to energy budget, which are widely used in recent literature.Albedo is a ratio between incidents to the reflected shortwave radiation.This study used MODIS daily composite product MCD43A3 to obtain albedo.This product consists of MODIS bands 1-7 and three broad bands, which cover visible (0.3-0.7 μm), near-infrared (0.7-5.0 μm), and shortwave (0.3-5.0 μm) spectral regions.It contains directional hemispherical reflectance (direct or black sky albedo) and bi-hemispherical reflectance (diffuse or white sky albedo).Blue sky albedo is the overall reflectance of shortwave radiation and is a function of white sky and black sky albedo [41], [42].Actual albedo lies between black sky and white sky albedo.A recent study has reported an ignorable difference and a high correlation between these two types of albedo [5].Therefore, this study obtains blue sky albedo by taking the arithmetic mean of black and white sky albedo.
from MODIS is retrieved using a split-window algorithm with extensive accuracy up to 0.5K [43].MOD11A2 (from Terra) and MYD11A2 (from Aqua) MODIS LST 8-day products at 1-km spatial resolution are used in this study.These products are the composite of MODIS level-3 daily LST products (Terra: MOD11_A1 and Aqua: MYD11_A1).Values for each pixel are generated when at least two out of eight retrievals are available after cloud masking.Available LST observations from Terra and Aqua MODIS were merged using the Simplified Merge Scheme, i.e., "use an average of overlapping pixels, or available a single-pixel," to increase the spatiotemporal coverage [44], [45].
Emissivity, a dimensionless quantity, represents the ratio of thermal energy radiating from a surface to that radiated from a black body.MOD11A2 product provides emissivity from various narrow bands in the middle and thermal infrared spectrum from MODIS bands 29 (8.4-8.7 μm), 31 (10.78-11.28μm), and 32 (11.77-12.22μm) [46].MODIS emissivity product has a high calibration accuracy over different land covers [43].
ET is a phenomenon that is associated with energy consumption in the evaporation and transpiration process.MOD16A2, an 8-day MODIS product with 500-m spatial resolution, provides LE by integrating many parameters such as albedo, land cover, leaf area index, and a fraction of photosynthetically active radiation with meteorological data.Derivation of LE is based on the modified Penman-Monteith equation [47].The product was validated against AmeriFlux tower observations and showed a mean absolute error of 0.33 mm•day −1 .The unavailability of LE data over barren lands is a core limitation of this product; however, it is still widely used in energy budget studies [21], [24], [35].
NDVI is a useful indicator of vegetation cover as well as land cover classification.MOD13A3 product provides monthly vegetation index NDVI, based on atmospherically corrected surface reflectance product (MOD09) [20], [48].In this study, monthly 1-km NDVI was used to identify vegetation cover.
To better understand the energy budget pattern, knowledge of underlying land cover is mandatory.For this purpose, MODIS annual land cover product (MCD12Q1) at 500-m resolution is used.It includes five land cover classification schemes from which International Geosphere-Biosphere Program (IGBP) scheme was used in this study.IGBP classifies land cover into 17 different classes.All the aforementioned MODIS products were downloaded from NASA's Level-1 and Atmosphere Archive and Distribution System and Distributed Active Archive Center. 2 All parameters from satellite datasets that are used in this study are summarized in Table I.

C. Data Processing and Analysis
Remote sensing data from different sensors need some certain preprocessing for calibration.Since all the scientific data sets used in this study have different spatial and temporal resolutions (see Table I), all data sets were resampled to 1000-m spatial resolution using the bilinear resampling method, which determines the value of a pixel using a weighted distance average of four neighboring pixels.After the resampling, all data sets were averaged at monthly temporal resolution.By merging similar land cover categories, 17 land cover classes of the IGBP scheme were reclassified as forest, shrubs, grasslands, crops, urban, barren land, and water.Next, annual land cover images from 2001 to 2018 were integrated to make a final land cover map shown in Fig. 1.Only pixels with the same value for at least 15 years (more than 80% of time duration) were assigned to that particular land cover class.
After the preprocessing, R N and SH were computed.As discussed in Section I, R N is composed of net shortwave and net longwave radiation, defined as downward minus upward flux.In this study, these fluxes were computed using the following equations: where R S N , R L N , R S↓ , and R L↓ are net shortwave, net longwave, downward shortwave, and download longwave radiation, respectively, and α, ε, σ, and T are albedo, surface emissivity, Stephan Boltzmann constant (5.67 , and LST respectively.Thus, R N is written as G in (1) can be ignored due to its negligibly small value, also, it is balanced to zero over an annual cycle [29].Moreover, the energy absorbed by the ground during the day is normally released into the atmosphere at night.Hence, G on the diurnal cycle remains negligible [18].Thus, the monthly temporal resolution, the basic time unit in the current study, provides the necessary liberty to ignore G.All other parameters are observed from satellite data, and thus SH can be expressed as ( It is worth mentioning that ( 5) is based on ideal energy conservation law, and any minor turbulent energy flux due to friction, ground absorption (G), or any other uncertainties were ignored.Thus, it was assumed that total incoming and outgoing energy are equal, and (5) covers energy transformation between the atmosphere and land surface (within a finite region, i.e., Pakistan).
Before the trend analysis of R N , SH, and LE on a spatiotemporal, temporal and seasonal scale, first, their correlation was established with albedo, NDVI, and LST.For this purpose, Pearson's correlation (R) was performed on spatial images and at the pixel level.For the spatiotemporal analysis, simple linear regression was applied to the mean annual spatial images of all parameters and found the slope value.For temporal and seasonal variations, spatial averages of each geographical region of Pakistan (NH, IP, and BP) were analyzed.The slope was used to observe the temporal trend over each geographical region.It is mentioned above that observations of LE, and hence SH, are unavailable over barren lands due to the limitation of MODIS ET product.For this reason, the regional averages of LE, SH, R N , albedo, NDVI, and LST have different spatial coverages.To overcome this problem, spatial images of R N , albedo, NDVI, and LST were masked according to the spatial images of LE and SH.By doing so, all the observed parameters have the same spatial coverage across all regions of Pakistan.
Bowen ratio was used to find the nonradiative energy distribution.Bowen ratio (β) is useful for analyzing energy distribution and can be defined as β = SH / LE [49], [50].Higher values of β correspond to more energy consumption in raising temperature, whereas lower values correspond to more moisture content being evaporated [29].Bowen ratio explains whether the climatic trend is toward SH or LE.

A. Correlation Analysis
Pearson's correlational analysis was performed to establish a relationship between the parameters that serve as causes (LST, albedo, and NDVI) and the parameters that are their effects (R N , LE, and SH).Fig. 2 shows the spatial correlation, whereas Fig. 3 displays the interpixel correlation.Table II is the numerical representation of Fig. 3 and presents the slope as well as R between the parameters.All values are 99% significant at a 95% confidence level.Here, these correlations are analyzed based on underlying land cover and geographical regions as shown in Fig. 1.
Correlation between R N and albedo shows dominant-negative relation in all spatial locations of Pakistan irrespective of the underlying land cover.A highly negative slope of −130 and R of −0.71 indicates that R N is negatively related to albedo.R N shows the same negative relation with LST for the vegetated land cover but where the land cover is barren land (in BP) or high-altitude forest area (NH); both parameters are positively related.Thus, in Pakistan, R N shows a weak positive correlation.The relationship pattern of R N with NDVI is different, as it shows only a minor negative relation on barren lands, but a strong positive correlation on the vegetation land cover in IP and NH.

TABLE II PEARSON'S R AND LINEAR REGRESSION SLOPE BETWEEN ENERGY FLUXES WITH LST, ALBEDO, AND NDVI
Values after ± are standard errors in the slope estimation.
The higher NDVI means dense vegetation cover, which resulted in a higher R N .
The limitation of the LE product used in this study is more visible in the spatial analysis.LE values over barren lands of BP are not available and the analysis was based on other land covers.The most prominent positive relation of LE is with NDVI which shows an overall strong positive correlation with 0.38 R and 33.83 slope value.This positive correlation enables NDVI to be an important parameter to analyze evaporative cooling which is caused by higher LE values [21].LE shows a negative correlation with albedo over vegetation land cover while a positive correlation with the shrubs and grasslands in upper regions of BP.The same pattern extends to the LE and LST relation with the only difference that over forests of the NH region, LE shows a prominent positive relation.
In this study, SH is obtained from (5), and thus the limitation of LE extends to SH and the observations are unavailable over barren lands.SH shows a positive correlation with albedo and LST over vegetated lands.A negative correlation was observed over the shrubs between SH with the albedo and LST.SH shows a strong positive correlation with NDVI over shrubs and a negative correlation over vegetated lands.Overall SH demonstrates a strong negative correlation of −0.81 R with albedo while a prominent positive correlation of 0.62 R with LST.A positive impact of LST on SH is well established in recent literature [51].These results from Table II summarize all these correlations as R N and SH have inverse relations with albedo and linear relations with NDVI and LST.LE has an inverse relation with LST while linear relation with albedo and NDVI.

B. Spatiotemporal Variation of the Surface Energy Budget
Fig. 4 shows the slope for the spatiotemporal variations of all the studied parameters.Negative values (orange and yellow) describe decreasing trend while positive slope values (green and blue) describe an increasing trend from 2001 to 2018.Blank areas (white) are the regions for which LE and SH data were not available.Linear regression slope was computed on mean annual images of each parameter for spatiotemporal analysis.
For albedo, a dominant decreasing trend is observed over IP and upper regions of BP where the land covers are vegetation, shrubs and grasslands.However, barren lands of south BP and high elevated NH show a moderate increasing trend.The reverse pattern is observed for NDVI which shows a moderate increase over the entire country except for some sparse regions of NH and south BP.The most prominent increase is observed for LST, which shows an increase over barren lands of south BP, east IP, and NH regions, whereas a moderate increasing trend is observed for the vegetation cover and a moderate decreasing trend is observed for grasslands of northwest IP.
R N shows the mixed spatiotemporal trends.A moderate increasing trend is observed over barren land of BP except for the few coastal regions of south BP.In IP, the regions adjacent to the Indus River show a prominent increasing trend with the 0.5 to 1 slope value.Contrarily, northeast IP shows a moderate decreasing slope.NH shows an interesting well-mixed trend, with adjacent increasing and decreasing pixels.A possible reason for such a mixed trend is the hill shade of the mountainous region.The side facing the sun behaves differently from the other side to incoming solar radiation.
As discussed in the previous section, the MODIS LE product does not provide data over barren lands; therefore, observations are not available for southwest BP, northeast NH, and a few regions of east IP.All the regions with available LE data show increasing trends except a moderate decrease over a narrow strip of south BP where the land cover is shrubs and grasslands.These regions also show a decrease in NDVI suggesting the reduction of vegetation cover of these pixels that translated into the reduction of LE.NH shows the most prominent increasing trend for LE.The high-altitude conifer trees are the distinguished land cover in this area.
SH data are also available for the regions with the LE observations.A noteworthy trend is observed over NH where LE shows a prominent increasing trend while SH shows a decreasing trend of the same intensity.This opposite trend is also observed over the narrow strip of BP where SH shows a moderate increasing trend.The upper regions of BP are the area over which both LE and SH show the same increasing trends of moderate density.
Analysis of spatiotemporal trends of energy fluxes based on biophysical parameters further explains the link between these factors.Over NH, which belongs to higher latitudes, albedo proved a key factor to shape the trend of R N , but as the latitude decreased and underlying land cover changed from forests to crops, LST dominates the R N behavior.Over low latitude barren lands, its albedo again dominates R N trends.It can be concluded that a brighter surface has strong effects on R N while for dark surfaces like vegetation; LST has strong effects on R N .For LE trends, NDVI is the key factor since both show almost the same trends over the whole of Pakistan.NDVI is also considered to estimate the evaporative cooling caused by LE [21].Table II suggests the negative relation between LE and LST which is also endorsed in Fig. 4 as, over the lower NH where the LST shows a strong decrease, LE shows a significant increase.using the regional average of the corresponding parameter.For this purpose, first R N , albedo, NDVI, and LST were masked according to the spatial coverage of LE and SH to remove the potential mismatch between these parameters due to different spatial coverage.NH is the high elevated region of Pakistan with an abundance of conifer trees and high mountains, whereas a few peaks are covered with snow.After eliminating the barren mountains from the spatial image, the remaining area is full of trees and hence shows a lesser mean annual albedo as compared to the other regions of Pakistan.This is the only region that shows a minor increasing trend with a 1.03 × 10 −4 slope value, while IP and BP show a decreasing trend.Pakistan accumulatively shows a minor decreasing trend with a slope of −3.94 × 10 −4 .Decreasing albedo corresponds to the increase of forests and vegetation cover.Literature shows the inverse relation of albedo with radiative forcing [17], [23].This fact is verified by NDVI temporal trends which show an increasing trend over all regions of Pakistan.IP shows the highest NDVI as the crops are the major land cover in this region.Accumulatively NDVI shows an increasing trend with a 0.03 slope value.LST shows the symmetrical increasing trend over all the regions and accumulatively shows the slope value of 0.005 over the whole of Pakistan.

C. Temporal Variation of the Surface Energy Budget
Increasing LST agrees with recent literature which reported increased LST over various regions of Pakistan [32], [33].There is a significant difference in mean annual LST over different regions of Pakistan.The variety of land cover and altitudinal differences are the major causes of LST behavior [52].High altitude NH containing snow cover peaks shows the least mean annual LST while the barren BP and heavily populated IP show the maximum LST.R N shows the only decreasing trend over NH with the slope values of −0.15 while BP and IP show the increasing temporal trends with the slope values of 0.18 and 0.11, respectively.Accumulatively, R N over the whole of Pakistan shows an increasing temporal trend with the slope value of 0.054.Temporal trends of R N are exactly opposite to albedo over all regions of Pakistan.This is well explained in the lite of the results presented in Section IV-A where it was observed that albedo and R N have a significant negative correlation.
LE shows significantly higher values and a prominent increasing trend with the slope value of 0.53 over NH.Thick coniferous trees along with wet air over glaciers are the main reasons for higher LE over this region.Along with NH, IP and BP also show the increasing temporal trend for LE with slope values of 0.18 and 0.11, respectively.LE over IP mainly corresponds to the widely spread crop fields [29], [30].Accumulatively LE shows the increasing temporal trend over Pakistan with a slope of 0.24.
SH shows the most significant decreasing trends over NH with the slope value of −0.79, the highest trend value observed in this study.Also, SH has the least mean annual value in this region.This is because LE is dominant over NH.Significant increasing LE leads toward a similar decreasing trend of SH.As discussed above, albedo and SH have a significant opposite correlation thus increasing albedo trend over NH resulted in the decreasing SH temporal trend.A minor decreasing trend of SH is also observed over IP and BP with the slope values of −0.06 and −0.03, respectively.SH shows a decreasing trend over Pakistan with a slope value of −0.21.
Bowen ratio is a useful method to observe energy distribution for nonradiative energy fluxes [53].Fig. 7 shows the temporal variations of the Bowen ratio over NH, IP, BP, and overall Pakistan.The temporal trend was analyzed using the slope values.From Fig. 6, it is observed that the Bowen ratio shows a decreasing trend overall regions and the whole of Pakistan, concluding that the nonradiative energy budget is shifting toward LE.It is endorsed by temporal trend analysis where LE showed a continuous increase while SH showed a continuous decrease over all regions of Pakistan.Studies showed the inverse temporal trends between LE and SH [4], [24].Increased vegetation cover that was indicated by increasing NDVI in the temporal analysis might be a major cause of this change as NDVI is directly related to LE [21].

D. Seasonal Variation of the Surface Energy Budget
The seasonal variation of the energy budget was calculated using spatial averages of mean monthly values for each parameter from January 2001 to December 2018.Fig. 7 shows the seasonal variations of R N , LE, SH, albedo, NDVI, and LST over NH, IP, BP, and the whole of Pakistan.Table IV presents the maximum and minimum monthly values of each parameter over all regions.The masking of R N , albedo, NDVI, and LST according to LE and SH was also performed for seasonal analysis.Albedo behaves differently over NH than in other regions.It shows a peak in February and the least monthly averages in July.Snow accumulation over glaciers after the winter season caused the maximum while snow melting in summer caused the minimum albedo over NH.IP has the maximum albedo in June when most crops have been harvested and the bare soil remains in the crop fields.Bare soil has a higher albedo than vegetation.Minimum albedo is observed in February when the Rabi crops are fully grown over the region.In BP, October shows  Numerical values of the maximum and minimum observation of each parameter and their occurrence month are all given in Table IV.These values are the monthly spatial average of the corresponding geographical region.

V. CONCLUSION
This study investigated the variations in land surface energy budget in Pakistan using remote sensing data.Multiple parameters were analyzed from MODIS and CERES products, amongst them, LST, Albedo, and NDVI were considered as influencing factors upon energy fluxes.Mean monthly values were computed for each dataset to analyze spatiotemporal, temporal, and seasonal variation patterns.Moreover, variations within three different geographic regions of Pakistan (NH, IP, and BP) were separately discussed.
The correlation analysis showed that R N and SH negatively correlate with albedo while LST and NDVI have a positive correlation.Moreover, LE has a negative correlation with LST and a positive correlation with albedo and NDVI.It was also observed that NDVI is the key parameter for LE variations while albedo and LST are the key parameters for R N and SH variations respectively.These findings are in agreement with recent literature [7], [29], [51], [53].It is found that due to the increase in albedo, R N showed a decreasing trend over NH.However other regions of Pakistan showed a nominal increase.Therefore, the positive albedo effect was counterbalanced by the negative effect of LST.Thus, the overall R N did not vary significantly.Despite a nominal increase in LST and NDVI, all regions of Pakistan showed a decrease in SH.This decrease is associated with increasing LE, which showed an overall increase in all regions of Pakistan.Previous literature [51] suggested that LE dominates the climate impacts of biophysical changes in irrigated agricultural areas like Pakistan.However, spatiotemporal analyses showed that NDVI is not the only controlling factor, as peak values of LE were observed at higher latitudes where snow melting evaporates surface water content.The Bowen ratio elaborated that the shifting of nonradiative energy flux toward LE was most prominent in BP.
Seasonal analysis showed that R N over NH has the extreme behavior, i.e., least observations in winter and maximum in summer.The LE trend over the whole of Pakistan is dominated by IP, as this region has the majority of crop fields which contributes the most to ET during the fully growing crop season.SH behaves similarly to LST in the seasonal cycle and shows a single peak in summer.It was also concluded that ignoring the influencing parameters (albedo, LST, and NDVI) could not comprehensively explain the variations in energy fluxes (R N , LE, and SH).The novelty of this research is the use of long-term remote sensing data to analyze the impacts of surface properties on the energy budget in Pakistan.Meteorological or observational data only provide limited spatial and temporal availability, in contrast, remote sensing can provide large-scale continuous data through which climate change phenomena can be observed at a broader level.This work opens the way for future long-term energy budget trends and the impact assessment of different parameters in Pakistan.

Fig. 1 .
Fig. 1.Study area map showing location, topography, geographical regions, and land cover classes of Pakistan.

Fig. 2 .
Fig. 2. Spatial correlation between LST, albedo and NDVI with energy fluxes using Pearson's R. The land cover map of Pakistan is shown in the inset.

Fig. 3 .
Fig. 3. Pearson's correlation of surface energy fluxes with (a,d and g) Albedo, (b,e and h) NDVI and (c,f and i) LST.236 251 points (pixels) are used to establish each relation.

Fig. 5 .
Fig. 5. Temporal analysis of (a) R N , (b) albedo, (c) LE, (d) NDVI, (e) SH, and (f) LST from 2001 to 2018.Dotted lines are the trend lines that show the slope of the corresponding parameter.

Fig. 5
Fig. 5 shows the temporal variations of R N , LE, SH, albedo, NDVI, and LST from 2001 to 2018, using spatial averages of NH, IP, BP, and complete Pakistan (PK).The slope is used to analyze the temporal trends.Table III presents the slope values, including each parameter's standard error (SE) over each region.The mean annual value of each parameter was computed

Fig. 6 .
Fig. 6.Temporal variation of Bowen ratio over NH, IP, BP, and PK.Red lines are the slope lines while the blue dotted lines are the mean value of the Bowen ratio over that region.

Fig. 7 .
Fig. 7. Annual variations of (a) R N , (b) albedo, (c) LE, (d) NDVI, (e) SH, and (f) LST over NH, IP, BP, and Pakistan (PK).Each value is the spatial average of the corresponding parameter over the specific region.
the maximum albedo when there is dry autumn over the region while the growing season of February shows the least albedo.Accumulatively, NH albedo observations are dominated over the whole of Pakistan.NDVI over IP has two significant peaks in February and August/September.The first peak belongs to fully grown Rabi crops and the second peak belongs to the native vegetation (naturally grown) grown after the monsoon.This temporal pattern of NDVI over IP has a core impact on the NDVI pattern over Pakistan which shows the same peaks with lesser density.NDVI over NH shows only one peak in august, resulting from growing native vegetation after the snow melting season.LST has the most obvious seasonal pattern overall geographical locations of Pakistan, showing the hot summer and cool winters.LST over NH shows the least monthly values across the year.The seasonal pattern of R N resembles LST and shows the peak values in summer while the least values in winter.R N over NH shows extreme behavior, i.e., the least values are observed in winter, whereas highest values are observed in winter months, compared with other geographical regions of Pakistan.As far as the seasonal pattern is concerned, SH exhibits the same behavior as R N .The seasonal pattern of LE resembles the seasonal pattern of NDVI.LE shows a single peak over NH in August when there is peak native vegetation.Over IP and the whole of Pakistan, LE shows two peaks: the first in February corresponds to the fully grown Rabi crops, and the second peak is in August which corresponds to the growth of native vegetation after the wet monsoon season.

TABLE I SUMMARY
OF SATELLITE DATA USED IN THIS STUDY Table III presents the slope values, including each parameter's standard error (SE) over each region.The mean annual value of each parameter was computed

TABLE III SUMMARY
OF THE TEMPORAL TRENDS OF ALL STUDIED PARAMETERS FROM 2001 TO 2018 OVER NH, IP, BP, AND PK

TABLE IV MINIMUM
AND MAXIMUM VALUES OF ALBEDO, NDVI, LST, RN, LE, AND SH OVER NH, IP, BP, AND THE WHOLE OF PAKISTAN