Seismically Induced Snow Avalanches: The Central Italy Case

Snow avalanches are critical events due to the sudden instability of snow in mountain regions. Seismically induced snow avalanches, due to a particular type of trigger perturbing snow-covered slopes, are occasional events in geodynamically active mountain regions around the world, which can lead to large anomalous avalanches. The seismic trigger of an avalanche is either simultaneous or delayed. In Italy, the Alps, and particularly the Apennines, feature large areas expected to be affected by seismically induced snow avalanches. In this work, the avalanches’ distribution in Central Italy (i.e., Abruzzo region) and its connection with moderate–high seismicity were investigated. With regional and historical analyses, we investigated the connection between avalanches and seismicity. We revealed two main events of seismically induced avalanches (April 6–9, 2009, western Gran Sasso-Assergi area, and January 18, 2017, Gran Sasso massif) characterized by an anomalous distribution of avalanches. A specific and detailed analysis focused on the northeast Mt. Siella slope, where well-documented and constrained events occurred on January 18, 2017 in the Gran Sasso area, affecting a tourist facility. This involves geomorphological and photogeological multitemporal investigations. Moreover, the Rigopiano avalanche was analyzed through 3D numerical modeling and stability analysis of the snow to investigate the possible seismic trigger. We found the avalanche was consistent with delayed seismically induced avalanche cases. This combination contributed to a more general assessment of the seismically induced avalanches in the Central Apennines area (i.e., Abruzzo region) and provided new examples of seismically induced avalanches at both the regional and detailed scales, which are largely consistent with the world distribution.

