Inventory of glaciers and glacial lakes of the Central Karakoram National Park (CKNP – Pakistan)

ABSTRACT This study presents a map reporting valuable information on the cryosphere of the Central Karakoram National Park (CKNP, the largest protected area of Pakistan and the highest park in the world). All the information is provided considering the CKNP as a whole, and in detail by dividing it into five basins (i.e. Shigar, Hunza, Shyok, Upper Indus, and Gilgit). The glacier inventory reports 608 ice bodies covering 3680 km2 (∼35% of the CKNP area), with a total glacier volume of ca. 532 km3. In addition, we modeled the meltwater from glacier ice ablation over the period 23 July to 9 August 2011. The total melt amount is ca. 1.5 km3. Finally, we considered glacial lakes (202 water-bodies, covering 4 km2). For these latter glacier features, we also analyzed their potentially dangerous conditions and two lakes were found having such conditions.


Introduction
The Central Karakoram National Park (CKNP) is the largest protected area of Pakistan (about 10,000 km 2 ), established in 1993. It is located in Northern Pakistan in the main glaciated region of the Central Karakoram. This area is situated in High Mountain Asia (HMA), which represents the largest glaciated region outside the Arctic and the Antarctic, the so-called Third Pole, covering an area of more than 100,000 km 2 (Gardner et al., 2013) and hosting about 40,000 km 2 of ice bodies (glaciers, glacierets, and perennial ice surfaces). The CKNP is the highest park in the world, as it is characterized by extremes of altitudes that range from 2000 m a.s.l. to over 8000 m a.s.l., including K2 (8611 m a.s.l.), the second highest peak in the world.
Although other glacier inventories covering the Karakoram region are available (Randolph Glacier Inventory, see Arendt et al., 2014; ICIMOD glacier inventory, see Bajracharya & Shrestha, 2011; GAM-DAM glacier inventory, see Nuimura et al., 2015; World Glacier Monitoring Service-WGMS glacier inventory), this work focuses on the specific area of the CKNP only, providing a high resolution, detailed inventory of glaciers. Moreover, compared to the previous CKNP glacier inventory by Minora et al. (2016), we considered the current park border and partitioned it into its five main catchments (i.e. Shigar, Hunza, Shyok, Upper Indus and Gilgit), which allowed us to obtain a different and more complete picture of the actual CKNP glaciation. In addition, we provide further information, such as supraglacial debris occurrence and thickness, the total freshwater stored in CKNP glaciers (i.e. glacier volume), the amount of freshwater released by glacier ice melt, and a glacial lake inventory with particular attention to potentially dangerous glacial lakes (PDGLs).

Methods
To produce the main map, we considered (i) the glacier boundaries in 2010 developed during the compilation of the CKNP glacier inventory, (ii) supraglacial debris cover and thickness in 2010 and 2011, respectively, (iii) magnitude and rate of glacier ice ablation during an 18-day period in summer 2011, and (iv) the occurrence of glacial lakes in 2013 and their features, which enabled us to classify PDGLs. In order to calibrate and validate our calculations, we coupled remote-sensing investigations with field observations collected during an expedition on the Baltoro Glacier (the widest ice body of the CKNP, 62-km long, with extensive debris cover) in summer 2011.

Glacier inventory
For the compilation of the CKNP glacier inventory, we followed the approach reported in Minora et al. (2016). To detect glaciers, mark their boundaries and calculate their area, remote-sensing investigations were applied. More precisely, Level 1T Landsat Thematic Mapper (TM) scenes from 2010 at 30 m spatial resolution were processed and analyzed (Table 1). We produced false-color images via a 543 band combination (referring to red, green, and blue channels, respectively). The presence of snow or debris could make it difficult to obtain a correct interpretation of the glacier perimeter (Collier et al., 2015;Paul et al., 2009). To detect and exclude steep (identified from breaks in slope), snow covered rock-walls in the accumulation area we used contour lines from the Shuttle Radar Topography Mission 3 DEM (SRTM3, CGIAR-CSI, 2012) (following Nuimura et al., 2015). In addition, we crosschecked the position of the actual glacier border under debris with other Landsat images (Table 1) and high-resolution images from Google Earth © .
Finally, when dealing with the production of glacier inventories through satellite images, inaccuracies may occur due to classification errors. A detailed investigation about the impacts of different sources of error is reported by Minora et al. (2016); based on this analysis, we evaluated the errors affecting our results. In particular, regarding the georeferencing error, the true geolocation is not too critical for our analysis because our Landsat data were processed in the same way by the USGS. We assessed the linear resolution error from the uncertainty due to the sources (Citterio et al., 2007;Vögtle & Schilling, 1999) considering that debris-covered pixels have an error of 45 m, three times that of clean ice (15 m). To limit the error dependent on specific scene conditions (i.e. presence of seasonal snow, cloud cover, shadows, and debris), we selected the scenes displaying minimum snow and cloud cover over the glaciers and we also used images from different sources (i.e. Landsat and Google Earth) and dates (see Table 1). This approach enabled us to cross-check the actual glacier limits and to minimize any possible interpretation error. As regards the error depending on the operator's misinterpretation, the manual digitization was carried out by an expert and a second check on the final mapping was performed.

