Estimating the volume of the 1978 Rissa quick clay landslide in Central Norway using historical aerial imagery

1 Introduction

Quick clay is a fine-grained sediment with high porosity and water content that is highly sensitive to disturbances such as vibration, loading, or bank erosion and may cause retrogressive landslides [1]. Quick clay is found in Scandinavia and Canada, but it is also found in northern parts of Russia, Alaska, New Zealand, Greenland, and Japan [1,2,3,4]. Landslides in clay have been a major contributor in shaping the landscape following the last major deglaciation [5].

In Norway, highly sensitive clays are found in areas where land has been submerged by isostatic depression due to glaciations, making up approximately 1,000 km² of land in total [6]. As the Scandinavian ice sheet retreated, large amounts of clay were deposited in the transgressing sea. The deposits are characterised by random orientation of particles due to the flocculation in saltwater, upheld by the positive charges in the salt. Following post-glacial isostatic rebound, the clay deposits were elevated above sea level and exposed to leaching of salt from groundwater flow and rainfall, compromising the internal structure of the clays. The weakened material structure makes the clay prone to liquidation following ground disturbance which causes landslides with high velocities and large run-out distances. Such landslides can occur on very gently inclined slopes, the quick clay landslides in Norway were reported on slopes as low as 4–5° [3,7]. Deposits are related to the post-glacial marine limit and are most prominent in the Trøndelag region and south-eastern Norway [8]. Depending on the location, the post-glacial marine limit was up to 220 m higher than the mean sea level in present days in Norway [9].

Major efforts have been made to map and characterise quick clay landslides [8,10]. For example, L’Heureux prepared an inventory of 37 well-documented quick clay landslides in Norway and reported that at least 1,150 people have died as a consequence of landslides in clay (as of 2012) [11]. Approximately 110,000 people live on quick clay in Norway. The occurrence of quick clay landslides larger than 50,000 m³ has nearly doubled over the last two decades [6]. One of the largest and most well-known quick clay landslides occurred on April 29, 1978 near Rissa, a small community close to the city of Trondheim, which killed one person and destroyed several buildings [12,13]. Apart from the Rissa slide, there are several other examples of destructive quick clay landslides in Norway, such as the 1893 Verdal slide [14] killing 116 people, the 1959 Sokkelvik landslide [15], which resulted in nine casualties, destructed farms and infrastructure, and initiated a tsunami wave, or the 2009 Kattmarka slide [16] which destroyed a highway and several buildings. The 2015 Skjeggestad slide resulted in the collapse of the Skjeggestad bridge and caused damages worth several million dollars [3,9]. A quick clay landslide occurred at Kråknes in Alta, northern Norway on June 3, 2020. The event destroyed several houses but did not lead to fatalities. Another retrogressive failure happened two days later and destroyed the old E6 road [17]. An even more recent event took place in the town of Ask, in Gjerdrum municipality, on December 30, 2020. Several buildings were destroyed by the landslide, killing ten people and leading to over 1,000 people being evacuated [6]. It has been described as the most deadly landslide in modern Norwegian history [18].

Due to the liquid nature of quick clay combined with human activities (e.g. agricultural activity), traces of past slides are often difficult to detect or are not visible in the landscape anymore [5]. Thus, the availability of historical remote sensing images is of great value to be able to better map, characterise, and understand past quick clay slides. Interpretation of digital elevation models (DEMs) can complement the analysis of quick clay slides, by either using hillshade models to visualise the terrain or by determining changes in surface elevation over time. DEMs can be derived from a variety of remote sensing sources including aerial or terrestrial laser scanning, stereo satellite imagery, or aerial photography. Aerial images have been used for several decades for preparing inventories of landslides in sensitive clays by visual image interpretation, but historical aerial images have seldom been used for DEM creation to calculate the volume of clay slides. Advanced and efficient image analysis techniques have opened a new era for performing retrospective analysis based on historical images [19].

