Modelling of mobility of Rissa landslide and following tsunami

Landslides in sensitive clays pose a major threat to life, property and the environment. The lack of warning signs and the extreme mobility of quick clays increase the risk. When occurring along fjords, lakes or large rivers, the landslides can generate destructive tsunamis. Historical records show that 45% of the shoreline landslides in marine sediments in Norway triggered tsunamis with run-up height of 1 to 15 m. To improve the management of the hazard and risk due to landslides in sensitive materials, the profession needs to develop tools for modelling the mobility of landslides and their tsunamigenic potential. The paper back-calculates the 1978 Rissa landslide, one of the largest quick clay landslides in Norway. Because of the quantity of data available, the Rissa landslide provides a unique benchmark to study the mobility of the quick clay, and the tsunami following the landslide. The paper describes a new approach to model landslide mobility and tsunami run-up, and presents the results of the analyses and a comparison of the results with the observations. The back-calculated mobility (runout distance, flow velocity and debris thickness) and wave run-up on the lakeshore agreed well with the measurements.


Introduction
Landslides in sensitive clays pose a major socio-economic threat to population, infrastructure, property, and the environment because of their retrogressive characteristics and extreme mobility.A recent example is the fatal quick clay landslide in Gjerdrum in Norway on the early morning of December 30, 2020: the landslide caused 10 fatalities, destroyed over 31 dwellings, and forced the evacuation of over 1,000 residents in pandemic times (which exacerbated infections), and resulted in chaos in the roads, services (e.g., sewers, fresh water pipelines) and the ecosystem downstream of the landslide.The cost for mitigation work around the Gjerdrum landslide will be over 16 M€, excluding the rebuilding of the local infrastructure, re-localization of houses and the environmental consequences.The Gjerdrum landslide is only one of many recent catastrophic landslides in Norway.Fig. 1 shows a few of the quick clay landslides in the past 12 years, including landslide in a fjord due to road construction (Kattmarka), landslide destroying a highway bridge (Skjeggestad), destruction of roads and houses (Kattmarka and Sørum).The Lyngen and Sørum landslides illustrate the extreme mobility of quick clay.
About 5,000 km 2 of Norway is covered by soft marine deposits, whereof 20% consist of highly sensitive or quick clay.These areas attract human settlement because they provide gently inclined and fertile land in otherwise rough mountainous terrain.Currently in Norway, over 110,000 people live on approximately 2,300 quick clay zones.Over 150 persons have perished in quick clay landslides in Norway over the last century.A review of the landslide database shows that over 85% of the quick clay landslides since year 2000 were triggered by human activity, most often in combination with natural erosion and/or unfavourable groundwater conditions (i.e., snowmelt and/or intense precipitation), both of which are aggravated by climate change.
The stability of slopes in quick clay depends on the material properties of the clay, imposed shear stresses and external factors, and changes to any of these elements will also impact the stability and thereby the risk.Even though quick clay challenges are well known to all stakeholders in the building, construction and transportation sectors in Norway, the frequency of quick clay landslides with volumes greater than 50,000 m 3 has nearly doubled (from 0.7 to 1.3 events per year) over the last two decades in Norway (L'Heureux et al., 2014;2018).Fig. 2 shows the frequency of recent large quick clay landslides (denoted QCL in figure) with volume >50,000 m 3 since 1970.There is also a large number of smaller landslides that disrupt everyday life, pose a risk to the population and have wide environmental impact, as well as lead the population to question the geotechnical profession ability to cope with landslide hazards.The reasons for this increase are: i) an increase in infrastructure development, ii) worsening slope stability over time due to natural processes such as erosion, intense rainfall and floods, and iii) probably a misunderstanding of the existing hazard and risk maps drawn up for quick clay areas.
When occurring along fjords, lakes or large rivers, the larger landslides can generate destructive flood waves or tsunamis.In August 1905 for example, glaciolacustrine clay deposits slid into the Thompson River at Spences Bridge in southwestern British Columbia (Evans, 2001).The landslide generated a 5-6 m high wave that rushed more than 1.5 km upstream, destroying twenty buildings and drowning 15 people.Three years later, in 1908, a landslide occurred suddenly on the Lièvre River in western Québec (Evans, 2001).The landslide created a displacement wave that overwhelmed part of the village of Notre-Dame-de-la-Salette, killing 27 people.Historical records show that 45% of the shoreline landslides in marine sediments in Norway triggered tsunamis with runup height of 1 to 15 m (L' Heureux et al., 2013a;b).Some examples are the 1930 Orkdalsfjorden landslide (15 m wave, L 'Heureux et al., 2014b), the 1959 landslide at Sokkelvik (up to 6 m wave, L 'Heureux et al., 2017) and the 1978 landslide in Rissa (up to 6.8 m high wave, NGI, 1978;Gregersen, 1981).The tsunamigenic hazard posed by landslides in nearshore areas underlain by marine clay deposits is usually not included in the hazard and risk maps today.Yet the generated tsunami can also have disastrous consequences.
To develop future hazard and risk maps and apply appropriate risk assessment methods, there is a need for an improved understanding and modelling of (1) the mobility of landslides (where mobility is runout distance, flow velocity and thickness of the debris with time), and (2) the tsunamigenic potential of larger landslides along bodies of water.
To further the profession's ability to deal with these hazards, the paper focuses on modelling the largest quick clay landslide in Norway in the last century, the 1978 Rissa landslide.As the Rissa landslide rapidly moved into Lake Botnen, it generated a tsunami wave that impacted the village of Leira 5 km from the landslide across Lake Botnen (NGI, 1978;Gregersen, 1981;L'Heureux et al., 2012).To model the landslide mobility and the tsunami generation, new formulations were developed and are used in this paper.The Rissa landslide was used for a validation of the calculation models because it has eye-witness reports and amateur videos of the landslide 1 , detailed bathymetry of the landslide, a mapping of the tsunami run-up depths on the day after the landslide, geophysical and geotechnical data from the site at the time of the landslide and some 20 years later.The Rissa landslide provides thus a unique benchmark allowing to study the failure, the mobility of the quick clay, and the tsunami following the landslide.
The paper is organized in six parts: (a) a summary of the Rissa Fig. 1.A few of the recent devastating quick clay landslides in Norway. 1 The reference list gives the YouTube link to the Rissa landslide video.
Z. Liu et al. landslide and available data; (b) a short review of current models that analyses the mobility of landslides; (c) a description of the new constitutive models and analysis methods for the prediction of (1) quick clay landslide mobility and (2) the run-up height of the tsunami generated by a landslide, both with more details in Appendix A; (d) an overview of the analyses done and a comparison of the results with observations; (e) a discussion of the results, their uncertainty and how the new approach can be implemented in practice; and (f) summary and conclusions.

Geological setting of Rissa landslide
Rissa is a small county located 25 km northwest of Trondheim in mid-Norway, resting along the shore of Trondheimsfjord.The village of Rissa is on the northwest shore of Lake Botnen, a narrow and long and brackish inlet 1 km by 5 km connected to the fjord.Fig. 3 shows a map of the Rissa agglomeration, Lake Botnen, the surrounding farms, the village of Leira at the northern end of Lake Botnen, and the location of the landslide in 1978.
Following the ice retreat, the sea level rose slower than the local land surface, resulting in an apparent marine regression and the emergence of the old sea bottom.Today, the lowlands along the shores of Lake Botnen are relatively flat and are almost entirely covered by thick glacio-marine and marine deposits, locally overlain by littoral deposits (Reite, 1987).During their emergence in the Holocene, the clays became exposed to a flux of fresh groundwater, which gradually leached the high salinity porewater from the sediment.The leaching led to an increase in interparticle repulsive forces (Rosenquist, 1953;Torrance, 1974;1983) and caused the high sensitivity of the marine sediments.The leached clay, or quick clay, has "high" peak shear strength and very low remoulded shear strength.When subjected to stresses that exceed the peak shear strength, the structure of the clay collapses.In Norway (NGF, 2011), a sensitive, brittle clay is defined as a clay with a remoulded shear strength less than 2 kPa; a quick-clay is defined as a clay with a remoulded shear strength less than 0.5 kPa (also in NGF, 2011).
Thick quick clay deposits are found along the southern shore of Lake Botnen.Earlier morphological analyses of Lake Botnen revealed a relatively flat sloping lake floor locally with steep shoreline slopes (up to 35 • ) (L'Heureux et al., 2012).The basin reaches a maximum water depth of 38 m in the central part of the lake.Large pockmarks (up to m wide and 5 m deep) and several mass-wasting deposits were found resting on the lake floor.
To investigate the landslide deposit and the soil profiles underwater, high resolution bathymetric surveys of Lake Botnen were acquired in 2010.Positioning within ±1 m was ensured with differential GPS.A grid of seismic reflection data was collected along 15 seismic lines using a 3.5 kHz parametric sub-bottom profiler.Two-way travel time was converted to water depth and sediment thickness assuming a constant sound velocity of 1470 m/s.The Geological Survey of Norway also carried out 2D resistivity measurements onshore and around Lake Botnen in 2009, 2010 and 2011.

The Rissa landslide
The 1978 landslide took place at the south-western corner of Lake Botnen (Fig. 3).Thirty-nine persons narrowly escaped from the landslide, and one person died.Gregersen (1981) described the landslide as a two-stage process.At first, a small "initial" landslide was triggered due to excavation to expand the barn of Fessöya farm.The excavated soil, stockpiled along the lakeshore, caused this initial landslide (Fig. 4a): to 90 m of the shoreline slid out into the lake, including half of the recently placed fill.The scarp of the landslide was at that time 5 to 6 m high and extended 15 to 25 m inland.
The landslide then developed retrogressively in the south-western direction over the next 20 to 40 min, based on eye-witness reports.This is the Stage 1 landslide (Fig. 4b).The sediments completely liquefied and the debris literally poured into the lake.In this first stage, the Z. Liu et al. landslide area had the shape of a long and narrow pit open towards the lake.The length of the sliding area was 450 m covering an area of 25 to 30,000 m 2 (or 6-8 % of the final slide area (Gregersen, 1981).
The second stage, the main landslide, started almost immediately after the retrogressive sliding had reached one of the roads.At this point, large flakes of dry crust (150 × 200 m) started moving towards the lake, not through the existing gate opening, but in the direction of the terrain slope (Fig. 4c, d, e).The velocity was initially moderate (Flake A in Fig. 4c) with velocity of about 10 to 20 km/h.The velocity increased to 30-40 km/h as Flake B moved (Fig. 4e).In the amateur video of the sliding (see additional material), houses and farms can be seen floating on the sliding masses.A series of smaller and retrogressive landslides followed over a short period of time.The sliding process propagated backwards up to the rocky mountainside where it stopped (Fig. 4f).The main landslide lasted approximately 5 min and covered 92 to 94% of the total slide area (0.33 km 2 ).The total volume of mobilized sediment was 5 to 6 x10 6 m 3 (Gregersen, 1981;L'Heureux et al., 2012).

Morphology of the Rissa landslide
Fig. 5 presents a morphological interpretation of the landslide.The final width of the failure scar above the shoreline was 400 m.The south eastern boundary of the landslide followed the base of the mountainside and the total distance from the lake to the head scarp was 1400 m.The deposits poured into the lake and still cover today an area of 760 km 2 , or 20% of the lake floor.The maximum length of the accumulation zone is 1200 m.Seismic data showed an average debris thickness of 6 to 8 m.This gives a total landslide volume in Lake Botnen of about 4 × 10 6 m 3 .Most of the flow took place through a main gate west of a morphological high (Fig. 5).This high, interpreted as bedrock from seismic data, forced the main landslide to open a new gate to the west (denoted as 'Main gate' on Fig. 5).Two types of mass transport deposits were identified on the shaded relief image in the lower Fig. 5: 1) The first type consisted of well-constrained sediment lobes of fine grained sediments and containing randomly deposited small blocks/ hummocks (less than 20 m in diameter).The lobes showed many flow structures, whereas the blocks show a semi-circular to triangular footprint.The distal portion of the lobes was characterized by transverse undulations (i.e., normal to the landslide direction), with arcuate ridges up to 100-200 m long.Such ridges are usually compressional and are mainly associated with flow deceleration combined with continued sediment flow piling up (Posamentier and Kolla, 2003).A 400-m long and up to 1 m high longitudinal ridge was formed on the basin floor in the vicinity of the lobes.From the illuminated swath bathymetry image, part of the sediment lobes can be traced back to the open gate west of the morphological high.The character of the lobes suggests that they were deposited in Stage 2 when the clay and debris had liquefied.The longitudinal ridge may result from simultaneous flows with different velocity, combined with eddies generated to the north of the morphological high.Large swirls were witnessed in the lake water during the event.Some of the hummocks on the basin floor are probably anthropogenic deposits (e. g., houses) that sank during or several days after the landslide.
2) The second type of mass transport deposits was found in the frontal part of the landslide deposit (Fig. 5).The blocks show rectangular and elongated forms and are typically 80 to 200 m long and 30 to 100 m wide.Some of these blocks travelled up to 1150 m from the shoreline.The blocks rise today 1 to 2 m above the lake floor and are oriented transverse to the flow direction.The grouping of the blocks at the outer rim of the deposit suggests that these features are related.Rafted blocks have been frequently observed in translational landslides (Mulder and Cochonat, 1996) and often are associated with debris flows (Ilstad et al., 2004).As such, the larger blocks appear to be remnants of the large flake sediments that rafted towards the lake early in the landslide.Seismic data also showed that some of the blocks were covered by thin sediment plumes in the final stages of the landslide.

Tsunami wave across Lake Botnen and impact
The Rissa landslide on the south of Lake Botnen generated a tsunami that damaged the village of Leira at the north end of the lake.The day after the landslide, NGI made observations of the damage and run-up heights and run-up distances in Leira.Fig. 6 presents the original record of the observations of run-up heights and the lower sketch defines run-up height and run-up distance as used in this paper.The vertical height above Lake Botnen for the uppermost wet point in the inundation zone after the large amount of water that the tsunami pushed onto the  (Gregersen, 1981).
shore is called Run-up height, or the maximum vertical height onshore above the lake level, and the horizontal distance reached on the ground is the run-up distance.In the inundation zone the water have high velocities and may lead to much larger damaging forces here than for the tsunami wave in the lake.During the impact the tsunami can destroy or damage all in its path.
The measured run-up heights on the shore were between 1.8 and 3.2 m on the northwest side of the lake and between 1.5 and 2.2 on the southeast side.The maximum run-up height was 6.8 m near the landslide and 3.2 m in the village of Leira.On the basis of eye-witness accounts, the tsunami wave took 6 to 7 min to cross the lake from Fessöya to Flyta, just close to the village of Leira.The tsunami had therefore an average speed of 11 to 13 m/s which agrees with the modelled arrival times The run-up height was recorded as 2.8 m at the location "Naust" (boathouse), 1400 m away from the landslide, 1.8 to 2.0 m north of Straumen, 2500 m from the landslide, and l.5 m close to Sjølia, about 2.5 to 3 km from the landslide.Close to the dwellings in Leira, the run-up height was recorded as 3.2 m.
At the boathouse ("Naust"), the wave was described as a grey-black mass of water arriving as if overturning, and it stood "like a sea wreck around a reef" at the 5 m high islet next to the houses.Approximately the same description was used at Straumen.The wave came rumbling forward over the ground and scared the people living by the lake.They fled their houses when they saw the wave.At Leira, on the other hand, the wave came as a surprise, the people did not have the time to run away.The wave was described as a large wave breaking and hitting the houses and the infield.It moved and destroyed, among others, the sawmill buildings (Fig. 7).In other locations, garages were deformed, houses were heavily damaged or destroyed.One person travelling on the road was moved out in the field at Flyta.The run-up height observation made on site in 1978 is penned in a dark line on the lower photograph in Fig. 7.

Geotechnical properties of Rissa sensitive clay
An extensive soil investigation program was carried out in the Rissa landslide area in 1978 immediately following the landslide, including soundings, vane shear tests and undisturbed sampling within and Fig. 6.Map of tsunami run-up heights measured around Lake Botnen (NGI, 1978), and definitions of run-up height and distance; Fessöya is where the landslide started; "Leirskred" means 'clay landslide').
outside the landslide area (Gregersen, 1981).In 2009, NGI carried out a site investigation westwards of the Rissa landslide, as part of a new road development planned by the Norwegian Public Road Administration (NGI, 2009).The investigation consisted of 35 rotary pressure soundings, 13 piston samples (54 and 72 mm), 20 piezocone tests (CPTU).Six piezometers were installed.The undrained shear strength of the undisturbed was measured in the laboratory with the fall cone and anisotropically-consolidated triaxial.Oedemeter tests were also conducted.The geotechnical laboratory tests followed state-of-the-art practice described in Sandbaekken et al. (1986) and Berre (1998).
At the site, the soil profile typically consisted of a 4-m thick crust and 11 m of sensitive clay overlying non-sensitive clay.Based on the 2009 investigations, the sensitive Rissa clay had a water content of 35%, plasticity index (PI) less than 10% and liquidity index (LI) above 1.5.The oedometer and piezocone tests suggested a slightly aged clay, with apparent overconsolidation ratio of 1.4.The sensitivity varied from 15 to 100 between 4 and 15 m depth, based on fall cone tests in the laboratory.
Fig. 8 shows the soil layering, the result of a typical piezocone test and the shear strength values measured in the laboratory and interpreted from the piezocone test.The figure also illustrates the profile of undrained shear strength for a normally consolidated clay with normalized shear strength ratios of s u /σ ′ v0 (undrained shear strength s u normalized with effective vertical stress in situ σ ′ v0 ) between 0.2 and 0.3.The value of 0.3 is typical for the undrained shear strength in triaxial compression (s uC ) of a normally consolidated Norwegian clay.The undrained shear strength to use for the sliding in Rissa on a quasihorizontal plane can be assimilated to the strength in simple shear with s u /σ ′ v0 ratio (normally consolidated) between 0.2 and 0.3.For a low plasticity clay like the Rissa clay, the anisotropy is expected to be high (Ladd et al., 1977;Thakur et al., 2017).For the Rissa clay, the normalized undrained shear strength ratio in simple shear (s uDSS ) for an overconsolidation ratio of 1.0 was taken as 0.24.The normalized undrained shear strength ratio in simple shear will be somewhat higher for an overconsolidation of 1.4.
The undrained shear strength from the piezocone test (s uCPTU ) is for the slightly overconsolidated Rissa clay in situ.This undrained shear strength, corresponding to the shear strength in triaxial compression, was obtained from Paniagua-López et al., 2019: (1) with α = 0.28-0.32,exponent m = 0.7 and k = 0.45.OCR is the over- , where q t is the corrected cone resistance, σ v0 is the total overburden stress and σ ′ v0 is the effective overburden stress.The correlation was established using the well-known SHANSEP relationship first published by Ladd et al. (1977) and the most recent CPTU correlations (Paniagua-López et al., 2019).
The data in Fig. 8 also shows that the fall cone test gives much too low peak shear strengths.However, the measured remoulded shear strength (s ur ) from the fall cone test is believed to be reliable and was between 0.5 and 1 kPa over the entire layer of quick clay.
For the mobility analyses, the shear strength of the clay model needs to be depth-averaged (this is one of the limitations of the model, for the time being).A single value of undrained shear strength needs to be determined over the depth where the sliding took place.The sensitive clay at the location of the landslides extends between depths of 4 and 15 m.
For the remoulded shear strength, a value of s ur of 0.7 kPa was used in the base case analysis.For the peak undrained shear strength, a value of 30 kPa was selected as representative of the peak undrained shear strength in triaxial compression of the sliding mass for the following reasons: (1) based on experience with other landslide in highly sensitive clays in Norway, the failure is believed to have taken place between depths of 5 and 10 m; (2) at depths of 5 to 7 m, the peak undrained shear strength was nearly constant (Fig. 8), with the triaxial compression tests indicating values between 25 and 32 kPa.The samples tested were probably somewhat disturbed; (3) the piezocone tests indicated undrained shear strength values in triaxial compression of 32 and 50 kPa between depths of 5 and 10 m, with a depth-averaged value of 42 kPa in this depth interval; (4) to account for anisotropy under simple shear conditions, a 30% anisotropy reduction was applied (NIFS, 2014), resulting in a depth-average undrained shear strength (s uDSS ) of 30 kPa in the 5 to 10 m depth interval.This peak undrained shear strength value was used as the initial undrained shear strength for all mobility simulations described in this paper.
For the selection of the undrained shear strength values, more weight was given to the results of the more recent (2009) site investigations than the 1978 investigations because of the higher quality of the laboratory tests, samples and the availability of CPTU tests.In Gregersen (1981), a value of undrained shear strength of 20 kPa was mentioned.However, this estimate was based on very little information and on fall cone tests which tend to give too low peak undrained shear strengths.Analyses of the mobility of the Rissa landslide (Liu, 2017;2018) have also been done earlier with peak undrained shear strength of 20 kPa.

Current models to analyse the mobility of landslides
The models to analyse the mobility of landslides fall into two categories: empirical, often statistical, models that rely on past observations of landslides, and analytical models.Hybrid "semi-empirical" numerical models also exist.McDougall (2017) presented an overview of mobility or runout distance estimation models and discussed the challenges associated with the models, including the need for improved guidance in the selection of the input parameters.In many cases, one relies heavily on empiricism because there are no straightforward constitutive laws that can describe adequately the dynamics of landslides in the numerical models (Pastor et al., 2012).Until recently, due to the complexity in describing constitutive laws that can capture the post failure softening, the estimation of landslide runout for quick clays has been based on engineering judgement.

Empirical models
The most practical empirical models are based on geometric correlations and until recently did not apply to landslide in quick clays.Two approaches are: (1) an inverse correlation between landslide volume and angle of reach (the angle of the line connecting the crest of the sliding source with the toe of the deposit) (e.g.Scheidegger 1973;Corominas 1996;Hunter and Fell 2003); and (2) a correlation based on Galileo scaling laws between landslide volume and the area covered by the deposit (e.g.Hungr 1995;Iverson et al., 1998).The empirical methods are simple.Several empirical correlations exist for estimating runout distance.For example, Rickenmann (1999) proposed, based on a worldwide dataset of 232 debris-flow events, that the maximum runout distance (L u , in m) is linked with the vertical drop (H D , in m) and the debris-flow volume (V, in m 3 ).The parameter H D is the vertical distance between the centre of gravity of the soil mass involved in the landslide and the centre of gravity of the slide debris on the downstream side: L u = 1.9 V 0.16 H D 0.83 .Corominas (1996), based on 52 debris flows, debris slides and debris avalanches suggested L u = 1.03 V − 0.105 H D (with a coefficient of determination r 2 = 0.763).Here the dimensions are the same as for Rickenmann (1999) with L u in m, V in m 3 and H D in m.However, due to the complexity of the debris flow process, including debris properties, mechanism of motion, topography condition, the simple empirical relationships based on the historical landslide data lead to considerable uncertainty in the prediction of runout distance.
McDougall (2017), on the other hand, suggested that the empirical (statistical) models can help establish limits of confidence for quantitative risk assessment.

Analytical models
Numerical models can provide more information than empirical models because they can be used to estimate velocity of debris flow, debris thickness and impact pressures, and if relevant, the tsunamigenic hazard if the landslide quickly moves into a body of water.With the newer software, the models usually provide visualization of the flowing process, which is important for validation of the model and understanding the complex dynamic process.
In continuum dynamics, the equations of motion use one of two frames of reference: Eulerian or Lagrangian.A Eulerian reference frame is fixed in space, while a Lagrangian reference frame moves with the Fig. 8. Soil profile, typical piezocone test result and undrained shear strength versus depth from laboratory tests and interpreted piezocone tests (q t : corrected cone resistance; u 2 : measured pore pressure by cone; s u : undrained shear strength; σ ′ 0 : in situ effective overburden stress).
Z. Liu et al. flow.To do the mass balance calculations, a discretization of the equations is usually done with a mesh (structured or unstructured) or with a meshless approach.In a meshless scheme, in lieu of a mesh, balance is determined from the spatial distribution of a number of moving reference masses (known as particles).
The majority of the numerical models are continuum models based on hydrodynamic modelling methods, including landslide-specific modifications to account for the effects of entrainment, internal stresses and rheology.McDougall (2017) listed the more common numerical models currently available.The list is updated in Table 1.Mathematical and numerical models of flow slides can be classified on two main groups: 1) Full 3D formulations where the equations of balance of mass and linear momentum are solved on a domain which changes with time.
The position of the interface between the soil and the air is tracked using different algorithms or functions.2) If some assumptions are made about the vertical structure of the flow, it is possible to integrate the balance equations on depth, arriving to the so-called "depth-integrated'' equations or "shallow water equations'' in coastal and hydraulic engineering.
Most of the continuum models are based on depth-averaged shallow flow equations adapted to simulate the flow of earth materials.Different computational methods solve the equations of motion.In the past years, there has been increasing interest in Riemann solvers, where a Riemann problem is solved at the boundaries between cells or elements at every time step.