Supraglacial debris occurrence and thickness
Supraglacial debris can originate from landslides from the steep rock-walls surrounding the glaciers, rock falls, and debris-laden snow avalanches, and it influences the glacier system by modulating the production of freshwater (Collier et al., 2015). In fact, supraglacial debris cover whenever thicker than a threshold called 'critical thickness' (sensu Mattson, Gardner, & Young, 1993) reduces magnitude and rates of buried ice melt with respect to bare ice melt at the same elevation Diolaiuti, D'Agata, Meazza, Zanutta, & Smiraglia, 2009;Mihalcea et al., 2006;Mihalcea, Brock, et al., 2008;Mihalcea, Mayer, et al., 2008). The critical thickness corresponds to the value for which buried ice melt rates equal those of bare ice located at the same altitude (equal to about 0.05 m on Baltoro Glacier according to Mihalcea et al., 2006). Therefore, in addition to the occurrence of supraglacial debris, which highlights the separation of the debris-free and debris-covered zones of each glacier, a map of the thickness of supraglacial debris over the whole glaciated area of the CKNP was generated.
Regarding the presence of supraglacial debris, a supervised maximum likelihood (SML) classification was applied on Landsat false-color composite (FCC, bands 543) images (Table 1). The SML algorithm assumes that values in each spectral band from Landsat TM are normally distributed and calculates the probability that a given image pixel is debris-covered or debris-free based on the values of all spectral bands (Minora et al., 2016). Similarly to the glacier boundaries, the actual presence of debris was cross-checked with different sources (other Landsat images, Table 1, and Google Earth © ).
To estimate supraglacial debris thickness, we followed the method described by Minora et al. (2015) where the input data are (i) debris thickness measured in the field in some selected representative debris-covered glacier areas and (ii) satellite-derived surface temperatures at the same sites. We selected Landsat images from the same period as field measurements undertaken on Baltoro Glacier (July-August 2011, Table 1). In particular, an SML classification was applied to the Landsat FCC image (bands 543). This approach involved training the classification algorithm with a number of sample areas from sites where the correct classification output (i.e. presence or absence of debris on the glacier surface) was known a-priori. Thus, we produced a map describing the supraglacial debris thickness for the year 2011 using the glacier boundary of 2010. This choice is supported by the stability of glacier area in the CKNP, found from the analysis of 2001-2010 data (Minora et al., 2016).

