- Email: [email protected]
Change detection of glaciers and snow cover and temperature using remote sensing and GIS: A case study of the Upper Indus Basin, Pakistan
Remote Sensing Applications: Society and Environment 18 (2020) 100308
Contents lists available at ScienceDirect
Remote Sensing Applications: Society and Environment journal homepage: http://www.elsevier.com/locate/rsase
Change detection of glaciers and snow cover and temperature using remote sensing and GIS: A case study of the Upper Indus Basin, Pakistan Abdul Jabbar a, b, Arsalan Ahmed Othman c, *, Broder Merkel a, Syed E. Hasan d a
Department of Hydrogeology, Institute of Geology, Technische Universit€ at Bergakademie Freiberg, Gustav-Zeuner Str. 12, 09599, Freiberg, Germany Baugrundbüro Dr. Matthias Mokosch, Dresdner Str. 39, 01683, Nossen, Germany c Iraq Geological Survey, Sulaymaniyah Office, Sulaymaniyah, Iraq d Department of Earth & Environmental Sciences, University of Missouri, Kansas City, MO, USA b
A R T I C L E I N F O
A B S T R A C T
Keywords: Snow cover Upper Indus Basin Himalaya NDSI Land surface temperature
Glaciers and snow cover in the Himalayas, Karakoram, and Hindu-Kush (HKH) mountain ranges are expected to decrease, resulting in serious consequences in terms of water availability for downstream population. This study is an attempt to detect temporal trend of snow cover in the Upper Indus Basin (UIB), in northern Pakistan, with respect to temperature change, spatial variation in land surface temperature, and to evaluate the accuracy of Landsat imagery for detecting land surface temperature. Landsat images, acquired for the period 1993 to 2019 and the Normalized Difference Snow Index (NDSI), were used to map the snow cover. Maximum snow cover of about 20,675 km2 was observed in April 2017 followed by 1996, and 1999, when it was about 18,800 km2; the minimum of about 14,500 km2 was recorded in April 1993, 2011 and 2018. The results show that there has been a very slight increase in snow cover area with time. The maximum elevation of snow-covered area, extracted from the Digital Elevation Model of Shuttle Radar Topography Mission (SRTM DEM) in the study area, is 7889 m while the minimum elevation is 991 m. Freezing temperature occurs above 3000 m elevation, indicating that these areas are covered by snow and glaciers for much of the year. Temperatures, of 20 � C to 25 � C, were present only at a few places–at the very high elevation of above 6000 m. We found inconsistencies in the Landsat temperature data, which were corrected using in-situ daily temperature data. Validation of Landsat extracted temperatures before and after using the in-situ daily temperature data improved the accuracy and removed the inconsistencies. We recommend use of in-situ diurnal temperature data to correct Landsat extracted temperature data.
1. Introduction Recent decades have witnessed significant impacts from climate change. Adverse effects of climate change have been quite pronounced on people, wildlife, and ecosystems, including the cryosphere. Since the 18500 s, receding of many highland glaciers all over the world have been observed and an exponential increase of receding during the last 30 years (Ren et al., 2006; UNEP, 2008; Zemp et al., 2006). It is estimated that fresh water from snow and glaciers provide fresh water supply to one-half of world’s population. Receding glaciers and shrinking ice sheets have adversely affected global fresh water supplies. The Himalayan ranges represent the most glaciated region in the world except polar ice sheets and Greenland (IPCC, 2007; Young and
Hewitt, 1990), but precise data on their characteristics and dynamics are lacking (Dyurgerov and Meier, 2005). These ranges not only play a vital role in providing water supply to about 1.3 billion people living in downstream basins of the 10 major rivers in Asia, but also in the global hydrological cycle, atmospheric circulation, and ecosystems services (Bajracharya and Shrestha, 2011). Recent research conducted by Gar €€ delle et al. (2012) and Ka ab et al. (2012) support the fact that while majority of glaciers are receding in the Himalayan range, mass gain or equilibrium condition has occurred in the Karakoram mountain glaciers as a result of favorable climatic conditions in the past 10 years (Gardelle €b et al., 2012). No need to repeat this reference in the et al., 2012; K€ aa same sentence; delete]. An increasing trend in winter precipitation (Fowler and Archer, 2006) and decreasing pattern in summer
* Corresponding author. E-mail addresses: [email protected] (A. Jabbar), [email protected] (A.A. Othman), [email protected] (B. Merkel), [email protected] (S.E. Hasan). https://doi.org/10.1016/j.rsase.2020.100308 Received 19 April 2019; Received in revised form 26 January 2020; Accepted 19 March 2020 Available online 24 March 2020 2352-9385/© 2020 Elsevier B.V. All rights reserved.
A. Jabbar et al.
Remote Sensing Applications: Society and Environment 18 (2020) 100308
Fig. 1. Study area, located in the extreme north of Pakistan, borders with China to the north, India to the east, and Afghanistan to the north-west.
temperature (Schmidt and Nüsser, 2009) has been observed. In the Karakoram, over 80% of the ice cover are located between 4000 m and 5500 m elevations and the extensive debris cover at the ablation zone reduces the glacier exposure which makes them insensitive to changing climate (Hewitt, 2011). Shrinking and thinning of glaciers is common in most of the Hima layan, Karakoram, and Hindu-Kush (HKH) regions, leading to formation of glacial lakes. Continued melting could lead to Glacial Lake Outbursts Floods (GLOFs) (Bolch et al., 2008; Zimmermann et al., 1986). Several studies in the HKH region have reported a large number of GLOFs due to glacier recession (Ashraf et al., 2012; Bajracharya et al., 2007; Zhang et al., 2015). Surging of a large number of glaciers have been reported in the Karakoram since 1980 (Copland et al., 2011; Hewitt, 2007), and the rate has doubled since 1990, which is likely due to positive mass budget (Copland et al., 2011). The study was designed to: (1) monitor the temporal trends of snow cover for the last 27 years in the Upper Indus Basin (UIB); and (2) evaluate the accuracy of Landsat imagery for detecting land surface temperature.
Batura (59.8 km); and Hisper (53.1 km) (Atif et al., 2015). The mean elevation of UIB is 4750 m most of which (about 60%) is above 4500 m elevation and only 12% is above 5500 m, which is the abode of per manent snow and glaciers (Tahir, 2011). The mighty Indus River along with all its major tributaries including Kabul, Jehlum, Chenab, Ravi, and Sutlej rivers originate in these mountainous regions; these are the main source of water to about 200 million people living in the Indus valley (Immerzeel et al., 2010). Ac cording to Archer and Fowler (2004) and Hewitt et al. (1989), the Indus River receives more than 80% of its annual flow from less than 20% of the extensively glaciated and snow covered zone above 3500 m eleva tion. Recent modelling carried out by Yu et al. (2013) suggests that the total contribution of the UIB to the Indus River system is about 135.68 km3, of which ~18% (23.44 km3) comes from the glacier melt, and ~82% comes from the winter snowpack. €ppen-Geiger classification, the climate of Upper According to the Ko Indus Basin is characterized as cold desert climate (BWk) with warm dry summer and cold winter with low precipitation. This climate type can be found in temperate regions that lie in the rain shadow of high mountains (Peel et al., 2007). As the influence of monsoon winds lessens north-westwards, about 90% of the UIB falls in the Himalayan rain shadow zone (Immerzeel et al., 2009). The annual precipitation of the arid valley floor at the UIB is 150 mm–200 mm, whereas in the Gilgit River valley at 4400 m altitude, precipitation increases to 600 mm; and, according to glaciological studies, the accumulation at 5500 m reaches 1500 mm–2000 mm (Fowler and Archer, 2005). There is a huge inequality in the climate setting of the UIB compared to the eastern Himalaya (Fowler and Archer, 2006; Young and Hewitt, 1990). The climate diversity across the entire Himalayas can be perceived from the fact that the eastern part receives more than 3000 mm and western part gets less than 300 mm precipitation annually that has direct and indirect effects on the rivers’ discharge (Immerzeel et al., 2009). The spatial and seasonal snow accumulation also varies widely in the region. The eastern and central Himalayas receive maximum snowfall from summer monsoon (Ageta and Higuchi, 1984), while in the northwest (Karakoram and UIB), westerly circulations winter most
2. Study area The UIB is located in the extreme north of Pakistan along the Chinese border (Fig. 1); its catchment extends upstream from Tarbela Dam reservoir and covers an area variously estimated to be about 206,000 km2 (Tahir, 2011); 220,000 km2 (Yu et al., 2013); and 200,677 km2 (Immerzeel et al., 2009). The UIB is home to some of the world’s tallest mountain ranges in the Himalayas, Karakoram, and Hindu-Kush (HKH) ranges with several peaks attaining heights of >7000 m. Bishop and Shroder (2010) report that out of the total 30 highest summits in the world, 13 are located in northern Pakistan. The world’s second tallest peak K2 (8611 m), the ninth highest, Nanga Parbat (6125 m) and Tirich Mir (7690 m) are all located here. Some well-known glaciers in the HKH Pakistan (in other words UIB) are: Siachen (75 km)–longest in the Kar akoram and the world’s second longest glaciers outside the Polar Re gions; Biafo (67.9 km) third longest in the world; Baltoro (62.1 km); 2
A. Jabbar et al.
Remote Sensing Applications: Society and Environment 18 (2020) 100308
winter snow accumulation comes from westerly circulating winds (Alford et al., 2016; Fowler and Archer, 2005; Quincey et al., 2011). Hewitt (2011) claims that the Himalayan glaciers, especially the central Karakoram, gain an extensive amount of snow from avalanches moving down the adjacent steep slopes. The topography of UIB shows marked variability: from the lowest elevation of less than 1000 m in the plain areas, near the Tarbela and Mangla dam reservoirs, to over 8000 m at a number of peaks (Tahir, 2011; Yu et al., 2013), and, according to Immerzeel et al. (2009), the lowest elevation in the UIB is 335 m and the highest is 8238 m. The lower Indus basin irrigation system is mainly dependent on the melt water from glaciers and snow cover that provide uninterrupted supply round the year, which has turned it into the world’s largest integrated irrigation system (Fowler and Archer, 2005). The inhabitants in the lower Indus Basin are mainly dependent on snow and glaciers melt water for agriculture purposes (Archer, 2003).
the region of interest. The cloud cover in some areas of the images taken in 1996, 1999 and 2015 was manually digitized and the shape files were created in Arc Map. Then, all the shape files were combined in one file through Boolean operation comprising all of the cloud cover area. Next, a union was created between the cloud cover file and the intersect area that resulted in a single file having a well-defined area for all scenes excluding cloud cover. This cloud cover shapefile was converted from vector to raster format due to the images in raster. Then, this raster cloud cover file and the raster images containing the threshold value (separate snow and non-snow area) were stacked together with each image one by one. Finally, a new file containing three values was prepared: -1 showing cloud cover area excluded from calculation, 0 representing non-snow area, and 1 displaying snow cover.
3. Methodology
The Normalized Difference Snow Index (NDSI) algorithm was used for mapping snow covered area which was calculated from Equation (1):
3.4. Normalized Difference Snow Index (NDSI) threshold selection
3.1. Scene selection and acquisition
NDSI ¼
Landsat images with 30 m pixel resolution captured in the month of April for the period 1993 till 2019 were searched in the United States Geological Survey (USGS) archive (http://glovis.usgs.gov/) and were received free of cost from the Earth Resources Observation and Science (EROS) center. The repeat period of Landsat TM (16 days) and its highaltitude position make the acquisition of cloud free images extremely difficult for a given time period. April was selected due to the avail ability of best-suitable images on the basis of image quality and cloud cover. Only images with excellent quality and less than 10% cloud cover were selected for this study. Data from Landsat sensors TM, ETMþ and OLI/TIRS were used in this study. In addition, eleven Operational Land Imager (OLI) scenes have also been used for different months between 2017 and 2018. Details are listed in Appendix A. The Landsat ETM þ data was missing since May 2003 due to scan line corrector failure (SLC off); also, in some other years owing to maximum cloud cover the data was not deemed suitable for inclusion in this study. The missing gap, to some extent, was filled by acquiring Landsat 5 TM images. The benefi cial characteristic of Landsat images from the EROS center is that they are already in (Level 1T format) geo-referenced, and orthorectified (without topographic distortion). The region of interest is positioned in zone 43 north and the path and row of the image are 149 and 35 respectively. The projected coordinate system for all images is UTM (Universal Transverse Mercator, WGS 1984 datum).
ðTM2 TM5 Þ ðTM2 þ TM5 Þ
(1)
The Atmospheric and Topographic Correction (ATCOR) model was employed to remove the atmospheric effects, sun illumination, and sensor viewing geometry from the images. In order to perform atmo spheric correction, the ATCOR model needs description of the compo nents in the atmospheric profile and data about the image of the area e.g. the exact time of image acquisition day, month, year; sensor type, calibration file, solar zenith and azimuth. The calibration file for at mospheric correction was created from the metadata file provided along with the image.
Where TM2 represents the reflectance in the band 2 (0.52 μm–0.6 μm) and TM5 refers to reflectance in band 5 (1.55 μm–1.75 μm) of the Landsat 5 and 7. In Landsat 8, band 3 and band 6 are used for NDSI because of the similar spectral bandwidths to the band 2 and band 5 of Landsat 5 and 7. Burns and Nolin (2014) state that atmospherically uncorrected data illustrate a conspicuous change in calculating the debris free glacier area and underestimates the glacier debris free area if the similar threshold is applied. For this reason, all scenes were atmo spherically corrected before selecting the threshold value. To get better results, the trial-and-error threshold method has been used instead of the automatic threshold method of Otsu, (1979), which is based on a one-dimensional (1D) histogram considering the spatial correlation between pixels. Therefore, its anti-noise capability is not satisfactory (Xiao et al., 2019). According to Hall et al. (1998), NDSI threshold value of �0.4 is an optimal value to map a snow pixel. Like wise, several investigators (Kulkarni et al., 2010; Xiao et al., 2001) used the threshold value of �0.4 while Burns and Nolin (2014) considered threshold value of �0.42 to be the best for debris free glacier areas after comparing the result of Landsat images to the IKONOS-2 and QuickBird high resolution imagery. Racoviteanu et al. (2008) report that NDSI is a very effective method for complex topography with highly shadowed regions whereas, Burns and Nolin (2014) applied the ratio technique TM4/TM5 and TM3/TM4 in examining the glacier area fluctuation in the Cordillera Blanca, Peru, but finally decided for NDSI because of good detection in shadows. Similarly (Krishna, 2005; Racoviteanu et al., 2008), also utilized the NDSI in the Himalayas for obtaining optimum results for shadowed areas. A number of threshold values were found in the literature (e.g., �0.42, �0.5, and �0.6) and all of them were tested; and results were carefully compared and analyzed visually to ascertain the impact on mapping the snow cover in shadowed regions. Based on these considerations, we selected the threshold value, � 0.42 (Burns and Nolin, 2014) because it enabled us to accurately determine snow pixels in shadow regions.
3.3. Cloud area masking
3.5. Classification of the area
Some images have missing data on the border, and others have an undefined area at the sides (all images did not cover exactly the same area), so frames in the form of shape files in ArcMap 10.4 were created for all images. Then an intersect area was developed for all scenes in order to get a well-defined area and to remove the error on the border of the images. Furthermore, a subset with the aid of the intersect area was created for all images to obtain exactly the same areas for all images in
The selected threshold value � 0.42 of all scenes was added together via raster calculator function in ArcMap that provided a combined file containing numbers ranging from 1 to 12. Where 1 exhibits the cloud cover region (excluded from the calculation), 0 displayed non-snow area and 1 to 12 (because there are 12 images) showed snow covered region. A very simple formula (the file containing the addition of all threshold value ¼ 12) was used in raster function that generated the snow cover
3.2. Atmospheric correction
3
A. Jabbar et al.
Remote Sensing Applications: Society and Environment 18 (2020) 100308
Fig. 2. Workflow procedures for snow/ice cover mapping and temperature extraction.
area (pixels) present in all images from 1993 to 2019 and was named as permanent (or minimum) snow area. Then, the function (the file con taining the addition of all threshold value > 0) & (the file containing the addition of all threshold value � 11) calculated the non-permanent (maximum) snow area that was not present in all images. Next, the files “permanent” (minimum) and “non-permanent” (maximum) were added together and the results presented as the total snow cover area in UIB since 1993 till 2019. Additionally, the positive and negative snow cover was calculated by performing the following function for each image i.e. (the threshold values file of each year 1994 – permanent (minimum) snow cover) ¼ Positive area and; the threshold values file of each year snow covers i.e. (non-permanent (maximum) snow cover – the threshold value of 1993) ¼ Negative area. Thus, positive area refers to the snow cover pixels that are present in 1993 images besides the permanent (minimum) snow area, whereas negative area refers to the snow cover area that is not present in 1993 images but existed in another images.
(Equation (2)). The atmospheric correction was performed using ATCOR Model for the thermal bands of each image before converting to radiance and radiance to surface temperature. The procedure given in NASA’s Landsat 7 Science Data Users Handbook (NASA, 2008), was used to estimate the surface temperature. A model was created in ERDAS IMAGINE version 16 and all the thermal bands were processed through this model. � � LMAXλ LMINλ Lλ ¼ (2) * ðQCAL QCALMINÞ þ LMINλ QCALMAX QCALMIN Where, Lλ ¼ spectral radiance at the sensor’s aperture in watts/(m2 * ster *
μm)
LMAXλ ¼ maximum detected spectral radiance in watts/(m2 * ster * μm) LMINλ ¼ minimum detected spectral radiance in watts/(m2 * ster * μm) QCAL ¼ the quantized calibrated pixel value in digital number (DN) QCALMAX ¼ the maximum quantized calibrated pixel value (255) QCALMIN ¼ the minimum quantized calibrated pixel value (1)
3.6. Land surface temperature estimation In order to compute temperature from the Landsat thermal bands, the digital number (DN) value of the images were converted to radiance and then the brightness temperature was estimated from the radiance 4
A. Jabbar et al.
Remote Sensing Applications: Society and Environment 18 (2020) 100308
Fig. 3. snow cover, non-snow cover, positive and negative areas during the years 1996 and 2015. Clouds cover areas are highlighted in yellow circle. (For inter pretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)
LMAXλ and LMINλ are the spectral radiances in digital numbers ranging from 1 to 255 for each band having different values and are given in the header file of the image. The radiance was converted to brightness temperature with Equation (3). K2 � T¼ � K1 ln Lλ þ 1
273:15
T ¼ effective at-satellite temperature in Kelvin K2 ¼ calibration constant 2 in Kelvin K1 ¼ calibration constant 1 in watts/(m2 * ster * μm) Lλ ¼ spectral radiance at the sensor’s aperture in watts/(m2 * ster * μm)
(3)
K2 and K1 are the pre-launch calibration constants provided in the header file of the image. Temperature records of two meteorological stations were obtained for validation. Point shape files for Gilgit and Skardu meteorological stations were created to validate the temperature
Where:
Fig. 4. Change of snow cover with cloud cover area for the period 1993 to 2019. A statistically significant trend is not perceptible. 5
A. Jabbar et al.
Remote Sensing Applications: Society and Environment 18 (2020) 100308
Fig. 5. Total snow and ice cover without cloud cover in the Upper Indus Basin in April for the period 1993 to 2019. A statistically significant trend is not perceptible.
The lowest snow cover area of about 14,500 km2 was recorded in April 1993 and 2011. Results show that there has been a very slight increase in the snow cover area with time. However, the non-parametric Mann Kendall (MK) test has been applied for the snow cover time series. The MK test shows that the 2-sided p-value is > 0.05, therefore, the trend is statistically not significant. Images showing change in snow cover for the years 1993, 1994, 1996, 1999, 2000, 2002, 2011, 2015, 2017, 2018, and 2019 are shown in Appendix D. Fig. 7 shows the positive and negative snow and ice cover areas in April for the period 1993 to 2019 (see heading 3.5 classification of the area, for positive and negative area). The year 1993 has the highest negative snow area of about 6300 km2 followed 5500 km2 for the year 2011. The year 2017 has the highest positive snow area of about 8000 km2. The years 1996, 1999, 2002, and 2003 have almost similar positive snow cover areas of about 5000 km2 and also nearly equal negative snow cover areas of about 2000 km2 with minor differences. Moreover, in 2000 and 2015, the positive and negative snow areas were almost equal, i.e. about 3500 km2 (Appendix B). Fig. 8 shows the change of snow cover for the period March 26, 2017 to November 24, 2018. The cloud cover in three images was masked to get only the snow cover area. The highest snow cover of about 6844 km2 was found in April 2017. The lowest snow cover area of about 20,676 km2 was recorded in August 2017.
of that particular pixel. One of these stations (better correlated with the temperatures extracted from Landsat data) was used to correct the bias of the temperature extracted from Landsat data. The other one was used to validate the corrected data. Following surface temperature estimation, the SRTM (1 arc or 30 m resolution) Version-3 DEM was downloaded from the USGS earth ex plorer website (https://earthexplorer.usgs.gov/). Elevation values for the study area, extracted from the SRTM DEM, gave maximum elevation of up to 7889 m and the minimum 991 m. Fig. 2 explains the workflow used for snow cover mapping and land surface temperature estimation. 4. Results 4.1. Results before exclusion of cloud cover area The above-mentioned procedure was performed once for all images, and layout were also created. However, as mentioned before, there were some cloud cover in three images that resulted in an error during the calculation of positive and negative area (Appendix B). For example, the highlighted (red) area in yellow circle is not a negative area but a cloud cover captured in 1996 image, that is represented as negative area and also the cloud cover area in 2015 image is shown as a positive area which is an error (see heading 3.5 classification of the area, for expla nation of positive and negative area). Finally, it was decided to exclude the cloud cover area and repeat the entire procedure for all images and re-check the results. Fig. 3 shows the cloud cover areas highlighted in a yellow circle. This area was estimated as positive snow cover area in one image and negative in another image. This error was removed by excluding/ masking the cloud cover areas from calculation. Fig. 4 shows the change of snow cover including cloud cover in April during the last 27 years. The cloud cover in three images was masked from the study area to get only the snow cover area. Figs. 5 and 6 represent the total snow and ice cover without cloud in the Upper Indus Basin in April for the period 1993 to 2019 (Appendix C). The highest snow cover of approximately 20,675 km2 was found in April 2017 fol lowed by 1996, and 1999 where the snow cover was about 18,800 km2.
5. Land surface temperature Fig. 9A shows the temperature extracted from the Landsat image. The dark brown color represents higher temperature and the dark blue lower temperature. Very low temperatures ( 20 � C to 25 � C) occur at very high altitudes of 6000 m to 7000 m. High temperatures reach 30 � C–40 � C at the valleys floor between 1000 m and 2500 m elevations. The very low and the very high temperature areas comprise a very small part of the total area–very low temperature of 30 � C is recorded in only 8 pixels, and very high temperature of 47 � C is recorded in only 4 pixels. The very low temperatures seem to be realistic because it is found at a very high altitude. Higher temperature, on the other hand, is exagger ated and contains an error. Fig. 9B displays the topographic variation in 6
A. Jabbar et al.
Remote Sensing Applications: Society and Environment 18 (2020) 100308
Fig. 6. Cloud cover, snow cover, non-snow cover, positive and negative area during different years.
the UIB region. The SRTM DEM of the Himalayas still has voids or data gaps. There is a very small area (white spots), on the eastern side of the study area where elevation data is missing. The maximum and minimum temperatures recorded at 8:30 a.m. on April 29, 2003 were used to verify temperature variations in the study area to display temperature changes in relation to the elevation (Ap pendix F). The altitude of meteorological observation station at Gilgit is 1460 m and at Skardu is 2317 m (Appendix F). Fig. 10 shows the maximum insitu temperature at the Gilgit and Skardu meteorological observation stations, along with temperature extracted from the Landsat imagery for six days (April 17, 1996; April 26, 1999; April 12, 2000; April 29, 2003; April 27, 2011; and April 22, 2015). The result does not match. In Gilgit, the average difference with the in-situ data is 4.33 � C. The Landsat re sults on an average are about 4.33 � C lower than the in-situ data. While, in Skardu, the average difference is 8.32 � C and the Landsat extracted results are about 8.32 � C higher than the in-situ temperature data
(Fig. 6). The coordinates of Gilgit and Skardu Metrological stations are shown in (Fig. 5). Some differences in temperature could be attributed to the difference in time the Landsat image was captured, and recording of the in-situ data. The maximum temperature data was used to validate the results due to lack of hourly data. The Landsat data was captured at about 9:00 a.m. which was close to the maximum temperature. So, this could be one of the reasons for the difference in both results (Gilgit and Skardu meteorological stations). In case of Gilgit, the results look satisfactory because the Landsat data is captured 3 h earlier than the maximum insitu temperature. So, the Landsat temperature may be close to the insitu temperature data at that time. Fig. 11 shows the maximum, mean and minimum temperature at the Gilgit and Skardu stations and Fig. 12 the sum of monthly precipitation data for April, taken from Pakistan Meteorological Department for the period 1990 to 2015. Fig. 11 shows the April monthly maximum, mean and minimum temperature at Gilgit and Skardu meteorological stations. There are 7
A. Jabbar et al.
Remote Sensing Applications: Society and Environment 18 (2020) 100308
Fig. 7. Positive and negative snow and ice cover area.
Fig. 8. Change of snow cover from March 26, 2017 to November 24, 2018.
large variations from year to year but the highest temperautre at both stations was recorded in 2007. Temperature at Gilgit exhibits a slight increase for maximum, mean and minimum temperautre while nearly stable conditions are observed in all cases at Skardu. However, the re sults are not significant statistically in both cases (Gilgit and Skardu). Significant variation in precipitation can be noted in Fig. 12 that shows the sum of April precipitation at Gilgit and Skardu meteorological stations for the period 1990 to 2015. High precipitation sums occur after
about 4–5 years (peaks) at both stations. In addition, the precipitation at Skardu area for most of the years is more than at Gilgit area. The highest precipitation recorded at Skardu in 2006 and at Gilgit in 1999 was about 145 mm and 90 mm, respectively. Overall, there are very small but statistically insignificant trends in precipitation at Gilgit and Skardu station.
8
A. Jabbar et al.
Remote Sensing Applications: Society and Environment 18 (2020) 100308
Fig. 9. Landsat temperature and elevation with meteorological stations location.
Fig. 10. Correlation of in-situ temperature and Landsat derived temperature for Gilgit and Skardu meteorological observatory.
6. Discussion and conclusions
acquisition of Landsat scenes even in a fixed week is impossible due to cloud cover, so some flexibility is needed for using suitable scenes. Considering these factors, the Landsat scenes during the month of April were acquired to examine the trend of snow cover for the period 1993 to 2019. Gathering reliable information regarding glacier/snow cover thick ness was not possible due to lack of data on high resolution and precise DEM, which are a very significant factors for investigating changing conditions. The SRTM DEM is not reliable for snow cover thickness investigation due to 9 m–16 m and greater vertical error, especially in high mountain region. Temperature extracted from Landsat was validated using the in-situ daily temperature data. Some differences are due to differences in time of recording Landsat data and the in-situ data. Maximum temperature data was used to validate the results due to lack of hourly data and
The Himalayas, the world’s highest mountain chain, is shrouded in cloud cover throughout the year and the chances to acquire images without cloud cover are slim. Several scenes for the study area at the required time period were found to have huge cloud cover, over 40%, which makes them inappropriate for time series analysis of snow cover monitoring. In this study, we used images having low cloud cover of about 10%. The cloud cover zones were masked from the study area that reduced the total area. The Landsat data coverage for the Upper Indus Basin, Pakistan since the 1980s–1990s is very poor, which limits the detailed investigation of glaciers and snow cover condition back in time. The common practice in remote sensing for time series studies of snow and ice is to select images on a fixed date of a particular month within each year. However, the 9
Remote Sensing Applications: Society and Environment 18 (2020) 100308
A. Jabbar et al.
Fig. 11. Monthly April maximum, mean and minimum temperature at Gilgit and Skardu meteorological stations (Source: Pakistan Meteorological Department).
Fig. 12. Variation of the sum of April precipitation over time at Gilgit and Skardu meteorological stations (Source: Pakistan Meteorological Department).
The largest snow cover of about 20,675 km2 was found in April 2017, followed by 1996 and 1999 when the snow cover was about 18,800 km2, while the lowest snow cover of about 14,500 km2 was recorded in April 1993, 2011, and 2018. The result shows that the snow cover in the month of April has slightly declined during 1993–2015. The year 1993 has the largest negative snow area, about 6300 km2 followed by the year 2011, when it was 5500 km2. The years 1996, 1999, 2002, and 2003 have almost similar positive snow area of about 5000 km2 and also nearly equal negative snow area of about 2000 km2, with minor dif ferences. Moreover, in 2000 and 2015, the positive and negative snow area (see heading classification of the area, for explanation of positive and negative area) was identical, about 3500 km2 (Appendix B). According to the SRTM DEM, the study area has a maximum eleva tion of 7889 m while the minimum elevation is 991 m. The maximum temperature reached about 42 � C while the minimum temperature dropped to approximately 25 � C. Temperatures of 0 � C and lower occurred above 3000 m elevation and decreased with increasing elevation; however, at an altitude of 7000 m all regions had temperature below 0 � C. The areas above 3000 m had temperature below 0 � C, which indicate that these areas are covered with snow and glaciers. As far as we could ascertain, there are no studies on estimating the temperature of snow cover with Landsat in the Upper Indus Basin area that we could use
because the Landsat data was captured at about 9:00 a.m., which was close to the maximum temperature. The lack of hourly climate data for the time when the satellite captured the data inhibited comparison of the satellite imagery to the in-situ data at the same time. Thus, in the case of Gilgit the average departure was 4.33 � C lower than the in-situ obser vation; and at Skardu the average departure was about 8.32 � C higher than the in-site temperature data. In both cases the in-situ temperature data and temperature data extracted from Landsat had a time difference of about 3 h. In case of Gilgit, the results are acceptable because the Landsat data was captured 3 h earlier than the maximum in-situ tem perature (more accurate than Skardu, where R2 is > 0.73). So, the Landsat extracted temperature may be close to the in-situ temperature at that time. The very low temperature values seem satisfactory because it is observed at a very high altitude. On the other hand, the higher tem peratures are exaggerated and are erroneous. Temperature at Gilgit exhibits a slight increase for maximum, mean and minimum temper autre while in the case of Skardu almost stable condition is observed in all cases. Therefore, we used the equation plotted in Fig. 10 of Gilgit station to correct the bias in the temperatures extracted from Landsat. We used the Skardu station to validate the bias corrected temperatures (extracted from Landsat data), which reduced the average error at Skardu station from 8.32 � C to 3.25 � C. 10
A. Jabbar et al.
Remote Sensing Applications: Society and Environment 18 (2020) 100308
to validate or compare the results. €€ According to Gardelle et al. (2012); Hewitt (2005); Ka ab et al. (2012), the glaciers in the Karakoram mountain ranges have increased in volume and have remained stable during the past 10 years due to cli matic conditions; Hewitt (2011) reports that more than 80% of ice cover is situated between 4000 m and 5000 m. Thus, the claim of these in vestigators seems to be true for higher elevation with freezing temperature. Topography of the Upper Indus Basin varies greatly (7889 m at the highest to 991 m at the lowest points in this study) whereas the tem perature and precipitation data are recorded at the very low Gilgit and Skardu meteorological observatories at the valleys floor, 1460 m and 2317 m, respectively. Therefore, they cannot be valid for the whole region with such a huge variation in topography and temperature. However, satellite data can be an ideal tool for evaluating temperature at a certain elevation, provided the Landsat or any other remote sensing data is precisely validated using in-situ temperature data because tem perature can be accurately estimated from the linear regression line for each elevation as shown in Fig. 10. In-situ temperature data recorded at the same time when the Landsat image was captured should be compared for validation and accuracy. Precise data about the glaciers/snow cover area and temperature esti mation would assist in understanding future climate-based changes. Further studies should be performed to find out the relation between snow cover estimation and retrieval of land surface temperature and to determine snow cover changes in relation to changes in temperature or precipitation.
relation to the development, writing, and publication of the article. All co-authors have seen and agree with the contents of the manuscript and there is no financial interest to report. We certify that the submission is an original work and has not been published nor has it currently been submitted for publication to another scientific journal, or being considered for publication elsewhere. Declaration of competing interest The authors declare no conflict of interest. CRediT authorship contribution statement Abdul Jabbar: Data curation, Formal analysis, Investigation, Methodology, Writing - original draft, Writing - review & editing. Arsalan Ahmed Othman: Conceptualization, Formal analysis, Investi gation, Methodology, Project administration, Supervision, Writing original draft, Writing - review & editing. Broder Merkel: Conceptu alization, Project administration, Supervision, Writing - review & edit ing. Syed E. Hasan: Writing - review & editing. Acknowledgements We are thankful to the United States Geological Survey Earth Re sources Observation and Science data center for providing Landsat im agery; the GEOSYSTEMS (Hexagon Geospatial) GmbH, Germany for ERDAS IMAGINE plus ATCOR software License; and the Pakistan Meteorological Department for providing climate data. This project would not have been possible without their valuable support.
Ethical statement The authors declare that all ethical practices have been followed in
Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi.org/10.1016/j.rsase.2020.100308. Appendix A. Landsat raw data Scene Type
Date
ID
Quality
Cloud Cover
TM5 TM5 TM5 TM5 TM5 TM7 TM7 TM5 OLI OLI OLI OLI OLI OLI OLI OLI OLI OLI OLI OLI OLI
1993/04/25 1994/04/28 1996/04/17 1999/04/26 2000/04/12 2002/04/10 2003/04/29 2011/04/27 2015/04/22 2017/03/26 2017/04/11 2017/08/01 2017/11/05 2017/12/07 2018/01/24 2018/02/09 2018/04/30 2018/08/04 2018/10/23 2018/11/24 2019/04/01
LT51490351993115ISP00 LT51490351994118ISP00 LT51490351996108ISP00 LT51490351999116AAA01 LT51490352000103AAA02 LE71490352002100SGS00 LE71490352003119EDC00 LT51490352011117KHC00 LC81490352015112LGN00 LC81490352017085LGN00 LC81490352017101LGN00 LC81490352017213LGN00 LC81490352017309LGN00 LC81490352017341LGN01 LC81490352018024LGN00 LC81490352018040LGN00 LC81490352018120LGN00 LC81490352018216LGN00 LC81490352018296LGN00 LC81490352018328LGN00 LC81490352019091LGN00
9 9 7 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9
1.1% 2.62% 08% 10% 10% 1.96% 0.98% 2.85% 14.69% 6.8% 1.06% 3.6% 3.28% 9.18% 3.7% 5.26% 4.36% 2.57% 7.38% 2.96% 6.61%
11
A. Jabbar et al.
Remote Sensing Applications: Society and Environment 18 (2020) 100308
Appendix B. Positive and negative snow cover area in each image Date
Positive area in pixels
Positive area in km2
Negative area in pixels
Negative area in km2
1993/04/25 1994/04/28 1996/04/17 1999/04/26 2000/04/12 2002/04/10 2003/04/29 2011/04/27 2015/04/22 2017/04/11 2018/04/30 2019/04/01
982,734 4,600,504 6,080,856 5,997,955 3,793,011 5,565,773 5,840,748 1,881,453 4,339,627 8,072,562 1,659,742 5,981,757
884.46 4140.45 5472.77 5398.16 3413.71 5009.20 5256.67 1693.31 3905.66 7265.31 1493.77 5383.58
7,028,045 3,410,275 1,929,923 2,012,824 4,217,768 2,445,006 2,170,031 6,129,326 3,671,152 36,243 634,583 222,059
6325.24 3069.25 1736.93 1811.54 3795.99 2200.51 1953.03 5516.39 3304.04 32.62 571.13 199.85
Appendix C. Total snow cover area in each image Date
Snow area in pixels
Snow area in km2
1993/04/25 1994/04/28 1996/04/17 1999/04/26 2000/04/12 2002/04/10 2003/04/29 2011/04/27 2015/04/22 2017/03/26 2017/04/11 2017/08/01 2017/11/05 2017/12/07 2018/01/24 2018/02/09 2018/04/30 2018/08/04 2018/10/23 2018/11/24 2019/04/01
15,883,117 19,500,887 20,981,239 20,898,338 18,693,394 20,466,156 20,741,131 16,781,836 19,240,010 22,022,869 22,972,945 7,603,917 12,155,484 14,955,918 21,066,522 19,699,663 16,560,125 8,003,115 17,525,150 21,991,987 20,882,140
14,294.81 17,550.81 18,883.12 18,808.50 16,824.05 18,419.54 18,667.02 15,103.65 17,316.00 19,820.58 20,675.65 6843.53 10,939.94 13,460.33 18,959.87 17,729.70 14,904.11 7202.80 15,772.64 19,792.79 18,793.93
Appendix DCloud cover, snow cover, non-snow cover, positive and negative area in April for all years
12
A. Jabbar et al.
Remote Sensing Applications: Society and Environment 18 (2020) 100308
13
A. Jabbar et al.
Remote Sensing Applications: Society and Environment 18 (2020) 100308
. (continued). 14
A. Jabbar et al.
Remote Sensing Applications: Society and Environment 18 (2020) 100308
Appendix E. Cloud cover, snow cover, and non-snow cover, areas in April for all years
15
A. Jabbar et al.
Remote Sensing Applications: Society and Environment 18 (2020) 100308
. (continued).
16
A. Jabbar et al.
Remote Sensing Applications: Society and Environment 18 (2020) 100308
Appendix F. Daily maximum temperature, in-situ data at Skardu and Gilgit observation stations along with temperature extracted from Landsat images for that date Date
Gilgit in-situ daily temperature
Gilgit/Landsat extracted temperature
Skardu in-situ daily temperature
Skardu/Landsat extracted temperature
1993.04.25 1994.04.28 1996.04.17 1999.04.26 2000.04.12 2002.04.10 2003.04.29 2011.04.27 2015.04.22
100* 100* 100* 30.2 29.4 26.2 29.5 34.5 26
25.4 24.11 100* 22.81 25.4 22.32 25.36 30.42 23.51
100* 100* 22.6 22.6 22.1 20.6 22.8 28 19.2
36.04 34.06 25.82 26.67 34.06 30.75 32.65 39.16 21.57
100* ¼ data not available.
Appendix G. Coordinates of the climate observation stations Observatory Name
Latitude
Longitude
Elevation
Gilgit Skardu
35� 180 35� 550
75� 410 74� 200
2317 m 1460 m
Source: Pakistan Meteorological Department.
17
A. Jabbar et al.
Remote Sensing Applications: Society and Environment 18 (2020) 100308
Appendix H. Modified Landsat temperature for the years 1999, 2000, 2002, 2003, 2011, and 2015
18
A. Jabbar et al.
Remote Sensing Applications: Society and Environment 18 (2020) 100308
References
Hewitt, K., Wake, C.P., Young, G.J., David, C., 1989. Hydrological investigations at Biafo glacier, Karakoram range, Himalaya: an important source of water for the Indus River. Ann. Glaciol. 13, 103–108. Immerzeel, W.W., Droogers, P., de Jong, S.M., Bierkens, M.F.P., 2009. Large-scale monitoring of snow cover and runoff simulation in Himalayan River Basins using remote sensing. Remote Sens. Environ. 113, 40–49. https://doi.org/10.1016/j. rse.2008.08.010. Immerzeel, W., van Beek, L.P.H., Bierkens, M.F.P., 2010. Climate change will affect the Asian water towers. Science 328, 1382–1385. https://doi.org/10.1126/ science.1183188. IPCC, 2007. IPCC Fourth Assessment Report (AR4). Ipcc, 1, p. 976 doi: ISSN: 02767783. €b, A., Berthier, E., Nuth, C., et al., 2012. Contrasting patterns of early twenty-firstK€ aa century glacier mass change in the Himalayas. Nature 488, 495–498. https://doi. org/10.1038/nature11324. Krishna, A.P., 2005. Snow and glacier cover assessment in the high mountains of Sikkim Himalaya. Hydrol. Process. 19, 2375–2383. https://doi.org/10.1002/hyp.5890. Kulkarni, A.V., Rathore, B.P., Singh, S.K., 2010. Distribution of seasonal snow cover in central and western Himalaya. Ann. Glaciol. 51, 123–128. https://doi.org/10.3189/ 172756410791386445. NASA, 2008. Landsat 7 Science Data Users Handbook. Science 80, 186. Otsu, N., 1979. A threshold selection method from gray-level histograms. IEEE Trans. Syst. Man Cybern. 9, 62–66. https://doi.org/10.1109/TSMC.1979.4310076. Peel, B.L., Finlayson, B.L., McMahon, T.A., 2007. Updated world map of the KoppenGeiger climate classification.pdf. Hydrol. Earth Syst. Sci. 11, 1633–1644. https:// doi.org/10.5194/hess-11-1633-2007. Quincey, D.J., Braun, M., Glasser, N.F., et al., 2011. Karakoram glacier surge dynamics. Geophys. Res. Lett. 38 https://doi.org/10.1029/2011GL049004. Racoviteanu, A.E., Williams, M.W., Barry, R.G., 2008. Optical remote sensing of glacier characteristics: a review with focus on the Himalaya. Sensors 8, 3355–3383. https:// doi.org/10.3390/s8053355. Ren, J., Jing, Z., Pu, J., Qin, X., 2006. Glacier variations and climate change in the central Himalaya over the past few decades. In: Annals of Glaciology, pp. 218–222. Schmidt, S., Nüsser, M., 2009. Fluctuations of Raikot glacier during the past 70 years: a case study from the Nanga Parbat massif, northern Pakistan. J. Glaciol. 55, 949–959. https://doi.org/10.3189/002214309790794878. Tahir, A.A., 2011. Impact of climate change on the snow covers and glaciers in the Upper Indus River Basin and its consequences on the water reservoirs (Tarbela Dam) – Pakistan. Thesis. UNEP, 2008. Global Glacier Changes : Facts and Figures, WGMS, 88. Xiao, X., Shen, Z., Qin, X., 2001. Assessing the potential of VEGETATION sensor data for mapping snow and ice cover: a Normalized Difference Snow and Ice Index. Int. J. Rem. Sens. 22, 2479–2487. https://doi.org/10.1080/01431160119766. Xiao, L., Ouyang, H., Fan, C., 2019. An improved Otsu method for threshold segmentation based on set mapping and trapezoid region intercept histogram. Optik (Stuttg) 196, 163106. https://doi.org/10.1016/j.ijleo.2019.163106. Young, G.J., Hewitt, K., 1990. Hydrology research in the Upper Indus Basin , Karakoram Himalaya , Pakistan. Hydrol. Mt Areas 139–152. Yu, W., Yang, Y., Savitsky, A., et al., 2013. The Indus Basin of Pakistan: Impacts of Climate Risks on Water and Agriculture. World Bank, Washington, DC. Zemp, M., Haeberli, W., Hoelzle, M., Paul, F., 2006. Alpine glaciers to disappear within decades? Geophys. Res. Lett. 33 https://doi.org/10.1029/2006GL026319. Zhang, G., Yao, T., Xie, H., et al., 2015. An inventory of glacial lakes in the Third Pole region and their changes in response to global warming. Global Planet. Change 131, 148–157. https://doi.org/10.1016/j.gloplacha.2015.05.013. Zimmermann, M., Bichsel, M., Kienholz, H., 1986. Mountain hazards mapping in the Khumbu Himal, Nepal, with prototype map, scale 1:50 000. Mt. Res. Dev. 6, 29–40. https://doi.org/10.2307/3673338.
Hewitt, K., 2005. The Karakoram anomaly? Glacier expansion and the ‘elevation effect,’ Karakoram Himalaya. Mt. Res. Dev. 25, 332–340 doi: 10.1659/0276-4741(2005) 025[0332:TKAGEA]2.0.CO;2. Ageta, Y., Higuchi, K., 1984. Estimation of mass balance components of a summeraccumulation type glacier in the Nepal Himalaya. Geogr. Ann. Phys. Geogr. 66, 249–255. https://doi.org/10.2307/520698. Alford, D.L., Archer, D.R., Bookhagen, B., et al., 2016. Monitoring of Glaciers , Climate , and Runoff in the Hindu Kush-Himalaya Mountains, vol. 1. World Bank Group, Washington, D.C., p. 152 Archer, D., 2003. Contrasting hydrological regimes in the Upper Indus Basin. J. Hydrol. 274, 198–210. https://doi.org/10.1016/S0022-1694(02)00414-6. Archer, D.R., Fowler, H.J., 2004. Spatial and temporal variations in precipitation in the Upper Indus Basin, global teleconnections and hydrological implications. Hydrol. Earth Syst. Sci. 8, 47–61. https://doi.org/10.5194/hess-8-47-2004. Ashraf, A., Naz, R., Roohi, R., 2012. Glacial lake outburst flood hazards in Hindukush, Karakoram and Himalayan Ranges of Pakistan: implications and risk analysis. Geomatics, Nat. Hazards Risk 3, 113–132. https://doi.org/10.1080/ 19475705.2011.615344. Atif, I., Mahboob, M.A., Iqbal, J., 2015. Snow cover area change assessment in 2003 and 2013 using MODIS data of the Upper Indus Basin, Pakistan. J. Himal. Earth Sci. 48, 117–128. Bajracharya, S., Shrestha, B., 2011. The status of glaciers in the Hindu Kush–Himalayan region. Mt. Res. Dev. 33, 114–115. Bajracharya, B., Shrestha, A.B., Rajbhandari, L., 2007. Glacial Lake outburst Floods in the Sagarmatha region. Mt. Res. Dev. 27, 336–344. https://doi.org/10.1659/mrd.0783. Bishop, M.P.J., Shroder, J.F., 2010. Satellite image Atlas of glaciers of the world— glaciers of Asia— GLACIERS OF Pakistan—an overview. US Geol Surv Prof Pap, 1386–F–4. Bolch, T., Buchroithner, M.F., Peters, J., et al., 2008. Identification of glacier motion and potentially dangerous glacial lakes in the Mt. Everest region/Nepal using spaceborne imagery. Nat. Hazards Earth Syst. Sci. 8, 1329–1340. https://doi.org/10.5194/ nhess-8-1329-2008. Burns, P., Nolin, A., 2014. Using atmospherically-corrected Landsat imagery to measure glacier area change in the Cordillera Blanca, Peru from 1987 to 2010. Remote Sens. Environ. 140, 165–178. https://doi.org/10.1016/j.rse.2013.08.026. Copland, L., Sylvestre, T., Bishop, M.P., et al., 2011. Expanded and recently increased glacier surging in the Karakoram. Arctic Antarct. Alpine Res. 43, 503–516. https:// doi.org/10.1657/1938-4246-43.4.503. Dyurgerov, M.B., Meier, M.F., 2005. Glaciers and the Changing Earth System: A 2004 Snapshot. Instaar. Fowler, H.J., Archer, D.R., 2005. Hydro-climatological variability in the Upper Indus Basin and implications for water resources. Iahs 295, 131–138. Fowler, H.J., Archer, D.R., 2006. Conflicting signals of climatic change in the Upper Indus Basin. J. Clim. 19, 4276–4293. https://doi.org/10.1175/JCLI3860.1. Gardelle, J., Berthier, E., Arnaud, Y., 2012. Slight mass gain of Karakoram glaciers in the early twenty-first century. Nat. Geosci. 5, 322–325. https://doi.org/10.1038/ ngeo1450. Hall, D.K., Foster, J.L., Verbyla, D.L., et al., 1998. Assessment of snow-cover mapping accuracy in a variety of vegetation-cover densities in central Alaska. Remote Sens. Environ. 66, 129–137. https://doi.org/10.1016/S0034-4257(98)00051-0. Hewitt, K., 2007. Tributary glacier surges: an exceptional concentration at Panmah glacier, Karakoram Himalaya. J. Glaciol. 53, 181–188. https://doi.org/10.3189/ 172756507782202829. Hewitt, K., 2011. Glacier change, concentration, and elevation effects in the Karakoram Himalaya, Upper Indus Basin. Mt. Res. Dev. 31, 188–200. https://doi.org/10.1659/ MRD-JOURNAL-D-11-00020.1.
19
Contents lists available at ScienceDirect
Remote Sensing Applications: Society and Environment journal homepage: http://www.elsevier.com/locate/rsase
Change detection of glaciers and snow cover and temperature using remote sensing and GIS: A case study of the Upper Indus Basin, Pakistan Abdul Jabbar a, b, Arsalan Ahmed Othman c, *, Broder Merkel a, Syed E. Hasan d a
Department of Hydrogeology, Institute of Geology, Technische Universit€ at Bergakademie Freiberg, Gustav-Zeuner Str. 12, 09599, Freiberg, Germany Baugrundbüro Dr. Matthias Mokosch, Dresdner Str. 39, 01683, Nossen, Germany c Iraq Geological Survey, Sulaymaniyah Office, Sulaymaniyah, Iraq d Department of Earth & Environmental Sciences, University of Missouri, Kansas City, MO, USA b
A R T I C L E I N F O
A B S T R A C T
Keywords: Snow cover Upper Indus Basin Himalaya NDSI Land surface temperature
Glaciers and snow cover in the Himalayas, Karakoram, and Hindu-Kush (HKH) mountain ranges are expected to decrease, resulting in serious consequences in terms of water availability for downstream population. This study is an attempt to detect temporal trend of snow cover in the Upper Indus Basin (UIB), in northern Pakistan, with respect to temperature change, spatial variation in land surface temperature, and to evaluate the accuracy of Landsat imagery for detecting land surface temperature. Landsat images, acquired for the period 1993 to 2019 and the Normalized Difference Snow Index (NDSI), were used to map the snow cover. Maximum snow cover of about 20,675 km2 was observed in April 2017 followed by 1996, and 1999, when it was about 18,800 km2; the minimum of about 14,500 km2 was recorded in April 1993, 2011 and 2018. The results show that there has been a very slight increase in snow cover area with time. The maximum elevation of snow-covered area, extracted from the Digital Elevation Model of Shuttle Radar Topography Mission (SRTM DEM) in the study area, is 7889 m while the minimum elevation is 991 m. Freezing temperature occurs above 3000 m elevation, indicating that these areas are covered by snow and glaciers for much of the year. Temperatures, of 20 � C to 25 � C, were present only at a few places–at the very high elevation of above 6000 m. We found inconsistencies in the Landsat temperature data, which were corrected using in-situ daily temperature data. Validation of Landsat extracted temperatures before and after using the in-situ daily temperature data improved the accuracy and removed the inconsistencies. We recommend use of in-situ diurnal temperature data to correct Landsat extracted temperature data.
1. Introduction Recent decades have witnessed significant impacts from climate change. Adverse effects of climate change have been quite pronounced on people, wildlife, and ecosystems, including the cryosphere. Since the 18500 s, receding of many highland glaciers all over the world have been observed and an exponential increase of receding during the last 30 years (Ren et al., 2006; UNEP, 2008; Zemp et al., 2006). It is estimated that fresh water from snow and glaciers provide fresh water supply to one-half of world’s population. Receding glaciers and shrinking ice sheets have adversely affected global fresh water supplies. The Himalayan ranges represent the most glaciated region in the world except polar ice sheets and Greenland (IPCC, 2007; Young and
Hewitt, 1990), but precise data on their characteristics and dynamics are lacking (Dyurgerov and Meier, 2005). These ranges not only play a vital role in providing water supply to about 1.3 billion people living in downstream basins of the 10 major rivers in Asia, but also in the global hydrological cycle, atmospheric circulation, and ecosystems services (Bajracharya and Shrestha, 2011). Recent research conducted by Gar €€ delle et al. (2012) and Ka ab et al. (2012) support the fact that while majority of glaciers are receding in the Himalayan range, mass gain or equilibrium condition has occurred in the Karakoram mountain glaciers as a result of favorable climatic conditions in the past 10 years (Gardelle €b et al., 2012). No need to repeat this reference in the et al., 2012; K€ aa same sentence; delete]. An increasing trend in winter precipitation (Fowler and Archer, 2006) and decreasing pattern in summer
* Corresponding author. E-mail addresses: [email protected] (A. Jabbar), [email protected] (A.A. Othman), [email protected] (B. Merkel), [email protected] (S.E. Hasan). https://doi.org/10.1016/j.rsase.2020.100308 Received 19 April 2019; Received in revised form 26 January 2020; Accepted 19 March 2020 Available online 24 March 2020 2352-9385/© 2020 Elsevier B.V. All rights reserved.
A. Jabbar et al.
Remote Sensing Applications: Society and Environment 18 (2020) 100308
Fig. 1. Study area, located in the extreme north of Pakistan, borders with China to the north, India to the east, and Afghanistan to the north-west.
temperature (Schmidt and Nüsser, 2009) has been observed. In the Karakoram, over 80% of the ice cover are located between 4000 m and 5500 m elevations and the extensive debris cover at the ablation zone reduces the glacier exposure which makes them insensitive to changing climate (Hewitt, 2011). Shrinking and thinning of glaciers is common in most of the Hima layan, Karakoram, and Hindu-Kush (HKH) regions, leading to formation of glacial lakes. Continued melting could lead to Glacial Lake Outbursts Floods (GLOFs) (Bolch et al., 2008; Zimmermann et al., 1986). Several studies in the HKH region have reported a large number of GLOFs due to glacier recession (Ashraf et al., 2012; Bajracharya et al., 2007; Zhang et al., 2015). Surging of a large number of glaciers have been reported in the Karakoram since 1980 (Copland et al., 2011; Hewitt, 2007), and the rate has doubled since 1990, which is likely due to positive mass budget (Copland et al., 2011). The study was designed to: (1) monitor the temporal trends of snow cover for the last 27 years in the Upper Indus Basin (UIB); and (2) evaluate the accuracy of Landsat imagery for detecting land surface temperature.
Batura (59.8 km); and Hisper (53.1 km) (Atif et al., 2015). The mean elevation of UIB is 4750 m most of which (about 60%) is above 4500 m elevation and only 12% is above 5500 m, which is the abode of per manent snow and glaciers (Tahir, 2011). The mighty Indus River along with all its major tributaries including Kabul, Jehlum, Chenab, Ravi, and Sutlej rivers originate in these mountainous regions; these are the main source of water to about 200 million people living in the Indus valley (Immerzeel et al., 2010). Ac cording to Archer and Fowler (2004) and Hewitt et al. (1989), the Indus River receives more than 80% of its annual flow from less than 20% of the extensively glaciated and snow covered zone above 3500 m eleva tion. Recent modelling carried out by Yu et al. (2013) suggests that the total contribution of the UIB to the Indus River system is about 135.68 km3, of which ~18% (23.44 km3) comes from the glacier melt, and ~82% comes from the winter snowpack. €ppen-Geiger classification, the climate of Upper According to the Ko Indus Basin is characterized as cold desert climate (BWk) with warm dry summer and cold winter with low precipitation. This climate type can be found in temperate regions that lie in the rain shadow of high mountains (Peel et al., 2007). As the influence of monsoon winds lessens north-westwards, about 90% of the UIB falls in the Himalayan rain shadow zone (Immerzeel et al., 2009). The annual precipitation of the arid valley floor at the UIB is 150 mm–200 mm, whereas in the Gilgit River valley at 4400 m altitude, precipitation increases to 600 mm; and, according to glaciological studies, the accumulation at 5500 m reaches 1500 mm–2000 mm (Fowler and Archer, 2005). There is a huge inequality in the climate setting of the UIB compared to the eastern Himalaya (Fowler and Archer, 2006; Young and Hewitt, 1990). The climate diversity across the entire Himalayas can be perceived from the fact that the eastern part receives more than 3000 mm and western part gets less than 300 mm precipitation annually that has direct and indirect effects on the rivers’ discharge (Immerzeel et al., 2009). The spatial and seasonal snow accumulation also varies widely in the region. The eastern and central Himalayas receive maximum snowfall from summer monsoon (Ageta and Higuchi, 1984), while in the northwest (Karakoram and UIB), westerly circulations winter most
2. Study area The UIB is located in the extreme north of Pakistan along the Chinese border (Fig. 1); its catchment extends upstream from Tarbela Dam reservoir and covers an area variously estimated to be about 206,000 km2 (Tahir, 2011); 220,000 km2 (Yu et al., 2013); and 200,677 km2 (Immerzeel et al., 2009). The UIB is home to some of the world’s tallest mountain ranges in the Himalayas, Karakoram, and Hindu-Kush (HKH) ranges with several peaks attaining heights of >7000 m. Bishop and Shroder (2010) report that out of the total 30 highest summits in the world, 13 are located in northern Pakistan. The world’s second tallest peak K2 (8611 m), the ninth highest, Nanga Parbat (6125 m) and Tirich Mir (7690 m) are all located here. Some well-known glaciers in the HKH Pakistan (in other words UIB) are: Siachen (75 km)–longest in the Kar akoram and the world’s second longest glaciers outside the Polar Re gions; Biafo (67.9 km) third longest in the world; Baltoro (62.1 km); 2
A. Jabbar et al.
Remote Sensing Applications: Society and Environment 18 (2020) 100308
winter snow accumulation comes from westerly circulating winds (Alford et al., 2016; Fowler and Archer, 2005; Quincey et al., 2011). Hewitt (2011) claims that the Himalayan glaciers, especially the central Karakoram, gain an extensive amount of snow from avalanches moving down the adjacent steep slopes. The topography of UIB shows marked variability: from the lowest elevation of less than 1000 m in the plain areas, near the Tarbela and Mangla dam reservoirs, to over 8000 m at a number of peaks (Tahir, 2011; Yu et al., 2013), and, according to Immerzeel et al. (2009), the lowest elevation in the UIB is 335 m and the highest is 8238 m. The lower Indus basin irrigation system is mainly dependent on the melt water from glaciers and snow cover that provide uninterrupted supply round the year, which has turned it into the world’s largest integrated irrigation system (Fowler and Archer, 2005). The inhabitants in the lower Indus Basin are mainly dependent on snow and glaciers melt water for agriculture purposes (Archer, 2003).
the region of interest. The cloud cover in some areas of the images taken in 1996, 1999 and 2015 was manually digitized and the shape files were created in Arc Map. Then, all the shape files were combined in one file through Boolean operation comprising all of the cloud cover area. Next, a union was created between the cloud cover file and the intersect area that resulted in a single file having a well-defined area for all scenes excluding cloud cover. This cloud cover shapefile was converted from vector to raster format due to the images in raster. Then, this raster cloud cover file and the raster images containing the threshold value (separate snow and non-snow area) were stacked together with each image one by one. Finally, a new file containing three values was prepared: -1 showing cloud cover area excluded from calculation, 0 representing non-snow area, and 1 displaying snow cover.
3. Methodology
The Normalized Difference Snow Index (NDSI) algorithm was used for mapping snow covered area which was calculated from Equation (1):
3.4. Normalized Difference Snow Index (NDSI) threshold selection
3.1. Scene selection and acquisition
NDSI ¼
Landsat images with 30 m pixel resolution captured in the month of April for the period 1993 till 2019 were searched in the United States Geological Survey (USGS) archive (http://glovis.usgs.gov/) and were received free of cost from the Earth Resources Observation and Science (EROS) center. The repeat period of Landsat TM (16 days) and its highaltitude position make the acquisition of cloud free images extremely difficult for a given time period. April was selected due to the avail ability of best-suitable images on the basis of image quality and cloud cover. Only images with excellent quality and less than 10% cloud cover were selected for this study. Data from Landsat sensors TM, ETMþ and OLI/TIRS were used in this study. In addition, eleven Operational Land Imager (OLI) scenes have also been used for different months between 2017 and 2018. Details are listed in Appendix A. The Landsat ETM þ data was missing since May 2003 due to scan line corrector failure (SLC off); also, in some other years owing to maximum cloud cover the data was not deemed suitable for inclusion in this study. The missing gap, to some extent, was filled by acquiring Landsat 5 TM images. The benefi cial characteristic of Landsat images from the EROS center is that they are already in (Level 1T format) geo-referenced, and orthorectified (without topographic distortion). The region of interest is positioned in zone 43 north and the path and row of the image are 149 and 35 respectively. The projected coordinate system for all images is UTM (Universal Transverse Mercator, WGS 1984 datum).
ðTM2 TM5 Þ ðTM2 þ TM5 Þ
(1)
The Atmospheric and Topographic Correction (ATCOR) model was employed to remove the atmospheric effects, sun illumination, and sensor viewing geometry from the images. In order to perform atmo spheric correction, the ATCOR model needs description of the compo nents in the atmospheric profile and data about the image of the area e.g. the exact time of image acquisition day, month, year; sensor type, calibration file, solar zenith and azimuth. The calibration file for at mospheric correction was created from the metadata file provided along with the image.
Where TM2 represents the reflectance in the band 2 (0.52 μm–0.6 μm) and TM5 refers to reflectance in band 5 (1.55 μm–1.75 μm) of the Landsat 5 and 7. In Landsat 8, band 3 and band 6 are used for NDSI because of the similar spectral bandwidths to the band 2 and band 5 of Landsat 5 and 7. Burns and Nolin (2014) state that atmospherically uncorrected data illustrate a conspicuous change in calculating the debris free glacier area and underestimates the glacier debris free area if the similar threshold is applied. For this reason, all scenes were atmo spherically corrected before selecting the threshold value. To get better results, the trial-and-error threshold method has been used instead of the automatic threshold method of Otsu, (1979), which is based on a one-dimensional (1D) histogram considering the spatial correlation between pixels. Therefore, its anti-noise capability is not satisfactory (Xiao et al., 2019). According to Hall et al. (1998), NDSI threshold value of �0.4 is an optimal value to map a snow pixel. Like wise, several investigators (Kulkarni et al., 2010; Xiao et al., 2001) used the threshold value of �0.4 while Burns and Nolin (2014) considered threshold value of �0.42 to be the best for debris free glacier areas after comparing the result of Landsat images to the IKONOS-2 and QuickBird high resolution imagery. Racoviteanu et al. (2008) report that NDSI is a very effective method for complex topography with highly shadowed regions whereas, Burns and Nolin (2014) applied the ratio technique TM4/TM5 and TM3/TM4 in examining the glacier area fluctuation in the Cordillera Blanca, Peru, but finally decided for NDSI because of good detection in shadows. Similarly (Krishna, 2005; Racoviteanu et al., 2008), also utilized the NDSI in the Himalayas for obtaining optimum results for shadowed areas. A number of threshold values were found in the literature (e.g., �0.42, �0.5, and �0.6) and all of them were tested; and results were carefully compared and analyzed visually to ascertain the impact on mapping the snow cover in shadowed regions. Based on these considerations, we selected the threshold value, � 0.42 (Burns and Nolin, 2014) because it enabled us to accurately determine snow pixels in shadow regions.
3.3. Cloud area masking
3.5. Classification of the area
Some images have missing data on the border, and others have an undefined area at the sides (all images did not cover exactly the same area), so frames in the form of shape files in ArcMap 10.4 were created for all images. Then an intersect area was developed for all scenes in order to get a well-defined area and to remove the error on the border of the images. Furthermore, a subset with the aid of the intersect area was created for all images to obtain exactly the same areas for all images in
The selected threshold value � 0.42 of all scenes was added together via raster calculator function in ArcMap that provided a combined file containing numbers ranging from 1 to 12. Where 1 exhibits the cloud cover region (excluded from the calculation), 0 displayed non-snow area and 1 to 12 (because there are 12 images) showed snow covered region. A very simple formula (the file containing the addition of all threshold value ¼ 12) was used in raster function that generated the snow cover
3.2. Atmospheric correction
3
A. Jabbar et al.
Remote Sensing Applications: Society and Environment 18 (2020) 100308
Fig. 2. Workflow procedures for snow/ice cover mapping and temperature extraction.
area (pixels) present in all images from 1993 to 2019 and was named as permanent (or minimum) snow area. Then, the function (the file con taining the addition of all threshold value > 0) & (the file containing the addition of all threshold value � 11) calculated the non-permanent (maximum) snow area that was not present in all images. Next, the files “permanent” (minimum) and “non-permanent” (maximum) were added together and the results presented as the total snow cover area in UIB since 1993 till 2019. Additionally, the positive and negative snow cover was calculated by performing the following function for each image i.e. (the threshold values file of each year 1994 – permanent (minimum) snow cover) ¼ Positive area and; the threshold values file of each year snow covers i.e. (non-permanent (maximum) snow cover – the threshold value of 1993) ¼ Negative area. Thus, positive area refers to the snow cover pixels that are present in 1993 images besides the permanent (minimum) snow area, whereas negative area refers to the snow cover area that is not present in 1993 images but existed in another images.
(Equation (2)). The atmospheric correction was performed using ATCOR Model for the thermal bands of each image before converting to radiance and radiance to surface temperature. The procedure given in NASA’s Landsat 7 Science Data Users Handbook (NASA, 2008), was used to estimate the surface temperature. A model was created in ERDAS IMAGINE version 16 and all the thermal bands were processed through this model. � � LMAXλ LMINλ Lλ ¼ (2) * ðQCAL QCALMINÞ þ LMINλ QCALMAX QCALMIN Where, Lλ ¼ spectral radiance at the sensor’s aperture in watts/(m2 * ster *
μm)
LMAXλ ¼ maximum detected spectral radiance in watts/(m2 * ster * μm) LMINλ ¼ minimum detected spectral radiance in watts/(m2 * ster * μm) QCAL ¼ the quantized calibrated pixel value in digital number (DN) QCALMAX ¼ the maximum quantized calibrated pixel value (255) QCALMIN ¼ the minimum quantized calibrated pixel value (1)
3.6. Land surface temperature estimation In order to compute temperature from the Landsat thermal bands, the digital number (DN) value of the images were converted to radiance and then the brightness temperature was estimated from the radiance 4
A. Jabbar et al.
Remote Sensing Applications: Society and Environment 18 (2020) 100308
Fig. 3. snow cover, non-snow cover, positive and negative areas during the years 1996 and 2015. Clouds cover areas are highlighted in yellow circle. (For inter pretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)
LMAXλ and LMINλ are the spectral radiances in digital numbers ranging from 1 to 255 for each band having different values and are given in the header file of the image. The radiance was converted to brightness temperature with Equation (3). K2 � T¼ � K1 ln Lλ þ 1
273:15
T ¼ effective at-satellite temperature in Kelvin K2 ¼ calibration constant 2 in Kelvin K1 ¼ calibration constant 1 in watts/(m2 * ster * μm) Lλ ¼ spectral radiance at the sensor’s aperture in watts/(m2 * ster * μm)
(3)
K2 and K1 are the pre-launch calibration constants provided in the header file of the image. Temperature records of two meteorological stations were obtained for validation. Point shape files for Gilgit and Skardu meteorological stations were created to validate the temperature
Where:
Fig. 4. Change of snow cover with cloud cover area for the period 1993 to 2019. A statistically significant trend is not perceptible. 5
A. Jabbar et al.
Remote Sensing Applications: Society and Environment 18 (2020) 100308
Fig. 5. Total snow and ice cover without cloud cover in the Upper Indus Basin in April for the period 1993 to 2019. A statistically significant trend is not perceptible.
The lowest snow cover area of about 14,500 km2 was recorded in April 1993 and 2011. Results show that there has been a very slight increase in the snow cover area with time. However, the non-parametric Mann Kendall (MK) test has been applied for the snow cover time series. The MK test shows that the 2-sided p-value is > 0.05, therefore, the trend is statistically not significant. Images showing change in snow cover for the years 1993, 1994, 1996, 1999, 2000, 2002, 2011, 2015, 2017, 2018, and 2019 are shown in Appendix D. Fig. 7 shows the positive and negative snow and ice cover areas in April for the period 1993 to 2019 (see heading 3.5 classification of the area, for positive and negative area). The year 1993 has the highest negative snow area of about 6300 km2 followed 5500 km2 for the year 2011. The year 2017 has the highest positive snow area of about 8000 km2. The years 1996, 1999, 2002, and 2003 have almost similar positive snow cover areas of about 5000 km2 and also nearly equal negative snow cover areas of about 2000 km2 with minor differences. Moreover, in 2000 and 2015, the positive and negative snow areas were almost equal, i.e. about 3500 km2 (Appendix B). Fig. 8 shows the change of snow cover for the period March 26, 2017 to November 24, 2018. The cloud cover in three images was masked to get only the snow cover area. The highest snow cover of about 6844 km2 was found in April 2017. The lowest snow cover area of about 20,676 km2 was recorded in August 2017.
of that particular pixel. One of these stations (better correlated with the temperatures extracted from Landsat data) was used to correct the bias of the temperature extracted from Landsat data. The other one was used to validate the corrected data. Following surface temperature estimation, the SRTM (1 arc or 30 m resolution) Version-3 DEM was downloaded from the USGS earth ex plorer website (https://earthexplorer.usgs.gov/). Elevation values for the study area, extracted from the SRTM DEM, gave maximum elevation of up to 7889 m and the minimum 991 m. Fig. 2 explains the workflow used for snow cover mapping and land surface temperature estimation. 4. Results 4.1. Results before exclusion of cloud cover area The above-mentioned procedure was performed once for all images, and layout were also created. However, as mentioned before, there were some cloud cover in three images that resulted in an error during the calculation of positive and negative area (Appendix B). For example, the highlighted (red) area in yellow circle is not a negative area but a cloud cover captured in 1996 image, that is represented as negative area and also the cloud cover area in 2015 image is shown as a positive area which is an error (see heading 3.5 classification of the area, for expla nation of positive and negative area). Finally, it was decided to exclude the cloud cover area and repeat the entire procedure for all images and re-check the results. Fig. 3 shows the cloud cover areas highlighted in a yellow circle. This area was estimated as positive snow cover area in one image and negative in another image. This error was removed by excluding/ masking the cloud cover areas from calculation. Fig. 4 shows the change of snow cover including cloud cover in April during the last 27 years. The cloud cover in three images was masked from the study area to get only the snow cover area. Figs. 5 and 6 represent the total snow and ice cover without cloud in the Upper Indus Basin in April for the period 1993 to 2019 (Appendix C). The highest snow cover of approximately 20,675 km2 was found in April 2017 fol lowed by 1996, and 1999 where the snow cover was about 18,800 km2.
5. Land surface temperature Fig. 9A shows the temperature extracted from the Landsat image. The dark brown color represents higher temperature and the dark blue lower temperature. Very low temperatures ( 20 � C to 25 � C) occur at very high altitudes of 6000 m to 7000 m. High temperatures reach 30 � C–40 � C at the valleys floor between 1000 m and 2500 m elevations. The very low and the very high temperature areas comprise a very small part of the total area–very low temperature of 30 � C is recorded in only 8 pixels, and very high temperature of 47 � C is recorded in only 4 pixels. The very low temperatures seem to be realistic because it is found at a very high altitude. Higher temperature, on the other hand, is exagger ated and contains an error. Fig. 9B displays the topographic variation in 6
A. Jabbar et al.
Remote Sensing Applications: Society and Environment 18 (2020) 100308
Fig. 6. Cloud cover, snow cover, non-snow cover, positive and negative area during different years.
the UIB region. The SRTM DEM of the Himalayas still has voids or data gaps. There is a very small area (white spots), on the eastern side of the study area where elevation data is missing. The maximum and minimum temperatures recorded at 8:30 a.m. on April 29, 2003 were used to verify temperature variations in the study area to display temperature changes in relation to the elevation (Ap pendix F). The altitude of meteorological observation station at Gilgit is 1460 m and at Skardu is 2317 m (Appendix F). Fig. 10 shows the maximum insitu temperature at the Gilgit and Skardu meteorological observation stations, along with temperature extracted from the Landsat imagery for six days (April 17, 1996; April 26, 1999; April 12, 2000; April 29, 2003; April 27, 2011; and April 22, 2015). The result does not match. In Gilgit, the average difference with the in-situ data is 4.33 � C. The Landsat re sults on an average are about 4.33 � C lower than the in-situ data. While, in Skardu, the average difference is 8.32 � C and the Landsat extracted results are about 8.32 � C higher than the in-situ temperature data
(Fig. 6). The coordinates of Gilgit and Skardu Metrological stations are shown in (Fig. 5). Some differences in temperature could be attributed to the difference in time the Landsat image was captured, and recording of the in-situ data. The maximum temperature data was used to validate the results due to lack of hourly data. The Landsat data was captured at about 9:00 a.m. which was close to the maximum temperature. So, this could be one of the reasons for the difference in both results (Gilgit and Skardu meteorological stations). In case of Gilgit, the results look satisfactory because the Landsat data is captured 3 h earlier than the maximum insitu temperature. So, the Landsat temperature may be close to the insitu temperature data at that time. Fig. 11 shows the maximum, mean and minimum temperature at the Gilgit and Skardu stations and Fig. 12 the sum of monthly precipitation data for April, taken from Pakistan Meteorological Department for the period 1990 to 2015. Fig. 11 shows the April monthly maximum, mean and minimum temperature at Gilgit and Skardu meteorological stations. There are 7
A. Jabbar et al.
Remote Sensing Applications: Society and Environment 18 (2020) 100308
Fig. 7. Positive and negative snow and ice cover area.
Fig. 8. Change of snow cover from March 26, 2017 to November 24, 2018.
large variations from year to year but the highest temperautre at both stations was recorded in 2007. Temperature at Gilgit exhibits a slight increase for maximum, mean and minimum temperautre while nearly stable conditions are observed in all cases at Skardu. However, the re sults are not significant statistically in both cases (Gilgit and Skardu). Significant variation in precipitation can be noted in Fig. 12 that shows the sum of April precipitation at Gilgit and Skardu meteorological stations for the period 1990 to 2015. High precipitation sums occur after
about 4–5 years (peaks) at both stations. In addition, the precipitation at Skardu area for most of the years is more than at Gilgit area. The highest precipitation recorded at Skardu in 2006 and at Gilgit in 1999 was about 145 mm and 90 mm, respectively. Overall, there are very small but statistically insignificant trends in precipitation at Gilgit and Skardu station.
8
A. Jabbar et al.
Remote Sensing Applications: Society and Environment 18 (2020) 100308
Fig. 9. Landsat temperature and elevation with meteorological stations location.
Fig. 10. Correlation of in-situ temperature and Landsat derived temperature for Gilgit and Skardu meteorological observatory.
6. Discussion and conclusions
acquisition of Landsat scenes even in a fixed week is impossible due to cloud cover, so some flexibility is needed for using suitable scenes. Considering these factors, the Landsat scenes during the month of April were acquired to examine the trend of snow cover for the period 1993 to 2019. Gathering reliable information regarding glacier/snow cover thick ness was not possible due to lack of data on high resolution and precise DEM, which are a very significant factors for investigating changing conditions. The SRTM DEM is not reliable for snow cover thickness investigation due to 9 m–16 m and greater vertical error, especially in high mountain region. Temperature extracted from Landsat was validated using the in-situ daily temperature data. Some differences are due to differences in time of recording Landsat data and the in-situ data. Maximum temperature data was used to validate the results due to lack of hourly data and
The Himalayas, the world’s highest mountain chain, is shrouded in cloud cover throughout the year and the chances to acquire images without cloud cover are slim. Several scenes for the study area at the required time period were found to have huge cloud cover, over 40%, which makes them inappropriate for time series analysis of snow cover monitoring. In this study, we used images having low cloud cover of about 10%. The cloud cover zones were masked from the study area that reduced the total area. The Landsat data coverage for the Upper Indus Basin, Pakistan since the 1980s–1990s is very poor, which limits the detailed investigation of glaciers and snow cover condition back in time. The common practice in remote sensing for time series studies of snow and ice is to select images on a fixed date of a particular month within each year. However, the 9
Remote Sensing Applications: Society and Environment 18 (2020) 100308
A. Jabbar et al.
Fig. 11. Monthly April maximum, mean and minimum temperature at Gilgit and Skardu meteorological stations (Source: Pakistan Meteorological Department).
Fig. 12. Variation of the sum of April precipitation over time at Gilgit and Skardu meteorological stations (Source: Pakistan Meteorological Department).
The largest snow cover of about 20,675 km2 was found in April 2017, followed by 1996 and 1999 when the snow cover was about 18,800 km2, while the lowest snow cover of about 14,500 km2 was recorded in April 1993, 2011, and 2018. The result shows that the snow cover in the month of April has slightly declined during 1993–2015. The year 1993 has the largest negative snow area, about 6300 km2 followed by the year 2011, when it was 5500 km2. The years 1996, 1999, 2002, and 2003 have almost similar positive snow area of about 5000 km2 and also nearly equal negative snow area of about 2000 km2, with minor dif ferences. Moreover, in 2000 and 2015, the positive and negative snow area (see heading classification of the area, for explanation of positive and negative area) was identical, about 3500 km2 (Appendix B). According to the SRTM DEM, the study area has a maximum eleva tion of 7889 m while the minimum elevation is 991 m. The maximum temperature reached about 42 � C while the minimum temperature dropped to approximately 25 � C. Temperatures of 0 � C and lower occurred above 3000 m elevation and decreased with increasing elevation; however, at an altitude of 7000 m all regions had temperature below 0 � C. The areas above 3000 m had temperature below 0 � C, which indicate that these areas are covered with snow and glaciers. As far as we could ascertain, there are no studies on estimating the temperature of snow cover with Landsat in the Upper Indus Basin area that we could use
because the Landsat data was captured at about 9:00 a.m., which was close to the maximum temperature. The lack of hourly climate data for the time when the satellite captured the data inhibited comparison of the satellite imagery to the in-situ data at the same time. Thus, in the case of Gilgit the average departure was 4.33 � C lower than the in-situ obser vation; and at Skardu the average departure was about 8.32 � C higher than the in-site temperature data. In both cases the in-situ temperature data and temperature data extracted from Landsat had a time difference of about 3 h. In case of Gilgit, the results are acceptable because the Landsat data was captured 3 h earlier than the maximum in-situ tem perature (more accurate than Skardu, where R2 is > 0.73). So, the Landsat extracted temperature may be close to the in-situ temperature at that time. The very low temperature values seem satisfactory because it is observed at a very high altitude. On the other hand, the higher tem peratures are exaggerated and are erroneous. Temperature at Gilgit exhibits a slight increase for maximum, mean and minimum temper autre while in the case of Skardu almost stable condition is observed in all cases. Therefore, we used the equation plotted in Fig. 10 of Gilgit station to correct the bias in the temperatures extracted from Landsat. We used the Skardu station to validate the bias corrected temperatures (extracted from Landsat data), which reduced the average error at Skardu station from 8.32 � C to 3.25 � C. 10
A. Jabbar et al.
Remote Sensing Applications: Society and Environment 18 (2020) 100308
to validate or compare the results. €€ According to Gardelle et al. (2012); Hewitt (2005); Ka ab et al. (2012), the glaciers in the Karakoram mountain ranges have increased in volume and have remained stable during the past 10 years due to cli matic conditions; Hewitt (2011) reports that more than 80% of ice cover is situated between 4000 m and 5000 m. Thus, the claim of these in vestigators seems to be true for higher elevation with freezing temperature. Topography of the Upper Indus Basin varies greatly (7889 m at the highest to 991 m at the lowest points in this study) whereas the tem perature and precipitation data are recorded at the very low Gilgit and Skardu meteorological observatories at the valleys floor, 1460 m and 2317 m, respectively. Therefore, they cannot be valid for the whole region with such a huge variation in topography and temperature. However, satellite data can be an ideal tool for evaluating temperature at a certain elevation, provided the Landsat or any other remote sensing data is precisely validated using in-situ temperature data because tem perature can be accurately estimated from the linear regression line for each elevation as shown in Fig. 10. In-situ temperature data recorded at the same time when the Landsat image was captured should be compared for validation and accuracy. Precise data about the glaciers/snow cover area and temperature esti mation would assist in understanding future climate-based changes. Further studies should be performed to find out the relation between snow cover estimation and retrieval of land surface temperature and to determine snow cover changes in relation to changes in temperature or precipitation.
relation to the development, writing, and publication of the article. All co-authors have seen and agree with the contents of the manuscript and there is no financial interest to report. We certify that the submission is an original work and has not been published nor has it currently been submitted for publication to another scientific journal, or being considered for publication elsewhere. Declaration of competing interest The authors declare no conflict of interest. CRediT authorship contribution statement Abdul Jabbar: Data curation, Formal analysis, Investigation, Methodology, Writing - original draft, Writing - review & editing. Arsalan Ahmed Othman: Conceptualization, Formal analysis, Investi gation, Methodology, Project administration, Supervision, Writing original draft, Writing - review & editing. Broder Merkel: Conceptu alization, Project administration, Supervision, Writing - review & edit ing. Syed E. Hasan: Writing - review & editing. Acknowledgements We are thankful to the United States Geological Survey Earth Re sources Observation and Science data center for providing Landsat im agery; the GEOSYSTEMS (Hexagon Geospatial) GmbH, Germany for ERDAS IMAGINE plus ATCOR software License; and the Pakistan Meteorological Department for providing climate data. This project would not have been possible without their valuable support.
Ethical statement The authors declare that all ethical practices have been followed in
Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi.org/10.1016/j.rsase.2020.100308. Appendix A. Landsat raw data Scene Type
Date
ID
Quality
Cloud Cover
TM5 TM5 TM5 TM5 TM5 TM7 TM7 TM5 OLI OLI OLI OLI OLI OLI OLI OLI OLI OLI OLI OLI OLI
1993/04/25 1994/04/28 1996/04/17 1999/04/26 2000/04/12 2002/04/10 2003/04/29 2011/04/27 2015/04/22 2017/03/26 2017/04/11 2017/08/01 2017/11/05 2017/12/07 2018/01/24 2018/02/09 2018/04/30 2018/08/04 2018/10/23 2018/11/24 2019/04/01
LT51490351993115ISP00 LT51490351994118ISP00 LT51490351996108ISP00 LT51490351999116AAA01 LT51490352000103AAA02 LE71490352002100SGS00 LE71490352003119EDC00 LT51490352011117KHC00 LC81490352015112LGN00 LC81490352017085LGN00 LC81490352017101LGN00 LC81490352017213LGN00 LC81490352017309LGN00 LC81490352017341LGN01 LC81490352018024LGN00 LC81490352018040LGN00 LC81490352018120LGN00 LC81490352018216LGN00 LC81490352018296LGN00 LC81490352018328LGN00 LC81490352019091LGN00
9 9 7 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9
1.1% 2.62% 08% 10% 10% 1.96% 0.98% 2.85% 14.69% 6.8% 1.06% 3.6% 3.28% 9.18% 3.7% 5.26% 4.36% 2.57% 7.38% 2.96% 6.61%
11
A. Jabbar et al.
Remote Sensing Applications: Society and Environment 18 (2020) 100308
Appendix B. Positive and negative snow cover area in each image Date
Positive area in pixels
Positive area in km2
Negative area in pixels
Negative area in km2
1993/04/25 1994/04/28 1996/04/17 1999/04/26 2000/04/12 2002/04/10 2003/04/29 2011/04/27 2015/04/22 2017/04/11 2018/04/30 2019/04/01
982,734 4,600,504 6,080,856 5,997,955 3,793,011 5,565,773 5,840,748 1,881,453 4,339,627 8,072,562 1,659,742 5,981,757
884.46 4140.45 5472.77 5398.16 3413.71 5009.20 5256.67 1693.31 3905.66 7265.31 1493.77 5383.58
7,028,045 3,410,275 1,929,923 2,012,824 4,217,768 2,445,006 2,170,031 6,129,326 3,671,152 36,243 634,583 222,059
6325.24 3069.25 1736.93 1811.54 3795.99 2200.51 1953.03 5516.39 3304.04 32.62 571.13 199.85
Appendix C. Total snow cover area in each image Date
Snow area in pixels
Snow area in km2
1993/04/25 1994/04/28 1996/04/17 1999/04/26 2000/04/12 2002/04/10 2003/04/29 2011/04/27 2015/04/22 2017/03/26 2017/04/11 2017/08/01 2017/11/05 2017/12/07 2018/01/24 2018/02/09 2018/04/30 2018/08/04 2018/10/23 2018/11/24 2019/04/01
15,883,117 19,500,887 20,981,239 20,898,338 18,693,394 20,466,156 20,741,131 16,781,836 19,240,010 22,022,869 22,972,945 7,603,917 12,155,484 14,955,918 21,066,522 19,699,663 16,560,125 8,003,115 17,525,150 21,991,987 20,882,140
14,294.81 17,550.81 18,883.12 18,808.50 16,824.05 18,419.54 18,667.02 15,103.65 17,316.00 19,820.58 20,675.65 6843.53 10,939.94 13,460.33 18,959.87 17,729.70 14,904.11 7202.80 15,772.64 19,792.79 18,793.93
Appendix DCloud cover, snow cover, non-snow cover, positive and negative area in April for all years
12
A. Jabbar et al.
Remote Sensing Applications: Society and Environment 18 (2020) 100308
13
A. Jabbar et al.
Remote Sensing Applications: Society and Environment 18 (2020) 100308
. (continued). 14
A. Jabbar et al.
Remote Sensing Applications: Society and Environment 18 (2020) 100308
Appendix E. Cloud cover, snow cover, and non-snow cover, areas in April for all years
15
A. Jabbar et al.
Remote Sensing Applications: Society and Environment 18 (2020) 100308
. (continued).
16
A. Jabbar et al.
Remote Sensing Applications: Society and Environment 18 (2020) 100308
Appendix F. Daily maximum temperature, in-situ data at Skardu and Gilgit observation stations along with temperature extracted from Landsat images for that date Date
Gilgit in-situ daily temperature
Gilgit/Landsat extracted temperature
Skardu in-situ daily temperature
Skardu/Landsat extracted temperature
1993.04.25 1994.04.28 1996.04.17 1999.04.26 2000.04.12 2002.04.10 2003.04.29 2011.04.27 2015.04.22
100* 100* 100* 30.2 29.4 26.2 29.5 34.5 26
25.4 24.11 100* 22.81 25.4 22.32 25.36 30.42 23.51
100* 100* 22.6 22.6 22.1 20.6 22.8 28 19.2
36.04 34.06 25.82 26.67 34.06 30.75 32.65 39.16 21.57
100* ¼ data not available.
Appendix G. Coordinates of the climate observation stations Observatory Name
Latitude
Longitude
Elevation
Gilgit Skardu
35� 180 35� 550
75� 410 74� 200
2317 m 1460 m
Source: Pakistan Meteorological Department.
17
A. Jabbar et al.
Remote Sensing Applications: Society and Environment 18 (2020) 100308
Appendix H. Modified Landsat temperature for the years 1999, 2000, 2002, 2003, 2011, and 2015
18
A. Jabbar et al.
Remote Sensing Applications: Society and Environment 18 (2020) 100308
References
Hewitt, K., Wake, C.P., Young, G.J., David, C., 1989. Hydrological investigations at Biafo glacier, Karakoram range, Himalaya: an important source of water for the Indus River. Ann. Glaciol. 13, 103–108. Immerzeel, W.W., Droogers, P., de Jong, S.M., Bierkens, M.F.P., 2009. Large-scale monitoring of snow cover and runoff simulation in Himalayan River Basins using remote sensing. Remote Sens. Environ. 113, 40–49. https://doi.org/10.1016/j. rse.2008.08.010. Immerzeel, W., van Beek, L.P.H., Bierkens, M.F.P., 2010. Climate change will affect the Asian water towers. Science 328, 1382–1385. https://doi.org/10.1126/ science.1183188. IPCC, 2007. IPCC Fourth Assessment Report (AR4). Ipcc, 1, p. 976 doi: ISSN: 02767783. €b, A., Berthier, E., Nuth, C., et al., 2012. Contrasting patterns of early twenty-firstK€ aa century glacier mass change in the Himalayas. Nature 488, 495–498. https://doi. org/10.1038/nature11324. Krishna, A.P., 2005. Snow and glacier cover assessment in the high mountains of Sikkim Himalaya. Hydrol. Process. 19, 2375–2383. https://doi.org/10.1002/hyp.5890. Kulkarni, A.V., Rathore, B.P., Singh, S.K., 2010. Distribution of seasonal snow cover in central and western Himalaya. Ann. Glaciol. 51, 123–128. https://doi.org/10.3189/ 172756410791386445. NASA, 2008. Landsat 7 Science Data Users Handbook. Science 80, 186. Otsu, N., 1979. A threshold selection method from gray-level histograms. IEEE Trans. Syst. Man Cybern. 9, 62–66. https://doi.org/10.1109/TSMC.1979.4310076. Peel, B.L., Finlayson, B.L., McMahon, T.A., 2007. Updated world map of the KoppenGeiger climate classification.pdf. Hydrol. Earth Syst. Sci. 11, 1633–1644. https:// doi.org/10.5194/hess-11-1633-2007. Quincey, D.J., Braun, M., Glasser, N.F., et al., 2011. Karakoram glacier surge dynamics. Geophys. Res. Lett. 38 https://doi.org/10.1029/2011GL049004. Racoviteanu, A.E., Williams, M.W., Barry, R.G., 2008. Optical remote sensing of glacier characteristics: a review with focus on the Himalaya. Sensors 8, 3355–3383. https:// doi.org/10.3390/s8053355. Ren, J., Jing, Z., Pu, J., Qin, X., 2006. Glacier variations and climate change in the central Himalaya over the past few decades. In: Annals of Glaciology, pp. 218–222. Schmidt, S., Nüsser, M., 2009. Fluctuations of Raikot glacier during the past 70 years: a case study from the Nanga Parbat massif, northern Pakistan. J. Glaciol. 55, 949–959. https://doi.org/10.3189/002214309790794878. Tahir, A.A., 2011. Impact of climate change on the snow covers and glaciers in the Upper Indus River Basin and its consequences on the water reservoirs (Tarbela Dam) – Pakistan. Thesis. UNEP, 2008. Global Glacier Changes : Facts and Figures, WGMS, 88. Xiao, X., Shen, Z., Qin, X., 2001. Assessing the potential of VEGETATION sensor data for mapping snow and ice cover: a Normalized Difference Snow and Ice Index. Int. J. Rem. Sens. 22, 2479–2487. https://doi.org/10.1080/01431160119766. Xiao, L., Ouyang, H., Fan, C., 2019. An improved Otsu method for threshold segmentation based on set mapping and trapezoid region intercept histogram. Optik (Stuttg) 196, 163106. https://doi.org/10.1016/j.ijleo.2019.163106. Young, G.J., Hewitt, K., 1990. Hydrology research in the Upper Indus Basin , Karakoram Himalaya , Pakistan. Hydrol. Mt Areas 139–152. Yu, W., Yang, Y., Savitsky, A., et al., 2013. The Indus Basin of Pakistan: Impacts of Climate Risks on Water and Agriculture. World Bank, Washington, DC. Zemp, M., Haeberli, W., Hoelzle, M., Paul, F., 2006. Alpine glaciers to disappear within decades? Geophys. Res. Lett. 33 https://doi.org/10.1029/2006GL026319. Zhang, G., Yao, T., Xie, H., et al., 2015. An inventory of glacial lakes in the Third Pole region and their changes in response to global warming. Global Planet. Change 131, 148–157. https://doi.org/10.1016/j.gloplacha.2015.05.013. Zimmermann, M., Bichsel, M., Kienholz, H., 1986. Mountain hazards mapping in the Khumbu Himal, Nepal, with prototype map, scale 1:50 000. Mt. Res. Dev. 6, 29–40. https://doi.org/10.2307/3673338.
Hewitt, K., 2005. The Karakoram anomaly? Glacier expansion and the ‘elevation effect,’ Karakoram Himalaya. Mt. Res. Dev. 25, 332–340 doi: 10.1659/0276-4741(2005) 025[0332:TKAGEA]2.0.CO;2. Ageta, Y., Higuchi, K., 1984. Estimation of mass balance components of a summeraccumulation type glacier in the Nepal Himalaya. Geogr. Ann. Phys. Geogr. 66, 249–255. https://doi.org/10.2307/520698. Alford, D.L., Archer, D.R., Bookhagen, B., et al., 2016. Monitoring of Glaciers , Climate , and Runoff in the Hindu Kush-Himalaya Mountains, vol. 1. World Bank Group, Washington, D.C., p. 152 Archer, D., 2003. Contrasting hydrological regimes in the Upper Indus Basin. J. Hydrol. 274, 198–210. https://doi.org/10.1016/S0022-1694(02)00414-6. Archer, D.R., Fowler, H.J., 2004. Spatial and temporal variations in precipitation in the Upper Indus Basin, global teleconnections and hydrological implications. Hydrol. Earth Syst. Sci. 8, 47–61. https://doi.org/10.5194/hess-8-47-2004. Ashraf, A., Naz, R., Roohi, R., 2012. Glacial lake outburst flood hazards in Hindukush, Karakoram and Himalayan Ranges of Pakistan: implications and risk analysis. Geomatics, Nat. Hazards Risk 3, 113–132. https://doi.org/10.1080/ 19475705.2011.615344. Atif, I., Mahboob, M.A., Iqbal, J., 2015. Snow cover area change assessment in 2003 and 2013 using MODIS data of the Upper Indus Basin, Pakistan. J. Himal. Earth Sci. 48, 117–128. Bajracharya, S., Shrestha, B., 2011. The status of glaciers in the Hindu Kush–Himalayan region. Mt. Res. Dev. 33, 114–115. Bajracharya, B., Shrestha, A.B., Rajbhandari, L., 2007. Glacial Lake outburst Floods in the Sagarmatha region. Mt. Res. Dev. 27, 336–344. https://doi.org/10.1659/mrd.0783. Bishop, M.P.J., Shroder, J.F., 2010. Satellite image Atlas of glaciers of the world— glaciers of Asia— GLACIERS OF Pakistan—an overview. US Geol Surv Prof Pap, 1386–F–4. Bolch, T., Buchroithner, M.F., Peters, J., et al., 2008. Identification of glacier motion and potentially dangerous glacial lakes in the Mt. Everest region/Nepal using spaceborne imagery. Nat. Hazards Earth Syst. Sci. 8, 1329–1340. https://doi.org/10.5194/ nhess-8-1329-2008. Burns, P., Nolin, A., 2014. Using atmospherically-corrected Landsat imagery to measure glacier area change in the Cordillera Blanca, Peru from 1987 to 2010. Remote Sens. Environ. 140, 165–178. https://doi.org/10.1016/j.rse.2013.08.026. Copland, L., Sylvestre, T., Bishop, M.P., et al., 2011. Expanded and recently increased glacier surging in the Karakoram. Arctic Antarct. Alpine Res. 43, 503–516. https:// doi.org/10.1657/1938-4246-43.4.503. Dyurgerov, M.B., Meier, M.F., 2005. Glaciers and the Changing Earth System: A 2004 Snapshot. Instaar. Fowler, H.J., Archer, D.R., 2005. Hydro-climatological variability in the Upper Indus Basin and implications for water resources. Iahs 295, 131–138. Fowler, H.J., Archer, D.R., 2006. Conflicting signals of climatic change in the Upper Indus Basin. J. Clim. 19, 4276–4293. https://doi.org/10.1175/JCLI3860.1. Gardelle, J., Berthier, E., Arnaud, Y., 2012. Slight mass gain of Karakoram glaciers in the early twenty-first century. Nat. Geosci. 5, 322–325. https://doi.org/10.1038/ ngeo1450. Hall, D.K., Foster, J.L., Verbyla, D.L., et al., 1998. Assessment of snow-cover mapping accuracy in a variety of vegetation-cover densities in central Alaska. Remote Sens. Environ. 66, 129–137. https://doi.org/10.1016/S0034-4257(98)00051-0. Hewitt, K., 2007. Tributary glacier surges: an exceptional concentration at Panmah glacier, Karakoram Himalaya. J. Glaciol. 53, 181–188. https://doi.org/10.3189/ 172756507782202829. Hewitt, K., 2011. Glacier change, concentration, and elevation effects in the Karakoram Himalaya, Upper Indus Basin. Mt. Res. Dev. 31, 188–200. https://doi.org/10.1659/ MRD-JOURNAL-D-11-00020.1.
19