Modelling of landslide runout in sensitive clays
Empirical models: The post-failure mobility of quick clay landslides is complex.L'Heureux (2012) studied the morphological and geotechnical properties of 30 Norwegian quick clay slides and their runout distance.Key factors controlling the slide evolution were the potential energy available for remoulding, the topography of the release area and of the outflow gate, the relative thickness of sensitive clay and overlying non-sensitive soil and the rheological properties of the clay (e.g.Thakur and Degago, 2013).Tavenas (1984) concluded that large landslides in eastern Canada occurred for liquidity indices higher than 1.2 or remoulded shear strengths less than 1 kPa.L'Heureux (2012) reconfirmed this for the 30 historical quick clay landslides in Norway.The mobility of the Norwegian landslides increased with the volume of sediment per unit width.For a given volume, the runout distances in Norwegian landslides appeared longer than those observed in eastern Canada.Locat et al. (2008), L'Heureux (2012) and NIFS (2013) proposed correlations to estimate the upper bound of the runout distance for Canadian and Norwegian sensitive clays, using also the average width of the landslide, W ave (in the formulations, all quantities are in m or m 3 ): For Canadian sensitive clays, the runout distance was: L u = 1.3 (V/ W ave ) 0.73 For Norwegian sensitive clays, the runout distance was: L u = 9 (V/ W ave ) 0.73   At the same time, Locat et al. (2008) suggested a maximum value for the runout distance for Canadian landslides of L u = 8.8 L 0.8 .Strand et al. (2017), based on 51 Norwegian historical landslides, suggested an empirical model to assess the runout distance and surficial extent of landslides in sensitive clays (NIFS, 2016;Dolva and Petkovic, 2017).Retrogression distance and the total volume of landslides were found to be positively correlated with runout distance (Fig. 9).The retrogression distance L and 'runout' distance L u depend on the landslide type (rotational slide, flake slide or flow slides).In Norway, the retrogression distance L, when estimated empirically, is taken as 15H to