Norway possesses an extensive archive of historical aerial images dating back to the 1930s. The Norwegian Mapping Authorities (Kartverket) are currently in the process of scanning the archive, with much of the collection already digitised and available at a nominal cost [20]. This opens a wealth of new possibilities for the generation of high-resolution DEMs in order to study changes in the landscape over considerable time periods.

Historical aerial photography is traditionally processed using conventional photogrammetric software such as PCI Geomatica or ERDAS Imagine. Such software packages typically require information on the interior orientation of the images (focal length, radial distortion, and image coordinates of fiducial marks) as well as sufficient tie points and ground control points (GCPs) to calculate the relative and exterior orientations. Generally, at least three GCPs must be located on each image that is processed, with each GCP containing both the image coordinates of the point as well as the X, Y, and Z coordinates from a reference dataset (for example a high-resolution aerial image, topographic map, or ground survey data). As such, GCPs should be located on stable terrain that can be identified in both the imagery being processed and the reference dataset. The location and spatial distribution of GCPs are crucial in solving the exterior orientation, estimating camera positions, and ultimately producing high-quality DEMs. The selection of suitable GCPs, however, can be a problem when working with historical imagery that covers large areas of change, for example glaciers, or landslides, when large portions of one of several images are unstable and therefore unsuitable for collecting GCPs [21].

Structure from Motion Multi-View Stereo (SfM-MVS) photogrammetry is a method for processing photogrammetric data, which compared to conventional photogrammetry requires less user input, less pre-processing such as georeferencing, and is capable of automatically solving the camera positions and geometry. SfM-MVS is typically used to process data from unmanned aerial vehicles (UAVs) where one survey can contain several hundred images [22,23].

SfM-MVS offers an alternative workflow for processing historical imagery [24,25,26,27]. SfM-MVS is more automated when compared to conventional photogrammetry, and while the processing principles remain the same, the sequence of processing is different. The raw aerial photographs are matched together based on the identification of key points within the images. If the same key point is found in several images, it is classed as a tie point and is used to solve the relative orientation. A bundle adjustment can then be performed, before a dense point cloud is generated, which can be gridded and saved as a DEM. An important distinction between the SfM-MVS and conventional photogrammetry methods relates to how GCPs are handled. SfM-MVS only requires a minimum of three GCPs per project as opposed to per image, and the GCPs can be added directly onto individual images or onto the orthorectified mosaic of all images. This reduces the amount of manual work required in processing sets of imagery and also mitigates the aforementioned problem of finding GCPs on images with little stable terrain. SfM-MVS is also less dependent on the initial camera calibration parameters, with the software able to estimate much of the inner orientation.

As such, SfM-MVS has increased the usability of archives of historical imagery [24,25], and has permitted the extraction of DEMs from otherwise complicated aerial surveys [26,28,29]. SfM-MVS has been used to process historical imagery to investigate glacier changes [25,29], measure rockfall volumes [30], study the morphology of volcanoes [31], quantify rates of cliff erosion [32], and assess river flood-plain systems [33,34].

Differencing of multi-temporal DEMs allows changes in surface elevation and volume to be quantified [35,36]. This is useful within many fields of geoscience such as glaciology [37,38,39,40], landslides [41,42,43,44,45], and fluvial and coastal geomorphology [33,46]. DEM-derived volume changes are also useful for calibrating and validating glacier mass-balance models [47], characterising landslide dams [48,49], and modelling landslide-induced tsunamis [50]. As a prerequisite to calculating elevation changes, it is necessary to first co-register the DEMs. This is typically done by examining elevation biases over stable terrain (i.e. terrain not affected by mass movements, glacial, fluvial, or other processes) where one would expect to see no changes in elevation over time. Several methods exist for co-registering DEMs and/or point clouds [51,52,53,54], although a comparison by Paul et al. [55] found that the most computationally effective and accurate method was that by Nuth and Kääb [53]. This method compares slope normalised elevation differences over stable ground with the aspect in order to derive the necessary shifts in X, Y, and Z dimensions. Once a pair of DEMs is co-registered, then changes in surface elevation and volume can be extracted.

