Lake Evolution, Hydrodynamic Outburst Flood Modeling and Sensitivity Analysis in the Central Himalaya: A Case Study

: Climate change has led to the formation of numerous high-altitude lakes of glacial origin in the Himalaya. Safed Lake is one of the largest glacial lakes, located at an elevation 4882 m a.s.l. in the state of Uttarakhand, central Himalaya, India. A temporal analysis of the lake surface using satellite imagery shows that the lake has grown more than double its size from 0.10 km 2 to 0.23 km 2 over the past 50 years. In this study, we performed a hazard assessment of the lake using 1D and 2D hydrodynamic modeling. We identiﬁed the potential glacial lake outburst ﬂood (GLOF) triggering factors and evaluated the impact of a moraine breach event of the lake on the nearest village located 16.2 km downstream of the lake. A series of dynamic simulations were performed for di ﬀ erent scenario-models based on varied breach depths, breach widths and time of moraine failure. In a worst-case GLOF scenario where breach depth reached up to 60 m, hydrodynamic routing of the breach hydrograph along the given channel revealed inundation depth up to 5 m and ﬂow velocities up to 3.2 m s − 1 at Milam village. Considering the ﬂat geometry of the frontal moraine, hazard assessment of the lake was performed by for di ﬀ erent breach incision depths (30 and 15 m). In addition, the study incorporated a series of hydrodynamic routing to understand the sensitivity of GLOF to di ﬀ erent model input parameters and terrain conditions. The sensitivity of the initial GLOF hydrograph to breach formation time ( T f ) was evaluated by considering di ﬀ erent hypothetical breach scenarios with a varied time of failure. Increases of 11.5% and 22% in the peak ﬂooding were recorded when the moraine failure time was decreased by 15 and 30 min respectively. The two-dimensional sensitivity revealed ﬂow velocity (m s − 1 ) to be more sensitive to change in Manning’s N when compared to the inundation depth (m). Changes of 10.7% and 0.5% in the mean ﬂow velocity (in m s − 1 ) and ﬂow depth (in m) were recorded when dN was 0.01. The ﬂow velocity was more sensitive to the slope and the top-width of the channel when compared to the inundation depths. A regression of ﬂow velocity versus slope gives a correlation coe ﬃ cient of 0.76. GLOF ﬂow hydraulics are sensitive to changes in terrain elevation, where ﬂow depth and velocity vary in a similar manner.


Introduction
Facing climate change-induced decline of glaciers associated with the formation and evolution of glacial lakes in the Hindu Kush-Karakoram-Himalaya (HKKH) region (e.g., [1,2]), glacial lake

Remotely Sensed Data and the Growth Assessment of Safed Lake
A growth assessment of Safed Lake was performed using multitemporal satellite images. We used Corona (KH4) (5.6 m), Landsat TM, Landsat ETM+ and OLI/TIRS (30 m) to map the lake's extent from 1968 to 2018 (years 1968, 1994, 2000 and 2018). The images used in the study were selected during the end of a hydrological year, towards the end of September to avoid seasonal or inter-annual variation. Terrain data for GLOF modeling was obtained from Advanced Spaceborne Thermal Emission and Reflectance Radiometer (ASTER) global digital elevation model (GDEM; version 2). ASTER Digital Elevation Model (DEM) is a freely available model (https://earthexplorer.usgs.gov/) that provides elevation information between 83° N and 83° S with a spatial resolution of 30 m and vertical accuracy of ≈10-20 m for hilly terrains [20,21]. ASTER GDEM has certain limitations related to the spatial resolution; e.g., underestimation of the valley floor elevations [22]. It is, however, the best freely available DEM for the remote Indian Himalaya. The land use land cover (LULC)

Remotely Sensed Data and the Growth Assessment of Safed Lake
A growth assessment of Safed Lake was performed using multitemporal satellite images. We used Corona (KH4) (5.6 m), Landsat TM, Landsat ETM+ and OLI/TIRS (30 m) to map the lake's extent from 1968 to 2018 (years 1968, 1994, 2000 and 2018). The images used in the study were selected during the end of a hydrological year, towards the end of September to avoid seasonal or inter-annual variation. Terrain data for GLOF modeling was obtained from Advanced Spaceborne Thermal Emission and Reflectance Radiometer (ASTER) global digital elevation model (GDEM; version 2). ASTER Digital Elevation Model (DEM) is a freely available model (https://earthexplorer.usgs.gov/) that provides elevation information between 83° N and 83° S with a spatial resolution of 30 m and vertical accuracy of ≈10-20 m for hilly terrains [20,21]. ASTER GDEM has certain limitations related to the spatial resolution; e.g., underestimation of the valley floor elevations [22]. It is, however, the best freely available DEM for the remote Indian Himalaya. The land use land cover (LULC) Figure 2. (a) The source zone of potential mass movements; (b) calving zone at the snout of the glacier. All images: Google Earth [19].