20H
, where H is the slope height.
Semi-analytical methods: Turmel et al. (2017a, b;2020) proposed a new method using the 'destructuration index' concept to estimate the energy available as kinetic energy for runout.The approach to predict landslide runout flow slides in sensitive clays has three components: (1) develop an energy reduction factor specifically for the location under study; (2) select the geotechnical and rheological data for the site where a prediction has to be made; and (3) do a numerical 3D modelling of the landslide runout.Most studies on landslide mobility use apparent rheological properties to describe the behaviour of the landslides.This works for back-analyses, but is of little help for the mapping of present or future landslide hazard and risk.Turmel and colleagues advocated the use of rheological data acquired with a rheometer to predict in advance the post-failure behaviour of a flow slide in sensitive clays.The rheological parameters, measured in the laboratory, are used to reproduce the runout distance in landslides that have already occurred.To ensure a best fit between measured and calculated runout characteristics, the energy present in the system is lowered by the use of an "Energy Reduction Factor" specific for each clay.
The numerical model used for the analysis is the multi-phase Inter-FOAM module of the OpenFOAM software 2 .The InterFOAM module tracks the interface between the different fluids.The solver uses the volume of fluid method to resolve Navier-Stokes equations over a finitevolume mesh.Considering the volume to be composed of only air and remoulded clay, an α value, corresponding to the proportion of the clay present in each element, will then be defined: a value α = 0 is given for  2 OpenFOAM (open Field And Manipulation) is a set of C++ modules used to build solvers to simulate specific problems in engineering mechanics (Weller et al., 1998).Multi-phase InterFOAM is a solver, in three dimensions, for two or more incompressible fluids.
Z. Liu et al. the elements composed of air, α = 1 for the elements composed of remolded clay, and an α value between 0 and 1 for elements located at the interface.The multi-phase InterFOAM solves equivalent mass conservation equation for flows involving one fluid.Turmel and colleagues showed that the consideration of the energy reduction factor made it possible to use site-specific geotechnical and rheological data rather than back-calculated values providing apparent rheological parameters.They also suggested that this new approach improved the accuracy of hazard zonation for flow slides in sensitive clay and provides realistic predictions of the runout behaviour.
Analytical methods: Sensitive clays, in particular, those with high sensitivity or quick clays evolve from a solid state to a viscous liquid during remoulding.The transition from a solid-to a fluid-like behaviour leads to extreme material deformation and rapidly changing free surface.Computational modelling of the post-failure stage of landslides in a sensitive clay thus requires not only an appropriate constitutive model to describe both the solid and fluid behaviour but also a robust numerical technique that can handle dramatic changes in the geometry.
Zhang, Sloan and Oñate (2018) gave an excellent review of the more recent models.A recent contribution was made by Wang et al. (2013) who used a large deformation finite element method with mesh regeneration to study the sliding of a marine sensitive clay deposit on a rigid surface.Dey et al. (2015;2016) implemented a strain-softening model into the coupled Eulerian-Lagrangian (CEL) approach in ABAQUS, and applied it to investigate the failure development in a deposit with a thin layer of sensitive clays.Wang et al. (2016) proposed an implicit material point formulation for modelling the retrogressive failure of a sensitive clay slope.In these studies, rate-independent continuum models were adopted.The neglect of the rheological property of sensitive clays, i.e., the viscosity, affects significantly the prediction of landslides in sensitive clays, particularly for clays with high sensitivity that behave more like liquids after they are fully remoulded.A Lagrangian formulation of elasto-visco-plasticity was proposed in Zhang et al. (2017) for analysing large deformation problems related to sensitive clays.Zhang et al. used an elasto-visco-plastic model, which is a mixture of Tresca model with strain weakening and the classical Bingham model, to describe the behaviour of sensitive clays.To handle mesh distortion and free-surface changes, the model was implemented into the framework of the particle finite element method (Oñate et al., 2011;Zhang et al., 2013).This approach was capable of simulating the progressive failure of sensitive clays (Zhang et al., 2017).Zhang, Sloan and Oñate, (2018) adopted the computational framework of Zhang et al. (2017) for a numerical investigation of the retrogressive landslide in sensitive clays with emphasis on the mechanism of multiple retrogressive failure modes, the kinematics of the sliding mass, and the resulting deposition (i.e., run-out distance and retrogression distance).The role of the clay sensitivity and viscosity on the failure and run-out distance was investigated.Locat et al. (2013;2015) modelled the progressive failure of the Sainte-Monique landslide in quick clay using the BIFURC program and a strain softening model (Jostad and Andresen, 2002), and obtained the Z. Liu et al. evolution of shear strength reduction in the horizontal shear band and the retrogression distance.Tran and Solowski (2019) investigated the retrogressive post-failure behaviour using the Material Point Method.The initial failure was introduced by an artificial steep cut, and the postfailure mechanism was found to be influenced by the uniformly oriented grid net whereby they recommended the anti-locking technique and mesh refinement to improve the issue.Shan et al. (2020) revisited the Sainte-Monique landslide, using a large deformation finite element (LDFE) method, with more general initiation history and advanced remeshing techniques, aiming at understanding the retrogressive failure mechanism, kinematics and controlling factors of the landslide.Issler et al. (2012;2014) did back-analyses of three Norwegian landslides in quick clay [the Byneset landslide in 2012 (0.3 × 10 6 m 3 ), the Finneidfjord landslide in 1996 (10 6 m 3 ) and the Rissa landslide in 1978 (5-6 × 10 6 m 3 )] with three existing numerical flow models to compare how well they performed with quick clay.The three models were: (i) the quasi-2D visco-plastic model BING (Imran et al., 2001a), (ii) the DAN3D-2009 version (McDougall andHungr, 2004) and (iii) MassMov2D (Begueria et al., 2009).The latter two use quasi-3D (visco-) plastic bed-friction laws.Buoyancy effects were considered by scaling the digital elevation model in the vertical dimension for DAN3D and MassMov2D and for the BING analysis in the subaqueous part of the Rissa slide path.Issler et al. (2012) concluded that: (1) the BING model could reproduce the observed runout distances of the Rissa slide only when using a remoulded yield strength much larger than the value measured in the laboratory and the velocities tended to be too high; (2) the DAN3D analyses predicted well the Byneset runout distance and the run-up height in adjoining rivers, but overpredicted runout distance and velocity for the Finneidfjord case, and the analyses did not match the observed landslide shape; and (3) the MassMov2D analyses correctly predicted slide mass running upstream in tributaries at Byneset, but overestimated the runout distance when the measured remoulded shear strength was used.It also predicted no deposits in the upper reaches of the slide path, whereas such had been observed.