In this study, we apply SfM-MVS photogrammetry on two sets of historical aerial photographs over the Rissa quick clay landslide, in order to generate pre- and post-landslide DEMs from 1964 and 1978, respectively. Based on this we derive the first geodetic volume estimate of the 1978 Rissa landslide and assess topographic and planimetric dimensions of the landslide.

2 Study area

The Rissa landslide occurred on April 29, 1978 on the southern end of the Botn inland fjord approximately 25 km northwest of the city of Trondheim in central Norway (Figure 1). Botn is connected to the sea through the river Straumen, which during high tides brings sea water into the lake, making it brackish. The maximum depth of the lake is measured to be ∼36 m [56]. Most of the superficial deposits around the lake are mapped to be marine deposits (mainly clay and silt), and the marine limit is estimated to be 160 m a.s.l. [57].

Figure 1

Location of the study area within central Norway (a) and within Trøndelag county (b). The setting of the Rissa landslide and the spatial extent of the aerial photography (c). Background sources are natural earth data (a), Sentinel-2 2019 mosaic (b), and 2018 Norge i bilder imagery (c). The imagery for (b) and (c) were accessed through ArcGIS online.

The Rissa slide was the largest quick clay landslide to occur in Norway in the last century and covered 0.33 km² with volume estimates ranging from 5 to 6 million m³ of quick clay [12,13]. The slide consisted of two stages. The initial stage has been linked to excavation work on a farm along the shoreline and covered an area of 25–30,000 m². A second larger slide was then triggered which propagated to the mountainside, covering an area of 330,000 m² in total [12,13]. One person was killed during this event, and in total seven farms and five houses were destroyed by the landslide [13]. The landslide also caused a tsunami, which hit the village Leira ∼5 km northeast of the slide area and caused damage to houses and infrastructure [6,12]. Apart from the Rissa landslide area, several stretches around the lake are mapped as potential quick clay zones [58,59]. Today, the landslide scar has been restored to farmland.

3 Data

Two sets of archived aerial photographs from 1964 (pre-landslide – ten images with a scale of 1:15,000) and 1978 (post-landslide – six images with a scale of 1:6,000) were delivered digitally from the Norwegian Mapping Authorities. The 1964 images were part of a routine acquisition, while the 1978 images were taken specifically to capture the landslide. An orthorectified 10 cm resolution Red-Green-Blue (RGB) orthophotomosaic from 2018 was downloaded through the Norge i bilder portal, which was used as a reference dataset along with a 25 cm gridded airborne LiDAR digital terrain model (DTM) from 2018 accessed through the Norwegian Høydedata service. The data used in this study are summarised in Table 1. The 1978 imagery was taken at a lower flying height specifically targeting the landslide. As such it has a finer scale than the 1964 aerial photographs but has a reduced spatial coverage. This meant that the majority of the stable terrain used to co-register the DEMs was found on the west of the landslide, while a smaller area was found on the east.

4 Methods

4.1 Photogrammetric processing, co-registration, and volume estimation