Seismically induced snow avalanches feature a particular type of trigger that induces a perturbation on snow-covered slopes. This type has been analyzed with several different approaches at the regional and local scales and with analogue and numerical modeling (Hackman, 1968;Singh and Ganju, 2002;Podolskiy et al., 2010b;Pérez-Guillén et al., 2014;Wang et al., 2015;Bulajić et al., 2018;Parshad et al., 2019;Puzrin et al., 2019). Regional analysis using a statistical approach allowed the inference that snow avalanches can be triggered by an earthquake of Mw > 2.0 (magnitude moment), up to a distance of more than 600 km, both during seismic shaking, and with a certain delay (Singh and Ganju, 2002;Podolskiy, 2010;Podolskiy et al., 2010a;Puzrin et al., 2019).
The delay in the seismic trigger of an avalanche is one of the most complex aspects to consider. The non-contemporaneity of the cause and effect could be understood as a lack of connection between the two events. A possible temporal delay from an earthquake to an induced avalanche was documented in the area of western Himalaya ( Figure 1A) (Singh and Ganju, 2002). Several cases were reported that highlighted the delay times of some avalanches from the seismic shaking (e.g., in Western Himalaya: March 4, 1996 Lagongama, delay 350 min;January 26, 1999 Chadan, 49 min;January 16, 2000 Drass Sector, 7 min; January 21, 2000 Jawala area, 82 min; and other cased less constrained; for more details see table 1 of Singh andGanju, 2002 andPuzrin et al., 2019). The release mechanisms were recently summarized by Schweizer et al. (2016); the main mechanism of the delay in earthquake-induced avalanches release was also investigated using a new analytical model that produces dynamics coherent with the Rigopiano case (Puzrin et al., 2019). The rupture mechanism occurs due to the interaction between creep phenomena within the mass about to slip, strain softening, and the strain-rate sensitivity of the snow, which guide the growth of a basal shear fracture, as explained by Puzrin et al. (2019) (Figures 1B-G). A broad range of realistic parameters was used for three potentially earthquake-triggered delayed avalanches (Gaume et al., 2013(Gaume et al., , 2018. The Alps and the Apennines in Italy, according to the global distribution reported by Podolskiy (2010), show large areas expected to be affected by seismically induced snow avalanches ( Figure 2A). The eastern side of central Italy (i.e., the Abruzzo Region) features the highest elevation mountain area in the entire Apennines chain, is frequently affected by deep snow cover and snow avalanches, and is also affected by strong seismicity. In the last decade, the 2009 L'Aquila earthquake (max. Mw 6.1) and the 2016-2017 central Italy seismic crisis (max. Mw 6.5) affected the area and, incidentally, had strong seismic shocks in winter and early spring (i.e., April 6-9, 2009, Mw 5.4-6.1, and January 18, 2017, the last Mw > 5 sequence of the 2016-2017 seismic crisis, Figure 2B (INGV, 2016; Chiaraluce et al., 2017;Luzi et al., 2017)). In mid-January 2017, the northeast Apennines flank and the entire Gran Sasso area was affected by a heavy snow precipitation event (snow cover greater than 3 m deep in the mountain area and 1-1.5 m in the hilly piedmont) (Di Lena and Giuliano, 2017;. After this combination of seismic and meteorological events, from January to February 2017, several avalanches, landslides, and flooding events occurred from the NE Apennines chain to its piedmont area and the Adriatic coast (Calista et al., 2019). These events caused severe damage, fatalities, and evacuees, and were declared an official natural disaster at the regional government level (Regione Romeo et al., 2018). For this reason, we focused on the investigation of these events and specifically on the general connection between avalanches and seismicity in central Italy.
Previous studies on these events mostly focused on the geology of the earthquakes and seismically induced surface deformation and landslides. Specific investigations particularly concerning avalanche features, field analyses, numerical simulations (Bocchiola et al., 2018;Frigo et al., 2018), and forensic applications , as well as analysis of the triggering conditions (Puzrin et al., 2019), also focused on the NE slope of Mt. Siella, where well-documented avalanches occurred on January 18, 2017 after four main seismic events (e.g., the catastrophic Rigopiano event, which affected a tourist facility). Other studies examined the landslides triggered in the same period in NE the Abruzzo area (Solari et al., 2018;Calista et al., 2019 and references therein). However, an overall analysis of the snow avalanches in the last decades and their connection with the seismicity, as well as an overall geomorphological detailed investigation integrated with numerical modeling of the Mt. Siella and Rigopiano case, are still lacking.
We attempted to fill this gap by analyzing the avalanches' distribution in the Central Apennines and its connection with moderate-to-high seismicity. The study was based on an integrated analysis, including snow depth data, seismicity data, and regional avalanche distribution analysis. The avalanche distribution was investigated and compared with the nivological features, and the seismic events occurred in winter. The analysis revealed two main events in April 2009 in the western Gran Sasso (Assergi) and January 2017 in the entire Gran Sasso massif. A specific detailedscale analysis then focused on the avalanches on the NE slope of Mt. Siella (eastern Gran Sasso massif). The study was performed through a geomorphological field investigation and photogeological analysis through multitemporal and multiscale air photo interpretation that allowed the definition of the main features and the temporal distribution of the avalanches. Geomorphological and photogeological data were enhanced with numerical modeling for investigating the possible seismic triggers. The detachment was investigated using stability analysis and 3D modeling in static (non-earthquake) and dynamic (earthquake) conditions using FLAC3D code (Itasca, 2019).
The main aim was to contribute 1) to the comprehension of the overall distribution and features of the avalanches in central Italy and the connection to the seismic events in April 2009 and January 2017, and possible previous events, in the framework of the global distribution of seismically induced avalanches, and 2) to the understanding of the 2017 Mt. Siella avalanches as delayed seismically induced avalanche in the framework of central Italy and the world.  Table 1 and Figure 2b). For each site and each year, the maximum snow depth was derived and plotted. The record period is 1935(ISPRA, 2018. (B) The average snow depth vs. snow gauges elevation in the Gran Sasso area. For each site and each year, the average maximum snow cover was derived and plotted. (C) The number of avalanches in the Abruzzo region from 1957 to 2017 (Regione Abruzzo, 2014), compared to the maximum snow cover recorded at all the snow gauges (Table 1) and to the main seismic sequences or moderate earthquakes (Mw > 4) occurring in winter.
Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 599611 5 area ( Figure 2B). The chain is composed of thrust ridges, made up of thick calcareous sequences with different thickness and erodibility, pertaining to Mesozoic-Cenozoic pre-orogenic paleogeographic domains. The Neogene deformation of these sequences caused the emplacement of the main mountain ridges including the major Gran Sasso ridge (Vezzani et al., 2010;Pace et al., 2014 and references therein). Pleistocene-Holocene extensional tectonics, mostly along NW-SE-to NNW-SSEoriented faults, defined the present landscape of the mountain area and the main intermontane basins, and is responsible for the recent and historical seismicity ( Figure 2B). The morphology of the chain is asymmetric, with its highest peaks (Gran Sasso, 2,912 m a.s.l.; Maiella, 2,793 m a.s.l.) on the NE side. A steep mountainside with large escarpments features the NE front of the Apennines, which heavily affects the climatic-meteorological features. The landscape shows a series of ridges trending from NW-SE to N-S (heights of 2,000-3,000 m), separated by narrow valleys parallel to the ridges or by broad intermontane depressions (elevation 250-1,000 m). The Gran Sasso ridge is the highest in the Apennines, with several peaks above 2,500 m a.s.l. It is arched-shaped, trending from W-E to N-S, and drops down to the piedmont area to an elevation below 1,000 m, defining a large and steep mountain escarpment. The mountain area is mostly affected by mass wasting, gravityinduced (e.g., rock falls, topples, flows, etc.), and fluvial-related processes (e.g., debris flows, alluvial fans, etc.). Ancient glacial processes are preserved as relict landforms (Miccadei et al., 2017). The area is largely characterized by deep snow cover during winter up to several meters thick. This has induced snow related processes, as recurrent avalanches documented for the last 60 years by the Historical Avalanches Map of Abruzzo ( Figure 2B) (Regione Abruzzo, 2014). In April 2009 and January 2017, avalanches occurred mostly across the NE side of the chain and the piedmont area.
The piedmont and coastal area feature a low relief area (heights ∼ 100-800 m) characterized by shaly-arenaceous rocks pertaining to Neogene sequences (Bigi et al., 1997) arranged in a cuesta and mesa relief incised by fluvial valleys and wide alluvial plains (D'Alessandro et al., 2003;Miccadei et al., 2017). The area is affected mainly by mass wasting, gravity-induced, and fluvial processes frequently triggered by heavy rainfall and seismic events and occasionally by heavy snowfall and rapid snowmelt (Piacentini et al., 2018;Calista et al., 2019;Carabella et al., 2019;Calista et al., 2020).

Climate, Nivological, and Meteorological Features
The central Italy and Abruzzo climate generally ranges from Mediterranean to temperate. However, it is affected by the orography of the Apennines chain and its eastern front, changing from a Mediterranean type with maritime influence along the coasts and the hilly piedmont area, to more temperate and continental-like in the inland (Peel et al., 2007;Di Lena et al., 2012). The average temperature is 16-18°C along the coast to 8-10°C in the mountain areas (with minimum values of 0-5°C, at high elevation). The winter temperature (average January values) is 8-10°C in the hilly coastal area to 0-2°C in the inland areas (with minimum values of approximately −5 to −10°C at high elevations). The rainfall distribution is also controlled by topography. The highest average annual rainfall (up to 1,500-2,000 mm/year) occurs along the main ridges and drop down to 600-800 mm/year in the piedmont hills and coastal areas. The rainfall distribution is frequently characterized by occasional heavy rainfall (>100 mm/day and 30-40 mm/h) (Chiaudani et al., 2013;Piacentini et al., 2018). In recent decades, the region has been affected by flash flood events and some heavy snowstorms, which were induced by heavy rainfall ranging from 60 to 100 mm in a few hours to >200 mm per day, and by snowfall up to >1 m/day (e.g., January 2003, April 2004and October 2007, March 2011, September 2012, December 2013, February-March 2015, and January-February 2017Piacentini et al., 2018;Calista et al., 2019).
The historical nivological data were analyzed specifically for the Gran Sasso and L'Aquila area, in which most of the avalanches occurred ( Figure 2B). The snow cover depth recorded by nine snow gauges was analyzed (1935-2003 data sets; Figure 3A). The maximum snow depth recorded ranges from 0 (from many years, at specific locations, at lower elevations) to a maximum of 200-210 cm (Castel del Monte, early 1970s). For many years in the time series, the maximum snow depth was far below 100 cm. The average maximum snow depth ranges from 17 to 66 cm for the entire time series  with a good elevation/cover correlation showing a gradient of ∼50 cm/ 1,000 m ( Figure 3B). However, according to the elevation distribution of the stations (Table 1), the recorded snow depth underestimates the snow depth expected in the high-elevation mountain ridges and flanks, which have been affected by snow avalanches. At high elevation and for the last decade, the snow depth was investigated using three additional snow gauges in ski areas and ranged from 100 cm to as much as 400-550 cm (e.g., 2014) with average cover of 150-318 cm (e.g., Campo Imperatore, Prati di Tivo, and Gran Sasso areas, stars in Figure 3B; nonofficial data from https://www.skiinfo.it). This means that snow  Figure 2B). cover >200-300 cm deep is typical on high-elevation mountain slopes (e.g., >1,500-2,000 m) and can induce avalanches development especially after strong snowstorms or rapid temperature increase. The yearly number of avalanches ( Figure 3C) has ranged from 0 to 70 in the last decades, and more than 40 avalanches were recorded only in the 8 years pre-2017 (1969, 1984, 1991, 1999, 2003, 2005, 2009, and 2015) with a poor direct correlation with the snow thickness. In this distribution, 2017 features a distinct anomaly. A total of 115 avalanches were reported, of which 105 occurred within a 3-day time span, on 18-20 January 2017. More specifically, 76 avalanches were reported on 18 January; this means that in a single day, more avalanches occurred than in any other entire year in the time series , indicating a large anomalous distribution. Moreover, the diagram shows that this anomaly is no correlation between this high number of avalanches and the snow depth both in the official station and in the high elevation ones (where snow depth above 2.5 m is common; Figure 3C).

MATERIALS AND METHODS
The study was carried out through an integrated multidisciplinary regional (Snow Avalanche vs. Earthquake Distribution in Central Italy) analysis combining: 1) regional meteorological and climatic analysis, 2) seismicity analysis, 3) avalanches distribution analysis, and 4) comparison between seismicity and avalanches in space and time (e.g., the cases of April 2009 and January 2017). Of the many avalanches, the Mt. Siella area was well documented, eyewitnessed, and constrained and was investigated in detail (The Mt. Siella Avalanches) combining: 1) morphometry of orography, 2) photogeological analysis through multitemporal and multiscale air photo interpretation, 3) geomorphological field mapping, and 4) numerical modeling and stability analysis in the detachment area. The multidisciplinary approach included different integrated investigations, from field-based to numerical (Miccadei et al., 2017;D'Alessandro et al., 2006).
The general meteorological and climatic analysis was based on the data from the weather stations provided by the Hydrographic Service of the Abruzzo region and specifically the snow depth data from 12 snow gauges in the Gran Sasso area ( Table 1)    Concerning the Mt. Siella area, the morphometric investigation was based on a digital elevation model (DEM) of the study area, derived from 1:5,000 topographic maps (DBTN of the Abruzzo Region from the Opendata service; Regione Abruzzo, 2017) and extracted through ArcGis software (ArcGis 10.6 Esri © , USA). The entire Rigopiano-Monte Siella slope was investigated using geomorphological field mapping at a 1:5,000 scale, focusing on the analysis and the comprehension of the surficial deposits and landform distributions.
The multitemporal-multiscale photogeological analysis performed on erial images allowed us to investigate the eastern Mt. Siella slope and the Rigopiano plain area, and to verify the evidence of geomorphological processes (with particular reference to pre-2017 avalanches). The survey was performed on eight series of erial images (erial photos and orthophotos, Table 3) in the 1945 to 2017 timespan, with scales ranging from 1:5,000 to 1:55,000. They were acquired from the Cartographic Portal and Opendata Portal of the Abruzzo Region (Abruzzo Region, 2017), the National Geoportal, from BING map services, or provided by the Geographic Military Institute of Italy (IGMI). The analysis was based on the criteria and methods of photogeology, defining the main photogeological parameters (i.e., tone, texture, shape, vegetation, and drainage) on each of the images analyzed and through the three-dimensional analysis by stereoscope on pairs of photographic images. For each of the images the identification and classification of features led to the erial photo interpretation. Through the deduction stage on the classified features, the photogeological analysis provided indications and evidence of geomorphological elements linked to instability. Specifically, avalanches and gravity-induced and water-related landforms, were investigated from the slope of Mt. Siella to the Rigopiano plain, as well as evidence of the anthropogenic changes.
The numerical analysis specifically focused on the detachment area of the Rigopiano avalanche ( Figure 4), which showed a possible connection with seismic shaking and was poorly constrained by the post-avalanche investigations with respect to the transit and accumulation area evident from field survey and erial photo analysis. The analysis included a combination of both 3D numerical modeling using FLAC3D software (Itasca, 2019) and stability analysis of the snow cover, which attempt to explain a different phase of the behavior of the snow layer. In the 3D modeling, we intended to verify the occurrence of important deformations in the snow cover with or without seismic shaking (investigating the factors controlling the avalanche trigger). With the desired stability index, we investigated the possible achievement of the ultimate conditions before the failure, considering gravity and seismic ground acceleration.
FLAC3D is a numerical code mostly used for landslide investigations and slope modeling. In this work, the software features were exploited for modeling the deformation occurring within a thick snow cover on a steep slope both under static and dynamic conditions to investigate the detachment area and the role of the seismic shaking in the avalanche trigger. The nivological data (Bocchiola et al., 2018;NGI, 2018) do not provide an accurate definition of the snow cover and profile in the detachment area. For this reason, a plausible overall cover of 3.0 m was assumed over the entire slope. The snow and bedrock volumes were discretized by hexahedral meshes (bricks) parallel to the topographic surface; the bricks were 4 × 4 m wide and varied in height from 3 m at the base of the model to 0.5 m near the topographic surface ( Figure 4A). The modeling area focused on the eastern Mt. Siella slope from the top to the mid-slope ( Figure 4B).
The static numerical analysis was performed to investigate the slope before the avalanche and the occurrence of seismic shaking. The dynamic analysis was performed for assessing the influence of seismic shaking on the trigger of the avalanche. For the analysis, the four Mw > 5 shocks occurring on January 18, 2017 were considered combined in direct sequence with each other ( Table 2). The accelerograms were acquired from the ITACA database for the Capestrano accelerometric station (Luzi et al., 2017), which is the closest to the study site (located ∼18 km south of the study area). From the input accelerograms (properly scaled according to Bindi et al. (2011)), the EW component was selected as representative of the most effective shaking on the slope, which is dipping from W to E. the EW component of the accelerogram was applied at the base of the model in the x direction ( Figure 4A). The stability analysis focused on the snow cover to control the fracture propagation in the slab avalanche detachment area (Heierli et al., 2008;Reuter and Schweizer, 2018). Fracture propagation is based on the traction tangential strain transfer in the snow slab from the shear plane to the detachment line, which produces the yield stress. The shear stress (τ D ) is derived as the snow weight component parallel to the slope acting on the shear plane. In slab avalanches, one or more weak layers are assumed between thick snow layers with cohesion and strength values higher than in the weak layer (Schweizer et al., 2003;Heierli et al., 2008;van Herwijnen and Heierli, 2009). The trigger of an avalanche must be investigated, defining the rupture within the weak layer through the stability index (S) as the ratio of the shear strength of the snow cover (τ R due to cohesion and friction) along the weak layer over the shear stress (τ D ) along the shear plane ( Figure 4C) (Supplementary Appendix B for details on the stability index calculation): When S > 1, the system is stable; when S 1, the system is in critical stability conditions and failure is incipient; when S < 1, the system is unstable and failure occurs. Failure occurs within a weak layer (W) between two snow layers (Sl and Su) with more cohesion or between a snow layer (Su) and the bedrock (B) if the lower snow layer (Sl) is not present (Schweizer et al., 2003). The weight of the snow layer above the weak layer and the sliding surface, for the calculation of S, depends on the snow cover (h) due to the snow fall over time and on the layering of the snow profile. Without detailed snow profiles in the detachment area, the analysis is parametric. When additional forces occur (i.e., additional snow, wind, seismic acceleration), the calculation is completed by additional factors. The wind component is not considered in this case, since there are no detailed data and no strong winds were reported for the avalanche's day. The seismic shaking induces inertial forces where Δτ s is the additional shear stress due to seismic shaking. The denominator can be calculated as (Pérez-Guillén et al., 2014): In the latter formula, the resulting shear stress is variable with time as a function of the seismic ground acceleration a(t) (i.e., the component parallel to the slope). The worst conditions within the weak layer and along the sliding surface occur at the maximum ground acceleration value (i.e., the peak ground acceleration, PGA). The calculation of τ R is based on the Abe (2004) diagram, considering the shear strength of the graupel snow expected to be present in the study area at the moment of the avalanche's occurrence.

SNOW AVALANCHE VS. EARTHQUAKE DISTRIBUTION IN CENTRAL ITALY
In recent decades, some combined events were revealed with moderate-to-strong earthquakes affecting mountain areas in winter/early spring with moderate-to-deep snow cover ( Figure 3C). From the historical dataset, the avalanches documented as possibly connected to seismic triggers are restricted to two main events: the April 2009 event in the Assergi area and the January 18, 2017 event in the entire Gran Sasso massif.

April 2009 in the Western Gran Sasso and Assergi Area
In April 2009, the Abruzzo region was struck by the second strongest earthquake in the Central Apennines. Specifically. On April 6, 2009, the L'Aquila area was hit by the main earthquake (Mw 6.1 at 2:32 UTC+1, Figure 2B). The main event was followed by four more (Mw > 5.0, Table 2). The highest part of the mountain area (i.e., >1,500-1,800 m a.s.l.) still had moderate snow cover, 0.5-1 m, preserved in mid-springtime with occasional avalanches of wet snow (Meteomont, 2017).
The April 2009 earthquake, due to the serious damage to the L'Aquila area, was poorly investigated in term of avalanche distribution. Most of the geomorphological effects were, instead, related to co-seismic fractures and scarps, seismically induced landslides, ground collapse, and local liquefaction (Boncio et al., 2010;Miccadei et al., 2013). However, some data were derived from erial photo interpretation and eyewitnesses, particularly in the Assergi-Pizzo Cefalone area. The avalanches occurred in connection with the main seismic shocks from April 6-9, 2009 and at a distance from 10 to 20 km from the main epicenters (Table 4) in the western Gran Sasso area. More specifically, the Corno Piccolo and Prati di Tivo areas were affected by avalanches on April 6, 2009; the Assergi area was affected by avalanches (mixed snow and debris) during the shocks of April 6-7, 2009 (evidence derived from orthophotos taken on April 8, 2009; Figure 5). In the internal Gran Sasso area, other events were documented but without precise timing ( Table 4).

January 2017 in the Entire Gran Sasso Area
In January 2017, the eastern part of the Central Apennines and the hilly piedmont was affected by intense snowfall and rainfall. This event was part of a major cold meteorological event that affected central-eastern to southern Europe from 5 to January 21, 2017 (Di Lena and Giuliano, 2017;Calista et al., 2019). An intense snowfall event affected the Abruzzo area starting on 6 and 7 January (temperatures −10 to −13°C) with a snow cover of ∼1 m. A thermal collapse occurred on 11 January with temperatures down to −20°C, even at intermediate elevations (1,000-1,500 m a.s.l.). After a few days, mainly on 16-18 January, another perturbation affected central Italy, and specifically the Abruzzo, with a temperature drop and a heavy snowstorm (Di Lena and Giuliano, 2017;Calista et al., 2019). This combination of events created a >3 m snow cover along the eastern Abruzzo Apennines ridges, characterized mostly by graupel snow. In the hilly piedmont, the snow cover (wet snow) was up to 1.5 m and in the coastal area, it was >0.5 m. Finally, a temperature increase and a further perturbation with rainfalls induced a progressive rapid snowmelt.
On January 18, 2017, a moderately intense seismic event occurred during the main snowstorm, with more than 400 events in 24 h and 11 main events (Mw > 4) from 10:25 to 21: 00 (UTC+1), with four strong shocks (Mw > 5) from 10:25 to 14: 33 ( Figure 2B; Table 2). In the Gran Sasso mountain area, these events induced a PGA up to 5.74 ms −2 (INGV, 2017). Specifically,   Several avalanches occurred in connection with the first main seismic event (Mw 5.1, 10:25) and most occurred in the period of the four main (Mw > 5) events; three more occurred after the last main (Mw > 5) event, including the Rigopiano avalanche, which occurred 135 min after the last Mw > 5 event and 32 min after an Mw 4.3 event ( Figure 6B). Finally, for the documented avalanches, the analysis of the distance from the triggering earthquakes produced values of 9.4 to 50.2 km, which are mostly consistent with the seismically induced avalanches known in the literature (Podolskiy et al., 2010a;Podolskiy, 2010). This distribution again confirms the anomaly in the distribution of avalanches on January 18, 2017 and a clear connection between avalanches and seismic events, as well as delayed avalanches with respect to the main seismic events.

THE MT. SIELLA AVALANCHES
The snowfall event occurred in the area (see January 2017 in the Entire Gran Sasso Area) created, in the Mt. Siella area, a thick snow cover, estimated up to >3 m thick (252 cm were documented at the lowest elevation on 20 January, Meteomont, 2017). The cover was composed of a lower layer (∼0.7-1 m) of more compacted and stronger snow (due to the snowfall occurred on 6-7 January) and an upper layer of softer, not compacted snow (due to the 16-18 January snowstorm). In these conditions, on the NE slope of Mt. Siella-Vado di Siella-Mt. San Vito (along the eastern part of the Gran Sasso ridge, Figure 6), three avalanches occurred on the afternoon of January 18, 2017 as large slab avalanches. One occurred at 16:48 UTC+1 (from local chronicles or 16:42-16:43 as recently defined geophysically by Braun et al., 2020) in the Rigopiano valley (also known as Bruciata valley), 135 min after the last main seismic shock of the four on 18 January (Mw 5.0 at 14:33) and 32 min after a moderate seismic event (Mw 4.3 at 16:16). It was one of the most significant, most damaging, and best constrained among those that occurred during the January 2017 event. It was confined within the Rigopiano valley and spread out along the Rigopiano plain, whereas the two additional avalanches on the NE slope of Mt. Siella and in the Cupa Valley affected the slope without reaching the plain (Figure 8).

Orography and Hydrography
The Mt. Siella-Vado di Siella-Mt. San Vito ridge is roughly NW-SE with an N-S segment in the Rigopiano valley area. From The hydrographic network and the catchment trend from WSW-ENE to SW-NE, with an overall parallel drainage network pattern consistent with the slope dip.

Photogeological Analysis
The analysis allowed us to outline evidence of small avalanches and/or combined debris flows before the 2017 avalanches in some images (1945, 1989, 1991, and 2017; Figure 9). The other images only showed evidence of tree canopy regrowth in previous avalanche areas (1954, 1982, 2007, and 2013).
In 1945, the photogeological interpretation provided evidence of geomorphological processes connected to small avalanches and/or debris flows of limited extent (area ∼8,800 m 2 , length 460 m) and localized in the middle part of the Rigopiano valley ( Figure 9A). The findings are not recent and can be likely attributed to processes occurred in the decade before the image was captured. In the lower part of the valley, at the junction with the Rigopiano plain, minor alluvial fans can be outlined, with evidence of minor debris accumulation, well entrenched in the valley and into an ancient inactive alluvial fan. From 1954 to 1982, the image analysis outlined progressive tree canopy growth in the Rigopiano valley and over the alluvial fans. In 1989, the image showed minor evidence of debris flows and/or small snow avalanches in the upper (area ∼18,200 m 2 , length 480 m) and middle (area ∼6,400 m 2 , length 220 m) parts of the Rigopiano valley. According to the freshness of the evidence, these processes were attributed to the year before the image (e.g., 1988-1989, and certainly not before 1982).
This evidence was better documented in high-resolution images from 1991 ( Figure 9B). These images confirmed the minor evidence of debris flows and/or small avalanches in the upper and middle parts of the Rigopiano valley. They outline two more areas affected by avalanches and mixed snow-debris avalanches in minor valleys (e.g., Cupa Valley) to the north of the Rigopiano valley (NE slope of Mt. Siella). These elements were not documented in the 1989 image and should have occurred between 1989 and 1991. Specifically, the Cupa Valley was  Figure 9D), which occurred in the Rigopiano valley and the minor northern valleys (e.g., Cupa Valley). The Rigopiano avalanche is large (area ∼191,500 m 2 , length 2,200 m), with a detachment area in the eastern slope of Mt. Siella, a transit area within the Rigopiano valley, and the accumulation area in the Rigopiano plain. This size is largely anomalous with respect to the avalanches documented from 1945 to before 2017.