Modelling of the Rissa landslide and tsunami run-up
In this paper, the model BingClaw was used to predict the mobility of the Rissa landslide and generate the input for the model BoussClaw to predict the run-up of the tsunami wave triggered by the landslide.

Landslide mobility
The BingClaw model is a visco-plastic model with Herschel-Buckley rheology, developed to predict the runout of landslide in sensitive clays (Kim et al., 2019).The model is an extension of the Bing model (Imran et al., 2001a(Imran et al., , 2001b) ) in Eulerian coordinates in two horizontal dimensions.BingClaw has the following capabilities: (1) it includes remoulding and locally varying yield strength, which allows for a gradual mass release; (2) it is a quasi 3D analysis; (3) it can account for hydrodrag; and (4) it can account for added mass, describing the inertia of the surrounding fluid that needs to be accelerated if the slide accelerates.
Rheology law.Locat and Demers (1988) were the first to suggest that the Herschel-Bulkley rheology is suitable for sensitive clays.Fig. 10 (left) illustrates five different fluid rheology law currently in use, and where for the Herschel-Bulkley (No.4 in figure) the shear stress, from a non-zero state, increases non-linearly with increasing shearing rate.The Herschel-Bulkley "shear-thinning" type of behaviour was confirmed for Norwegian sensitive clays, as illustrated in the right of Fig. 10.Grue et al. (2017) showed with viscosity measurements in the laboratory that the shear stress (torque) vs shear rate (rotation speed) exhibited a Herschel-Bulkley behaviour.The shear stress grows sub-linearly with shear rate.This rheology was used in the BingClaw computation of the Rissa landslide.
For simple shear conditions, the Herschel-Bulkley rheological model can be described as (Kim et al., 2019): where γ is strain rate; and γr the reference strain rate at which the viscous contribution to the shear stress τ equals the contribution from the yield strength τ y ; γr can be expressed in terms of the dynamic consistency μ and the flow exponent n: γr = ( τ y /μ ) 1/n .For shear-thinning materials such as clay, the exponent n is between 0 and 1.A value of n = 1 describes a Bingham fluid.Exponent n can be greater than 1 in thickening fluids (Curve 1 on left in Fig. 10).
The reference strain rate, γr , and the flow exponent, n are kept constant during remoulding in the model.This is equivalent to keeping the shear-thinning behaviour unchanged, that the consistency μ decays at the same rate as the yield strength.This assumption on flow exponent and the reference shear rate remaining constant during remoulding was for convenience and in the absence of experimental data.
Fig. 11 illustrates the rheology model in BingClaw with constant velocity profile for a plug overlying a shear layer and a parabolic velocity profile for the shear layer.A Herschel-Bulkley material behaves as a solid where the shear stress is below the yield stress (or peak undrained shear strength) and flows as a power law fluid past this threshold.A feature of the free-surface, gravity-driven flow of a Herschel-Bulkley fluid is the emergence of regions of "plug flow" where the shear stress is below τ y (Fig. 11).Appendix A describes in more detail the model and the governing equations.
In the landslide simulations, the mass balance was integrated over the flow depth (Eq.(A1), Appendix A) and two separate momentum balance equations were integrated, one over the plug (Eq.( A2)) and one over the shear layer (Eq.( A3)) (Huang and Garcia (1997;1998); Imran et al., 2001a).The shear layer and plug layer thicknesses are computed from the mass and momentum equations Remoulding.De Blasio et al. (2005) proposed to approximate the remoulding process by reducing the yield strength (peak undrained shear strength) as a function of accumulated shear deformation: where τ y,0 and τ y , ∞ are the initial and remoulded yield stress (peak and remoulded undrained shear strength), γ is the accumulated shear deformation at the sliding bed and Γ is a dimensionless coefficient describing the rate of remoulding.The change in τ y , occurring as the clays remoulds, is assumed to depend on the accumulated shear deformation.The remoulding function in Eq. ( 3) was chosen because it was one of the simplest one-parameter functions with monotonic decrease (De Blasio et al., 2005).Fig. 12 illustrates the spatially-averaged strength of the Rissa clay as a function of the accumulated shear strain for four values of the remoulding factor Γ. The remoulding process is modelled to nonreversible, and the average yield stress decreases as deformation increases.Small values of the remoulding Γ imply that large accumulated shear strain γ is needed for remoulding.Kim et al. (2019) pointed out that the early phase of the sliding process and tsunami genesis is sensitive to the value of Γ. Kim et al. (2019) suggested that including remoulding in the model allowed for a gradual mass release, which can mimic the sliding behaviour, as captured by more sophisticated models such as in Gauer et al. (2005).A gradual release does not occur in simpler visco-plastic models without remoulding such as Bing (Imran et al., 2001a;b), Volcflow (Kelfoun et al., 2010), and Geoflow-SPH (Pastor et al., 2009).In BingClaw, large parts of the release area with gentle slope and small earth-pressure gradients will be initially stable for high values of τ y,0.If a slope failure and remoulding occur in one location, a new small area may become unstable and start to slide.BingClaw is not able to simulate the breakup into individual blocks.However, simulations of the Storegga slide made by Kim et al. (2019) suggested that the overall deformation pattern and the velocity of the landslide compared well with those predicted by Gauer et al. (2005).
Experimental study of the rate of remoulding coefficient, Γ. Sinding-Larsen (2019) did a study of the rate of remoulding coefficient, Γ, using the data from over 60 locations at the sites of 20 different Norwegian clays, including both sensitive and non-sensitive clays.He used multiple regression analysis to express the rate of remoulding coefficient, Γ, as a function of, among others, index parameters, overconsolidation ratio and undrained shear strength of the clay.Fig. 13 presents Sinding-Larsen's relationships between the rate of remoulding coefficient, Γ, and liquidity index, peak remoulded shear strength and remoulded shear strength.For the Rissa clay, the following characteristics are considered representative: liquidity index >1.5,Even if there is scatter in the data in Fig. 13, the following values of Γ are suggested for the Rissa clay: based on liquidity index and overconsolidation ratio, Γ could be between 0.06 and 0.08, and perhaps up to 0.1; based on the peak and remoulded shear strength, Γ could be between 0.02 and 0.09 (neglecting the four higher values of Γ on Fig. 13).A best estimate could perhaps be 0.08.A Γ value of 0.08 was used as a best estimate in the analyses, but parametric analyses will look at values of Γ between 0.02 and 0.1.
In a study of the remoulding due to cone and ball penetration, Einav and Randolph ( 2005) introduced a normalising parameter, γ P 95 (equal to 3/Γ), representing the plastic shear strain to achieve 95% reduction in shear strength.They suggested typical values of γ P 95 of 10 to 50 by considering the local shear strains within an assumed discontinuity layer across a laboratory sample.They suggested a value of Γ in the range of 0.06 to 0.3.For highly sensitive clays like the Rissa clay, it is expected that the Γ-value would place at the lower end of the range suggested by Einav and Randolph.
Numerical implementation.BingClaw combines a finite volume method with a finite difference method.It builds on the GeoClaw variant (Berger et al., 2011) of the Clawpack library for solving conservation law equations (Clawpack Development Team, 2015;Mandli et al., 2016).Clawpack uses the Eulerian approach on structured meshes, combined with shock-capturing finite volume methods and Riemann solvers.The plug layer thickness, shear layer thickness and plug and shear layer velocities are computed by the mass and the momentum conservation equations in Appendix A.
For each time step, the numerical scheme proceeds with one of two alternatives: 1) The earth pressure gradient combined with gravity is compared to the yield strength in each cell where the material is at rest.If the yield strength is larger than the driving forces in a cell, the cell is stable and there is no movement at their interface.2) If at least one of the cells is unstable or in motion (i.e. one of the cells deforms), the equations are solved as follows: a) the set of equations are solved first without friction terms for the time step from t n− 1 to t n.
At each cell interface, a Riemann problem is solved with the wave propagation algorithm of the finite volume method (LeVeque, 2002) and obtain the predictor step at time t n ; and (b) the friction forces are then applied using a Godunov fractional step method, as described in Kim et al., 2019.