Glacier volume and thickness
To assess the total fresh-water resource stored in CKNP glaciers, we followed the method introduced by Haeberli and Hoelzle (1995). Ice thickness and volume data were estimated from an indirect approach, which considers glacier geometry data recorded in the inventory (2010 data base). Taking into account all approaches to quantify glacier volume reported by Farinotti et al. (2017) and Frey et al. (2014) (i.e. minimization approaches, mass conserving approaches, shearstress-based approaches, velocity-based approaches, area-related thickness estimations, GlabTop2, HFmodel), the one we applied can be considered the most appropriate for the purpose, because information for example about ice flow dynamics is not available for all CKNP glaciers. This approach strongly depends on the parameterization of the average basal shear stress (τ), in particular the upper limit (τ max ), which is used for glaciers with elevation ranges of more than 1600 m. It was found to overestimate τ for Hindu-Kush glaciers with an elevation range (ΔH ) between about 500 and 2000 m, but to underestimate it for glaciers with ΔH > 2500 m (Frey et al., 2014). However, this method was widely applied (e.g. Baumann & Winkler, 2010) and it gave good results in analyzing glaciers from the New Zealand Alps and Norway, suggesting a wide applicability. Moreover, Hoelzle, Haeberli, Dischl, and Peschke (2003) applied such method to estimate changes and evolution of glaciers worldwide, which supports the use of this parameterization for CKNP glaciers. The geometry data required as input in the Haeberli and Hoelzle (1995) analytical approach are the glacier altitudinal range (i.e. the difference between glacier maximum and minimum elevation), the glacier maximum length (measured along the main flow line), and the area.
Average ice depth along the central flow line was estimated from average surface slope (derived from the ratio of altitude range and glacier maximum length) and a mean basal shear stress along the central flow line. The latter depends in a nonlinear way on the altitudinal range as a function of mass turnover (Driedger & Kennard, 1986;Haeberli, 1985;Haeberli & Hoelzle, 1995;Hoelzle et al., 2003).
In 1954, Ardito Desio promoted an expedition to the Baltoro Glacier with the aim of acquiring important geological and glaciological information. In particular, gravimetric surveys were carried out to assess the glacier depth (Marussi, 1964). These data permitted a comparison between modeled and measured ice thickness values on the Baltoro Glacier and thus a discussion of the reliability of our calculations and results.