Geomorphological Features
The bedrock lithologies of the Rigopiano valley consist of massive and layered Mesozoic-Cenozoic limestone and dolomites ( Figure 10). The Rigopiano plain features a Neogene conglomerate and sand-clay bedrock. These lithologies were deformed by folds, thrusts, and normal faults, and are highly jointed. The bedrock is covered by thin surficial deposits consisting of slope debris and colluvial deposits at the base of the slopes and within the valleys and the plain. At the Mt. Siella slope-Rigopiano plain junction, patches of ancient alluvial fan deposits (cemented breccias and consistent gravel-debris) are present.
The main landforms (nival, gravity-induced, and waterrelated) are summarized in the geomorphological scheme of the Rigopiano valley (Figures 10, 11A), which combines the results of the field mapping and the photogeological analysis. The upper slope of Mt. Siella-Vado di Siella was affected by erosion incisions and grooves with transport and accumulation of slope debris ( Figure 11B). In the upper part of the Rigopiano valley, a blanket of slope deposits was observed. The deposits feature a low thickness (not more than a few meters in alignment with the incisions). At the mouth of the Rigopiano valley in the plain area, an irregular convex shape was observed that consists of highly heterogeneous debris deposits (with centimetric to metric blocks) firmly cemented, indicating an ancient form and deposit. It is reincized from the terminal part of the Rigopiano valley with a scarp up to 5-8 m high. This is assumed to be the remnant of an ancient debris fan, re-incized and no longer active in the current geomorphological conditions, possibly referable to the Late Pleistocene. Minor debris flows and alluvial fan are also present in the junction between Mt. Siella slope and Rigopiano plain (as well as on the SW slope of Mt. Siella).
The nival landforms consist of large nivation hollows and avalanches areas (pre-2017 and 2017 avalanches) along the Mt. Siella-Vado di Siella ridge. The SW slope is mostly affected by avalanche paths and accumulation areas. More specifically, the NE side of the slope is characterized by the landforms related to the avalanches that occurred on January 18, 2017 ( Figures 11A-C). The detachment developed in the upper part of the slope with a 30°-40°gradient, consistent with the typical slope distribution of avalanche detachment (van Herwijnen and Heierli, 2009). The transit area develops from the upper to the lower part of the slope, confined in the Rigopiano valley and minor valleys. The minor avalanches (e.g., Cupa valley) feature transit areas within poorly incised valleys; the accumulation area is confined within the Mt. Siella slope and did not affect the underlying plain. The major Rigopiano avalanche features a sinuous plan shape confined within the valley; the section profile is concave (deeply incised in the valley); the transit channel is asymmetrical, particularly in the valley turns, where it is upraized on the outer side of the valley ( Figure 11B). The accumulation area develops from the outlet of the Rigopiano valley and spread out over the entire plain ( Figure 11A).
The upper part of the slope was investigated in detail to outline the evidence of the detachment area. It was defined through tree cover removal ( Figures 11B,C, 12). No evidence of debris FIGURE 10 | Geomorphological scheme of the Rigopiano area. The scheme outlines the scattered superficial deposits covering the east slope of Mt. Siella-Vado di Siella and the mostly exposed bedrock. The main nival (pink), gravity-induced (red), and water-related landforms (green) are mapped.
Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 599611 16 removal was found in the detachment area, suggesting that the snow slide involved the snowpack and poorly affected the slope debris deposits. However, the field survey allowed us to map new evidence of co-seismic deformation. Some scarps and fissures in the debris cover were documented. These features occur along the Rigopiano valley at 1,535 m a.s.l. in a roughly NNW-SSE orientation ( Figures 11D,E) and the NE slope of Mt. Siella at 1,465-1,510 m a.s.l. in NW-SE orientation ( Figures 11D,E). According to their distribution and pattern, these elements were interpreted as co-seismic fractures connected to the seismic shaking of the slope. This further confirmed the connection between this avalanche and the seismic trigger. For this reason, we investigated the detachment area of the Rigopiano avalanche using 3D numerical modeling and stability analysis. Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 599611