Tsunami wave generation and run-up
For modelling of tsunamis within the class of so-called "long waves", the pressure is assumed to be hydrostatic and standard shallow water models are appropriate.However, for "shorter waves" and wave propagation over longer distances, it is necessary to use dispersive tsunami models where the wave propagation speed depends on the wave spectrum (longer waves propagate faster than shorter ones).Since landslidegenerated tsunamis often include shorter wavelengths, especially during tsunami generation, hydrostatic wave models can overestimate the wave height, while dispersive models will be more accurate (Glimsdal et al., 2013).For inundation modelling, non-linear effects need to be included (e.g.steepening of wave fronts leading to wave breaking).
For the simulations of the tsunami triggered by the Rissa landslide, the Boussinesq solver called BoussClaw (Kim et al., 2017) was used.The solver models non-linear dispersive wave propagation, also accounting for inundation on dry land.The solver is an extension of the GeoClaw package.The wave dispersion in the tsunami propagation becomes important, Boussinesq-type models are preferred to nonlinear shallow water model.The BoussClaw model implements a hybrid of finite volume and finite difference methods to solve Boussinesq's equations, based on depth-averaged velocity and including enhanced dispersion properties.The BoussClaw is similar to the often used general purpose models such as Funwave-TVD and Coulwave-TVD (e.g., Kim et al., 2009;Shi et al., 2012), but is based on a different and simpler set of governing equations, as well as a slightly different numerical scheme.The Bouss-Claw solver is not as vulnerable to instabilities for nonlinearities in shallow water as other fully nonlinear Boussinesq models.More information on the tsunami model is available in Appendix A.
The governing equations for the BoussClaw solver is modified from the equations in Schäffer and Madsen, 1995 which reads: where f D is a Manning friction term.H(x,t) and u(x,t) are the total flow depth and depth averaged velocity of the water respectively, h(x) is the still water depth, η(x, t)is the surface elevation, and therefore H(x, t) = h(x) +η(x, t).Further, t is time, x is the spatial dimension, g is acceleration of gravity, and B is a dispersion parameter (typically 1/15 for dispersion relation similar to linear potential theory).The operator D is defined in terms of the dummy variable w according to: The differential equations include first (subscript "x" for spatial or "t" for time) and second order derivatives (subscript "xx").
For landslide generated waves, the time evolution of the landslide debris thickness is added to the still water depth, and the total flow depth is then adjusted to H(x, t) = h(x, t) + η(x, t).For the Rissa case, the BoussClaw model was fed with the computed landslide debris thickness as computed by BingClaw.
To validate the BoussClaw solver, Kim et al. ( 2017) compared the results with analytic solutions, solutions obtained by pre-existing models, and laboratory experiments.Even though the equations of BoussClaw are not fully nonlinear, they performed far better than standard Boussinesq equations with only linear dispersion terms.Kim et al. (2017) demonstrated that the point of wave breaking may be wrongly identified in many of the earlier Boussinesq models.They also found a good pre-breaking performance with the BoussClaw Boussinesq equations where only some nonlinearity was retained in the dispersion.They concluded that retaining full nonlinearity in Boussinesq shoaling/ run-up models could be relaxed, which helps reducing the stability problems experienced with fully nonlinear Boussinesq equations (Løvholt et al., 2013).
As mentioned, simulations of the Storegga slide made by Kim et al. (2019) suggested that the overall deformation pattern and the velocity of the landslide compared well with those predicted by Gauer et al. (2005).Deformation pattern (thickness and extent of debris) and the velocity of the landslide are of most significance for the prediction of a tsunami triggered by a landslide in a body of water.
The input to the tsunami simulations are the bathymetry and the time-dependent seafloor variations computed by the landslide modelled in BingClaw.

