Towards a rapid assessment of highway slope disasters by using multidisciplinary techniques

Mega-earthquakes and extreme climate events accompanied by intrinsic fragile geology lead to numerous landslides along mountain highways in Taiwan, causing enormous life and economic losses. In this study, a system for rapid slope disaster information integration and assessment is proposed with the aim of providing information on landslide occurrence, failure mechanisms, and subsequent landslide-affected areas to the highway authority rapidly. The functionality of the proposed system is deployed into three units: (1) geohazard rapid report (GeoPORT I), (2) multidisciplinary geological survey report (GeoPORT II), and (3) site-specific landslide simulation report (GeoPORT III). After landslide occurrence, the seismology-based monitoring network rapidly provides the initial slope disaster information, including preliminary location, event magnitude, earthquake activity, and source dynamics, within an hour. Within 3 days of the landslide, a multidisciplinary geological survey is conducted to collect high-precision topographical, geological, and remote-sensing data to determine the possible failure mechanism. After integrating the aforementioned information, a full-scale three-dimensional landslide simulation based on the discrete element method is performed within 10 days to reveal the failure process and to identify the areas potentially affected by subsequent disasters through scenario modeling. Overall, the proposed system can promptly provide comprehensive and objective information to relevant authorities after the event occurrence for hazard assessment. The proposed system was validated using a landslide event in the Central Cross-Island Highway of Taiwan.


Introduction
Rock slope instabilities along highways severely threaten the lives of road users and the safety of infrastructure (Sidle and Ziegler 2012;Bar and Barton 2017;Jaboyedoff et al. 2019;Sudmeier-Rieux et al. 2019). In the past two decades, mega-earthquakes and extreme climate events accompanied by intrinsic fragile geology have led to numerous landslides in the Central Cross-Island Highway of Taiwan, causing enormous life and economic losses (Hovius et al. 2000;Khazai and Sitar 2004;Dadson et al. 2004;Lin et al. 2006). On June 12, 2020, a destructive landslide completely buried more than 50 m of the road (Fig. 1). Such slope disasters frequently occur along mountain highways in Taiwan (Lin and Jeng 2000;Chen and Petley 2005;Lee et al. 2018); therefore, it is essential to provide useful information rapidly to the highway authority and the public for reducing maintenance cost and risk of traffic disruption.
The progress in slope stability analysis and failure prediction is closely associated with the development of new monitoring and assessment techniques, including satellite images, aerial images, light detection and ranging (LiDAR), interferometric synthetic aperture radar (InSAR), global navigation satellite system (GNSS), extensometers, crack gauges, inclinometers, tiltmeters, and seismometers (Burjánek et al. 2010;Wasowski and Bovenga 2014;Stähli et al. 2015; Barla and Antolini 2016;Dai et al. 2020;Rault et al. 2020). In the past decade, landslide seismology, which is a noninvasive approach, has been applied to comprehensively understand the physical process of mass movement (Chen et al. 2013;Pankow et al. 2014;Manconi et al. 2016;Chao et al. 2017). The failure mechanism of an unstable slope can be revealed through numerical approaches, including limit equilibrium analysis, finite element method, and discrete element method (Dawson et al. 1999;Collins and Znidarcic 2004;Eberhardt et al. 2004;Chiu and Weng 2019;Weng et al. 2019). Although some of the aforementioned techniques have been applied in specific cases only, a rapid, reliable, and integrated assessment system is highly needed to comprehensively monitor highway slope disasters and to provide crucial information, including the landslide scale, failure mechanism, influencing area, and subsequent potential failure zones for hazard mitigation strategies and emergency response implementation. To the best of our knowledge, there exists no above integrated system, and such a system is under development.
synthetic aperture radar (ATS-InSAR), the desk geological work, and the in situ geological survey is conducted to characterize the potential failure mechanism. GeoPORT II then can produce highprecision topographical, geological, and geomechanical models for subsequent analysis. The GeoPORT III report that is based on the full-scale 3D landslide simulation of discrete element method is provided to clarify the failure process and to identify the areas potentially affected by subsequent disasters through scenario modeling.
Overall, the proposed system can provide comprehensive (integrated disaster information) and objective (quantified analytical results) information in the context of the rapid assessment to highway slope instability. The detailed methods of three units are described below.