Numerical Modeling
The numerical modeling focused on the investigation of the snow deformations, under static and dynamic conditions using FLAC3D software (Itasca, 2019), in the Rigopiano avalanche detachment area from the Mt. Siella slope to the upper Rigopiano valley (Figures 4, 12). Accurate dynamic modeling was conducted to verify the connection with the trigger of the snow avalanches and the mapped co-seismic fractures/scarps. The static analyses outlined a critical area in the snow cover, whereas the dynamic analyses highlighted the neoformation of permanent deformations inside the snowpack accumulated during the main seismic events. Combining the obtained result with a specific stability analysis on the snow cover, with a congruent calculation, we determined the overall conditions of the snow on the slope.

Numerical Modeling of Snow (FLAC3D)
The first analysis was carried out under static conditions ( Figure 13A), which outlined that the overall slope is generally stable. Minor deformations were highlighted in the upper part of the slope (steepest part, with slope >35°). According to the modeling, these areas did not reach rupture conditions during the static analysis, as evidenced by the maximum unbalanced force graph ( Figure 13A). A second analysis was run under dynamic conditions. The analysis used as input an accelerogram derived from the inputs of the four Mw > 5.0 earthquakes on January 18, 2017 (Table 2, Figure 13B), combining the most significant part of the four main shocks (red boxes in Figure 13B, combined in Figure 13C). The accelerograms from the Capestrano station were scaled according to Bindi et al. (2011) obtaining a maximum acceleration of 0.17 m s −2 . This input configuration was assumed to be representative of the effects of the entire seismic sequence on the snow stability. The resulting accelerogram was applied at the base of the model along the x-axis. The simulation was run for a total of 65 s of input. The dynamic analysis ( Figure 13D) highlighted two critical areas of deformation within the modeled area. These critical areas were reached after 11 s of modeling run (i.e., at the end of the Mw 5.1 input; Figure 13D) and increased in terms of displacement amount for the entire input sequence. The first area was located in the upper part of the Mt. Siella slope, from 1925 to 1830 m a.s.l. (zone A in Figure 13D) in the steepest areas (>35°). It largely corresponds to the area affected by low deformation in the static simulation ( Figure 13A). The second area was located in the lower part of the Mt. Siella slope and the upper part of the Rigopiano valley from 1,540 to 1,430 m a.s.l. (zone B in Figure 13D). This second area showed the largest displacements and corresponded to the upper and narrow part of the Rigopiano valley, where the co-seismic features were also documented. The areas in a plastic state ( Figure 13E) were even larger than those highlighted by the displacement contour ( Figure 13D). Not considering the possible side effects, the critical areas highlighted here are larger but correspond to those shown by the displacement diagram of the upper part of the Mt. Siella slope and the upper part of the Rigopiano valley.