Input data
The topography of the area and of the different phases of the Rissa landslide was reconstructed based on the study by L' Heureux et al. (2012).To simulate the landslide propagation, a digital terrain model with 10 m square cells and high resolution bathymetry acquired in using a 250 kHz interferometric sonar system (GeoAcoustics) was used.A maximum release depth in Lake Botnen of 18 m and a released volume of about 4 × 10 6 m 3 were considered in the numerical back-analysis.
Fig. 4 illustrated the sequence of the Rissa landslide and Fig.
showed the extent of the landslide deposits at the bottom of Lake Botnen at the end of the landslide.As mentioned earlier, the landslide occurred in two stages.The first stage was relatively small, and the sliding was a flow-type slide with successive retrogression.The second stage, which is the main landslide, was a combination of flow-and flake-type landslides.The release of two large flakes, Flakes A and B, generated the largest waves.In the present study, only the second stage that triggered the tsunami, was simulated.
Fig. 14 presents the spatial distribution and debris thickness from the source area, as modelled by BingClaw.The volume of material used as input was estimated by interpreting the topography before and postfailure.The figure also superposes, in a red dotted contour, the extent of the landslide debris at the end of the sliding based on the interpretation of seismic images.
The best estimate for the three BingClaw quick clay parameters were selected in the previous sections as: Peak shear strength, τ y,0 = 30 kPa Remoulded shear strength, τ y,∞ = 0.7 kPa Rate of remoulding coefficient, Γ = 0.08 However, a few parametric analyses were also done to illustrate the effect of changing the remoulded shear strength and the rate of remoulding coefficient Γ .Table 2 lists the analyses discussed in this Fig.14.Spatial distribution and debris thickness from the source area, as modelled by BingClaw (volume of material used as input was estimated by interpreting the topography before and post-failure); the red dotted curve is the mapped landslide runout, the colours illustrate the debris thickness.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)paper.Analyses with other values of peak shear strength and rate of remoulding coefficients were presented earlier (Liu, 2017;2018).
Case 1, called the base case, with the parameter combination (τ y,0 , τ y,∞ , Γ) of (30 kPa, 0.7 kPa, 0.08) gave a very good match of the extent of the landslide, the velocity of the moving mass and a reasonable match for the debris thickness at the end of the landslide.