Remotely Sensed Data and the Growth Assessment of Safed Lake
A growth assessment of Safed Lake was performed using multitemporal satellite images. We used Corona (KH4) (5.6 m), Landsat TM, Landsat ETM+ and OLI/TIRS (30 m) to map the lake's extent from 1968 to 2018 (years 1968, 1994, 2000 and 2018). The images used in the study were selected during the end of a hydrological year, towards the end of September to avoid seasonal or inter-annual variation. Terrain data for GLOF modeling was obtained from Advanced Spaceborne Thermal Emission and Reflectance Radiometer (ASTER) global digital elevation model (GDEM; version 2). ASTER Digital Elevation Model (DEM) is a freely available model (https://earthexplorer.usgs.gov/) that provides elevation information between 83 • N and 83 • S with a spatial resolution of 30 m and vertical accuracy of ≈10-20 m for hilly terrains [20,21]. ASTER GDEM has certain limitations related to the spatial resolution; e.g., underestimation of the valley floor elevations [22]. It is, however, the best freely available DEM for the remote Indian Himalaya. The land use land cover (LULC) classification, to determine Manning's roughness coefficient (Manning's N) of the given terrain, was obtained using GlobCover (version 2) and cross-verified using high-resolution geo-referenced CNES/Airbus imagery tiles of Google Earth [19]. GlobCover is the global product of land cover maps created from the 300 m MERIS sensor onboard of the ENVISAT satellite [23].
The normalized difference water index (NDWI) [24] was employed to identify and map the lake's surface area. In order to check the accuracy of mapping, the lake's outlines were overlapped on high-resolution Google Earth imageries using GIS-based tools. The orthorectification of the historical image (Corona) has been carried out using ground control points (GCP). The Corona subset image of the lake with no cloud cover and minimum snow was used to calculate the total area for the year 1968. Further, we employed the rubber sheeting method to orthorectify the given image using a total of 50-70 GCP's from high-resolution geo-referenced CNES/Airbus imagery tiles [19], which reduced the uncertainty significantly to less than a pixel (in the total area calculations). The lake boundaries were delineated manually using the rectified Corona images. Corona images had previously been used successfully for mapping the glaciated terrain in Himalaya [11,25,26]. Since the present study aimed to evaluate the GLOF hazard of Safed Lake at a present scenario, we considered the latest extent (2018) of the lake for volume estimation and flood routing.

Hazard Assessment of Safed Lake
Due to the absence of bathymetric data for Safed Lake, we employed an area-based scaling relationship to calculate the total volume of the lake. Several empirical relationships are available to roughly estimate the total volume of glacial lakes based on the total surface area [24,[27][28][29][30][31]. The relation given by Huggel et al. [24] is among the most widely used and has been employed in several previous GLOF hazard studies, where no direct bathymetric measurements are available [14,[32][33][34][35][36][37]. The total water volume was calculated considering the total surface area in 2018. The equation [24] is given as: (1) where V is the total volume of the lake and A is the area of the lake. The total volume of the lake was calculated to be 4.34 × 10 6 m 3 . A difference of 7% (higher) was calculated when compared to the volume estimates obtained from the Huggel et al. [24] equation in the total lake volume using the equation of Cook and Quincey [27] (V = 0.1217 A 1.4129 ). Despite this is a rough estimate, we considered this volume estimation for the modeling of the worst-case GLOF scenario for Safed Lake. Apart from that, we also considered releases of 50% and 25% of the total lake volume (see below). One-dimensional modeling (HEC-RAS) was employed to calculate the breach hydrograph. The average breach width and failure time were calculated using the empirical relations given by Froehlich [38]: where B w is the breach width (in m), V w is the volume of the lake (in m 3 ), h b is the height of the breach (in m) and T f (in h) is the time taken for the breach to form (where distances B w and h b are fully reached; see Figure 3). Here, the volume calculated using Equation (1) is considered to derive the breach parameters. The h b is derived from DEM and is considered to be 60 m, assuming that the moraine breaches up its base, in a worst-case scenario. Based on the above equations, the breach width Equation (1) and breach formation time Equation (2) in a worst-case scenario are calculated to be 73.13 m and 0.21 h respectively (Scenario 1; Table 1). Water 2020, 12, x FOR PEER REVIEW 5 of 19 However, considering the flat geometry of the dam (a relatively large width of >400 m, and a distal face of the dam with a slope <25%, which is comparatively gentle compared to the maximum slope which moraine dams can reach), it is unlikely that the moraine incision would take place up to a depth of 60 m. Therefore, apart from the worst-case scenario, we considered two different scenarios based on varied breach depth (hb) (30 and 15 m) and total volume of water released in a potential GLOF event (Scenarios 2 and 3; Table 1). Assuming that the volume released in a potential GLOF event of Safed Lake varies linearly with the breach depth, we constructed the model in a way that half of the total volume and quarter of the total volume were released for breach depth of 30 and 15 m respectively (Scenarios 2 and 3). These GLOF scenarios allowed us to capture the wide variability of diverse GLOF triggers and their magnitudes, which we hoped would lead us to different process chains. The breach parameters for all the given models were calculated based on Froehlich [38]m and are given in Table 1. The time of failure (Tf) is an inverse function of the breach depth (hb), and therefore, it is longer for shallower depth [38]. The HEC-RAS 1D model [39] was employed to model a series of potential GLOF events of Safed Lake. Despite that this tool has not been widely employed in Uttarakhand previously, it has been successfully employed in GLOF modeling (both reconstructing past GLOFs and modeling potential future GLOFs) in different geographical settings and contexts (e.g., [40][41][42]). As such, we were convinced this tool would provide physically plausible modeling results in our case too.  However, considering the flat geometry of the dam (a relatively large width of >400 m, and a distal face of the dam with a slope <25%, which is comparatively gentle compared to the maximum slope which moraine dams can reach), it is unlikely that the moraine incision would take place up to a depth of 60 m. Therefore, apart from the worst-case scenario, we considered two different scenarios based on varied breach depth (h b ) (30 and 15 m) and total volume of water released in a potential GLOF event (Scenarios 2 and 3; Table 1). Assuming that the volume released in a potential GLOF event of Safed Lake varies linearly with the breach depth, we constructed the model in a way that half of the total volume and quarter of the total volume were released for breach depth of 30 and 15 m respectively (Scenarios 2 and 3). These GLOF scenarios allowed us to capture the wide variability of diverse GLOF triggers and their magnitudes, which we hoped would lead us to different process chains. The breach parameters for all the given models were calculated based on Froehlich [38] m and are given in Table 1. The time of failure (T f ) is an inverse function of the breach depth (h b ), and therefore, it is longer for shallower depth [38].
The HEC-RAS 1D model [39] was employed to model a series of potential GLOF events of Safed Lake. Despite that this tool has not been widely employed in Uttarakhand previously, it has been successfully employed in GLOF modeling (both reconstructing past GLOFs and modeling potential future GLOFs) in different geographical settings and contexts (e.g., [40][41][42]). As such, we were convinced this tool would provide physically plausible modeling results in our case too.
The initial breach hydrograph was constructed for each scenario based on the parameters given in Table 1. The worst-case GLOF scenario of the lake is revealed in Scenario 1 (B w = 73.13 m and T f = 0.21 h). Two-dimensional hydraulic routing of the breach hydrographs calculated for Scenario 1, Scenario 2 and Scenario 3 was performed for a distance of 16.2 km from the lake up to Milam village to evaluate its impact on the village. The HEC-RAS 2D (version 5.0.5; Brunner [43]) was used to rout the initial breach hydrographs along the main flow channel. The primary input of the model includes terrain data and boundary conditions for simulating an unsteady hydraulic process. ASTER DEM raster was used as a continuous terrain input. A flow area was initially defined over the given terrain, containing the area of interest; i.e., the lake to Milam village. A 2D mesh was constructed within the given flow area, with an individual cell dimension of 30 × 30 m. The land use/land cover (LULC) for the flow area was extracted from GlobCover [23]. Further, a value of Manning roughness was defined to each cell based on the LULC. An average Manning's N within the flow area was calculated to be 0.05, which was considered for routing. Spatially distributed hydraulic properties (flow velocity and inundation depth) were mapped along the channel for a distance of 16.2 km from the lake to the village.

Model Sensitivity to Input Parameters
A series of hydraulic simulations were performed to evaluate the sensitivity of the hydrodynamic model to the input parameters (Manning's N and breach parameters). Here, we assess the sensitivity of the model to the dam-breach parameters (breach width and breach formation time) in order to evaluate its effect on the initial breach hydrograph. A total of six hypothetical one-dimensional moraine breach simulations were performed with varied breach width (B w ) and time of failure (T f ). Three scenarios with T f = 0.50 h, T f = 0.75 h and T f = 1.0 h with a constant B w =100 m were considered to evaluate the sensitivity of the model to T f . Similarly, keeping T f constant, the model sensitivity to B w was evaluated by considering three scenarios with varied B w (50, 75 and 100 m). The peak discharge and time of peak in each scenario were compared to assess the sensitivity of the model to the input parameters. Moreover, to evaluate the two-dimensional sensitivity of the model to the above hypothetical breach hydrographs, the modeled hydrographs were routed along the flow channel from the lake to Milam village. A Manning's N of 0.05 was considered for dynamic routing (Section 3.3.2). The model sensitivity was evaluated by assessing the change in the flow hydraulics (inundation depth and flow velocity) for the different upstream boundary conditions (input hydrographs).
In order to evaluate the sensitivity of the model to Manning's N, a series of two-dimensional routings of the dam-breach hydrograph (obtained in a worst-case scenario, Scenario 1) was performed. This was achieved by the two-dimensional routing of the potential moraine-breach hydrograph for varied Manning's N (0.03, 0.04, 0.05, 0.06 and 0.07). The model sensitivity was evaluated by assessing the change in the flow depth and flow velocity (maximum and mean) for different values of Manning's N.

Sensitivity to Channel Characteristics and DEM Used
The sensitivity of the flow hydraulics of a modeled potential GLOF (Section 4.2) to the slope and top-width of the channel (GLOF inundation width) was evaluated. The slope was calculated for each 100 m elevation interval along the entire length of the given flow channel (from the lake to Milam village). The mean velocity in each 100 m elevation band was calculated from the obtained spatially-distributed 2D modeled results and plotted against its respective slope values to evaluate its relationship. Similarly, the modeled 2D flow depth was extracted for every elevation band, and the mean depth calculated for each band was plotted against the slope of the channel for the sensitivity evaluation.
The sensitivity of the modeled results to the top-width of the flow channel was evaluated by analyzing the variation in the hydraulic properties for different channel top-width. For this, the width of the modeled GLOF inundation in the given channel was measured across 20 cross-sections. Further, the top-width (in m) was plotted against the mean flow velocity (ms −1 ) and inundation depth (m) calculated along the different cross-sections. The variations of flow depth and velocity to the different top-widths of the channel were then evaluated.
Since DEM is a critical element for GLOF route modeling, we evaluated the sensitivity of the GLOF hydraulics to the terrain data (DEM) by comparing the flow hydraulics, performed over two different DEMs (ASTER and SRTM). We considered the worst-case GLOF scenario (Scenario 1; B w = 73.13, T f = 0.21) for flood routing over both ASTER DEM (Section 3.2) and SRTM DEM. We divided the flow channel (3200-4960 m a.s.l.) into elevation bands of 250 m each (7 bands). The terrain elevation values of each DEM (ASTER and SRTM) were extracted within the inundation boundary of each elevation band. The difference in the elevation was calculated using raster-based operations for the elevation bands (SRTM elevation -ASTER elevation). Similarly, we calculated the difference in the routed flow depth (m) and flow velocity (m/s) over each DEM (SRTM flow hydraulics -ASTER flow hydraulics). Further, we plotted the average flow depth difference (in m) and average flow velocity (in m/s) difference against the elevation difference (m) for the given elevation bands.

Hazard Assessment of Safed Lake
The occurrence of GLOF is controlled by two factors-dam stability and the possibility of a triggering event. As for dam stability, we consider moraine/periglacial debris dam of Safed Lake to be potentially unstable, and moreover possibly filled with buried ice cores (see Section 2); therefore, we consider the dam susceptible to breaching and failure, despite relatively flat geometry. As for the triggers, Safed Lake is located in the close proximity to steep slopes and we found mass movements from these slopes to be the most likely triggers of potential GLOF from this lake.
Froehlich's [23] equations were employed to model a series of potential GLOFs of Safed Lake. In the worst-case scenario (Scenario 1, see Table 1), the breach event was modeled considering a breach width (B w ) of 73.13 m and time of failure (T f ) of 0.21 h. The GLOF hydrograph produced a peak discharge of 8181 m 3 s −1 that was achieved within 6 min after the initial breach event. Figure 5 shows the breach hydrograph in a worst-case GLOF scenario of Safed Lake. Since we modeled the breach formation in a "Sinewave" progression and not in a "linear" progression, we expected a peak before the breach was fully developed. In addition, two different potential GLOF scenarios were modeled, for which breach parameters were calculated based on varied breach depth (h b ) and the volume of water (m 3 ) released (see Table 1). In Scenario 2 (h b = 30, B w = 51.33 m, T f = 0.27 h) the GLOF hydrograph releasing the half (1/2) the volume of the lake produced a peak of 5110 m 3 s −1 (peak attained in 12 min after the initial breach event). Similarly, a GLOF hydrograph was constructed considering a new set of breach parameters (h b = 15, B w = 35.97 m, T f = 0.34 h), releasing 1/4 of the lake volume (Scenario 3). A peak discharge of 1495 m 3 s −1 was attained 22 min after the initial moraine breach event.    Table 1).
In order to evaluate the potential impact of a Safed Lake GLOF on Milam village, twodimensional hydraulic routing of the breach hydrographs obtained in Scenario 1 (hb = 60 m, releasing full lake volume), Scenario 2 (hb = 30 m, releasing 1/2 the lake volume) and Scenario 3 (hb = 15 m,    Table 1).
In order to evaluate the potential impact of a Safed Lake GLOF on Milam village, twodimensional hydraulic routing of the breach hydrographs obtained in Scenario 1 (hb = 60 m, releasing full lake volume), Scenario 2 (hb = 30 m, releasing 1/2 the lake volume) and Scenario 3 (hb = 15 m, Figure 5. The potential moraine-breach hydrograph of Safed Lake. A peak discharge of 8181 m 3 s −1 was reached within 6 min after the initiation of the dam breach event (worst-case scenario: Scenario 1); breach hydrographs for varied breach depth (h b ) (30 m and 15 m) and total volume of water released in a potential glacial lake outburst flood (GLOF) event (Scenarios 2 and 3) (see Table 1).
In order to evaluate the potential impact of a Safed Lake GLOF on Milam village, two-dimensional hydraulic routing of the breach hydrographs obtained in Scenario 1 (h b = 60 m, releasing full lake   In Scenario 2, the maximum inundation depth and flow velocity at Milam village reach up to 2.3 m and 1.3 m s −1 respectively. The GLOF arrives in the village at 2 h after the breach event and the peak is reached within 15 min of the flood wave arrival. The total GLOF inundation reduces by 31% when compared to Scenario 1. In Scenario 3, the inundation depth reaches up to a maximum of 0.9 m at Milam village with flow velocities reaching up to 0.3 m s −1 . The peak in the flow hydraulics is attained in 6 h 35 min after the initial breach event of Safed Lake. The longer time taken here for the flood to reach Milam village can be attributed to the fact that the GLOF discharge is very low and the bed-resistance dominates the convective acceleration of the flow. The total inundation is reduced by 63% when compared to the worst-case scenario (Scenario 1). Table 2 summarizes the hydraulics of the GLOF hydrographs and the routed flow properties at Milam village. In Scenario 2, the maximum inundation depth and flow velocity at Milam village reach up to 2.3 m and 1.3 m s −1 respectively. The GLOF arrives in the village at 2 h after the breach event and the peak is reached within 15 min of the flood wave arrival. The total GLOF inundation reduces by 31% when compared to Scenario 1. In Scenario 3, the inundation depth reaches up to a maximum of 0.9 m at Milam village with flow velocities reaching up to 0.3 m s −1 . The peak in the flow hydraulics is attained in 6 h 35 min after the initial breach event of Safed Lake. The longer time taken here for the flood to reach Milam village can be attributed to the fact that the GLOF discharge is very low and the bed-resistance dominates the convective acceleration of the flow. The total inundation is reduced by 63% when compared to the worst-case scenario (Scenario 1). Table 2

Model Sensitivity to Input Parameters
The sensitivity of the hydrodynamic model to breach width (Bw), time of failure (Tf) and Manning's N were evaluated by a series of hypothetical-scenario models. In order to assess the influence of one-dimensional breach parameters on the breach hydrograph, six hypothetical scenarios with varied Bw and Tf values were considered. At first, three scenarios were considered with varied Tf (30, 45 and 60 min) and a constant Bw (100 m). The peak discharge changes by ±11.5% when dTf is 15 min. The time of peak discharge changes by ±4 to ±5 min when Tf is increased or decreased by 15 min (Figure 8a). Similarly, the remaining three scenarios were modeled by keeping Tf (30 min) constant and varying Bw (50, 75 and 100 m). The peak discharge decreased by 2.3% when Bw was reduced by 25 m. It further decreases by 14.4% when Bw was reduced by 50 m. The time of peak was delayed by 1 min and 4 min when Bw was reduced by 25 and 50 m respectively (Figure 8b). From the above GLOF scenario assessment, it is evident that GLOF hydrographs are more sensitive to the time of moraine failure compared to the width of the breach.
To evaluate the two-dimensional sensitivity of the model to the different boundary conditions (input breach hydrographs), the modeled hypothetical GLOF hydrographs with varied Tf were (Figure 8a) routed along the flow channel from the lake to Milam village. The inundation depth and flow velocity of the entire channel were evaluated. The maximum inundation depth decreases by 3.7% and 5.6% when Tf is increased by 15 = and 30 min respectively. Similarly, the mean flow depth decreases by 2.9% and 5.4% when Tf is increased by 15 and 30 min respectively ( Table 3). The evaluation of the maximum flow velocity along the given channel reveals a decrease of 2% and 3% when Tf is increased by 15 and 30 min respectively. Also, the mean velocity shows a decrease of 4%

Model Sensitivity to Input Parameters
The sensitivity of the hydrodynamic model to breach width (B w ), time of failure (T f ) and Manning's N were evaluated by a series of hypothetical-scenario models. In order to assess the influence of one-dimensional breach parameters on the breach hydrograph, six hypothetical scenarios with varied B w and T f values were considered. At first, three scenarios were considered with varied T f (30, 45 and 60 min) and a constant B w (100 m). The peak discharge changes by ±11.5% when dT f is 15 min. The time of peak discharge changes by ±4 to ±5 min when T f is increased or decreased by 15 min (Figure 8a). Similarly, the remaining three scenarios were modeled by keeping T f (30 min) constant and varying B w (50, 75 and 100 m). The peak discharge decreased by 2.3% when B w was reduced by 25 m. It further decreases by 14.4% when B w was reduced by 50 m. The time of peak was delayed by 1 min and 4 min when B w was reduced by 25 and 50 m respectively (Figure 8b). From the above GLOF scenario assessment, it is evident that GLOF hydrographs are more sensitive to the time of moraine failure compared to the width of the breach.
To evaluate the two-dimensional sensitivity of the model to the different boundary conditions (input breach hydrographs), the modeled hypothetical GLOF hydrographs with varied T f were (Figure 8a) routed along the flow channel from the lake to Milam village. The inundation depth and flow velocity of the entire channel were evaluated. The maximum inundation depth decreases by 3.7% and 5.6% when T f is increased by 15 and 30 min respectively. Similarly, the mean flow depth decreases by 2.9% and 5.4% when T f is increased by 15 and 30 min respectively ( Table 3). The evaluation of the maximum flow velocity along the given channel reveals a decrease of 2% and 3% when T f is increased by 15 and 30 min respectively. Also, the mean velocity shows a decrease of 4% and 7.4% when T f is increased by 15 and 30 min respectively ( Table 3). The hydraulic properties of a GLOF at any point downstream is influenced by the initial breach hydrograph. The results show that the maximum inundation depth and the mean velocity are slightly more sensitive compared to the mean depth and maximum velocity respectively.
Water 2020, 12, x FOR PEER REVIEW 11 of 19 and 7.4% when Tf is increased by 15 and 30 min respectively ( Table 3). The hydraulic properties of a GLOF at any point downstream is influenced by the initial breach hydrograph. The results show that the maximum inundation depth and the mean velocity are slightly more sensitive compared to the mean depth and maximum velocity respectively.  At Milam village, maximum inundation depth decreases by 4.2% and 8.5% when Tf is increased by 15 and 30 min respectively. Likewise, the maximum flow velocity decreases by 4.3% and 7.5% when Tf is increased by 15 and 30 min respectively ( Table 4). The temporal evaluation of the hydraulic parameters reveals a delay in the peak by 8 and 14 min when breach formation is increased by 15 and 30 min respectively. The timing of GLOF (peak) is sensitive to the initial flow hydrograph, which in turn is a function of the breach width (Bw) and the breach formation time (Tf). Bw and Tf are driven by various factors, such as the geometry and internal structure of the moraine dam and the type and magnitude of triggering events [44,45]. The sensitivity of the model to Manning's N was tested by performing a series of twodimensional routing of the dam-breach hydrograph (see Figure 5). Keeping the dam-breach parameters constant (as in Section 4.2), two-dimensional routing was performed for varied Manning's N (0.03, 0.04, 0.05, 0.06 and 0.07). The results were evaluated by plotting the maximum  At Milam village, maximum inundation depth decreases by 4.2% and 8.5% when T f is increased by 15 and 30 min respectively. Likewise, the maximum flow velocity decreases by 4.3% and 7.5% when T f is increased by 15 and 30 min respectively ( Table 4). The temporal evaluation of the hydraulic parameters reveals a delay in the peak by 8 and 14 min when breach formation is increased by 15 and 30 min respectively. The timing of GLOF (peak) is sensitive to the initial flow hydrograph, which in turn is a function of the breach width (B w ) and the breach formation time (T f ). B w and T f are driven by various factors, such as the geometry and internal structure of the moraine dam and the type and magnitude of triggering events [44,45]. The sensitivity of the model to Manning's N was tested by performing a series of two-dimensional routing of the dam-breach hydrograph (see Figure 5). Keeping the dam-breach parameters constant (as in Section 4.2), two-dimensional routing was performed for varied Manning's N (0.03, 0.04, 0.05, 0.06 and 0.07). The results were evaluated by plotting the maximum and the mean flow depth and velocity against the Manning's roughness coefficient (Figure 9). The maximum flow velocity decreases with the increase in the channel roughness showing a correlation coefficient of 0.91 (strong correlation) (Figure 9a). Similarly, the mean velocity versus Manning's N revealed a similar trend with a correlation coefficient of 0.93 (strong correlation) (Figure 9b). The primary motive was to evaluate which parameter (flow velocity or flow depth) is more sensitive to the channel roughness.
Water 2020, 12, x FOR PEER REVIEW 12 of 19 primary motive was to evaluate which parameter (flow velocity or flow depth) is more sensitive to the channel roughness. The results reveal that flow velocity is more sensitive to the channel roughness compared to the inundation depth, as both maximum and mean depth do not show significant changes with change in the channel roughness. The Pearson correlation coefficient for maximum depth versus Manning's N was calculated to be 0.83 (strong correlation), and that of mean depth versus Manning's N was 0.98 (strong correlation). Figures 10 and 11 show the spatial distributions of the maximum depth and flow velocity from Safed Lake to Milam village for different Manning's roughness coefficients (N) along the channel.  The results reveal that flow velocity is more sensitive to the channel roughness compared to the inundation depth, as both maximum and mean depth do not show significant changes with change in the channel roughness. The Pearson correlation coefficient for maximum depth versus Manning's N was calculated to be 0.83 (strong correlation), and that of mean depth versus Manning's N was 0.98 (strong correlation). Figures 10 and 11 show the spatial distributions of the maximum depth and flow velocity from Safed Lake to Milam village for different Manning's roughness coefficients (N) along the channel.

Sensitivity to Channel Characteristics
This section of the study evaluates the sensitivity of the flow hydraulics (flow velocity and inundation depth) to channel characteristics (slope, top-width and terrain elevation). The slope calculated for every 100 m elevation contour was plotted against the mean velocity calculated for each elevation band. A similar analysis was performed by considering the mean depth and the slope of the channel. It is evident that the flow velocity linearly varies with change in the slope of the given channel with a correlation coefficient of 0.76 (Figure 12a). However, flow depth is independent (R 2 = 0.14) of the channel slope, as it shows a random variation (weak correlation) (Figure 12b). Figure 9. Plot of (a) Maximum inundation depth/flow velocity versus Manning's N; (b) mean inundation depth/flow velocity versus Manning's N.  To evaluate the sensitivity of the hydraulic properties of a GLOF to the top-width of the flow channel, the flow velocity and depth are plotted against the top-width of the channel. The top-width was calculated along 20 cross-sections covering the entire flow channel, for which the mean flow depth and velocity were calculated. It is evident that the flow velocity and depth show similar trends when plotted against the respective top-width with correlation coefficients of 0.83 and 0.93 respectively ( Figure 13). Both the flow velocity and depth tend to decrease as the top-width of the channel increases.
The sensitivity of GLOF hydraulics to terrain elevation (DEM) was evaluated by plotting the (SRTM DEM -ASTER DEM) elevation difference (in m) versus flow depth difference (in m) for every 250 m elevation contour band (see Section 3.3b) (Figure 14a). A correlation coefficient of 0.46 was calculated. Similarly, the plot of elevation difference (in m) versus flow velocity difference (in m/s) shows a correlation coefficient of 0.54 (Figure 14b). It is evident that GLOF hydraulics is sensitive to terrain elevation changes, with flow velocity being slightly more sensitive when compared to the flow depth. Both these coefficients show that there is a positive correlation between the elevation difference and flow depth difference and flow velocity difference, respectively. Practically, this reveals that the larger the difference between the DEMs used (3833-3583 m.a.s.l. elevation band in our case), the larger difference in modeled results (flow depth and flow velocity) observed among different DEMs.

Sensitivity to Channel Characteristics
This section of the study evaluates the sensitivity of the flow hydraulics (flow velocity and inundation depth) to channel characteristics (slope, top-width and terrain elevation). The slope calculated for every 100 m elevation contour was plotted against the mean velocity calculated for each elevation band. A similar analysis was performed by considering the mean depth and the slope of the channel. It is evident that the flow velocity linearly varies with change in the slope of the given channel with a correlation coefficient of 0.76 (Figure 12a). However, flow depth is independent (R 2 = 0.14) of the channel slope, as it shows a random variation (weak correlation) (Figure 12b).

Sensitivity to Channel Characteristics
This section of the study evaluates the sensitivity of the flow hydraulics (flow velocity and inundation depth) to channel characteristics (slope, top-width and terrain elevation). The slope calculated for every 100 m elevation contour was plotted against the mean velocity calculated for each elevation band. A similar analysis was performed by considering the mean depth and the slope of the channel. It is evident that the flow velocity linearly varies with change in the slope of the given channel with a correlation coefficient of 0.76 (Figure 12a). However, flow depth is independent (R 2 = 0.14) of the channel slope, as it shows a random variation (weak correlation) (Figure 12b).  To evaluate the sensitivity of the hydraulic properties of a GLOF to the top-width of the flow channel, the flow velocity and depth are plotted against the top-width of the channel. The top-width was calculated along 20 cross-sections covering the entire flow channel, for which the mean flow depth and velocity were calculated. It is evident that the flow velocity and depth show similar trends when plotted against the respective top-width with correlation coefficients of 0.83 and 0.93 respectively ( Figure 13). Both the flow velocity and depth tend to decrease as the top-width of the channel increases.  (Figure 14b). It is evident that GLOF hydraulics is sensitive to terrain elevation changes, with flow velocity being slightly more sensitive when compared to the flow depth. Both these coefficients show that there is a positive correlation between the elevation difference and flow depth difference and flow velocity difference, respectively. Practically, this reveals that the larger the difference between the DEMs used (3833-3583 m.a.s.l. elevation band in our case), the larger difference in modeled results (flow depth and flow velocity) observed among different DEMs.

Discussion
The study evaluated a series of potential GLOF scenarios of Safed Lake and their impact on the nearest settlement Milam village, using 1D and 2D dynamic modeling. The potential breach event in a worst-case GLOF scenario (h b = 60 m) of Safed Lake produced a GLOF hydrograph with a peak discharge of 8181 m 3 s −1 . Hydrodynamic routing of the potential GLOF hydrograph along the given channel revealed inundation depth up to a maximum of 5 m and flow velocities up to 3.2 m s −1 at Milam village. The GLOF potentially inundates a part of the village, thereby affecting a total of 25 to 30 houses located at the site. The major part of the settlements in the village remains unaffected, as the flow channel takes a sharp turn towards the south and the flow is obstructed by the elevated terrain (see Figure 7b). Further, hazard assessments of the lake were carried out based on different breach depths (30 and 15 m). The maximum inundation depth at Milam village was reduced by 54% and 94% when h b = 30 m and h b = 15 m, respectively. Similarly, the flow velocity reduced by 59% and 90.6% when h b = 30 m and h b = 15 m respectively. In terms of hazard posed to Milam village, we further highlight the hazard of lateral erosion of the river banks, which may affect buildings not modeled to be flooded. This specific hazard also needs to be taken into consideration when designing flood protection and mitigation measured.
Hydrodynamic modeling of GLOFs involves a series of process-modeling, which includes dam breach formation and flood routing. A challenging task is to evaluate the accuracy of the modeled outputs, as physically-based hydraulic modeling is sensitive to the various inputs the model demands [46,47]. The process of moraine dam breach is complex and is dependent on many factors like height, width, and slope of the moraine [48][49][50]. The hydrograph obtained during a breaching process determines the flow hydraulics in the downstream region [47,50]. Moreover, the channel properties also affect GLOFs and their nature, by influencing the flow over a given terrain. The Manning's roughness coefficient has a significant effect on GLOFs, as it directly influences the peak discharge and flow velocity over a given channel [51][52][53][54]. Anacona et al. [50] performed a sensitivity analysis to evaluate the effect of Manning's N and breach hydrograph on GLOF. These parameters were reported to be the main controlling factors for a given GLOF event. Bajracharya et al. [51] stated the importance of Manning's N for accurate GLOF simulations. A 9% decrease in the peak discharge was observed as the Manning's N value was increased by 10% [14], whereas only 8% variation in the discharge was reported when Manning's N was varied by 30% [15]. This was attributed to the dominance of convective acceleration over bed-resistance, and the breach forming process too. Similar findings by Yochum et al. [54] suggest that hydraulic resistance is underestimated in high-gradient channels due to the dominance of flow acceleration on steep slopes. The sensitivity of the hydraulic models to simulate GLOF events varies from case to case depending upon the morphometric setting and input parameters.
It is therefore elementary to evaluate the sensitivity of these dynamic models and its simulated outputs to various hydraulic inputs and morphometric parameters.
Our study evaluated the sensitivity of the dynamic hydraulic model and its simulated outputs to various hydraulic inputs and morphometric parameters. For this, a series of hydrodynamic moraine-breach modeling and flood routing was performed to assess the sensitivity of flow hydraulics to different model input parameters and terrain characteristics. The sensitivity of the initial GLOF hydrograph to the total breach width (B w ) and breach formation time (T f ) was evaluated by considering six hypothetical moraine-breach scenarios with different combinations of B w and T f . It was evident that the GLOF hydrograph was more sensitive to the breach formation time (T f ) than the breach width (B w ). In order to evaluate the two-dimensional sensitivity of the model to the breach hydrographs, the modeled hypothetical GLOF hydrographs with varied T f were (Figure 8a) routed along the flow channel from the lake to Milam village. The analysis revealed that the initial breach hydrograph determines the flow hydraulics of a routed GLOF. The hydraulic properties of a GLOF at any point downstream is sensitive to the initial breach hydrograph.
The relationships between the flow hydraulics (flow velocity and inundation depth) and channel characteristics were evaluated. The flow velocity is linearly dependent on the slope of the given flow channel. However, the flow depths do not show any relation with the slope of the channel and tend to vary in a random manner. In addition, the flow hydraulics were analyzed with respect to the top-width of the channel. Both flow velocity and inundation depth reveal a similar trend of decrease with an increase in the top-width of the channel. The choice of DEM is also crucial in the dynamic routing of GLOFs, as it affects flow hydraulics (both flow depth and velocity) ( Figure 14).

Conclusions
The Himalaya hosts numerous glacial lakes and has experienced disastrous GLOFs in the past. Climate change-driven glacial retreat is the major reason behind the formation of such lakes. In this work, we analyzed the dynamics and evolution of Safed Lake (Uttarakhand, Central Himalaya), employed HEC-RAS to model several GLOF scenarios and evaluated the sensitivity of this model to changes in selected input parameters and DEMs used. Safed Lake is one of the largest glacial lakes, and underwent rapid glacier retreat-driven growth in past decades. We showed that the lake surface has grown from 0.10 km 2 (1968) to 0.23 km 2 (2018), and we estimated current lake volume to be 4.34 × 10 6 m 3 , with the potential for future growth. The modeled flow hydraulics of a potential GLOF event revealed the part of Milam village, located 16.2 km downstream the lake, to be at great risk in a worst-case GLOF scenario, considering the dam incision till the base of the moraine (60 m breach depth), releasing the entire volume of the lake. However, the impact of a potential GLOF significantly reduces when the breach depth is reduced to half (30 m) and quarter (15 m). GLOF resulting from 15 m breach depth presents no risk of flooding to the village. Due to the fact the moraine of Safed Lake has a flat geometry, it is highly unlikely that the breach depth would reach its maximum. However, site-specific investigation and monitoring of the lake, moraine dam and its structure, and a detailed investigation of the magnitude of possible GLOF triggers, are highly recommended for advanced hazard and risk assessment. Further, steps for flood protection and water management can be implemented.