Two sets of historical aerial photographs were ingested into Agisoft MetaShape 1.6.2. The same procedure was followed for both image sets. The camera frame information was masked out for each image, the fiducial markers automatically detected, and the camera calibration (fiducial marker location in image coordinates and the radial lens distortion) were entered into the software. The remainder of the aerial photography processing was mostly automatic with the processing parameters set to “Ultra High” for each step within the software. Points of noticeable image contrast, or “key points” were identified and matched between each set of images to solve the relative orientations. Any point with a reconstruction uncertainty of 15 or a reprojection error greater than 1 pixel was removed. This resulted in a dense cloud with a point density of 5 and 13 points per m² for the 1964 and 1978 datasets, respectively. The resulting point cloud was then meshed before being converted into a gridded DEM. The aerial photographs were then orthorectified and mosaicked together. Using the 2018 aerial photography and LiDAR digital surface model (DSM), several GCPs (30 for the 1964 imagery and 19 for the 1978 imagery) were added onto the orthophotomosaic to solve the exterior orientation. In total, the GCPs had root mean square errors (RMSEs) of 6.9 and 3.3 pixels for the 1964 and 1978 datasets, respectively. The orthomosaicked images were exported at 0.2 m resolution and the DEMs at 0.5 m resolution. The DEMs were filtered to produce DTMs as opposed to DSMs using the DSM2DTM tool within Catalyst Professional setting an object size to 30 m and a gradient threshold of 30°. Bumps and pits were removed from the DTMs using a 7 × 7 moving window filter with a gradient threshold of 5°. Finally, a 5 × 5 median filter was used to smooth each DTM. The filtering was successful at removing buildings and sparse vegetation, although some patches of forest remained.

The two DTMs were co-registered using the method of Nuth and Kääb [53] applied on stable pixels, i.e. those not affected by the landslide, but also excluding water, and areas of poor image contrast where DTM values are most likely erroneous, such as shadows and dense forest. The co-registration was repeated six times, until the improvement of the co-registration was less than 2%, as suggested by Nuth and Kääb [53]. The surface elevation change was smoothed using a 3 × 3 median filter. Lastly, the mean elevation change as well as the volume change within the slide area was calculated.

A combination of the mosaicked images, the hillshade model from the post-landslide DTM, and the surface elevation changes was used to guide the manual digitisation of the landslide area. The top of the scarp was especially noticeable in both the hillshade model and the orthophotomosaic from 1978 (Figure 2). As we did not have sufficient post-landslide elevation data over the portion of the lake affected by the landslide, we set these pixels to the mean elevation of the shoreline (0.6 m a.s.l.).

Figure 2

Orthophotomosaics from 1964 and 1978 for two locations that were affected by the Rissa landslide. The yellow line indicates the landslide outline.

5 Results and discussion

The Rissa landslide was typically 220–300 m wide and extended approximately 1 km in southwest–northeast direction. Several farm buildings have been either displaced or destroyed in the event (Figure 2). One area of the shoreline of the lake retreated ∼150 m inland. In total, the landslide covered an area of 364,940 m² (0.36 km²) as measured on the aerial photography, including the land that fell into the lake in 1978.

The Rissa landslide, as measured by the difference between the two DTMs, had a geodetic volume of 2.53 ± 0.52 × 10⁶ m³. Up to 20 m were vertically excavated at the eastern side of the landslide, falling to 4–6 m at the edge of the lake (Figure 3). The flanks of the landslide were very prominent on all sides, with typically a change in elevation of ∼5 m. The mean slope for the area that failed was 11.9° in 1964. The majority of surface elevation change happened in the centre of the slide area with some areas dropping by 20 m. The surface dropped by an average of 6.93 ± 0.82 m.

Figure 3

(a) Surface elevation changes between 1964 and 1978. Two profiles (X–X′, Y–Y′) that transects the landslide are shown in (b) and (c). Profile Y–Y′ overlays the area of greatest surface elevation changes, where the surface dropped by up to 20 m. The result from transect Z is shown in Figure 5.

In the delineation of the landslide area we included the north western part close to the shoreline of the lake since obvious changes were visible in the 1978 photographs and the corresponding DTM. This area, however, shows a gain in elevation (<3 m), indicating that landslide material was deposited (Figure 3). This is in line with the findings from Gregersen [13], who reported that during the second stage of the sliding, the mass started moving in the direction of the terrain slope, and not towards the lake.

5.1 Uncertainty of DTMs and DTM differencing