Runout distance:
The results of the landslide simulation in terms of runout distance are presented for Case 1 in Fig. 15 at times 1, 3, 10 and 20 min after the start of the second stage of the Rissa landslide.The results show that the runout was completed after 20 min.There is a very good agreement between the simulated runout distance and the observed runout distance.There is also a general good agreement in the shape of the areal distribution of the final debris.
Velocity: Fig. 16 presents the evolution of the velocity field with time for the Case 1 simulation.The contour of the final landslide extent is also shown in the figure in the red dotted curve.A maximum velocity of 10 m/s was reached about 1 min after the start of the second stage of the Rissa landslide.The moving mass reached a velocity of about 12 m/s about 3 min after the start of the second stage.The landslide movement stopped at about 19 min.The velocity was very small at that time, as illustrated in the Fig. 16d.
Debris thickness: Fig. 17 compares the thickness of the debris along the cross-section shown on Fig. 17.For comparison, the seismic line of the lake bottom for the same cross-section, obtained in 2010, is used.A distinct reflector shows the buried lake bottom of 1978 below the debris material.The lower bound of the debris computed by BingClaw is shown with the red line on Fig. 17.The predicted debris thickness extends from that red line up to the lake bottom observed from the geophysical surveys.The debris were spread from distance 500 m to 1700 m along the track on the horizontal axis on Fig. 17.
There is a good agreement on the debris thickness from 450 to 600 m and from 1000 to 1700 m on the seismic line.The BingClaw prediction is less satisfactory from distance 600 to 1000 along the track, where the debris thickness is underpredicted by 30 to 50%.The discrepancy can be due to several factors, including the facts that geophysical reflectors do not provide exact measures of depth, and that there are inaccuracies in the estimated initial debris volume.One other important reason for the discrepancy is that the landslide occurred in 1978 and the seismic data   are from 2010, nearly 30 years after the slide.Settlement of the lake bottom must have occurred as the soft seafloor consolidated under the accumulated debris.The settlement could well be 1 to 2 m over 30 years under the weight of the average 6 to 8 m of debris measured (the debris was up to 10 m in the middle of the landslide).The consolidation of the clay is also probably the cause for the concave shape of the lake bottom in 2010, as shown by the reflector.The modelling in BingClaw is for the situation before any settlement under the debris had occurred.Another source of discrepancy between calculated and observed debris thickness is that the input bathymetry for the BingClaw landslide mobility analyses was the bathymetry from the multibeam measurements in 2010 and not the bathymetry of 1978 before the landslide occurred.The latter is not known, but one can safely assume that the bottom of Lake Botnen in the debris area was lower than the sea bottom profile used as input in BingClaw.Since the lake bottom is not absolutely flat, this will lead in local differences between the predicted and measured landslide debris thicknesses.

Effect of rate of remoulding coefficient on landslide mobility
The rate of remoulding coefficient Γ has a strong influence on the rate of remoulding of the sensitive clay and therefore the runout distance and the landslide velocity.Fig. 18 presents the maximum velocity of the Rissa landslide as a function of time, as computed by BingClaw for four rates of remoulding coefficient, Γ, 0.08, 0.1, 0.05 and 0.02 (Cases 1 to 4).The maximum velocity is the maximum value of the landslide movement at a given time, over the entire flow domain.The velocity for each calculation starts from 0. The velocity reaches a maximum velocity of 16 m/s after just a few seconds because Stage 1 of the Rissa landslide was set, in the analyses, as the initiation of the landslide flow subject to remoulding.In Fig. 18, the curves for the four calculations are superposed for the first few seconds.
Fig. 18 shows that the landslide movement stops completely after about 20 min, but very little movement occurs after 15 min, for Γ ≥ 0.08 (Cases 1 and 2).The landslide debris are still mobilizing after 20 min for Γ = 0.02 (Case 4).The velocity is zero at 20 min for Γ = 0.05.Liu et al. (2018) presented additional studies on the runout distance as a function of the variation of the rate of remoulding coefficient Γ.
Fig. 19 presents the landslide propagation for the cases where the remoulding factor Γ is 0.1 and 0.02.Case 2 has the parameter combination (τ y,0 , τ y,∞ , Γ) set to (30 kPa, 0.7 kPa, 0.1) and Case 4 the combination (τ y,0 , τ y,∞ , Γ) set to (30 kPa, 0.7 kPa, 0.02).At 20 min, the extent of the landslide debris in Case 1 is slightly off the observations.For Case 4 with significantly lower Γ, and thus less remoulding for a given shear rate (as illustrated in Fig. 12), only one-half of the total landslide debris are transported to the bottom of Lake Botnen at time 20 min.

Effect of remoulded shear strength, τ y,∞ , on landslide mobility
The remoulded strength impacts the runout distance.Fig. 20 presents the extent of the landslide debris for remoulded shear strength, τ y,∞ , of kPa and 0.5 kPa (Cases 5 and 6).Neither predict the observed landslide extent as well as the analysis with remoulded shear strength, τ y,∞ , of 0.7 kPa.For a τ y,∞ of 1 kPa, the landslide extent was too short.The movement stopped completely at about 16 min: For a τ y,∞ of 0.5 kPa, the extent of the debris was slightly overpredicted, but still reasonably good.The movement stopped completely at about 22 min.

Input data
The high resolution bathymetry and topography obtained from the geophysical surveys in 2010 and the evolution of the landslide debris thickness with time from the BingClaw analyses of the Rissa landslide were the input of the Rissa tsunami analyses with the BoussClaw software.Two analyses were done for the Cases 1 and 5 in the landslide runout analysis (Table 2).The simulations were performed on grids with 10 m resolution.Each simulation was run for 1200 s (20 min) to ensure that maximum run-up values were calculated also in the northern part of the lake.

Simulation results
Fig. 21 presents the map of the predicted maximum wave heights and maximum run-up heights for the base case (Case 1).The predicted runup heights and can be compared the run-up height measured in and shown on the left side of Fig. 21.Table 3 compares the numerical values of the predicted run-up heights around Lake Botnen due to the Rissa landslide at the location of the measured run-up heights in 1978.The agreement is reasonably good.The back-calculation of the tsunami run-up heights using the landslide simulation of the base case (Case giving the best match of the landslide observations) as input, gives results that are close to the observed run-up heights, except for a few locations, most of them close to the area where the landslide started.
The changes in the predicted run-up heights as the remoulded shear strength is increased to 1 kPa are not very important, except near the locations where the landslide started.

Discussion
The results of the simulations of the landslide runout and the run-up heights of the tsunami following the 1978 Rissa landslide are encouraging.More work and further validations of the models and the parameters with other landslides are required, although no landslide provide as much information as what has been gathered over the years for the Rissa landslide.

Uncertainties
The new models gave a reasonably good prediction of the runout distance and maximum velocity over the flow range for the Rissa landslide, and the run-up heights of the tsunami on the shore of Lake Botnen.The debris thickness was underpredicted in some areas, but there are sources of error, including the settlement of the bottom of Lake Botnen that has occurred under the weight of the landslide debris since and could not be accounted for, and the fact that the input bathymetry for the numerical analyses is not exactly the same at that in 1978 before the landslide occurred.
There are uncertainties in the analysis, including the input parameters of the sensitive clay, seismic observations, measurements in 1978, processes that may have occurred between 1978 and 2010, the most appropriate rheology law, the assumptions and approximations in the BingClaw and BoussClaw numerical analyses and the value of the rate of remoulding coefficient Γ.Nevertheless, the comparisons are good.Liu et al. (2017) did Monte Carlo analyses with BingClaw to estimate the uncertainty in the predicted mobility of the Rissa landslide.The input parameters to the BingClaw model were taken as random variables in these earlier Monte Carlo analyses, while the initial bathymetry was left as a deterministic parameter.The peak undrained shear strength of the clay had a coefficient of variation (CoV, equal to the ratio of the standard deviation to the mean) of 17%, the remoulded shear strength a CoV of 10% and the rate of remoulding coefficient Γ a CoV of 50%.The Monte Carlo analyses, after 1000 simulations, indicated the following variation in the response calculated by BingClaw: a CoV of 8% on the runout distance, 8% on the maximum velocity and 9% on the average debris thickness.