Seismic waveform inversion
Data preprocessing before the inversion involved removing the instrument response, mean, and linear trend, integrating the seismograms from velocity to displacement, rotating the horizontal records to the radial and transverse components for each station, and applying a bandpass filter with a frequency range of 0.4-0.8 Hz. The lower bound of filtering frequency should cover the duration of source process (toppled-mass impacting discussed in the section of geohazard rapid report (GeoPORT I)). The upper cutoff frequency of 0.8 Hz is chosen due to the inaccuracy in a simple one-dimensional (1D) velocity structure. With a priori knowledge of source location, a 1D velocity model (Chen and Shin 1998) and sinusoidal force-time function with a duration of 5 s were used to calculate synthetic seismograms. We adopted the general source inversion (GSI, Brown et al. 2015;Chao et al. 2016Chao et al. , 2017) assuming a single-force model to determine the optimal solutions for the sliding direction, dipping angle, and force magnitude.

Adaptive time-series interferometric synthetic aperture radar analysis
In recent years, the primary time-series analysis methods are mainly divided into two types: Persistent Scatter Interferometric SAR (PS-InSAR) and small baseline subset (SBAS). The permanent scatterers in the former image set can be inverted through the atmospheric phase to obtain the final high-precision surface deformation. The latter method mainly uses small-baseline image connections and adopts the interferogram after all-phase recovery for performing time-series deformation calculations. For heavily vegetated areas, Shih et al. (2019) mentioned that the traditional PS technique has difficulties in obtaining high PS density to promise the interpretation of landslide activity. Given the insufficient point distribution of PS-InSAR in processing nonurban areas, this study proposed the adaptive time-series interferometric synthetic aperture radar (ATS-InSAR) approach, which is based on PS-InSAR and extends the steps of processing distributed scatters. The processing flowchart is shown in Fig. 3. After conventionally calculating PS points by using the statistical verification method to identify homogeneous pixels and applying the multiprimary image, the small-baseline connection method was used to generate image pairs. Simultaneously, we used the statistically homogeneous pixel to filter analytical data and used the coherence as the weight. Temporal surface deformation, which is called points of the distributed scatter, was finally merged with PS points to form the final surface deformation map. For assessing landslide activity, the ATS-InSAR method provides denser ground deformation information than do traditional InSAR approaches.

Unmanned Aerial Vehicle (UAV) survey
On the basis of the results of aforementioned ATS-InSAR, we further conducted a field investigation and a UAV survey that focused on the active area with high precursory deformation inferred from ATS-InSAR. In the present study, a total of 1,366 aerial photographs with both vertical and oblique views were processed to obtain point cloud data by using ContextCapture (Bentley 2018). The point cloud Fig. 2 Flowchart of the rapid slope disaster information integration and assessment system data were georeferenced using UAV onboard GPS, two ground control points, and one checkpoint measured in the field through realtime kinematic positioning. Further processing of the point cloud data enabled orthophoto map production with pixel dimensions of 0.05 m and a digital elevation model of outcrop with similar resolution. Combining ATS-InSAR analysis and in situ descriptions provided multidisciplinary data for establishing a database for clarifying the failure mechanism of the landslide.

Numerical simulation
To assess the area affected by the landslide and its subsequent effect, we performed failure inversion analysis by using discrete element software 3DEC (Itasca Consulting Group Inc. 2016). The numerical simulation workflow proceeded as follows. First, a triangulated mesh representing the topography surface was constructed according to the pre-event and post-event topographic data. Here, the meshing process was conducted using Rhinoceros CAD software and CloudCompare (Dewez et al. 2016). Second, the observed discontinuity sets (DSs) investigated from the field survey were explicitly introduced in the numerical model. Simulation was performed and validated with the actual influenced area to reproduce the dynamic process of the landslide. Finally, a scenario modeling was also performed with a considering of the potential failure area to identify the impact of subsequent event.