Stability Analysis
For the stability analysis, without a detailed snow profile, a lower consistent layer (0.5-1.0 m thick) was assumed due to a first snowfall event around 6 and 7 January, ∼10 days before the avalanches. Above this layer, an upper less-consistent snow layer, related to the snowfall on 16-18 January 0-3 days before the avalanches, is defined with four different possible depths (0.5, FIGURE 12 | Geomorphological scheme of the Rigopiano avalanche detachment area (zoom of Figure 10). For the legend, see Figure 10; the red crosses mark the co-seismic fractures and scarps ( Figures 11D-G).
Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 599611 Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 599611 1.0, 1.5, and 2.0 m). The weak layer and initial sliding surface are assumed at the boundary between the upper and lower layers or at the soil snow boundary. In the analysis, a τ R ranging from 1.5 to 1.65 kPa was used, as derived for a snow density above 200 kg/m 3 according to Abe (2004). The stability index (S) was calculated for the four main seismic shocks on January 18, 2017 (

DISCUSSION
We discussed the regional distribution of the snow avalanches vs. earthquakes and compared the findings with the detailed features of the Mt. Siella avalanches to outline the role of seismicity in the triggering of the avalanches. Finally, the delayed trigger of the Mt.
Siella avalanches was discussed and compared with other documented cases.

Seismically Induced Snow Avalanches in Central Italy (Abruzzo Region)
The analysis of the snow avalanches vs. earthquakes distribution in central Italy (Abruzzo) documented two main events affected by seismically induced avalanches: the April 2009 event in the Assergi-Pizzo Cefalone area (western Gran Sasso, L'Aquila) and the January 2017 one (entire Gran Sasso area). On 6-9 April 2009, 24 snow avalanches in the western Gran Sasso out of 49 that occurred in the entire 2009 (in the Abruzzo mountain area) were connected, although without a precise timing, to the L'Aquila earthquake, the second strongest earthquake in the last decade in central Italy (five main shocks, Mw 6.1-5.4, from 6 to 9 April). Although the snow was present only at high elevation and the inventory of the avalanches was not systematic (heavy earthquake damages had just occurred), some events were documented through the analysis of post-earthquake erial photos and personal communications; for this reason, the number could be even underestimated. The 24 documented avalanches  (i.e., mixed snow-debris avalanches) occurred from high elevation slopes (∼1700-2,150 m a.s.l.), with an estimated snow depth of 50-100 cm, through steep slopes and channels, down to gentle slope bases (elevation ∼1,000-1,450 m a.s.l.). The distribution in terms of avalanches number vs. snow thickness is moderately anomalous, even if compared with the yearly distribution (in the lower part of the graph, green dot, Figure 14A), indicating a low snow depth for the number of avalanches that had occurred. This is confirmed by the connections of avalanches with seismic shocks. The distance of the avalanches to the epicenters was ∼10-20 km for Mw 5.4-6.1, entirely falling within the range of the distance vs. magnitude of triggering earthquakes documented at a global level (Podolskiy et al., 2010a;Podolskiy, 2010).
In January 2017, an anomalous number of avalanches occurred in the entire Gran Sasso area. A total of 76 were documented on 18 January, 14 on 19 January, and 15 on 20 January. Only 10 more were recorded throughout the rest of the year. The avalanches occurred from slopes at elevation ∼1,200-2,400 m a.s.l., with an estimated snow depth of >250-300 cm, proceeding from steep slopes, channels, and deep valleys, down to gentle slope bases or valley floors (elevation ∼700-1,450 m a.s.l.). In this case, the distribution in terms of avalanches number vs. snow thickness was largely anomalous, considering both the average snow depth recorded by the snow gauges and the maximum snow depth recorded at high-elevation ski areas. A total of 115 avalanches (105 of which occurred on 18-20 January) were never recorded in the historical 1957-2017 dataset and this value is two/three times the maximum previous yearly values. The values are on the extreme right of the avalanches number vs. snow thickness graph considering both average and high-elevation snow depth (red dots, Figure 14A), indicating high snow depth but a very high number of avalanches. The avalanches were well documented (Supplementary Appendix C), confirming the definite time connection to the four main seismic shocks (Mw 5.0-5.5) on 18 January, with synchronous or delayed release ( Figure 6B). During this event, the distance of the avalanches to the shocks' epicenter was ∼10-50 km for Mw 5.0-5.5 shocks and the distribution again falls within the ranges documented at the global level (Podolskiy et al., 2010a;Podolskiy, 2010).
Finally, the comparison between the avalanches documented on April 6-9, 2009 and January 18-20, 2017 and those recorded from 1957 to 2017 revealed how the 2009 and 2017 events were anomalous within the decadal distribution ( Figure 14A, Table 6, see Figure 3B for the yearly distribution). The general comparison between the yearly number of avalanches and the maximum recorded snow cover depth showed a poor connection of avalanche numbers with snow cover (scattered distribution in Figure 14). This could be due to several factors, mainly including the number of snow gauges and their low elevation, which may underestimate the snow cover on high-elevation slopes affected by avalanches ( Figures 3A,B). This is also confirmed by the high elevation stations where snow depth above 2,5 m are common and only in 2017 a large number of avalanches occurred even if the snow thickness was lower than previous years ( Figure 3C).
Moreover, on April 6-9, 2009, 24 events occurred in 4 days (2.56% of 938 avalanches recorded in 60 years). On 18-20 January 2017, 105 avalanches occurred in 3 days (11.19% of 938 recorded in 60 years), and the affected area was as high as ∼19% of the total area documented in the 60-years dataset (Table 6; Figure 6). In the April 2009 and definitely January 2017 events, the highly anomalous number of avalanches (red dots in Figure 14A) and related area ( Table 6) are well evidenced, strongly indicating an anomalous trigger related to seismic shaking.

The 2017 Mt. Siella Avalanches
At a detailed scale for the specific and well-constrained case of the eastern Mt. Siella and Rigopiano avalanches, the combination and integration of the field geomorphological investigations, remote photogeological analysis, and 3D numerical modeling allowed us to contribute 1) to the definition of the geomorphological features of the area, 2) to the definition of pre-existing evidence of avalanches, and 3) to the assessment of the geomorphological, nivological, and seismical factors contributing to the avalanches' triggering.
The eastern Mt. Siella slope is characterized by a combination of landforms: mainly nival and gravity-induced in the upper part and gravity-induced and water-related in the lower part (Figures 10, 12). The slope was affected before 2017 by snow  (minor avalanches) and gravity-induced processes from Mt. Siella on the upper part of the slope and the Rigopiano valley; the lower part of the slope and valleys were only affected by gravity and water-related processes, as confirmed by the air photo interpretation and photogeological analysis (Figure 9). In January 2017, the entire valley was affected by the processes related to the Rigopiano avalanches, to which gravity-induced (upper part of the valley) and water-related processes (lower part of the valley) were largely overprinted. The comparison of photogeological and geomorphological features outlined the connection between the avalanche and the seismic shaking through a combined geomorphological cross-section ( Figure 14B). The anomaly of the Rigopiano avalanche was also confirmed. The comparison of the evidence documented in images from 1945, 1989-1991) showed that the single Rigopiano valley 2017 avalanche was largely anomalous with respect to the previous events that occurred in the same valley. Specifically, the 2017 event area was 10.5 to 50 times bigger and 4.8 to 10.7 times longer than the events documented in the 1945 and 1989-1991 images ( Table 7). Moreover, this avalanche is anomalous also within the general area frequency distribution of all the avalanches documented in the 1957-2017 dataset, corresponding to the 94th percentile of the 938 avalanches dataset ( Figure 14C). The avalanche modeling ( Figure 14B) and the stability analysis showed that the detachment (eastern Mt. Siella slope and upper part of the Rigopiano valley) was triggered by seismic shaking combined with the considerable snow depth. The critical areas for displacement shown by the FLAC3D numerical modeling ( Figure 13D) correspond to the detachment area (Zone A). This is also consistent with the distribution of the new surface evidence of coseismic fractures and scarps documented in the geomorphological mapping. More specifically, the location of the co-seismic scarp observed in the field correspond to the lower part of the detachment area and the second critical area evidenced by the modeling (Zone B, with stability conditions worse than Zone A; Figure 14B). This suggests that the avalanche was triggered in the middle part of the slope (upper part of the Rigopiano valley), the snowpack was destabilized from middle part to the upper section of the slope of Mt. Siella and thus was able to release an exceptionally large avalanche.
Finally, the analysis of the Mt. Siella and Rigopiano avalanches outlined that the January 18, 2017 event was largely anomalous in terms of avalanche distribution and magnitude, even on a detailed scale, as shown for the regional distribution (Seismically Iinduced Ssnow Aavalanches in Ccentral Italy (Abruzzo Rregion)). This provides further confirmation that the 2017 event resulted from an anomalous combination of the heavy snowstorm (with high snow depth on the mid-high elevation slopes) and the seismic shaking induced by the January 18, 2017 earthquakes. This was confirmed even if the avalanches were delayed with respect to the closest main seismic shock.

