Anomalous winter-snow-amplified earthquake-induced disaster of the 2015 Langtang avalanche in Nepal

Coseismic avalanches and rockfalls, as well as their simultaneous air blast and muddy flow, which were induced by the 2015 Gorkha earthquake in Nepal, destroyed the village of Langtang. In order to reveal volume and structure of the deposit covering the village, as well as sequence of the multiple events, we conducted an intensive in situ observation in October 2015. Multitemporal digital elevation models created from photographs taken by helicopter and unmanned aerial vehicles reveal that the deposit volumes of the primary and succeeding events were 6.81± 1.54× 106 and 0.84± 0.92× 106 m3, respectively. Visual investigations of the deposit and witness statements of villagers suggest that the primary event was an avalanche composed mostly of snow, while the collapsed glacier ice could not be dominant source for the total mass. Succeeding events were multiple rockfalls which may have been triggered by aftershocks. From the initial deposit volume and the area of the upper catchment, we estimate an average snow depth of 1.82± 0.46 m in the source area. This is consistent with anomalously large snow depths (1.28–1.52 m) observed at a neighboring glacier (4800–5100 m a.s.l.), which accumulated over the course of four major snowfall events between October 2014 and the earthquake on 25 April 2015. Considering long-term observational data, probability density functions, and elevation gradients of precipitation, we conclude that this anomalous winter snow was an extreme event with a return interval of at least 100 years. The anomalous winter snowfall may have amplified the disastrous effects induced by the 2015 Gorkha earthquake in Nepal. Published by Copernicus Publications on behalf of the European Geosciences Union. 750 K. Fujita et al.: Anomalous winter-snow-amplified earthquake-induced disaster


Introduction
The avalanches and rockfalls in Langtang Valley, about 70 km north from Kathmandu, resulted in one of the most remarkable tragedies in Nepal triggered by the 25 April 2015 Gorkha earthquake. Coseismic avalanches and rockfalls, as well as their simultaneous air blast and muddy flow, killed more than 350 people and destroyed the village of Langtang (Kargel et al., 2016). Satellite imagery and oblique helicopter photographs suggested that the deposit covering the village consisted of abundant snow and ice extending over an area of about 0.75 km 2 . The extent of muddy deposits on the opposite slope (∼ 200 m running up from the valley floor) provided an estimated debris speed of 63 m s −1 at the valley floor. Based on a mean thickness (≥ 2 m) and density (2200 kg m −3 ) of the deposited material, the released gravitational potential energy was estimated as ≥ 3.2 × 10 13 J, although the thickness was highly uncertain (Kargel et al., 2016). Utilizing a set of high-resolution SPOT6/7 images, Lacroix (2016) estimated a deposit volume of 6.95 × 10 6 m 3 over the village, with a larger amount of avalanche material deposited above 5000 m a.s.l. (9.66 × 10 6 m 3 ), and a total avalanche volume of 14.38 × 10 6 m 3 . Differences between two digital elevation models (DEMs) created from images on 21 April 2014 and 10 May 2015 showed distribution of not only the deposit but also destroyed forest on the opposite slope ( Fig. 5a in Lacroix, 2016) with errors not exceeding 3 m over the unaffected terrain. Nevertheless, the contributions and sequences of avalanches (snow and ice) and rockfalls are not yet fully understood, though some of the event was described in Kargel et al. (2016).
In this study, we aim to unravel the sequence of the multiple events that led to this tragic event and quantify the contributions of snow, ice, and rock by using in situ observation from October 2015 in combination with pre-event data and aerial photography directly after the earthquake. We created DEMs and orthomosaiced images (orthomosaics) from oblique photographs taken by unmanned aerial vehicle (UAV) and multiple helicopter flights in May and June 2015 together with geographical information obtained by differential GPS (dGPS) survey. Using a satellite-based highresolution DEM as the pre-event control, we estimate temporal changes in deposit volumes over the village. In addition, we investigated detailed structures of the deposit during the in situ observation from which we estimate the volumetric contributions of snow avalanche, collapsed glacier ice, and rockfalls. During the in situ survey and in Kathmandu, we interviewed villagers who witnessed the events in the valley and their narratives provide additional insights to the Langtang village tragedy. By analyzing meteorological data recorded at multiple sites in the valley, we reveal that anomalous winter snowfall prior to the earthquake might have amplified the coseismic hazards.

Field observation
To survey the deposit covering the Langtang main village (Fig. 1), we carried out a field campaign in October 2015, during which we performed several flights with different UAVs in combination with dGPS surveys to create DEMs and orthomosaics. We investigated the boundary between debris and ice, rock presence in the ice, and debris density in the ice at a few dozens of ice cliffs distributed near the Langtang River. We also observed size, condition, and age of trees flattened by the air blast and muddy flow on the opposite slope. During the in situ observation and in Kathmandu, we interviewed 20 Langtang villagers to record their narratives of the events.

UAV campaign
We utilized a fixed-wing UAV (Hobbyking SkyWalker X-5), which had a wing span of 1.18 m and a weight of 1.4 kg including camera. This plane is capable of approximately 45minute flights at cruise speeds of 60 km h −1 with either manual or autopilot modes. The camera mounted on the UAV, a Ricoh GR, has a lens of fixed focal length of 18.3 mm (equivalent to 28 mm focal length of 35 mm film camera) and 14.2 megapixel sensor, i.e., 4352 by 3264 pixels, capturing both JPEG and RAW format images in the visible light range. We performed two flights on 23 and 24 October 2015 for image acquisitions. We launched the UAV from an upstream/eastern hill of the village (Fig. 2a). The desired image overlap was set to be greater than 70 % in both lateral and longitudinal directions with respect to the UAV flight path. In order to generate the desired ground resolution finer than 0.1 m, the flight altitude was set 120-170 m from the ground.