Both the standard deviation and mean elevation change over stable terrain can be used as a first estimate of the uncertainty of a DTM; however, it can significantly overestimate the error as spatial correlation is not considered [65,66]. This is especially true in our case, where elevation biases are typically <0.5 m on gently sloping terrain with high image contrast (such as farmland and bare soil), where the majority of the landslide occurred, and elevation biases are larger (5–10 m) over steep terrain, or in areas of low image contrast, such as forests or shadows. For comparisons between the historical DTMs and the LiDAR data, as well as between the two historical DTMs, the standard deviation over stable ground was <2 m (Table 2). Both the historical DTMs are positively skewed against the LiDAR DTM, while the 1964 DTM is strongly negatively skewed relative to the 1978 DTM. E ∆ V is most likely a better estimate of our volume change uncertainty, as it considers the hypsometry of the landslide, i.e. the majority of the surface elevation change occurred on gently sloping terrain with lower elevation biases. Our estimate of the vertical uncertainty (E∆_h) in our surface elevation change calculation is 0.8 m.

Table 2

Metrics for assessing the accuracy of the two DTMs produced as well as the elevation change between them. The accuracies for the 1964 and 1978 DTMs were calculated based on comparison with the 2018 LiDAR DTM on stable terrain. The inter-comparison of the 1964 and 1978 DTMs was based on stable ground

DEM of difference	Mean deviation (m)	SD (m)	Skewness (m)	RMSE (m)	Elevation change uncertainty, E∆h (m)
1964–2018	0.06	1.27	0.87	1.39	—
1978–2018	0.08	1.34	1.60	0.31	—
1964–1978	0.01	1.21	0.47	0.24	0.52

When comparing the DTMs to the 2018 LiDAR DTM, the 1964 and 1978 DTMs have mean elevation bias over stable ground of −0.06 and −0.08 m, with standard deviations of 1.27 and 1.34 m, respectively. The mean elevation changes in stable ground between the two DTMs is even lower (−0.01 m) with a standard deviation of 1.21 m. If the elevation biases are assumed to be distributed normally, then 95% of the stable ground elevation biases are <1.47 m (Figure 4). This indicates that the systematic errors relating to DEM co-registration are minimal in our analysis. A degree of uncertainty in the accuracy values might be associated with the circumstance that the reference LiDAR DTM was acquired several decades after the aerial photographs. Thus, even the area assumed to be stable ground may have undergone minor changes, such as building construction or vegetation change. This most likely explains the skewing of our results and indicates that our accuracy assessment is most likely too conservative. The accuracy and precision of the DTMs based on historical imagery is heavily dependent on the flying height, scanning spatial resolution (dots per inch), and radiometric resolution [67]. In our case, the photographs were digitised professionally by Kartverket and no scanning artefacts were visible, and the flying height allowed sub-metre DTMs to be produced, although in cases of poorer image quality we would not expect such precise results.

Figure 4

Plotted elevation biases calculated by comparison with the 2018 LiDAR DTM over stable ground for the 1964 and 1978 DTMs. The inter-comparison between the 1964 and 1978 DTMs is also shown.

5.2 Comparison with other studies

Our geodetic estimate of the Rissa landslide volume (2.53 ± 0.52 × 10⁶ m³) varies considerably from existing estimates. A volume of 4.4–5.8 × 10⁶ m³ was estimated by L’Heureux et al. [12], which was based on bathymetric data in order to examine the morphology of the landslide deposits and combine it with an estimate of the deposit depth based on seismic data. Gregersen [13] estimated a total landslide volume of 5–6 × 10⁶ m³ but did not specify how this number was determined. The areal extent of the landslide as determined by this study (0.36 km²) is however very similar to that demarked by Gregersen 0.33 km² [13]. A recent study modelled the Rissa landslide based on the post-event topography and bathymetry, and modelled the resulting tsunami [6]. Although the scope of their research was different from ours, part of their analysis involved estimating the thickness of the material that slid into the lake. They showed a maximum thickness of 18 m, which corresponds well with our maximum measured surface elevation change of 20 m.