Geological and topographical background
The practical application of the proposed system was executed on the destructive landslide on June 12, 2020. Here, we describe the topographical and geological backgrounds of the event. The landslide occurred on the left riverbank of Dajia River downstream of the Deji reservoir (Fig. 1a). The landslide blocked the middle and lower lanes and destabilized the upper lane of the highway, affecting an approximately 250-m-long and 50-m-wide area (Fig. 1a,c). The geology of the highway can be categorized as the Paileng Formation of the Oligocene period, which comprises massive metamorphic sandstone with occasional slate interlayers. The metamorphic sandstone is highly cemented and can be classified as strong rock (Fig. 1a). The terrain revealed that the slope height is approximately 250 m from the source area and has a steep slope angle of 50 degrees on average, which revealed that the slope failure type was slope toe sliding induced by overall toppling failure.

Camera-recorded video
The authorities reported that rocks started falling on the road approximately an hour before the landslide occurred; hence, they shut down the road in advance for safety. When the falling rocks became larger and more frequent, the authorities started recording the video from a safe distance (Movie 1). Figure 4 shows the landslide process obtained from the video recorded using the camera in a time series, indicating that the slope failure process was slope toe sliding-induced overall toppling failure. A total duration of entire event process is about 72 s. In the first 17 s, rock mass fell from the high slope and rolled toward the downslope ( Fig. 4a-e); then the block mass of source area toppled and hit the middle lane ( Fig. 4f-g). In the end of video, the rock mass broke into fragments and deposit on lanes and river (Figs. 1c and 4h).

