Engineering-Geological Analysis of a Subaerial Landslide in Taan Fiord, Alaska

On 17 October 2015, a large-scale subaerial landslide occurred in Taan Fiord, Alaska, which released about 50 Mm3 of rock. This entered the water body and triggered a tsunami with a runup of up to 193 m. This paper aims to simulate the possible formation of a weak layer in this mountainous slope until collapse, and to analyze the possible triggering factors of this landslide event from a geotechnical engineering perspective so that a deeper understanding of this large landslide event can be gained. We analyzed different remote-sensing datasets to characterize the evolution of the coastal landslide process. Based on the acquired remote-sensing data, Digital Elevation Models were derived, on which we employed a 2D limit equilibrium method in this study to calculate the safety factor and compare the location of the associated sliding surface with the most probable actual location at which this landslide occurred. The calculation results reflect the development process of this slope collapse. In this case study, past earthquakes, rainfall before this landslide event, and glacial melting at the toe may have influenced the stability of this slope. The glacial retreat is likely to be the most significant direct triggering factor for this slope failure. This research work illustrates the applicability of multi-temporal remote sensing data of slope morphology to constrain preliminary slope stability analyses, aiming to investigate large-scale landslide processes. This interdisciplinary approach confirms the effectiveness of the combination of aerial data acquisition and traditional slope stability analyses. This case study also demonstrates the significance of a climate change for landslide hazard assessment, and that the interaction of natural hazards in terms of multi-hazards cannot be ignored.