We digitised the profile and volumetric zones presented by Gregersen [13] in order to compare his results to our own (Figure 5 and Table 3). In both cases it is clear that our estimates are noticably smaller than those reported by Gregersen [13]. When we calculated the mean surface elevation change per volumetric zone reported by Gregersen [13], it is noticable that our results are smaller by 25.05–64.94% (Table 3).

Figure 5

Comparison between (a) the transect after Gregersen [13] across the landslide and (b) the surface elevation change derived from the comparison of the photogrammetric DEMs generated in this study. The house visible on “a” is left unfiltered on “b” for reference. Note that the terrain data shown in (b) is a DSM to include the buildings. The position of the lake in 1978 is shown with a blue dashed line.

Table 3

Comparison between the volumes reported by Gregersen [13] and the geodetic surface elevation change derived from this study

Volumetric range (derived from Gregersen [13] (m)	Mid-point of volumetric range (m)	Mean surface elevation change (this study) (m)	Percentage deviation (%)
0 to −5	−2.5	−0.88	−64.96
−5 to −10	−7.5	−4.71	−37.23
−10 to −15	−12.5	−9.37	−25.05
>(−15)	−15.0	−8.80	−34.70

This leads us to infer that the landslide was correctly delineated, yet the disparity in analysis could stem from the different methods for volume calculation. Our geodetic approach does not account for material that was deposited in the landslide area. We calculated the surface volume change based on the difference between the pre- and post-event DTMs, while the actual slip surface was at least partly covered by the deposited material and thus likely lower than the surface as represented by the post-event DTM.

Moreover, our geodetic volume estimation does not fully account for the parts of the landslide that occurred within Rissa-Botn lake which can partly explain the difference between our results and previous estimates. Our results show that at the northern front of the landslide – at the shore of the lake – there are elevation changes between −7 and −10 m. It can therefore be expected that a portion of the landslide failure occurred underwater.

In order to fully quantify the Rissa landslide, it is necessary to cover both the areas where clay was removed and where it was deposited. As such combined multi-temporal and high-resolution topography and bathymetry datasets would be ideal to get an overall picture of the total volume changes [68]. However, since the landslide entered the lake it is not possible to apply the same method for estimating the volume on land and below the water, i.e. analysing the depletion zone and the accumulation zone in one step. This is an additional reason why a direct comparison of our calculated volume, i.e. the removed material, with the existing volume estimation by L’Heureux et al. [12], i.e. the accumulated material, is difficult, even when accounting for an underestimation of our volume calculation. Furthermore, apart from potential technical and methodological limitations, additional factors such as reworking of the slide material during movement [41,69,70] and erosional and depositional processes that shape the deposited material [6] – the bathymetry data was acquired more than three decades after the landslide – may contribute to differences in volume estimations. We also believe that our remote sensing-based elevation changes could be useful in calibrating models of landslides and their resulting tsunamis (for example, Liu et al. [6]).

6 Conclusion

This study derived a geodetic volume of the 1978 Rissa landslide of 2.53 ± 0.52 × 10⁶ m³ with drops in surface elevation of up to 20 m, and a mean surface lowering of 6.93 ± 0.82 m. The actual landslide volume calculated is underestimated due to the challenges remaining in the computation of such landslide volumes. A large portion of the landslide failure occurred underwater, and additionally, we could not account for material deposited in the landslide area. However, our study provides an updated estimate of the landslide considering the stochastic and systematic uncertainties. Our chosen workflow was largely automated and the DTMs produced were both precise and accurate when compared to a reference dataset, which can help quantify and understand quick clay landslides in high-resolution. As such, we can recommend SfM-MVS for processing large datasets of historical imagery in order to quantify landscape changes over large spatial and temporal scales. Such high-resolution volume estimates of quick clay landslides are necessary to better comprehend their development, and model their impacts, such as the triggering of tsunami waves. Knowledge derived from the analysis of historical aerial images can fill critical information gaps and provide useful information for hazard and risk analysis.