Improvement to the models
The new numerical model to predict landslide mobility in sensitive clays is still under development.In addition, it is a challenge to find the most representative soil parameters to use in the numerical analyses, especially for the selection of the rate of remoulding coefficient, Γ.
Three improvements are suggested at this time: (1) the model allows today only one value for soil parameters τ y,0 , τ y,∞ , and Γ of the clay.It would be important to develop a formulation to enable varying the shear strength in the clay, by including two of three layers of sensitive clay in the model; (2) the modelling would also benefit from including a layer of non-sensitive topsoil riding above the sliding quick clay, as a layer of desiccated crust or sand overlying the sensitive clay is often encountered in Norway; and (3) although some work have been done to quantify the value of the rate of remoulding coefficient Γ, more work should be done to find an easy way to select this parameter for different clays.The parametric analyses in this paper and in Liu et al. (2017;2018) illustrated that a realistic description of the remoulding rate coefficient Γ is indispensable for predicting the landslide mobility of sensitive clays.
It is also envisioned that the BingClaw analysis model could also be modified and used for landslides in non-sensitive materials.This needs, however, to be further studied before recommendations can be made.

Implementation in practice
In a design situation, it is recommended to do Monte Carlo analysis to set bounds for the runout distance and debris thickness that could be expected.In design, one could then consider the computed mean runout distance ± one or two standard deviations around the predicted mean when doing a consequence analysis.The analyses in this paper look at only the hazard aspect of risk (landslide hazard and tsunami hazard, and the tsunami is a consequence of the landslide).For a complete risk mapping and risk assessment, one should include the consequences.In practice, it is recommended to combine the BingClaw analyses with both Monte Carlo simulations and a GIS application, and to develop hazard maps.
With respect to wave run-up height due to landslides occurring by a large body of water, it is very important to do a tsunami-type analysis as done herein with BoussClaw, to enable the prediction of the run-up heights at a distance, as the wave usually strikes as a surprise, causes significant damage and poses a risk to life.For the Rissa landslide, the  Hazard and risk maps should also be adjusted with time, because the risk is dynamic as the outer conditions change due to urbanisation, climate change, building and construction, new transportation corridors, etc. Methods to map this changing risk should be developed.With the increased availability of remote sensing, satellite images, and GIS applications, updating of risk assessments (either qualitative or quantitative) should become a requirement to prevent disasters such as the Rissa landslide and the recent Alta and Gjerdrum landslides.

Summary and conclusions
The runout distance of landslides is of great significance for the identification of the elements at risk, the consequence of a landslide, and for risk management.This is even more important for landslides in sensitive materials because of their extreme mobility and the difficulty of rescue operations after the landslide has occurred.
The paper presented a new numerical model in the computer software BingClaw and BoussClaw to predict the runout of landslide in quick clays and, in the case of the Rissa landslide on the border of Lake Botnen, the propagation and run-up heights of the generated tsunami across the lake.
The model, although still under development, already provides promising results.The new model gave a good back-calculation of the runout distance, maximum velocity over the flow domain and debris thickness for the Rissa landslide, and good back-calculations of the runup heights of the following tsunami on the shoreline of Lake Botnen.The challenge resides in finding the representative soil parameters to use in the numerical analyses and developing a more realistic modelling BingClaw of the soil layering and variable soil properties.
The calculation of landslide mobility (runout distance, flow velocity and debris thickness) and tsunami run-up heights, even if approximate, is essential for establishing maps of landslide hazard and risk in our communities.

Table 3
Comparison of predicted and measured run-up heights around Lake Botnen due to tsunami wave after Rissa landslide (locations are the same as in NGI (1978) shown in Fig. 6, and move clockwise, starting at the south end of the lake where the Rissa landslide started).where h 0 is a reference depth which is chosen as the maximum equilibrium depth.The numerical model, called BoussClaw, was used.BoussClaw is an extension of GeoClaw (Clawpack Development Team, 2015), and solves the Boussinesq-type equations derived by Schäffer and Madsen (1995).The extended model is formulated in two horizontal directions.The BoussClaw model employs a finite volume technique for the NonLinear Shallow Water (NLSW) part of the equations and a finite difference discretization in fractional steps.The GeoClaw software is designed to solve the nonlinear shallow water equations, and is a part of Clawpack (Clawpack Development Team, 2015) developed by LeVeque (1997), George (2008) and Berger et al. (2011).Schäffer and Madsen (1995) derived Boussinesq-type equations where addition of a higher order O(µ 4 ) term enabled optimization of linear dispersion properties.A value of B 2 = 0 from the Schäffer and Madsen formulation was selected.Section 4.2 in the main text presents the governing equations.The BoussClaw model solves the Boussinesq-type Eqs. ( 4) and (5) (in the main text) numerically with a hybrid combination of the finite volume and finite difference methods.There have been several studies of this type of hybrid schemes, e.g., Tissier et al. (2011;2012); Shi et al. (2012); Dutykh et al. (2013).

A.2.2. Governing equations
To facilitate a fractional step method, as outlined below, the hydrostatic terms of Eq. ( 5) were moved inside the (1 − D) operator, while balancing with extra terms in the function Ψ , to obtain Eqs.(A5) and (A6): where Eq. ( 5) in the main text and Eq.(A5) may be solved by a fractional step method as described in LeVeque ( 2002).First, it is observed that Eq. (A5) may be rewritten as: At the first stage of the hybrid schemes, Hu is integrated over a time step taking into account all hydrostatic terms (those inside the curly brackets in Eq. (A7), and omitting the source terms involving Ψ.When this is combined with the continuity Eq. ( 5), this corresponds to advancing the shallow water equations one time step forward.To achieve this, the high-order finite volume solver for the shallow water equations in GeoClaw is used.The Manning resistance term is also accounted for.In the final stage, the H value is retained, and Hu (essentially the momentum density) is integrated further from the two first stages by solving: Since the differential operator D contains spatial derivatives, a system of differential equations need to be solved for the spatial and time discretization.The second order-centered scheme is used for the spatial discretization, and a four-stage Runge-Kutta method is used for the tine integration.
Illustrations of the approach and verifications can be found in Kim et al. (2017).

Fig. 3 .
Fig. 3. Map of Rissa community and shaded relief image of Lake Botnen, with location of the Stage I landslide and Flakes A and B. The black arrows indicate direction of movement.

Fig. 4 .
Fig. 4. Rissa landslide sequence: a) initial slide; b) end of Stage 1; c) sliding of Flake A; d) after Flake A had slid in lake; e) sliding of Flake B (compared to d), part of the road had slid); f) Final scarp(Gregersen, 1981).

Fig. 5 .
Fig. 5. Shaded relief image of the Rissa landslide deposit.Top: From south with lobes, morphological high, ridge, hummocks and debris; Bottom: From north with houses taken by the landslide visible south of the ridge.

Fig. 7 .
Fig. 7. Destruction of Leira village after the Rissa landslide tsunami.Top and middle: Sawmill from the southwest and north; Bottom: Delimitation of wave run-up in bay south of Halvspannet with penned-in dark line.

Fig. 9 .
Fig. 9. (Top) Empirical model for landslides in sensitive clays based on 51 landslides in Norway: Travel distance vs landslide volume (left) and runout distance vs retrogression distance (right); (Bottom): Notation used in empirical correlations (after Strand et al., 2017).

Fig. 13 .
Fig. 13.Variation of rate of remoulding coefficient, Γ, as a function of liquidity index, peak undrained shear strength and remoulded undrained shear strength in clays in Norway.

Fig. 15 .
Fig. 15.Predicted Rissa landslide propagation versus time, plan view.Case 1 (red dotted curve is the mapped landslide extent).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 16 .Fig. 17 .
Fig. 16.Predicted Rissa landslide velocity versus time, Case 1 (red dotted curve is the mapped landslide extent).Note: scale for debris thickness differs slightly for panel b.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 18 .
Fig. 18.Maximum velocity of the Rissa landslide as a function of time as computed by BingClaw for four rates of remoulding coefficient, Γ (Γ = 0.08 is the base case); at time 0 and a few seconds later, the four curves are superposed.
(continued on next page) Z.Liu et al.

Table 2
Overview of the BingClaw simulations of the Rissa landslides in this paper.