Landslide in the Context of Natural Multi-Hazard and Recent Climate Change
Although they are a common natural hazard, the understanding of landslide initiation and failure processes under complex conditions still represents a challenge for potential hazard assessment. Landslides mobilize variable volumes of rock and debris, commonly varying from hundreds to millions cubic meters, which can have severe impacts on the surrounding environment. Landslides can lead to cascading natural hazards; for example, a landslide could initiate a tsunami when the failing slope is situated next to a waterbody [1][2][3]. Landslides, and other natural disasters triggered by landslides, have already caused a large number of casualties and property loss worldwide [4][5][6]. In many cases, the subaerial slopes involved in landslide processes are located in complex environments, where many factors affect slope stability. To properly understand a landslide destabilization and collapse process, geotechnical analyses are indispensable, in addition to related geological and geomorphological analyses.
Landslides are the culmination of slope instability processes, who evolve over time due to a variety of factors, including local hydrological conditions, regional geological In this study, we investigate the 17 October 2015 landslide in Taan Fiord, Alaska, from a geotechnical perspective by applying the limit equilibrium method. Previous studies on this landslide and related tsunamis in this region were mainly conducted from geological, geomorphological and geophysical perspectives [2,3,[34][35][36]. In the present study, the multi-temporal aerial radar data were used to reconstruct the slope morphology over time and to constrain the subsequent slope stability analyses via the limit equilibrium method. The applicability of remote-sensing data for landslide investigation has been confirmed in various studies [37][38][39][40][41][42][43]. The aim of this study is to gain a better understanding of the landslide process, and of which factors or external triggers may have induced the slope failure. Our findings help to clarify possible triggering mechanisms of this landslide by comparing our results and the observed sliding surface in an actual situation.

The 2015 Taan Fiord Tsunamigenic Landslide Event
The large-scale landslide occurred near Icy Bay, Alaska ( Figure 1). Taan Fiord, as an arm of Icy Bay, is located in the St. Elias Mountains of coastal southeastern Alaska [3]. In this region, the high mountainous terrain, active faults, frequent earthquakes [44,45] and accelerated glacier retreat all increase the potential for large-scale landslides.
On 17 October 2015, after a period of heavy rainfall, a massive landslide occurred at the terminus of Tyndall Glacier and induced a tsunami with a runup of up to 193 m [2], which was the fourth-highest tsunami ever recorded in the world [46]. The landslide released approximately 50 Mm 3 [47] of rock into the Taan Fiord. The length and width of this landslide are approximately 1630 m and 915 m, respectively. The slide area is 1.234 km 2 . The maximum thickness of the landslide body is 93 m (Figure 1), which means that the maximum depth of the sliding surface is around 100 m.
No one directly observed this landslide event. Although field surveys have been conducted and related research works have been performed [1][2][3]35,36,[47][48][49], the triggering mechanism of this landslide is still not elucidated. Seismic waves from an earthquake about 500 km away arrived shortly before this landslide, and a few seconds of mild shaking occurred [2], which could also be a possible trigger of this event. The precipitation in this area before this event (a month of above-average rain) was also higher than usual [2], which was also not beneficial to the stability of this slope. Glacier retreat and thinning due to climate change weakening the support of the toe of the slope was also observed in this area.
Previous studies related to this landslide event were conducted in 2015, including the geological and geomorphological conditions of this area. For example, Meigs et al. [34,50] investigated the ultra-landscape response and sediment yield following glacier retreat in this area, and described the scarps within the slope. After this landslide event, more research was conducted on landslide and landslide-generated tsunami phenomena in this area. For example, Higman et al. [2] conducted the field investigations and their field observations provided a benchmark for the landslide and tsunami hazards modelling. Haeussler et al. [35] focused on the submarine observations of bathymetry, submarine seismic profiling and surface change, and calculated the total volume of this landslide, including the volume that entered the fiord. Dufresne et al. [3] described the landslide mass and the landslide deposit using the debris on land, supra-glacial debris and submarine deposition. George et al. [1] presented a new methodology with the D-Claw model to compute the tsunami generation by the subaerial landslide. Gualtieri et al. [48] carried out a broadband analysis of the seismic signals by this landslide event, deduced the mass, trajectory and characteristics of the landslide dynamics, and also estimated the landslide volume to be about 55 million cubic meters. Williams et al. [49] analyzed the geomorphic and sediment accumulation history of Taan Fiord for the Tyndall Glacier retreat for over 50 years, and carried out a comparison of glacial and paraglacial denudation responses to rapid glacier retreat. Bloom et al. [36] investigated the characteristic modifications of the fan deltas in Taan Fiord, which were associated with the tsunami. Franco et al. [47] simulated a tsunami event in Taan Fiord with Flow-3D to verify the validity of the applied  numerical models in reproducing the wave dynamics. Although many scientists have  analyzed and studied the tsunami generated after the landslide event and its effects, only  rudimentary related geotechnical investigations are performed on the 2015 Taan Fiord  landslide (called the Taan landslide later in this paper).

Geological and Structural Setting
The slope is located on the St. Elias Mountains, which is the product of the collisional history of the Yakutat terrane and the North American continent, the world's highest coastal mountain range [3,51]. The topographical relief of this area has changed due to the collision between microplates and the influence of the external natural environment, e.g., the impact of the 2015 landslide event. The data collected with various remote-sensing methods were used to derive several Digital Elevation Models (DEM), from which essential information such as the topography and altitudes of different places in this area were collected [47]. Table 1 summarizes the different remote-sensing campaigns and the derived DEMs, which are used for analysis from 2000 to 2016. The DEM data for the Taan Fiord area are provided by the Elevation Portal of Alaska-Discrete Global Grid System (DGGS) and Haeussler et al. [35] for the years 2000-2014 and 2016, respectively. For the Taan Fiord in October 2015, just before the landslide occurred, a reconstructed, pre-collapse model based on the Interferometric synthetic aperture radar (InSAR) DEM of 2012 by George et al. [1] was presented. Since 2017, a new dataset (Arctic DEM AK V.2-2014, see Table 1) is available, and was used to update the landslide volume estimation just before the final collapse. The topographic surface after this catastrophic landslide was created by Haeussler et al. [35]. In comparison to the DEMs of 2012 and 2016, the Arctic DEM-2014 shows significant elevation and horizontal offsets [47,52]. To obviate this disagreement, a re-projection and a co-registration were applied for the Artic DEM-2014. The re-projection process is conducted by the "Coordinate transformation (Grid)", a tool in SAGA GIS [53] (Tool Libraries, Projection, Proj.4). The result provides a better re-adaption of the chosen digital model on the referred EPSG (in this case, 32607, WGS84/UTM zone 7N) and the related DEM of 2016 (from Haeussler et al. [35]). Initially, a manual co-registration with a mean offset of 25 m in a vertical direction was used. The result gave a good correspondence between the Arctic DEM-2014 and the one of 2016, with an accuracy of ±10 m in the vertical direction [47]. This accuracy was shown to be sufficient to identify very large mass movements.
However, some horizontal discrepancies were still present after manual correction. Thus, a co-registration process was performed for the DEMs of 2014 and 2016 to the DEM of 2012 with the method by Nuth and Kääb [54]. For this "demcoreg", which is a collection of python and shell scripts, was utilized [55]. For the co-registration process, some stable sections (e.g., sections that are non-glaciated, or without mass movements or rockfalls) were selected in this area. Those areas have been highlighted in purple in Figure 2a. For the results before and after the co-registration process, Figure 2b Table 2. It can be seen that the Root Mean Square Difference (RMSD) between the DEMs of 2014 and 2012 and between the DEMs of 2016 and 2014 reduced to 2.03 and 1.85, respectively, after the co-registration process. The statistics seem to be reasonable when considering high geomorphological dynamics and the high variability in the snowpack in this area over the year.  The geomorphological analysis is given using the open source software QGIS [56] and SAGA GIS. To be consistent with the analysis, some digital elevation models were downsampled to a resolution of 5 m with the method of bicubic spline interpolation that is available in SAGA GIS. Utilizing the raster calculator method implemented in QGIS, the difference in elevation for the landslide area from 2000 to 2016 can be obtained and is shown in Figure 3. In this period, especially after this catastrophic landslide, the elevation at different slope locations changed a lot. For the vertical displacement, a maximum negative value (mass loss) of 180 m and a positive value (material deposition) of 60 m are observed in this area. Before the year 2000, some facing scarps on the upper part of this slope were already present, as illustrated in the work of Meigs et al. [34]. Based on the aerial photographs and mapping [34], a mountain-scale landslide failure of the fiord wall was detected between 1983 and 1996. This documented landslide covered an area of approximately 0.4 km 2 , and a head scarp of about 200 meters appeared with a dip of around 60°. The resulting head scarp cross-cut the bedrock structure and struck subparallel to the fiord wall [34].
To analyse possible triggers of the 2015 Taan landslide, some studies refer to the glacial retreat [3,36,49], as well as to a large amount of precipitation before the occurrence of this landslide or frequent earthquakes in this area [2]. George et al. [1] also put forward that the mobilization of 2015 Taan landslide cannot be directly attributed to the glacial retreat, as the old slump was suggested to be a result of the instability caused by postglacial deglaciation and the consequent increase in valley wall height [34], while the height of the valley wall in this area did not significantly increase from 1996 to 2015.
Regarding the pre-evolution and formation process of the 2015 Taan landslide, Franco et al. [47] gave a four-stage interpretation of the landslide dynamics: (1) Before the year 2000, slope instability may be induced due to the fast ice loss and glacier retreat, which leads to the formation of internal antithetic discontinuities or local normal faults with a dip ranging between 60°and 80°against the slope; (2) Before the year 2012, the slope may have undergone some minor rotational displacements, such as a rigid body, rotating counter clockwise by about 20°to 30°; (3) After the year 2012, the orientation of these internal discontinuities gradually changed and shifted to a sub-vertical dip due to the rotation of the sliding mass. This led to a sudden significant vertical displacement during the period from 2012 to 2014; (4) In the years prior to the 2015 landslide event, the glacial retreat may lead to the final collapse due to its weakened buttressing effect. A schematic diagram that can be used to interpret the landslide dynamics prior to slope failure is shown in Figure 4. It is not known whether the rotation of the sliding mass was caused by deformation along the newly formed, internal discontinuities prior to 2012. The actual triggers that led to the 2015 landslide event are also unknown, as well as the causes of the newly formed discontinuities within the slope. Therefore, a detailed investigation of possible failure processes and triggers is presented in the study.
In the area around the location of the Taan landslide, there is a Chaix Hills Thrust fault along the east-west orientation. There are several known stratigraphic units for this fault, e.g., the Kulthieth Formation and the Yakataga Formation [3]. For the stratigraphy of Taan landslide site, the Kulthieth Formation dominates. This consists of interlayered sandstone, coals, mudstones, and conglomerates [58]. Due to the ongoing collision of the Yakutat microplate with North America, the surface compositions of the this slope are weakly lithified sedimentary rocks [59].
The geotechnical parameters (e.g., soil cohesion or friction angle) are unavailable for the slope, which provides a huge challenge for geotechnical analysis. Since the geological compositions of this area are dominated by sand-and mudstone, the friction angle should be correspondingly large, and their cohesion is closely related to the ratio of each component and the degree of weathering. As the surface of this slope is severely affected by weathering, its strength parameters should be close to those of weathered sandstone or mudstone. Based on their similar geological composition, the strength parameters of this slope could be estimated accordingly, by referring to the range of values taken from the other literature, i.e., the cohesive strength lies between a few tens and a few hundreds of kPa, and the friction angle is usually between 30°and 50°(and between 40°and 50°in many cases) [60][61][62][63][64][65][66][67][68][69][70][71][72][73][74]. The literature review of strength parameters with a similar geological composition can be seen in Table 3.

The Cross Sections
For this study, the elevation changes in this slope and the adjacent glacier are investigated by selecting representative cross-sections in the area (

Methods
The purpose of this case study is to investigate and verify the possible triggers and processes of this landslide using two-dimensional limit equilibrium analysis. To facilitate a comparison of the results, all calculation results in this paper were obtained with the Morgenstern-Price method.
The computational model is based on the slope profile of 2014 at the cross-section S-S', as shown in Figure 5. Since the geotechnical parameters and internal strength distribution of this slope are not available, a homogeneous material is assumed for the slope body (weak lithified stone). In the calculation model, as shown in Figure 6, the horizontal starting point of the cross-section S-S' in this study area is also the zero point of x-coordinate (distance is zero in Figure 5b), and the elevations are the y-axis coordinates. Since the altitude of the mountain does not vary much within a few hundred meters uphill from the starting point of the known cross-section, away from the Taan Fiord, which is at a height of approximately 900 m, the model is extended by 100 m along the horizontal and vertical directions, respectively, to reduce the artificial constraints.
The height of this slope is about 890 m, the horizontal length is about 1700 m, and the average slope is slightly less than 30°, with a steeper upper part and a gentler lower part of about 35°and 20°, respectively. The groundwater level in this slope is also unknown, which makes it difficult to analyze the effect of precipitation on slope stability. To account for the influence of groundwater on the results, different groundwater levels were adopted for seepage and no-seepage conditions. In Slide (a Rocscience software), the resulting groundwater table is obtained based on the steady-state groundwater calculation. Setting the height of the groundwater table at the left boundary (uphill-side, at the position with a horizontal coordinate of −100 m) of this model as the initial groundwater table, as shown in Figure 6, different steady-state tables were adopted in the calculation, marked by dashed lines of different colors, ranging from 800 m to 200 m. The range of values for the groundwater table in calculations is listed in Table 4. We emphasize that this is a wide range of possible groundwater tables, and which one leads to realistic results will be revealed later (see Section 4). In the calculation model, the parameter values of geotechnical materials are based on the literature review of cases with a similar geological composition. The density values of sandstones are usually between 2150 and 2650 kg/m 3 [76]. The values of 2350 kg/m 3 [2] and and 2700 kg/m 3 [48,77] are also adopted in the estimation of landslide mass. Based on rocks with a similar composition, age and burial history from the the Susitna and Cook Inlet basins, Alaska [78,79], the density of the landslide material is estimated as 2350 kg/m 3 . Combining the reported density values, we chose 2450 kg/m 3 . Thus, the unit weight for the geotechnical materials of the slope is taken as 24 kN/m 3 . The values of the strength parameters for the slope body will be discussed in combination with the analysis. The unit weight for the glacier is taken as 9 kN/m 3 , since it is composed entirely of ice. Based on the typical values for the Mohr-Coulomb parameters of ice [80,81], 4.5 MPa and 6°are adopted for the strength parameters of the glacier, i.e., the cohesion and the internal friction angle, respectively.

Preliminary Slope Stability Analysis
First, some preliminary calculations were performed to investigate the effect of precipitation on the position of possible sliding surfaces. The Mohr-Coulomb parameters were adopted for the landslide body in this section. Here, the combined, extreme values of the cohesion and friction angle of the involved materials were assumed to test their influence on the results. The cohesion is taken as 20 kPa or 100 kPa, while the friction angle is taken as 30°or 45°according to most of the strength parameters, which are defined by other authors for the same materials (Table 3). Therefore, there are four parameter combinations for the sliding mass ( Table 5). The four sets of parameters are denoted as p1, p2, p3 and p4. The results are shown in Figure 7, and the FoS values calculated for different groundwater levels are listed in Table 6.
In Figure 7, the red line represents the position of the assumed slip surface in the 2015 Taan landslide. The four subfigures show that the sliding surfaces appear only in the upper part of this slope, regardless of the high or low groundwater levels, when only the effect of precipitation on the slope stability is considered. This is not the most probable actual location of the sliding surface in the 2015 Taan landslide. This is also true when performing the calculation with a predefined discontinuity in only the upper part according to Figure 4. Moreover, since the calculated sliding surface remains well above the toe of this slope, the glacial retreat does not have a strong influence on the location of the failure surface.   The outcomes of the preliminary 2D limit equilibrium calculations revealed that a sliding surface close to the assumed shear zone in 2015 Taan landslide could not appear if only the influences of glacial retreat and different groundwater levels are considered. This suggests that the glacier retreat in the decades before 2015 may have had an impact on the stability of this slope during the formation of the internal discontinuities, but was perhaps not the direct, dominant factor in its generation throughout the whole slope body.
These seemingly unsuccessful preliminary attempts led us to think more deeply about the triggers of this landslide event and the evolution of highly likely discontinuities, which have a critical impact on slope stability. Based on preliminary calculations, various influencing factors and possible discontinuities are taken into account in the simulation and analysis of slope failure in the following sections.

The Formation of a Weak Basal Layer
As the slope is in the earthquake-prone zone, the influence of earthquakes on the stability of this slope cannot be ignored. Large-magnitude historical earthquakes led to active faults in this area, e.g., earthquakes with M8.2 in 1899 and seven earthquakes >M6 since 1958 [3]. Peak ground motions generated by nearby earthquakes are important in earthquake safety analysis for slopes, especially the peak ground acceleration. In general, the horizontal acceleration generated by earthquakes can have an impact on the stability of the slope, and even on the shape and location of the sliding surface when a landslide occurs, while the vertical acceleration is usually neglected [82] especially in the conventional pseudo-static slope stability assessment. From the information available on the United States Geological Survey (USGS) website, combined with the empirical formulae [83,84] for the horizontal peak ground acceleration (PGA) generated by earthquakes, several earthquakes near the landslide site were found to generate horizontal accelerations of 0.3 g or even larger at this region before the 2015 Taan landslide.
One new hypothesis for a possible interpretation of the historical formation process of the internal weak layers in the slope and the process of landslide occurrence is as follows. On the whole, the landslide formation could be broadly separated into two stages: (1) With earthquakes and precipitation as the possible main influencing factors, multiple destabilizations of this slope occurred and, accordingly, the sliding mass moved several times; however, there was no overall slope collapse. These motions caused shear deformations within the slope, which were localized in a relatively small shear band. The landslide body initially consisted of weathered rock, and the shear zone was most probably a granular material, which is produced by crushing the original rock mass. This crushing led to the disintegration of the material, resulting in a lower shear strength. Obviously, this decrease in shear strength was not strong enough to yield a final collapse; however, a zone of lower-strength material remained, which we call a weak layer. (2) Despite the presence of this weak layer, the slope did not fail (FoS > 1) due to the buttressing effect provided by the glacier at the toe of the slope. Before the 2015 Taan landslide, the glacier substantially retreated, causing a loss of support at the slope toe and triggering this large-scale landslide. To account for the assumed two stages of landslide formation, the below analyses follow these two stages: the formation of a weak layer, and the final slope collapse.
In the first stage of calculation, groundwater level and seismic intensity are investigated as influencing factors of the calculation results. As illustrated in above sections, different groundwater levels are adopted (Table 4). Similarly, different horizontal acceleration values for seismic loading (0.0 g, 0.1 g, 0.2 g, 0.3 g, 0.4 g and 0.5 g) were adopted in the calculations (Table 7). Table 7. The range of values for the horizontal seismic load coefficient in calculation.

The Influence of Groundwater Level
The two main investigated factors are the height of groundwater table and the seismic load. Their effects are correlated, and we begin with the influence of the groundwater table height for two exemplary seismic loads. The mechanical parameters of the basal weak layer and overlying sliding mass (e.g., unit weight, strength parameters) are listed in Table 8 as the input for the material properties in the caluclation model. As illustrated in the above section, the cohesion and friction angle of weathered rock with a mixture of sandstone and mudstone are usually between a few tens and hundreds of kPa, and 30°to 50°, respectively. As an example, the Mohr-Coulomb parameters for Set 1 were as follows: cohesion and internal friction angle were taken as 100 kPa and 45°, respectively. Since no significant differences were found in the results of the computations with parameter Set 1 and 2, only the results of parameter Set 1 are presented. Table 8. The material properties assumed in the limit equilibrium method (LEM) analyses. ince the location of the calculated sliding surface significantly differs from that assumed for the 2015 landslide when seismic loading is not considered, the results are compared for different groundwater levels under the horizontal seismic acceleration conditions of 0.3 g and 0.2 g, as shown in Figures 8 and 9. The corresponding FoS values are listed in Table 9. Other seismic accelerations are investigated in the following subsection.   For a seismic acceleration of 0.3 g, the most likely sliding surface occurs in the upper part of this slope when the groundwater table is very high (e.g., initial groundwater table is 700 m or 800 m) or very low (e.g., initial groundwater table less than or equal to 400 m). For cases with low groundwater levels (less than or equal to 400 m or without considering the groundwater), sliding surfaces appear at the same location, as they are not affected by the groundwater. Moreover, the sliding surfaces resulting from the very low groundwater table conditions are located slightly higher than those computed with very high groundwater levels. For an initial groundwater table at a moderate level, e.g., 600 m or 500 m, the sliding surface cuts through the entire slope from the top to the bottom, and is close to the position of the assumed Taan landslide slip surface (red solid line). It is also noted that the sliding surface for an initial groundwater table of 600 m is more close to the assumed slip surface than that of 500 m.
For the seismic acceleration of 0.2 g, the sliding surface locations resulting from high and low groundwater tables have similar characteristics to those of the 0.3 g case. However, for the medium groundwater table, the sliding surface location produced by the initial groundwater table of 600 m is slightly deeper than that of the 0.3 g condition, which is closer to the assumed slip surface. For an initial groundwater table less than or equal to 500 m, the slope did not fail (FoS > 1). The calculation results for the seismic acceleration condition of 0.3 g are, therefore, more realistic than that of 0.2 g.
Futher results of a large number of calculations, conducted by adopting different geotechnical parameters (cohesion: 50 kPa to 500 kPa, friction angle: 30°to 50°), show that the sliding zone under moderate groundwater level and seismic loading conditions is relatively close to the observed Taan landslide slip surface.

The Influence of Horizontal Seismic Acceleration
To investigate the influence of different seismic accelerations on the calculation results, two sets of geotechnical parameters and an initial groundwater table at 600 m are selected. Cohesions and friction angles in these two sets of Mohr-Coulomb parameters (Set 1 and Set 2 in Table 8) are 100 kPa and 45°, and 120 kPa and 50°, respectively. The resulting slip surfaces for these two sets under different seismic load conditions are shown in Figures 10 and 11. The corresponding FoS values for the two sets are listed in Table 10.   Since the upper part of this slope is relatively steep, the most likely sliding surface is located in the upper part of this slope when no earthquake is considered or when the earthquake acceleration is relatively small, and the FoS values are larger than 1 for these two parameter sets. When the seismic acceleration is no less than 0.2 g, the most likely sliding surface is located from the top to the toe of the slope. When the seismic acceleration is relatively large (≥0.3 g), the sliding surface is more close to the assumed Taan landslide slip surface than that with small seismic loading (≤0.2 g), which is also reflected in the previous analysis of the groundwater level effects.
By comparing the results for different parameter values, the potential sliding surface calculated under the seismic condition was found to be very close to the actual sliding surface location assumed for the Taan landslide when the Mohr-Coulomb parameters were within a certain range. For example, a well-fitted result could be obtained when the friction angle was around 45°(±5°), the cohesion was around 100 kPa (±20 kPa), the initial groundwater table was around 0.65 (±0.05) times the slope height and the horizontal seismic acceleration was at least 0.2 g (better for ≥0.3 g).
Since the obtained safety factor is only slightly less than 1 (usually larger than 0.8), this sliding only produces a relatively small displacement during the earthquake. The block mass above the shear zone slid for a certain distance, and the slope returned to a stable state after the earthquake, fitting the main assumption of any displacement estimation based on the Newmark method [85]. One can assume that several such movements occurred within a small distance in the decades before the 2015 Taan landslide. Due to the block mass movement above a localized shear zone, a weak layer developed within the slope, called a fault.

The Final Slope Collapse
In the 2D LEM analyses, the weak shear zone at the base of the landslide, formed during the first slope instability stage, was considered a discrete sliding surface. The sliding surface was modelled as a pre-defined discontinuity with the shear strength parameters of the weak layer. To calculate the second stage, the glacial retreat at the toe of this slope was focused on as a trigger of the Taan landslide. Since there was only one mild shaking, from a small earthquake far from this site, shortly before the Taan landslide in 2015, the seismic loading was not considered in the calculations in this phase. The two sets of Mohr-Coulomb parameters (Set 1 and Set 2; see Table 8) are still used for the analysis with two scenarios-with and without glaciers at the toe of this slope-as shown in Figure 12.
Based on the values taken in the previous subsection, the height of the groundwater table at the left boundary was taken as 600 m. The parameters of the weak layer were assumed to be different from the overlying sliding mass due to the shearing process that occurred in recent decades. To set the Mohr-Coulomb parameters at the weak layer position, the cohesion is usually small and the friction angle is reduced to some extent. For example, the strength of the material in the shear zones of deep seated landslides, measured by Strauhal et al. [86], is rather low.
Four sets of Mohr-Coulomb parameters (a1 and a2, b1 and b2) were chosen as the strength parameters of the weak layer, summarized in Table 8. These values represent the weakened material in the evolved shear zone. The results of the safety factors for these four cases are listed in Table 11.  The results clearly show that the safety factor obtained under the condition with the glacier at the toe was greater than 1, while the result under the condition without glacier was less than 1. This indicates that the glacier at the toe had a buttressing effect on the slope, and, during the period from 2014 to 2015, before the landslide occurred, the glacier at the toe melted and substantially retreated, which caused the buttressing effect at the toe to weaken, and the safety factor of this slope gradually decreased. The influence of the glacier, together with a combination of factors such as a larger than usual amount of rainfall in September and October in 2015, before the landslide, and a mild shaking due to the small earthquake in the distance, all contributed to the occurrence of this large landslide event.

Conclusions
In this work, a geotechnical study on the 2015 Taan Fiord landslide is provided, where stability analyses were carried out by employing 2D-LEM. Diverse possible external triggers for the landslide were investigated. The available DEMs were co-registered and used to better understand the morphological evolution of this coastal landslide. Representative cross-sections (of both slope and glacier body) were derived from DEMs for further numerical investigation. The main conclusions are as follows: (1) In the first stage of landslide formation, earthquakes are critical to the formation of a weak basal layer. The LEM analyses carried out on the 2015 Taan landslide qualitatively demonstrated that rainfall and seismic ground motions may have initiated the slope instability process.
(2) In the second stage of landslide formation, the final slope collapse resulted from glacial retreat at the toe of the slope. The weak basal layer formed in the first stage acted as the final sliding zone.
Remote-sensing data can be employed to constrain slope stability analyses to reconstruct the morphology of the slope when direct topographical data are lacking. This is especially true when multi-temporal aerial radar data are also available, as this paper demonstrates. In this interdisciplinary approach, the limit equilibrium method applied to sections of the DEMs is well-suited to identifying weak layers and possible formation processes. In this case, the weak layer likely formed through previous earthquake activities. The final triggering factor was most probably a combination of glacier retreat and precipitation.
Although the geotechnical parameters of the slope are unavailable in this study, plausible calculation results are obtained and the location of the calculated weak layer is close to the actual observed sliding surface location. This is achieved by parameter studies, in combination with reasonable reference values of similar conditions in other papers. This can also illustrate the applicability of the 2D limit equilibrium analysis in the analysis of large-scale landslides under complex conditions. Finally, this study effectively confirms the importance of an interdisciplinary approach that combines aerial radar data and traditional slope stability analyses. The implications of climate change for landslide hazard assessment and the interaction of physical processes in a multi-hazard context cannot be ignored.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.