dGPS survey and ground control points (GCPs)
In order to obtain precise positions of GCPs and validation data for generated DEMs, we conducted dGPS surveys around the village using dual-frequency carrier-phase GPS receivers (GEM-1 and 2, GNSS Technologies, Inc.; and R10, Nikon-Trimble Co., Ltd.). One receiver was set up at the Mundu campsite (Fig. 1a) as a base station and two other receivers were operated in kinematic mode with a logging interval of 1 s. The location of the base station was determined by the precise point positioning service operated by Natural Resources Canada (http://webapp.geod.nrcan. gc.ca/geod/tools-outils/ppp.php?locale=en, last accessed on 22 August 2016). We then post-processed the dGPS data with the RTKLIB software (http://www.rtklib.com/, last accessed on 22 August 2016), from which the dGPS points calculated with the fixed solution were used in the following analyses. For the evaluation of DEMs generated from multiple sources, we convert the dGPS point data into gridded DEM of 1 or 5 m resolution (depending on the targeted DEMs with which the comparison made) by the inverse distance weighted (IDW) interpolation method with a search radius of x/ √ 2 m (here x is the resolution). Grid cells with no dGPS point were excluded following the method proposed by Fujita et al. (2009. Tshering and Fujita (2016) tested alternative interpo-lation methods such as arithmetic averaging, spline, and kriging with various settings, and then they demonstrated that the height difference from the IDW method (as root-meansquare errors) did not exceed 0.3 m except for the spline method.
During the survey, we evenly distributed and surveyed eight orange fabrics of 1 × 1 m over the deposit-covered area (Fig. 2a). We additionally surveyed seven obvious rocks located outside the deposit-covered area as the GCPs.

Generation of the DEM and orthomosaic
We processed the UAV photographs (665 images, Table 1) into a DEM and orthomosaic using the structure from motion (SfM) workflow, which we implemented using the software package Agisoft PhotoScan Professional version 1.2.5 (Agisoft, 2016). Finally, we exported the DEM and orthomosaic with a resolution of 0.06 m. For estimating deposit volumes and comparisons with the dGPS data, we resampled the DEM to resolutions of 1 or 5 m, depending on the targeted DEMs used for comparison.

Aerial photographs
As our UAV survey was conducted 6 months after the Gorkha earthquake, some of the deposited snow and ice had likely melted away. Some coauthors of this study have conducted aerial photography from a helicopter on 7 and 10 May (D. F. Breashears) and on 1 June (H. Yagi). We processed these sets of photographs with the PhotoScan too. Unfortunately, no GPS measurements were acquired when the photographs were taken. In addition, the surface conditions drastically changed by the successive rockfalls and the subsequent melting. For georeferencing the imagery, we selected eight tie points over an area unaffected by massive deposit (off-debris area hereafter) from the UAV DEM/orthomosaic products (Fig. 2). Details of camera and lens are listed in Table 1. We processed these sets of helicopter photographs in Agisoft PhotoScan with the same parameter settings as the UAV DEM and then generated three sets of helicopter-based DEMs (heli-DEMs hereafter) and orthomosaics (7, 10 May and 1 June). For a better overall overview, we also processed a set of oblique aerial photographs taken in December 2007 and created a 3-D image of the Langtang village and Mount Langtang Lirung, which is a possible source area of the avalanche (Fig. 1b). The village location was manually delineated from the Advanced Land Observing Satellite (ALOS) Panchromatic Remote-sensing Instrument for Stereo Mapping (PRISM) image taken in December 2010 (background image of Fig. 1a) and the outlines of the debris/ice deposits are shown to indicate how the avalanche destroyed the village (Fig. 1b).

Pre-earthquake data
In order to estimate the deposit volumes and their changes over time, a pre-event DEM is required. We tested two DEMs (5 m resolution) derived from the 2.5 m resolution PRISM on board the ALOS. The first one is a DEM created with the ALOS PRISM image acquired on 3 December 2010 (A10DEM). We first automatically produced an ALOS stereo-model and DEM with rational polynomial coefficients (RPCs) data on the Leica Photogrammetric Suite (LPS) software. Such an automatically extracted DEM often contains many errors, especially in areas of highly irregular and abruptly changing topographies due to shadows and low contrast in the images. We therefore manually edited the point cloud by removing false mass points in the air or under the ground and by placing adequate and representative mass points exactly on the terrain surface assisted by the LPS Terrain Editor (Sawagaki et al., 2012). The second DEM is a commercial product named "ALOS World 3D" (AW3D), which is created by stacking and mosaicking using multitemporal ALOS digital surface models to achieve the target accuracy within 5 m in both horizontal and vertical (Tadono et al., 2014).
To minimize horizontal and vertical errors, we compare both ALOS DEMs of 5 m resolution with the UAV DEM over the off-debris area (Fig. 2). We horizontally shift the ALOS DEMs at a 1 m step along south-north and west-east directions, and then compared the result with the UAV DEM resampled from 1 m resolution into 5 m resolution. The elevation difference and its standard deviation of the ALOS DEMs against the UAV DEM are depicted in Fig. 3. The performance of the AW3D (the minimum standard deviation of 1.46 m, also listed in Table 2) is much better than that of the A10DEM (the minimum standard deviation of 2.63 m). The DEM based on single scene (A10DEM) has large errors in the shaded portion along the Langtang River, while stacking multiple images may reduce those errors (AW3D). Given this, we adopt the AW3D in the following analysis. Even in the AW3D, however, there are many void (no-data) areas at elevations higher than 5500 m a.s.l., where the glacier/snow surface lacks sufficient contrast to create stereo views. It is therefore difficult to evaluate the extent and amount of snow, ice, and rocks that detached in the avalanche source area. Figure 3. Distribution of elevation bias (dz, thick contour lines) and its standard deviation (color shadings with thin contour lines) of (a) AW3D and (b) A10DEM against the UAV DEM on 23 October, which are calculated over the off-debris area (Fig. 2a). ALOS DEMs are shifted at 1 m interval to both directions of west to east (de, horizontal axis) and of south to north (dn, vertical axis), and then mean elevations for 5 m cells (corresponding to the ALOS DEMs) are calculated from the UAV DEM, which is resampled to 5 m from the original 1 m resolution. Table 2. Accuracy evaluation of digital elevation models over offdebris area (Figs. 3 and 4). Average and standard deviations of elevation difference (unit is meter) are listed on the lower left and upper right in the table, respectively. Subtraction is performed as recent DEM (lower in the column or righter in the line) is subtracted by previous one. Cell size for the comparison with AW3D is 5 m while the others are 1 m. Comparison with dGPS data is shown in Fig. 4

Automatic weather stations
To understand the interaction between meteorology, glaciology, and hydrology, several automatic weather stations and precipitation gauges have been operated in the basin (e.g., Immerzeel et al., 2014;Ragettli et al., 2015; in addition to a long-term meteorological observation at Kyangjin, which was started during the mid-1980s (Fujita et al., 2006) (Fig. 1a). The region is monsoon dominated, and winter precipitation accounts for only 20-30 % of annual totals (Immerzeel et al., 2014;Shea et al., 2015). In order to measure both liquid and solid precipitation precisely at high elevation, advanced pluviometers were installed near Yala and Gangja La glaciers (YL_aws, YL_pv, and GL_pv). Along a slope towards Gangja La Glacier, two more tipping buckets with temperature sensors (not air but instrument tem-perature) were installed (GL_tb1 and GL_tb2). Air temperature and snow depth were also measured at the three sites (YL_aws, YL_pv, and GL_pv). All instruments, except for two tipping buckets, were damaged by a prior snow fall event or the air blast triggered by the Gorkha earthquake so that some parameters were not fully recorded. Precipitation, snow depth, and air temperature data are retrieved from the data loggers and analyzed to understand the meteorological conditions until and during the event.
2.8 Recurrence period of winter snow As described later, we find that the winter snowfall since November 2014 was anomalous. In order to evaluate how rare this winter snow event was, we calculate the recurrence period of the cumulative precipitation between 1 November to 25 April, which is a proxy for the total amount of snowfall. As the in situ precipitation data observed at Kyangjin ( Fig. 1a) were available for only 12 years, it was not possible to calculate a probability density function. We therefore adopt the Aphrodite precipitation data for the period 1978-2007 (29 winters), in which the observed data was utilized (Yatagai et al., 2012) and thus representative even without calibration. It has been reported previously that a vertical precipitation gradient exists along the slope between Kyangjin (3830 m a.s.l.) and Yala Glacier (5100 m a.s.l.) (Seko, 1987;Fujita et al., 1997;Immerzeel et al., 2014). As the station in Kyangjin malfunctioned in the winter of 2014-2015, we estimate cumulative winter precipitation from observations at Yala Glacier. Two precipitation gradients, 30.8 % km −1 in 2012/13 (Immerzeel et al., 2014) and 38.0 % km −1 in 1985/86 (Seko, 1987), which were observed precipitation ratios between Yala Glacier and Kyangjin in winter, were applied to the observations made near Yala Glacier.
To calculate the recurrence period of the 2014/2015 winter snowfall at Kyangjin, we tested major probability density functions such as Gumbel, generalized extreme value, lognormal, gamma, logistic, and normal. Among all functions, which were acceptable by the Kolmogorov-Smirnov test, the Gumbel distribution function showed the best representativeness for the long-term gridded precipitation.

Interviews of local villagers
We interviewed 20 villagers in the Langtang valley and in Kathmandu during our visit from October to November 2015. We asked where the interviewee was during the earthquake and the amount of time between the earthquake and the first avalanche that hit the Langtang village. In addition, we asked about the color, wetness, and surface condition of the primary deposit. We exclude hearsay statements to avoid fictional narratives, which may have resulted from the traumatic course of events. We also collected traditions with respect to the former huge earthquake in 15 January 1934 Figure 4. Evaluation of DEMs, expressed as (a) histogram and (b) box plot of elevation difference from dGPS data (dz, vertical axes). Evaluation was performed over the off-debris area (Fig. 2a), at which no/less debris was found, while "full domain" denotes the whole area at which both elevation data of UAV DEM and dGPS were obtained during our in situ observation in October 2015. The cell size of the DEMs is 1 m except for that of the AW3D (5 m). Also note the bias in the AW3D (left axes). Box, line, and circle in the box-and-whisker plot denote interquartile, median and mean, and 1 standard deviation, respectively.
(M w 8.0, epicenter at 26.86 • N, 86.59 • E; Singh and Gupta, 1980) to understand how the village has relocated and expanded. All villagers verbally consented to the interviews.

DEM evaluation
The dGPS survey method used in this study has a geodetic accuracy of ∼ 0.2 m both horizontally and vertically (Fujita et al., 2008. By analyzing the differences between all DEMs with the dGPS data in the off-debris area (excluding built-up area, Fig. 2), we assess the accuracy of the DEMs (dz, Fig. 4 and Table 2). The mean deviations of the UAV DEM (23 October) is ±0.26 m over the off-debris area and 0.34 m over the whole domain where elevation data are available. These errors are equivalent to or slightly greater than that from the dGPS method. Mean deviations of the other heli-DEMs over the off-debris area increase toward 0.46 m (1 June), 0.96 m (10 May), and 1.51 m (7 May), while deviation of the AW3D (1.21 m) is of the same magnitude as those of two DEMs in May. In addition, deviations among the DEMs, which are necessary to evaluate uncertainty in volume change, are rather large, ranging from 1.0 to 1.8 m ( Table 2).
The heli-DEM on 7 May shows biases ranging from −1.37 to −1.00 m against the other heli-/UAV DEMs (Table 2). Together with the positive bias of 1.54 m against the dGPS, this implies that the 7 May surface was higher than all other dates ( Fig. 4 and Table 2) and that the initial deposit was also present in the off-debris area. The spread of material outside the main debris deposit was also inferred from the bound-ary observed on the 7 May orthomosaic (Fig. 2d), although we selected obviously identifiable rocks as the tie points on all orthomosaics. The AW3D DEM suffers from a significant negative bias of −5.08 m against the dGPS. For estimating deposited volume, the bias of AW3D is therefore corrected by increasing the elevation of the AW3D by 4.81 m, which is an average of biases of three DEMs (10 May, 1 June, and 23 October) against the AW3D. Except for the AW3D and the heli-DEM of 7 May, biases among the other heli-/UAV DEMs are within 0.4 m and therefore these were not corrected.

Change in deposit volume
Based on the multitemporal DEMs, we evaluate thickness and volume changes of the deposit covering Langtang village ( Fig. 5 and Table 3). The average thickness and total volume of the deposit by the primary avalanche, estimated from elevation difference between the heli-DEM of 7 May and the AW3D over the deposit-covered area of 0.58 km 2 , are 11.31 ± 1.85 m and 6.55 ± 1.07 × 10 6 m 3 , respectively. Remarkable increases (> 30 m) of surface elevation (blue color in Fig. 5a), suggesting deposition of avalanche material, are found along a gorge (locally named Pääbe Chu) situated at the western side of the village and in the Langtang River bed situated south of the village. A decrease of surface elevation (red color in Fig. 5a) on the opposite slope (southwestern part) is likely to be the results from destruction of forest by the air blast and/or muddy flow. If the positive bias of the heli-DEM of 7 May found over the eastern off-debris area (0.25 km 2 , Figs. 2d and 5a) is taken into account as additional deposit, the total area and volume increase to 0.83 km 2 and 6.81 ± 1.54 × 10 6 m 3 , respectively.
Between 8 and 10 May 2016, satellite images indicated a secondary large mass movement over the village (Kargel et al., 2016). Our heli-DEMs show that surface over the deposit-covered area elevated by 1.45 ± 1.58 m as an average thickness, corresponding to additional volume of 0.84 ± 0.92 × 10 6 m 3 ( Fig. 5b and Table 3). In addition, Langtang villagers witnessed rockfall on 12 May, which could have been triggered by the largest aftershock of M w 7.3 (Kargel et al., 2016), although the volume of this third large event could not be quantified. However, the observed thickness of the debris mantle on ice cliffs (0.90 ± 0.87 m, locations shown in Fig. 2a) is thinner than the estimated thickness from the DEM difference (Figs. 2a and 6). This discrepancy may be due to the limited sample number and condition of the ice cliff exposures or the second mass movement also contained some ice and not only debris.
During the 3 weeks between 10 May and 1 June, and the roughly 5 months between 1 June and 23 October, the surface elevation lowered on average by 2.00 ± 1.02 and 2.98 ± 0.72 m, respectively, which is equivalent to volumes of 1.14 ± 0.58 and 1.64 ± 0.39 × 10 6 m 3 , respectively (Figs. 5c and d, and Table 3). Ice covered by a thick debris   (Table 2) are used to estimate the error in volume. For the first deposit (7 May -AW3D), the volume of 0.26 ± 0.47 × 10 6 m 3 over the eastern off-debris area (0.25 km 2 , Figs. 2d and 5a) is additionally taken into account, and thus the total volume of the primary avalanche is estimated to be 6.81 ± 1.54 × 10 6 m 3 (see also Sect. 3.2).

AW3D
7  mantle (> 0.5 m) is hard to melt due to the insulating effect of debris (Mattson et al., 1993). The observed pattern of surface lowering suggests that ice melting likely occurs at the bottom of the Pääbe Chu and of Langtang River (Figs. 5c and d) due to thermal erosion of the water flowing underneath.

Cheese boulder
Based on a suggestion and guided by two of the Langtang villagers, we surveyed a large boulder (Fig. 7), on which sculptures of cheese and bread were carved in 1998 as the landmark for a cheese factory developed in the center of the village (hereafter referred to as the "cheese boulder") ( Fig. 7a). The cheese boulder was found on a terrace between the village and Langtang River (Fig. 7b). We surveyed both the original and found locations by the dGPS.
The original location was difficult to identify as the scenery was completely changed by the deposit so we relied on guidance of the villagers. We took photographs of the boulder from multiple angles and positions and then processed them in Agisoft PhotoScan to estimate the volume of the boulder (Fig. 7c). We estimate density of the boulder by measuring weight and volume of 10 samples of small stones of the same material as the cheese boulder. We calculated the volume of the cheese boulder to be 1.67 ± 0.21 m 3 . With a measured density of 2.65 × 10 3 kg m −3 , a value of typical granite (Lillie, 1999), we determined the mass of the boulder to be 4.4 ± 0.2 × 10 3 kg. The villagers assert that the boulder was relocated by the blast because many farm animals, including horses and yaks, were also blown away. We therefore estimate initial speeds of the boulder at the original location though it is unknown whether the air blast or muddy flow moved the boulder. We assume three possible scenarios, which fit both original and found locations of the boulder: (1) no touchdown between both original and found locations (78 m s −1 , upper curve in Fig. 7d); (2) overshooting the neighboring 4.5 m height hill and then rolling down the slope (106 m s −1 , lower curve in Fig. 7d); (3) an eject angle of 45 • giving the minimum speed (58 m s −1 , not shown in the figure). It should be noted that the mass of boulder is not required to calculate the initial speeds. The estimated initial speeds are reasonably consistent with the speed of muddy flow (63 m s −1 ), which was estimated from run-up height on the opposite slope (Kargel et al., 2016). This boulder movement can be utilized for validation of avalanche simulations in further study.

Other in situ data
The extent of flattened forest on the opposite slope evidenced the remarkable impact of the air blast (Kargel et al., 2016). Field measurements of size, age, and genera of fallen trees provide more insight into this event (Fig. 8a). The mean and standard deviation of the diameter of 113 fallen trees are 0.16 ± 0.07 m (Fig. 9). Most of them were uprooted because of thin soil layer on the bedrock. The average age of these trees, which consist of Abies and Rhododendron species, is estimated to be around 35 years from tree ring counts. It suggests that no event of this magnitude has occurred in the recent decades. Although a 30-year mass balance record of Yala Glacier reconstructed from multiple sources suggests large annual accumulations (solid precipitation > 1200 mm year −1 , greater than +1.4σ ) in 1975, 1976, 1980, and 1981(Fujita et al., 2006, it is unknown whether the precipitation fell as winter snow, which could have led to similar destructive avalanches, or as summer rain, which would not have affected the forest. There is no other record or villagers' narrative suggesting any event to destroy trees in 1980 or earlier. We performed visual investigations of the boundary between snow/ice avalanche and rock deposits, the thickness of debris on the avalanche deposit, and other features of the ice deposit. The boundary between ice and debris deposits at 20 exposed ice cliffs distributed near the Langtang River was clearly distinguishable (Fig. 8b), suggesting that the primary snow and ice avalanche and subsequent rockfalls occurred at different timings. Visual investigation of the ice  tunnel and ice cliff over Langtang River also supports that these occurred as different events because no rocks greater than 1 m were found in the ice cliffs (Fig. 8c). To measure the density of sand contained in the ice deposit, we took six samples (weight ranging from 0.9 to 1.7 kg) from ice cliffs. On average the percentage of sand and silt in the ice was 3.9 ± 0.76 % in weight. Within the sand-rich avalanche deposit, we found rounded "clear ice balls", whose diameter ranged on the order of several centimeters (Fig. 8d and e), suggesting that this deposit could have different sources. We speculate that the dirty ice originated from entrained sand as a result of the snow avalanche, whereas the clear ice balls were derived from detached glacier ice that broke into small fragments without entraining dirty materials. Therefore, the ice cliffs we investigated were originally composed mostly of snow avalanche deposits, which was compressed, melted, and refrozen to approximate ice.

Meteorological data
In mid-October 2014 (14-15 October), a heavy snowfall event hit the valley during the passage of extratropical Cyclone Hudhud (Wang et al., 2015). Although this event was recorded in both the precipitation and snow depth records at all sites, the precipitation totals vary between sites while the snow depth data are more consistent at the different stations (Fig. 10a). Therefore, we calculate cumulative precipitation since 1 November by excluding records in October (Fig. 10b). Snow depth and cumulative precipitation show that, including Cyclone Hudhud, four major snowfall events occurred in mid-October, mid-December, early January, and early March (Fig. 10a and b). The instrument at Gangja La (GL_pv) stopped recording data on 2 January 2015, most likely because of an avalanche caused by a heavy snowfall event. However, this station supports the consistency of records among the meteorological sites and therefore it is included in the analysis. Precipitation records at the two Gangja La sites showed much less precipitation than those near Yala Glacier (Fig. 10b). This is partly because of the elevation gradient of precipitation (Immerzeel et al., 2014) but mainly because the tipping bucket tended to miss solid precipitation in winter. Comparing long-term averages of cumulative precipitation at Kyangjin (12 years, blue line with shading) and the Aphrodite precipitation data (29 years, purple line with shading; Yatagai et al., 2012), the amount of precipitation during the 2014/2015 winter was anomalously large.
Compared with the long-term average of air temperature at Kyangjin (3830 m a.s.l.), the winter air temperature recorded near Yala Glacier (4830 m a.s.l.) was 0.5 • C above average if a local and seasonal lapse rate of air temperature is taken into account (−5.8 • C km −1 ; Immerzeel et al., 2014). However, the air temperature near Yala Glacier (−7.6 • C for the 4 months from December to March) suggests that no snowmelt occurred above 5000 m a.s.l. throughout the winter (Fig. 10c).
Temperatures measured inside the tipping bucket sensors decreased after the Gorkha earthquake and associated avalanche (green arrow in Fig. 10d). Compared with two averages of diurnal change under normal conditions (blue lines with shadings), observed temperatures significantly dropped was synchronous with the timing of the earthquake and associated avalanche. The YL_pv record stopped because of the destruction due to another strong blast from a separate but coincident avalanche (orange line). It is unclear whether the observed temperature decreases were caused by air temperature cooling in the whole valley or by snow accretion to the instruments.

Recurrence period of the anomalous winter snow
We calculate winter snowfall recurrence periods to establish how rare the 2014-2015 winter snowfall was (Fig. 11). Although the recurrence period highly varies depending on various factors such as representativeness of the Aphrodite precipitation data and elevation gradient of precipitation, it suggests that the winter snowfall between November 2014 to 25 April 2015 was an extreme event which occurs only once during more than 100-500 years. Considering that the Aphrodite precipitation data are greater than the observed . Blue lines with shadings in (a) and (c) denote long-term means observed at Kyangjin (KJ_aws, 1988(KJ_aws, -1992(KJ_aws, , 2002(KJ_aws, -2008. Purple line with shading in (a) denotes long-term means of Aphrodite precipitation at a grind including Kyangjin (PRa, 1978(PRa, -2007. (d) Diurnal temperatures recorded by a sensor embedded in a tipping bucket at two Gangja La sites. Thick green lines and thin blue lines with shadings denote temperatures on 25 April 2015, the date when the Gorkha earthquake occurred, and means created from weeks before and after the earthquake (14 days in total but excluding 25 April). Also shown is air temperature recorded at Yala site (orange line), which was destroyed by wind blast or dense flow. Locations of the site are shown with abbreviations in Fig. 1a. Green arrow in (d) denotes the temperature drop after the earthquake. one at Kyangjin (Fig. 10b) and that Cyclone Hudhud supplied more snow in October, the recurrence period of the anomalous winter snow could be greater than those evaluated above. Figure 11. (a) Cumulative density function of cumulative precipitation from 1 November to 25 April (PR a , purple circles) and the Gumbel density functions (black line); (b) recurrence interval. Red and orange circles denote calibrated precipitation for the elevation of Kyangjin (3830 m a.s.l.) from two sites near Yala Glacier (Fig. 1a), assuming two elevation gradients for precipitation.

Sequence of the events
Together with witness statements of villagers, we reconstruct the sequence of events triggered by the Gorkha earthquake as follows: 1. Triggered by the main shock of the Gorkha earthquake, a predominantly snow avalanche, with an estimated ice volume of 6.81 ± 1.54 × 10 6 m 3 , hit the Langtang village and covered the area of 0.58 km 2 with an average thickness of 11.3 m (Figs. 2d, 5a, and 12, Table 3). East and north of the main deposits, debris with an average thickness of 1.5 m (Fig. 5a) was spread over the village area of 0.25 km 2 . Villagers witnessed that multiple avalanches hit the village within a short time frame and it is difficult to attribute volumes to individual avalanches.
2. Glacier ice detached by the main shock could have triggered the avalanche although the detached volume is much smaller than the estimated deposited volume (Sect. 3.4 and Fig. 8d). 3. The avalanche-deposited material also consisted of snow-entrained fine particle debris such as sand and silt, which is estimated to be about 3-5 % in weight density. The sand and silt were entrained during descent of the snow and ice over several thousands of meters (Sect. 3.4 and Fig. 8c) and most of villagers witnessed "black avalanches". Although some villagers observed rockfall after the first avalanche, witness statements about the slippery ice surface of the Langtang village and the photographs taken on 7 May (D. F. Breashears) indicate that rockfall was probably limited. This was already pointed out using the same oblique photographs and high-resolution satellite images (Kargel et al., 2016).
4. Because of large elevation difference of 2500 m between source area (∼ 5890 m a.s.l.; Fig. 13) and the village (∼ 3400 m a.s.l.), gravitational potential energy should melt snow and result in the "wet avalanche deposit" that most of villagers witnessed. If all released energy was used to melt snow, water content is esti-760 K. Fujita et al.: Anomalous winter-snow-amplified earthquake-induced disaster mated to be 7.3 % by following equation: where θ is water content in weight (%), l m is latent heat for ice melt (333.5 J kg −1 ), g is gravitational acceleration (9.8 m s −2 ), and h is elevation difference (assumed to be 2500 m). Even if the cold snow (−7.6 • C as the winter mean temperature) is assumed, the estimated water content is 6.7 %. It should be noted, however, that not all of the released gravitational potential energy would have been used to melt the deposited ice. Some amount of energy would have been released in the air blast and muddy flow. Considering the volume of ice deposit (6.81 ± 1.54 × 10 6 m 3 ) and ice density of 900 kg m −3 , this further implies that the energy of 1.50 ± 0.34 × 10 14 J was released, which is 4.7 times greater than the previous estimate (3.2 × 10 13 J by Kargel et al., 2016). Kargel et al. (2016) already pointed out, the released energy generated into a strong air blast and muddy flow, which swept away vulnerable nearby buildings and flattened mature forest on the opposite slope, which consisted of Abies and Rhododendron species with an average age of 35 years. Most of the trees were uprooted from the thin soil layer.

6.
The initial speed of the cheese boulder could have been between 58 and 106 m s −1 , depending on a number of assumptions (Fig. 5d). This estimate is reasonably consistent with the other estimation of the minimum speed of the muddy flow at the valley bottom (63 m s −1 ) based on the run-up height of muddy material on the opposite slope (Kargel et al., 2016) 7. Some villagers witnessed that the deposit covering the village also temporarily dammed the Langtang River but an ice tunnel formed quickly and the river water drained without causing an outburst flood.
8. Between 8 and 10 May, coseismic rockfalls occurred and the surface covering the Langtang village changed from ice and snow to rock debris. The event resulted in approximately 0.84 ± 0.92 × 10 6 m 3 of debris and rocks and elevated the deposit surface by 1.45 ± 1.58 m in average (Figs. 2c and 5b, Table 3). This is also evidenced by the photographs taken on 10 May (D. F. Breashears) as Kargel et al. (2016) described.
9. A villager who had been searching for the missing villagers witnessed that additional rockfalls, which were triggered by the largest aftershock on 12 May, hit and covered the village. However, its thickness and volume are not quantifiable.
10. During the succeeding 3 weeks between 10 May and 1 June, the surface elevation lowered by 2.00 ± 1.02 m, suggesting that 1.14 ± 0.58 × 10 6 m 3 of ice melted away (Figs. 2b and 4c, Table 3). If the additional rock debris deposited on 12 May is taken into account, by assuming average thickness of 1 to 2 m based on the elevated surface by the rockfall between 8 and 10 May as addressed above, the melted ice could be 1.5 to 2.0 times greater than the above estimation. During roughly 5 months, including monsoon season, the surface elevation lowered by a further 2.98 ± 0.72 m, suggesting that 1.64 ± 0.39 × 10 6 m 3 of ice melted away (Figs. 2a and 5d, Table 3). The greatest surface lowering occurred along the Pääbe Chu between Gomba and the village, as well as along Langtang River (Figs. 5c and d). Because ice covered with thick debris (> 0.5 m) melts slowly due to the insulation effect of the debris mantle (Mattson et al., 1993), the high melt rates are probably caused by thermal erosion of the water flowing underneath. The debris-covered area was slightly reduced to 0.54 km 2 during the monsoon, which was mainly due to change in surface condition on the opposite slope.
11. By the end of October 2015, 6 months after the main shock, about 60 % of the initial volume of deposited ice and snow remained over the village (6.81 ± 1.54 × 10 6 m 3 as initial volume and −2.73 ± 0.68 × 10 6 m 3 as melted volume).

Avalanche sources
A set of pre-and post-event images of SPOT6/7 successfully quantified the mass movement in the Langtang valley due to the earthquake-induced avalanche (Lacroix, 2016). The main deposit covering the Langtang village was estimated to be 6.95 × 10 6 m 3 , which is consistent with the volume evaluated in this study from different data sources (6.81 ± 1.54 × 10 6 m 3 ). At high elevations the imagery was of high quality (clear and with a lot of contrast) ( Fig. 5b in Lacroix, 2016) such that the estimated detached volume (14.38 × 10 6 m 3 ) seemed to be reliable. It is problematic, however, that the pre-event image was taken on 21 April 2014, 1 year earlier than the event, so that snowfall on and melting of the glacier surface during the 2014 monsoon season were fully ignored in the differentiation of two DEMs. Given the fact that the amount of snow prior to the event was anomalous, it is likely that the detached volume should have been greater than 14.38 × 10 6 m 3 . Another potential source of error is glacier dynamics: submergence in the accumulation zone and emergence in the ablation would affect the detached volume calculations (Cuffey and Paterson, 2010). In addition, change in density from snow to ice was not taken into account. Although other images available on Google Earth were carefully checked to confirm that no landslides had occurred during the 2014 monsoon season,  (Kargel et al., 2016). c Villagers witnessed that the rockfall was equivalent to or greater than those on 7-10 May. Figure 13. Hypsometry and box plot for the upper catchment above the Langtang village, from which the main avalanche might have occurred (orange polygon in Fig. 1a). Also shown is that limited above 5000 m a.s.l. Elevation data are from ASTER GDEM2 (Tachikawa et al., 2011). Box, line, and circle in the box-andwhisker plot denote interquartile, median and mean, and 1 standard deviation, respectively. none of the above issues were addressed in the source region of the main avalanche.
To further constrain the source of the avalanche material, we performed visual inspections of the pre-event (22 January 2015) and post-event (3 May 2015) satellite images on Google Earth and the 3-D image product created from the helicopter oblique photographs taken by D. F. Breashears on 7 May 2015 (not shown as figures). We found evidence for a collapsed block of hanging glacier ice and roughly estimated the volume of the ice block to be 66 × 10 3 m 3 (110 m × 20 m × 30 m) by assuming a glacier thickness of 30 m, which is similar to that observed for nearby Yala Glacier (Sugiyama et al., 2013). However, other possible areas for the avalanche source are either highly distorted in the satellite image (due to insufficient orthorectification of the post-event image of 3 May 2015) or unclear in the 3-D image because of the upward viewing angle of the photographs out of a helicopter and the reduced contrast on snow-covered areas (7 May 2015), which complicates the evaluation of the contributing volume of the collapsed glacier ice. Although the earthquake-triggered glacier collapse may have played a key role in initiating the entire avalanche event, we think that the collapsed glacier ice is not the dominant source for the avalanche. This hypothesis is further supported by the in situ observation of glacier-derived clear ice balls in the dirty ice deposit (Fig. 8d).
Meteorological data suggest that multiple snowfall events during the 2014/2015 winter led to an accumulated seasonal snowpack with a thickness of up to 1.5 m above 5000 m a.s.l. (Fig. 10a). Limiting the source area of the avalanche to areas higher than 5000 m a.s.l. (8.73 km 2 , Fig. 13), above which no snowmelt can be assumed in winter (Fig. 10c), the deposit over the Langtang village (6.81 ± 1.54 × 10 6 m 3 ) is equivalent to 0.78 ± 0.18 m of ice thickness over the entire source area. Using a density of the winter snow (385 ± 9 kg m −3 ) obtained from snow depth and cumulative precipitation at the Yala Glacier sites (Figs. 10a and b), the snow thickness and volume over the avalanche source area are estimated to be 1.82 ± 0.46 and 15.89 ± 4.05 × 10 6 m 3 , respectively. This snow thickness is reasonably consistent with the observed snow depths near Yala Glacier (1.28 m at 4830 m and 1.52 m at 5100 m, Fig. 10a). Conversely, if all of the winter snow came down during the avalanche, then this snow volume (12.22 ± 1.05 × 10 6 m 3 ) is the equivalent of 5.23 ± 0.44 × 10 6 m 3 deposited ice using the density conversion. The residue of 1.58 ± 1.10 × 10 6 m 3 can then potentially be attributed to the detached glacier ice, which is equivalent to 23 % of the mass of the primary avalanche (Table 4). It should be noted that volume and mass during entrainment and deposit during the avalanching are unknown. Another point of discussion is whether the entire slope became instable by the earthquake. Although the volume detached from the upper catchment (14.38 × 10 6 m 3 ) estimated by Lacroix (2016) seems consistent with the estimate in this study, this could be a mere coincidence if both gain and loss of mass on a typical Himalayan glacier in summer season (Ageta and Higuchi, 1984), the glacier dynamics (Cuffey and Paterson, 2010), and the winter snow are taken into account.

Anomalous winter snow
Although there is no doubt that the 2015 Gorkha earthquake triggered the primary snow avalanche that destroyed the Langtang village, the anomalous winter snow, which was an extreme event with a recurrence period of 100 to 500 years, likely amplified the disaster that took place. According to villagers statements, coinciding with the previous large earthquake (M w 8.3) recorded in 1934 (Singh and Gupta, 1980), ice avalanches hit the Numthang and Mundu areas, which were former sites of the main village (Fig. 1a). Since then, Langtang village expanded towards its present location, and a recent surge in tourism and trekking has accelerated this expansion since the 1990s. It is evident that the unlikely combination of an anomalously deep winter snowpack and largemagnitude earthquake contributed to the magnitude of the Langtang village disaster. It is still unclear whether such an amount of anomalous winter snow could cause such a huge avalanche by itself without an earthquake trigger. It should be noted that summer precipitation, which accounts for 70-80 % of annual precipitation totals (Immerzeel et al., 2014;Shea et al., 2015), could not result in such a large avalanche volume. The equilibrium line altitude, the elevation where accumulated snow should be melted away at the end of summer, fluctuates around 5300-5400 m a.s.l. and has risen in recent decades (Fujita and Nuimura, 2011).

Conclusions
We use multiple lines of evidence to re-construct the series of events leading to the devastating avalanche that occurred in Langtang Valley, Nepal, during the Gorkha earthquake. Through intensive in situ observation in October 2015, we evaluate the volume and structure of the deposit covering the Langtang village and the sequence of the multiple events. Multitemporal digital elevation models created from photographs taken by helicopter and unmanned aerial vehicle together with differential GPS survey reveal that the deposit volumes of primary and succeeding events were 6.81 ± 1.54 × 10 6 and 0.84 ± 0.92 × 10 6 m 3 , respectively. Visual investigations around ice cliffs exposed around Langtang River and witness statements of villagers suggest that the primary event was snow avalanche; the contribution of the collapsed glacier ice could not be the dominant source for the entire mass of avalanche, though it probably played a key role in initiating the entire event. Our observations also indicate that successive rockfalls, possibly triggered by aftershocks up to nearly 3 weeks after the initial earthquake, buried the initial deposit. This study would provide valuable constraints on avalanche dynamics models.
An average snow depth (1.82 ± 0.46 m) estimated from (a) the deposited volume, (b) the upper catchment area above 5000 m a.s.l., and (c) densities of snow and ice is reasonably consistent with snow depths (1.28-1.52 m) observed at nearby meteorological stations (4830-5100 m a.s.l.). Both amounts are anomalously large, as a result of four major snowfall events that occurred since October 2014. Considering the vertical gradients of winter precipitation (e.g., Collier and Immerzeel, 2015), the 1300 m difference in elevation between the lower village (where long-term observational data are available) and the high-elevation stations, as well as a range of probability density functions, we conclude that this anomalous winter snow was an extreme event, which occurs only once every 100 to 500 years. The anomalous snow conditions likely amplified the most remarkable disaster induced by the 2015 Gorkha earthquake in Nepal. To support the village and community reconstructions in Langtang Valley, an avalanche hazard map based on the facts revealed through this study will be provided to the villagers' community and local authorities.