Seismicity and rainfall data before event
Slope failure may be triggered by seismic excitations or rainfall. To verify the triggering factors, precipitation data recorded using a rain gauge station were collected, and earthquake forcing (F EQ ) was estimated from May 12 to June 13, 2020. The precipitation data recorded at the rain gauge station C1F9W0 (Deji rain gauge) operated by the Central Weather Bureau (CWB), which is approximately 4 km away from the landslide, were collected in this study for investigating the rainfall response on the slope (triangle shown in Fig. 1a). To clarify the earthquake response at the study site, we first collected the earthquake catalogs of Wu et al. (2008) and BATS (https:// bats. earth. sinica. edu. tw), and then, we established a linear relationship between the local magnitude (M L ) and moment magnitude (M w ) as shown in the below equation: For our study site, we collected the earthquake catalog of the CWB (https:// gdms. cwb. gov. tw/) from May 12 to June 13, 2020 (Fig. 5a). The M L of earthquake events was converted using Eq. (1), and a formula was applied for linking the seismic moment (M 0 ) and M w of Hanks and Kanamori (1979) to compute the logarithm of M 0 . We defined the earthquake forcing caused by a specific earthquake considering the epicentral distance (D) correction as follows: Hourly earthquake forcing is a summation of earthquake forcing per hour, which is used to discuss the earthquake response at the study site. Within 10 days before slope failure, no obvious indication was detected based on rainfall and earthquake seismicity at the study site ( Fig. 5b). An earthquake with a local magnitude (M L ) of 4.7 occurred 38 min before slope failure at 60.3 km away from the sloping site, resulting in a weak seismic response at the study area (Fig. 5c). Thus, additional observations and investigations, such as seismic constraints, geological survey, and numerical simulation, are required to advance the understanding of the slope failure mechanism.

Seismic records
The data from two seismic stations with an epicentral distance ranging between 665 and 916 m, namely, TWT and TDCB, were used in this study (Fig. 1a). Station TDCB is maintained by the Broadband b A rock fell to the road. c,d Rocks falling and rolling. e Rock mass detached from the slope, and toe sliding occurred. f,g Rock mass toppling and interacting with the slope. h Recording ended Array in Taiwan for Seismology (BATS, https:// bats. earth. sinica. edu. tw, last accessed July 2021) and is equipped with STS-1 seismometers and a Q330HR datalogger. The Central Weather Bureau (CWB) deployed station TWT, which is equipped with an S13 short-period seismometer and a force-balance accelerometer (https:// gdmsn. cwb. gov. tw, last accessed July 2021). Both broadband and short-period waveforms are used in seismic waveform inversion. Furthermore, station TWT accelerograms and simulated force history were compared (Fig. 5d) in the section of results and discussion.

Results and discussion
Geohazard rapid report To rapidly provide the magnitude, and dynamics of the event, we conducted a series of seismic analysis to initially determine event characteristics, especially the source volume and mass-movement process based on the waveforms recorded by two seismic stations (Fig. 1a). In our case, only two detected stations cannot offer the accurate location, which is demonstrated by Chang et al. (2021). Thus, we fixed the source location at the center of landslide area for the following seismic analysis. The seismograms and the video recorded using the camera (Movie 1) were analyzed to comprehensively extract information on the entire movement, which can be divided into three stages (each stage is shown in Fig. 6a). Stage I involved multiple fine granular mass falling or rolling, with weak power spectral density (PSD) at a frequency of 4-20 Hz for initiation. In Stage II, a column-shaped spectrogram was obtained (a period of 37-42 s shown in Fig. 6a) that was induced by the toppled mass affecting the ground, which generated relatively strong PSD, especially for the vertical component, and the massive rock mass impact led to a low frequency (< 1 Hz), which is consistent with the results of the field experiment of Huang et al. (2007). Stage III, the final stage, involved multiscale granular material interacting with the topographic surface; the spectrogram exhibited an up-gliding feature caused by fragmentation of the toppled-boulder mass or mobilization of the basal rock deposits. Relatively low-frequency (0.4-0.8 Hz) seismic waveforms generated by the toppled mass impact during Stage II, which can be observed clearly in spectrograms (Fig. 6a), provide an opportunity to study force properties, including the magnitude, direction, and dipping angle of impacting force. Thus, the general source inversion approach (GSI) was conducted using the single-force (SF) model. The aforementioned unknown parameters were determined through minimization of the difference between synthetic and Fig. 6 Time series of seismic signals, spectrograms, force history, and results of seismic inversion. a Seismograms and spectrograms for station TDCB and the simulated vertical force history. A horizontal solid black line indicates the waveform window used in GSI. Colors shown in spectrograms are of PSD. Vertical dashed lines mark specific time points for defining three stages. b Results of GSI inversion. Gray and red curves are observed and synthetic filtered (0.4-0.8 Hz) displacement seismograms, respectively. The arrow shows the horizontal force vector, and numbers indicate the absolute maximum force magnitude (F max ), sliding direction, and dipping angle. The station name, signal-tonoise ratio value, time shift, normalized cross-correlation coefficient (CC), and variance reduction (VR) are given at the top of each trace observed seismograms (Fig. 6b). SF inversion yields a maximum absolute force (F max ) of 1.285 × 10 8 N in the northwestward direction (330°), which is discussed in the site-specific landslide simulation section. A high dipping angle of 70° supported that the impacting force acted vertically primarily (Fig. 6b). Using an empirical equation of mass = 0.405 × F max derived by Chao et al. (2016), we found a mass of 21.7 million kg and a volume of 8,680 m 3 (assuming an average density of 2,500 kg/m 3 ). Furthermore, we simulated Wood-Anderson seismograms by using seismic records from the two seismic stations (TWT and TDCB) and estimated the average M L to be 0.79. Moreover, we used the peak amplitude in the horizontal envelope function of the closest station as the seismic amplitude at the source location (A 0 ). The horizontal envelope function was extracted from the filtered (1-8 Hz) root-mean-square (RMS) amplitude of the horizontal waveforms (N-S and E-W components). Based on two scaling equations linking local magnitude (M L ) with Fig. 7 The ATS-InSAR analysis with ESA Sentinel-1A images from December 3, 2019, to June 12, 2020. a Slope deformation distribution before the landslide. The symbol colors represent LOS cumulative displacement. b The two greatest deformation points located in the study area are highlighted, and the slope deformation became intensive from May 19, 2020 volume and the seismic amplitude at the source location (A 0 ) with volume established by Chang et al. (2021), the source volumes of 9,220 m 3 and 5,238 m 3 can be predicted, respectively. In summary, the seismic technique resulted in an event volume ranging from 5,238 to 9,220 m 3 based on the peak values of force and seismic records, which constrain the lower limit of the volume of the original collapse mass. Since the real-time stream of seismic waveforms was ready, aforementioned seismologically inferred source characteristics could be determined rapidly after the event, which are crucial parameters for performing numerical simulation and field investigation.

Multidisciplinary geological survey
To probe the slope activity before sliding, we conducted ATS-InSAR analysis immediately after detecting substantial slope failure. The time-series measurement showed that LOS displacement began to increase on May 19, instigating a landslide within 1 month (Fig. 7), which may be associated with abundant precipitation and strong earthquake forcing in late May (Fig. 5b). The results simply help us to better understand the slope activity preceding sliding and be a reference to discuss its possible failure mechanism.
The sliding volume of the landslide was determined through the following processes. The differences between the 2012 LiDAR 1-m digital terrain model (DTM) and 2020 digital elevation model at the source area were calculated to preliminary obtain a fallen volume. The result was then calibrated through a comparison of the deposition and entrainment volumes between digital elevation models of two periods. Finally, the sliding volume was determined by confirming with the results of the geological survey. The processing result indicates that the main landslide area resulted in a sliding rock mass with a volume of approximately 22,500 m 3 (Fig. 8a). Considering the complex landslide structure and site morphology, the measurement uncertainty of the landslide volume accounted for an error rate of up to 10%. Elevation profiles generated using high-resolution terrain data contributed to the identification and visualization of major changes and the deposition feature (Fig. 8b). The main landslide source located between the upper and middle lanes is clearly shown in profile A-A', as well as the deposition of landslide material in the slope toe. Profile B-B' shows that the source block was cut by lateral erosion gullies before the landslide. The analysis revealed that the greatest depth of the sliding mass was approximately 22 m and its location is identical to the observed displacement through ATS-SAR interferometry (Fig. 7).
Discontinuity orientation data were obtained using the UAVphotogrammetry approach and were validated using the field investigation from outcrops (Fig. 9a). Four DSs were identified in the source area (Table 1; Fig. 9b). The point cloud data analysis indicated that the mean orientations of DS1 and DS2 were 83°/225° and 41°/030°, respectively. They formed the main failure surface of the landslide. DS3 dipped into the slope with a high dip angle (80°/041°), inducing rock blocks to topple in the upper slope. On the basis of its long persistence and the orientation being approximately parallel to the orientation of erosion gullies at the site, DS4 appeared to be the bedding plane with an orientation of 62°/099°. In addition to discontinuity orientations, point cloud data were used to interpret the persistence (P) and spacing (S) of each DS by using CloudCompare (Dewez et al. 2016). The analysis allowed the consideration of nonpersistent discontinuities and realistic block size for further scenario modeling.
Kinematic analysis can show the potential failure mode for the landslide. On the basis of predisaster 1-m LiDAR DTM (Fig. 8), the slope orientation at 44°/047° was considered in the analysis. The friction angle along all discontinuities was assumed to be 20° due to surface condition discontinuity. Kinematic analysis results revealed that the main failure might be initiated by the plane slide along DS2 (Fig. 9c). Due to the initial removed blocks, the secondary toppling phenomenon occurred with the blocky rock mass cutting by the DS1 (Fig. 9d).

Site-specific landslide simulation
According to in situ measurements, the identified landslide area was 30-m wide and 75-m long (Fig. 10a). A 2012 LiDAR 1-m DTM was used for this study, which could satisfactorily reproduce the original topography feature for the proposed 3DEC mechanical model. The DTM data used for topography construction were obtained before the landslide event. Thus, topography description obtained from the point cloud of the LiDAR DTM was not affected by the landslide. In this study, the UAV survey was performed few days after the landslide event. The obtained UAV data could help us to identify the sliding area and estimate the volume loss due to the sliding event. The triangulated surface mesh, in STL or DXF format, was imported into 3DEC as a geometry object. An extrusion process, from a surface object to a solid object (i.e., block), using 3DEC macro language (FISH), was then executed to obtain the model geometry. The extrusion depth (35 m), from each single triangle mesh, should be deep enough compared with the real sliding depth (roughly 20 m) to allow the natural development of progressive Fig. 9 Discontinuity set analysis. a Discontinuity features within the source area were obtained from the unmanned aerial vehicle-based point cloud model. b The dip and dip direction of discontinuity sets were automatically acquired using the plane fitting algorithm in CloudCompare. The persistence and spacing of discontinuities were semi-automatically determined in point clouds. c Possible failure mode: plane sliding. d Possible failure mode: flexural toppling rock sliding depth in the model. The meshing process is summarized in Fig. 10. The dimension of the 3DEC model was 200 m along the slope strike, 550 m along the slope dip, and 370 m high (difference in the vertical elevation between the highest and lowest grid points), with approximately 82,000 rigid blocks. Once the 3DEC geometry model was prepared, the sliding area was marked based on the data collected from the geological survey. According to the in situ investigation, the main sliding area was located between the upper and middle lanes as green blocks in Fig. 10c. Moreover, site investigation showed that due to the main event, the region located immediately above the upper lane displaced along the slope dip. Therefore, a subsequent landslide might occur in that area. The scarp of the main event and two well-developed gullies provide a potential landslide area for subsequent failure. The potential area of the scenario event is represented with blocks in blue in Fig. 10c, located within the range of the two gullies (marked with a dotted red line in Fig. 10b).
Four DSs were identified from the exposed slope surface in the studied site (Fig. 11). The blocks representing the two events were further subdivided into smaller blocks based on those fracture sets by using 3DEC joint elements. Discontinuity persistence and discontinuity spacing were accounted for when proceeding with the block cut. Considering the computational efficiency, only blocks representing the two sliding areas could move. In view of the mechanical properties of metamorphic sandstone, the adopted numerical model is a rigid block-based model. Each individual block behaves as a rigid body and can interact with other blocks through contacts of discontinuity. The required contact parameters of discontinuity include normal and shear stiffness, and friction angle of Coulomb slip criterion. A high joint stiffness value of 1 10 3 MPa/m was used to reflect the collision force between blocks during the sliding. According to Barton (2002), the inter-block friction angle could be estimated from the Jr and Ja rating tables. The estimated rating of Jr and Ja were 1.5 and 6, respectively. Therefore, an initial friction angle of 21° was obtained. The weathering effect was then accounted for through a degradation process on the friction angle until the instability of the slope was triggered. In order to facilitate the tensile cracks before the landslide and could better reproduce the phenomena of a slope toppling, friction angles of the two sub-vertical joints sets were mainly degraded. The initial slope instability can be captured numerically when a drastic increase in kinetic energy is observed. Table 2 summarizes the material properties adopted in the numerical analysis.  The simulation of the main landslide event showed that small blocks at the slope toe initially detached sporadically and slid along the slope surface (Fig. 12a), which induced a series of tensile failures on the subvertical DSs (DS1 and DS3) in the middle slope. These nearly parallel failed discontinuities separated the slope into several rock pillars, and a curved failure surface formed in the lower slope and connected to the failed subvertical joints. Finally, a series of flexural toppling of the rock pillars occurred, and the main failure surface continued to propagate to the slope top. The simulated failure mechanism was consistent with the findings of the recorded video (Movie 1) and seismic waveform analysis (Fig. 6). Thereafter, the sliding blocks traveled along the concave part of the topography. They finally stopped and accumulated at the riverbed of Dajia River (Fig. 12b). Figure 12c provides a comparison the 3DEC model-estimated affected area and DTM data. The simulation could estimate the influenced area revealed by DTM data. However, there is discrepancy observed in influenced area between simulation and DTM. A part of discrepancy originated  due to vegetation. In reality, the speed of sliding blocks might be affected by the vegetation, which was not accounted for in the simulation model. Kinetic energy was dissipated through the frictional force and gravitational damping force at each block. Thus, almost all the sliding blocks eventually reached and stopped at the riverbed, resulting in a wider influenced area at that place. The main event may induce subsequent landslides in future; thus, we further simulated and evaluated the scenario event in the upper lane area. The sliding area of the scenario event in the upper lane is situated between two gullies that developed and connected in the upper slope. The scarp of the main event and two well-developed gullies provide a potential landslide area for subsequent failure. Unlike the mechanism of the main event, the scenario event exhibits translational sliding (Fig. 13a). Due to the loss of rock blocks after the main event, a daylight condition forms at the upper slope, and a more pronounced sliding surface develops along DS2. The block volume loss produced by the scenario event is less than that of the main event (blue part in Fig. 13b). Block volumes of 13,901 and 8,627 m 3 were lost due to the main and scenario events, respectively. Most sliding blocks would be eventually distributed on the existing deposit area of the main event. The subsequent landslide would enlarge the influenced area of the upper slope of the upper lane and may induce backward slope instability in the future (Fig. 13c). These areas should be recommended to promote the preparation of disaster prevention and mitigation.

Capability and limitation of proposed system
Overall, the information of slope failure can be temporally updated via our proposed system. As shown in Fig. 2, the unit 1 (GeoPORT I) of our proposed system is a pilot approach associated with the seismic monitoring, providing information about timing, location, collapsed volume, and source dynamics of event within a short time to the highway authority for rapid hazard assessment. However, uncertainties in the results inferred from seismic analysis mainly rely on the data quality of observed seismograms and the number of available seismic stations. To improve the station coverage  of seismic network, an additional broadband seismic stations deployed by the GeoPORT Working Group recently were included in our proposed system. To further discuss the reliability in the source dynamics of the GeoPORT I, an integration of the results of numerical analysis and seismological-inferred source dynamics (gray line shown in Fig. 6a) showed impact force evolution of the main event due to sliding blocks on the middle lane. After aligning the time point of the peak vertical impacting force of 33 × 10 6 N and with the peak acceleration inferred from seismogram ( Fig. 5d and Stage II shown in Fig. 6a), force duration history was consistent with the signal duration of landslide-induced waveform. However, there is a clear difference existed in the forces of toppled mass results from simulation and seismic waveform inversion. Several factors may contribute to this discrepancy. In waveform inversion, seismic record quality, imperfect station coverage, and velocity model may lead to uncertainty in the inverted force magnitude.
In practical application for the highway authority, rapidly estimating the source volume is a key to guide safety assessment for a purpose of resuming the traffic. Based on the multidisciplinary geological survey (GeoPORT II), the total volume composed of the sliding rock mass and subsequent fragmentary rockfalls, which is approximately 22,500 m 3 , can be directly computed. However, the GeoPORT II cannot report the source volume within a short time. For a purpose of emergency response, seismic constraints on event volume provide a potential solution to above gap, even though event volume estimated by using the peak values of seismologically determined force and seismic signals could only constrain the lower limit of the volume of the original collapse mass (a range of 5,238-9,220 m 3 ). Ideally, a fully automatic process of seismic analysis that can yield source information including the occurrence time, preliminary location, event volume and force mechanism rapidly, and manual spectrogram analysis and identification is required to extract the physical process. Then, with the knowledge of source location and volume, we could issue the emergency traffic regulation and public notice for reducing the potential impact due to traffic disruptions. Furthermore, seismology-determined force properties and physical process could also be provided soon after landsliding for designing the engineering-based slope protection.
In a case of seismic data lacking, the GeoPORT II and III could also be conducted if the event captured by video and/or eyewitness. However, the system cannot provide timely information without the GeoPORT I. Thus, a direct way to improve the performance of current system is to deploy the temporal seismic arrays along highways for rock slopes with high-potential failure (Chang et al. 2021). Most techniques used in the GeoPORT II are in situ, which may be limited in extreme weather conditions. Current GeoPORT system has been examined successfully by a destructive rock slope failure only. Rock slope failures with the variable source mechanisms occurred at the area with different geological setting will be used to test our proposed system.

Conclusions
An integration system of the seismic analysis, multidisciplinary geological survey, and numerical simulation, named GeoPORTs (Fig. 2), is proposed with the aim of providing information on landslide occurrence, landslide scale, failure mechanisms, and subsequent landslide-affected areas to the highway authority rapidly, which would be helpful in designing the mitigation structures or to identify similar locations prone to hazard. Since the real-time stream of seismic records is ready, GeoPORT I can be implemented automatically to provide a preliminary disaster evaluation as a trigger for the following reports of GeoPORT II and III. Although a discrepancy existed in source characteristics between seismic-inferred result and simulation was found, overall, the GeoPORTs provide comprehensive and objective information to relevant authorities for rapid hazard assessment and determining emergency mitigation strategies. More landslide cases along highways should be collected and analyzed in the future to determine the applicability of the proposed system in slope engineering practice.