Meltwater
Once glacier area and supraglacial debris occurrence and thickness were defined, we generated a map highlighting the magnitude of ice ablation and evaluated the derived meltwater amount. Unlike glacier volume, which represents the total water resource stored in CKNP glaciers, meltwater is the actual contribution to runoff and therefore represents the water supply available for civil use, hydropower production and farming. To estimate this daily water amount, we applied two distributed melt models (full details in Minora et al., 2015) to calculate ablation under debris-covered (M DC ) and debris-free (M DF ) conditions (Mihalcea, Mayer, et al., 2008;Pellicciotti et al., 2005). The M DC model depends on the energy available at the debris-ice interface estimated assuming a linear temperature gradient from the top of the debris layer to the ice surface under mean daily conditions (Mihalcea, Mayer, et al., 2008). The daily ice melt at each pixel representing debris-free ice (M DF ) was estimated by applying an enhanced T-index model where air temperature and incoming solar radiation were distributed with elevation gradient approaches (Gambelli, Senese, D'agata, Smiraglia, & Diolaiuti, 2014).
In particular, to model the ice melt amount in the whole CKNP glacier ablation area, we considered the following input data: (i) the 2010 glacier boundaries, (ii) the DEM describing the CKNP area, (iii) the supraglacial debris cover map of 2010, (iv) 2011 meteorological input data (daily mean air temperature and daily mean incoming solar radiation measured by the permanent automatic weather station installed at Askole, Minora et al., 2015), and (v) the 2011 supraglacial debris thicknesses, daily surface debris temperatures (computed from daily incoming solar radiation and debris thickness, Minora et al., 2015) and debris effective thermal resistance (evaluated from debris thickness, Minora et al., 2015).
The maximum errors in the air temperature and solar radiation models are ±1.3°C and ±125 W m −2 , respectively. During the 2011 ablation season, we collected 29 measurements on Baltoro Glacier (both debris-covered and debris-free conditions) at altitudes ranging from 3699 to 5200 m. We divided this dataset into two subgroups: one used to calibrate our melt models and the other to validate them. Considering the second field dataset used for validation, we found a mean error of +0.01 m w.e. (corresponding to 2%) and a RMSE (root mean square error) equal to 0.09 m w.e. (17%).

Glacial lake inventory
In the Main Map, we also show the presence of glacial lakes in the CKNP, derived from a general glacial lakes inventory developed by PARC (Pakistan Agricultural Research Council) and PMD (Pakistan Meteorological Department) in 2015 for the whole HKH (Hindu-Kush Himalaya) area. Eleven Landsat-8 scenes from 2013 were used for identification and characterization of glacial lakes in the HKH region. Table 1 reports the input datasets pertaining to the CKNP area used in the PARC glacial lake inventory. Field surveys were carried out in Hunza, Gilgit and Chitral basins during 2013 in order to assess risk of flood hazards and investigate glacial environments. The DEM was used to decide the criteria for discrimination of features and land-cover types in GIS and to provide better perspective viewing. Different spectral band combinations in False Color Composite (FCC) and individual spectral bands were used to study glacial lakes. Image enhancement using FCC of bands 2, 3, 5 (RGB) of Landsat-8 proved useful in observing glacial lakes and other glacial features in false color.
Glacial lakes originate from glacier activities and/or retreating processes of a glacier and can be classified as: (i) Glacial erosion lakes (water bodies formed in a depression after the glacier has retreated), (ii) Cirque lakes and (iii) Trough Valley lakes (two specific types of glacial erosion), (iv) Supraglacial lakes (in any position of the glacier surface), (v) Lateral moraine dammed and (vi) End moraine dammed lakes (derived from the retreating process of a glacier), and (vii) Blocking lakes (formed through glacier and other factors).
The criteria to identify the PDGLs are based on geomorphological, geotechnical characteristics, and records of past processes and events of the lake. To classify a lake as being potentially dangerous, we considered the lake physical conditions and features and its surroundings as discussed by Mool, Bajracharya, and Joshi (2001), Bajracharya, Mool, and Shrestha (2007), ICIMOD (2011) andPARC et al. (2015). In particular, the following conditions were considered: (i) presence of a group of closely spaced supraglacial lakes on glacier tongues, when they merge and form large lakes that may become potentially dangerous; (ii) the conditions of the damming material in moraine dammed lakes; (iii) the nature of the parent glaciers i.e. presence of a large parent glacier near the lake, debris cover in the glacier snout area and steep gradient at the snout; (iv) presence of crevasses, ponds on the glacier tongue, collapses of glacier masses on the tongue and ice blocks draining to the lake; and (v) physical conditions of the surrounding area like potential rockfalls, mass movements, hanging glaciers, and snow avalanches around the lake which can fall into the lake suddenly. Although a standard index to define a lake as a source of potential danger of out bursting is difficult to develop, because of high topographic and climatic variations in glaciated regions, some general rules can be formulated for specific vulnerable areas for effective risk mitigation. The physical characteristics of the lakes and associated glaciers like their extent, distance to each other, altitudinal/topographic variations, drainage and moraine conditions can be extracted carefully through image interpretation and analysis. Lake extent is one of the most commonly derived attributes, which is important because it is normally related to lake volume, in turn influencing likely GLOF volumes. Similarly, large lakes associated with hanging glaciers at higher altitudes are highly critical as they present outbursting potential in case of avalanche occurrence. The potentially dangerous lakes are generally located at the lower part of the ablation area of the glacier near the end moraine, and the parent glacier should be sufficiently large to create a potentially dangerous lake environment.

Glacier inventory
In the CKNP, there are 608 glaciers (among which some of the largest Karakoram glaciers: Baltoro, Biafo, and Hispar) with a mean size of 6 km 2 . Their total area in 2010 is 3680 ± 61 km 2 , ∼35% of the CKNP area. This area represents ∼24% of the glacier surface of the entire Karakoram Range within Pakistan (total area from Bajracharya & Shrestha, 2011).
The Shigar glaciated area is the widest of the CKNP basins, covering more than half of the whole glaciated surface of the park (i.e. 2308 km 2 , Table 2), and featuring the highest number of glaciers (i.e. 294 bodies, 48% of the total CKNP census, Table 2). In addition, four of the biggest CKNP ice bodies are located in this basin: namely Baltoro Glacier (604 km 2 ), Biafo Glacier (438 km 2 ), Chogo Lungma Glacier (265 km 2 ), and Panmah Glacier (264 km 2 ). Gilgit basin hosts the lowest number of glaciers (36, Table 2, corresponding to 6% of the whole CKNP glacier census) and the glaciated area is only 2% (84 km 2 , Table 2) of the total CKNP glaciation, the smallest compared to the other basins.
In the widest basin (i.e. Shigar), most glaciers (36% of all Shigar glaciers) feature an area lower than 0.5 km 2 , covering only 1% of the whole Shigar glaciated area. On the other hand, glaciers larger than 50 km 2 cover 71% of the whole Shigar glaciated area. Similar conditions are found in the other smaller basins.
The mean glacier terminus elevation in Shigar basin is found to be about 4450 m a.s.l. (in agreement with the other four basins), ranging from 2750 to 5750 m a.s.l. The mean glacier length, in line with the mean value of the whole CKNP, is 3.4 km and the maximum length (not only among glaciers in Shigar basin but also among all CKNP glaciers) is reached by the Biafo Glacier (63.7 km).
We compared our glacier outlines against those from the Randolph Glacier Inventory, version 5.0 (RGI, Arendt et al., 2014), another region-wide inventory. To make the comparison consistent, we selected only the glacier polygons that were mapped in both inventories. The comparison was made for the entire glacier area and for the accumulation area only, because only minor changes over time are expected to occur in the accumulation area. An elevation of 5200 m a.s.l. was used as the equilibrium line altitude (ELA; Minora et al., 2015). We found a total area of 3659 km 2 in our inventory (25% less than the RGI area of 4565.1 km 2 ); the difference in area between the RGI inventory and our results decreases if we consider only the area above the ELA: 1053 km 2 against 1223 km 2 in the RGI, corresponding to −16%. Our inventory thus shows a lower glacier area considering both the whole surface and the accumulation zones. This might be caused by different strategies of mapping the upper glacier limits in the different inventories. In particular, we used a slope criterion to exclude all the headwalls steeper than 60°from the upper glacier limit, while the RGI includes steep headwalls of the accumulation basins in the glacier outlines, thus leading to larger glacier areas, as also reported by Nuimura et al. (2015). In addition, the presence of seasonal snow cover and rock outcrops within glacier areas were considered in the source data of the RGI. These different approaches can partly explain the lower overall glacier area found in our inventory, compared to the RGI 5.

Supraglacial debris occurrence and thickness
The supraglacial debris cover was found to be 920 ± 59 km 2 in 2010, i.e. about 21% of the total ice-covered area. In general, 27% of the CKNP glaciers were found to be debris-covered. Therefore, if CKNP glaciers are divided into debris-free and debris-covered types, we can immediately recognize two patterns. Debris-covered glaciers are generally larger (Baltoro and Hispar glaciers belong to this group, see also the Main Map) and they reach the lowest elevations (even below 3000 m a.s.l.). In fact, supraglacial debris covers between 20% and 27% of glaciers in the size classes larger than 2 km 2 , with a maximum in the size class from 20 to 50 km 2 . Moreover, they are covered by debris almost  entirely up to about 4000 m a.s.l.: the maximum supraglacial debris cover is found at 4300 m a.s.l. Conversely, debris-free glaciers are in general smaller (see also the Main Map), and their termini are found on average at higher elevations (4500 m a.s.l., almost 700 m above the mean termini of debris-covered glaciers). This can be probably explained by reduced average melt rates in the lower tongue sectors of debris-covered glaciers. In fact, in these areas, we found a very thick supraglacial debris cover (>0.5 m, see the Main Map). This abundant presence of supraglacial debris reduces the melting processes affecting the underlying ice (Nicholson & Benn, 2012), allowing glaciers to survive at these lower elevations. In addition, larger glaciers generally extend to lower elevations due to their large ice flux, irrespective of their debris cover. Shyok basin hosts the highest number of debris-covered glaciers (62 ice bodies, Table 2) which make up almost the whole glaciated area (313 km 2 corresponding to 94% of Shyok glacier cover). Considering all CKNP glaciers, supraglacial debris thickness in 2011 (Main Map) is very high at the terminus (up to ∼3 m for example on the Baltoro Glacier) with a mean value of 0.2 m, and an overall rock debris volume of about 0.20 km 3 . The error in debris thickness estimation is ±0.1 m (Minora et al., 2015). The Gilgit basin is characterized by the highest mean thickness, equal to 0.3 m ( Table 2), but the maximum value (3 m) is found to occur over glaciers in the Hunza basin. In fact, all glaciers with mean debris thicknesses higher than 1 m belong to the Hunza basin. Sorting the debris thickness values according to elevation bands (500 m interval), the maximum mean thickness (1.5 m) is found between 2000 and 2500 m.

Glacier volume and thickness
The mean ice thickness was found ranging from more than 200 m (totally two glaciers: 285 and 213 m at Biafo and Baltoro Glaciers, respectively) to 5 m (only one glacier), with an area-weighted average of 145 m, in agreement with findings by Frey et al. (2014). Very small glaciers (i.e. with a surface area smaller than 0.1 km 2 ) are characterized by lower thickness values. Debris-covered glaciers have a mean ice thickness of 41 m (ranging from 9 to 213 m), which is thicker than that of debris-free glaciers (equal to 29 m, ranging from 5 to 285 m). The maximum ice thickness was found on the Biafo Glacier (1362 m) and deep ice thicknesses were also found on the Baltoro (1016 m), Braldu (984 m), and Hispar (906 m) Glaciers. Generally, higher ice thickness values were found in the ablation area compared to the accumulation zones (mean value of 53 m, ranging from 6 to 545 m). We compared data for the Baltoro Glacier with those acquired in 1954 during the well-known expedition led by Desio. On that occasion, gravimetric surveys gave a maximum glacier depth of about 900 m (Marussi, 1964), suggesting that our computations are reliable. In addition, Hewitt, Wake, Young, and David (1989) reported ice thickness of Biafo Glacier derived with a monopulse radar system. The maximum values range from 500 to 700 m at Sim Gang and Baintha to 1400 m at the equilibrium line. These findings are also in agreement with our results.
The total fresh-water resource stored in the CKNP glaciers was estimated as ca. 532.38 km 3 , of which 308.30 km 3 is in debris-covered glaciers and 224.07 km 3 in debris-free ones. Baltoro Glacier has the largest volume (129 km 3 ) and is the largest glacier (with an area of 604 km 2 ), even if Biafo Glacier has the highest mean ice thickness (and is the second largest glacier, with an area of 438 km 2 ). However, more than half of all CKNP glaciers (69%) contains a volume of water lower than 0.05 km 3 (Figure 1), accounting for about 1% of the total volume ( Figure 2). In particular, ice bodies such as glacierets (with an area of about 0.02 km 2 ) have the minimum volume, lower than 0.001 km 3 (Figure 3).

Meltwater
In the Main Map, we report the cumulated ice ablation in the time frame 23 July-9 August 2011 (i.e. 18 days) for which field ablation data were available (see also Minora et al., 2015). The error in the melt model is ±0.09 m w.e. in the considered period (Minora et al., 2015). The total ice melt from the CKNP was equal to 1.54 km 3 w.e. (0.32 km 3 from debris-covered parts of all glaciers and 1.22 km 3 from the debris-free parts), with a daily average of 0.09 km 3 w.e. d −1 . As expected, the highest contribution comes from glaciers located into the Shigar basin (0.92 km 3 w.e., Table 2).

Glacial lake inventory and PDGLs
The CKNP area hosts 202 glacial lakes (Main Map and  Table 2), corresponding to about 7% of the total 3044 glacial lakes reported for the HKH region (PARC et al.,

2015)
. This percentage increases to 23% if compared to the 887 glacial lakes identified in different river-basins of the Karakoram Range by Ashraf, Roohi, and Naz (2011). The park lakes feature a cumulative extent of 4 km 2 ( Table 2, about 3% of the total glacial lake area in the HKH). Detailed information about each glacial lake is available online at http://users.unimi.it/glaciol/ catasto_CKNP.pdf (Smiraglia & Diolaiuti, 2016). Considering the lake type (Table 3), supraglacial lakes are the most common in the CKNP, representing 69% of the total number, and they cover 2 km 2 . Blocked type lakes are the second most abundant accounting for 20% of the total number. The distribution of glacial lakes in CKNP by type shows a different picture with respect to the HKH general conditions. In fact, in the greater HKH region erosion lakes prevail (857 water bodies, 28% of the total number), followed by end moraine dammed lakes (791 water bodies, 26% of the total number) (PARC et al., 2015).
In most cases, larger lakes are a more likely source of GLOF hazards than smaller ones. Thus, we analyzed lakes with a surface area greater than 0.02 km 2 . The CKNP hosts 37 major lakes, corresponding to 18% of the glacial lakes. Most of these lakes (65%) have an area between 0.02-0.05 km 2 . Overall, 17 major lakes belong to the supraglacial type and 16 to the blocked type. In particular, only two PDGLs are found; both of them lie in the Gilgit catchment and are identified as supraglacial (Main Map and Table 2). These PDGLs have caused frequent flooding events in the recent past (Iturrizaga, 2005). In fact, the ephemeral lake developed on the surface of the Hinarchi Glacier has a history of multiple breaching in the Bagrot valley of Gilgit basin (PARC et al., 2015). The other supraglacial lake in the Gilgit basin is also growing rapidly due to the melting of the associated glacier (i.e. Gargo Glacier) in the Bagrot valley, posing a threat of outburst flood for downstream communities (PARC et al., 2015). Moreover, many other supraglacial lakes identified in the park area could develop into PDGLs. This suggests that lake monitoring should continue and early strategies for risk mitigation and disaster management should be developed. The information reported in this study can provide a basis for future monitoring of glacial lakes and GLOFs and for planning and prioritizing disaster mitigation efforts in the park. An automatic early warning system developed and implemented by Haemmig et al. (2014) at Kyagar glacier can provide a potential solution for monitoring active GLOF sites. However, regional GLOF monitoring by satellite observations is still required due to potential undiscovered or inactive GLOF sites (Chan, 2016).

Conclusions
As regards the CKNP as a whole, 608 glaciers are found with a total area of 3680 ± 61 km 2 , ∼35% of the CKNP area. The widest basin (with respect to the number of ice bodies, glacier extent, and ice volume) is the Shigar basin, where the largest glaciers are present (among which Baltoro Glacier), and the smallest one is the Gilgit basin. Finally, the highest number of debris-covered glaciers is located in the Shyok basin (62 glaciers).
For a period of 18 days in summer 2011, we quantified a total water volume of 1.54 km 3 derived from ice melting. Even though we considered a relatively short period, this water volume equals ∼11% of the reservoir capacity of the Tarbela Dam, similar to findings by Minora et al. (2015).
In addition to glacier information, we provided data on glacial lake occurrence, as these ephemeral water bodies can produce actual glacial risk conditions. While only two PDGLs were identified in the park territory, they are located in a highly vulnerable and fragile area and their recent flooding caused destructions of villages and communities.

Software
Esri ArcGIS was used to generate the final map layout of the CKNP reporting glacier boundaries, supraglacial debris occurrence and thickness, meltwater amount and position of glacial lakes, and PDGLs. Google Earth © was used for cross-checking the position of the actual glacier border especially under debris and the actual presence of supraglacial debris.

Acknowledgments
The Central Karakoram National Park (CKNP) Glacier Inventory and CKNP Glacial Lake Inventory were developed under the umbrella of a project managed by Ev-K2-CNR Pakistan and carried out through the cooperation of the Università degli Studi di Milano (Italy) and the Pakistan Meteorological Department. This inventory is an open access data base published in a book in 2016 (Smiraglia and Diolaiuti Editors) and also available online at: users.unimi.it/glaciol (pdf file also reporting tables, diagrams, and maps). The reviewers and the editor are acknowledged for their useful suggestions which improved the article and the maps. The authors are also grateful to Andrea Zerboni for checking and improving the first draft of the main map of this manuscript.

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
This research was developed in the context of the 'Social Economic Environment Development' (SEED) project which focused on the Central Karakoram National Park (CKNP) Gilgit Baltistan Region, Phase II, funded by the Government of Italy and the Government of Pakistan in the framework of the Pakistan-Italian Debt for development Swap Agreement (PIDSA see also http://openaid.esteri.it/en/ projects/initiative/008942/). The main aim of the Project is the promotion of an integrative development of the CKNP region through the support of the implementation and management of the CKNP, the improvement of local wellbeing and livelihood options, through poverty alleviation, community development, livelihood improvement and conservation by means of an integration of intrinsic scientific ecosystem management-oriented research, indigenous practices for natural resource management and ecotourism principles to support the development and implementation of the CKNP. The present study was also carried out by early career researchers supported by DARA (Department of Regional Affairs and Autonomies) of the Presidency of the Council of Ministers of the Italian Government through the Glacio-VAR project (P.I. G. Diolaiuti grant number: COLL_MIN15GDIOL_M).