Delay Between Earthquakes and Avalanche Trigger
The delay between earthquakes and avalanche trigger was explicitly investigated for the 2017 event, in which the time distribution of the avalanches was, at least in some cases, well documented. For the 2009 case, besides the avalanches being snow-debris mixed, their precise time distribution was poorly constrained. On January 18, 2017, the avalanche time distribution with respect to the main earthquakes showed that several avalanches occurred right after the first main (Mw > 5.0) seismic event (i.e., at 10: 25 UTC+1); other events occurred in the timespan during the main shocks, whereas other events were triggered with a delay from a few minutes to tens or hundreds of minutes after the main seismic events ( Figure 6B). For the latter, the Rigopiano avalanche is one of the temporally best-constrained avalanches on 18 January. Even if the snow slope was in poor stability condition, as suggested by the 3D numerical modeling and by the stability analysis, the avalanche was not triggered at the specific time of any of the main (Mw > 5.0) seismic events. It was delayed after the events, triggering 135 min after the last major seismic shock (Mw 5.0 at 14:33, distance ∼40 km) and 32 min after a moderate seismic shock (Mw 4.3 at 16:16, distance ∼44 km). The delays occurred in this event are largely consistent with those documented before in Western Himalaya, ranging from a few minutes to several hours (specifically, 7 min on January 16, 2000 and 36 min on February 12, 2000 in Drass area; 49 min on January 26, 1999 in Chadan; 82 min on January 21, 2000 in Jawala; 350 min on March 4, 1996 in Lagongama; several hours on February 1, 2000 and February 14, 1999 again in Drass area; Singh andGanju, 2002 andPuzrin et al., 2019). However, in Western Himalaya cases, the magnitude of the seismic events was moderate (Mw 2.7-4.1, comparable or lower than in the 2017 event of Central Italy) and the distances were larger than in our case studies, so as for the cases with the largest distance (300-500 km) the direct earthquakeavalanche connection appears unlikely if compared with the analysis of Podolskiy et al. (2010a).
The occurrence and mechanism of the delayed seismically induced avalanche are consistent with what has already been documented in literature in the Himalayan area (Singh and Ganju, 2002) and for the Rigopiano avalanche (Puzrin et al., 2019). In particular, the new analytical model produced dynamics coherent with the documented case of Rigopiano. The proposed model, which is particularly complex and supported by a well-defined mathematical demonstration, calculates the mechanical development of the delayed release of earthquake-induced avalanches on the basis of the following points (Puzrin et al., 2019): The earthquake induces an early shear fracture of length l 0 within an existing weak layer. Due to this fracture, the forces in the sliding slab assume a new distribution, and a Kelvin-Voigt constitutive model in the creeping mass and in the intact weak layer external of the shear fracture will activate. This Kelvin-Voigt creep decays, producing a new defined growth-triggering length l g . If l 0 > l g , the same decaying creep will produce sufficient dislocation in the intact weak layer to activate its viscoplastic behavior. At time t 0 , the shear fracture increases with the rate of its development, which depends on the visco-plastic rate-hardening of the weak layer. If the shear strength is rate-dependent (if not, the spread of the shear fracture would produce a catastrophic movement already at time t 0 ), the rate-hardening experiences a gradually accelerating evolution of the shear fracture until a serious strain-rate is reached in the weak layer. Once this critical strain-rate is reached, the shear strength in the weak layer depends on the visco-plastic rate-softening instead of the rate-hardening. At t f t 0 + Δt cr , the shear fracture reaches its critical extent l cr , which reveals the start of its catastrophic propagation mainly to avalanche release.
The results for Rigopiano are depicted in Figure 1C. A large range of realistic parameters were used for three potentially earthquake-triggered delayed avalanches. The parameter k that was introduced in Figure 1C represents the stress ratio (tanα/tan ϕ)-1, as derived by Gaume et al. (2013Gaume et al. ( , 2018.
The model shows that delayed avalanches caused by earthquakes are rare but possible, and could lead to significant damage, especially on wide and milder slopes. Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 599611 24