Tectonic deformation within fault damage zones can influence slope stability and landslide failure mechanisms due to rock mass strength effects and the presence of tectonic structures. Here, we used detailed site investigations to evaluate controls on deformation within the Half Moon Bay landslide complex, located ~1 km from the surface trace of the Hope fault in the South Island of New Zealand. During the 2016 Mw 7.8 Kaikōura earthquake, the slope experienced up to ~13 m of displacement and partially transitioned into a rock avalanche (with a volume of ~350,000 m3). Deep-seated deformation of the entire slope predated the 2016 earthquake. Results of geomorphological analysis, field mapping, geophysical surveys, slope displacement, and a 60-m-deep borehole in the incipient portion of the landslide indicated the presence of a subvertical tectonic fabric and intense fracturing and weathering of the rock mass, which gradually decrease with depth. Based on these results, we established a conceptual model wherein the landslide failure mechanism is a combination of flexural toppling along the subvertical structures coupled with joint-step-path sliding along preexisting, closely spaced discontinuities within the graywacke rock mass. Coseismic slope displacements revealed a large area of incipient failure behind the headscarp of the 2016 rock avalanche, which will likely result in further avalanching at the site. This case study demonstrates that inherited tectonic structures (combined with seismicity and weathering in an oversteepened coastal slope) play an important role in the evolution of hillslopes near active faults.

Coseismic landslides are among the most hazardous secondary effects of earthquakes and have caused widespread damage to infrastructure and loss of life (Keefer, 1984; Marano et al., 2009; Fan et al., 2019). Understanding the factors that influence coseismic slope failure is therefore critical for effective risk management.

In tectonically active regions, steep slopes within fault damage zones have previously been identified as particularly susceptible to coseismic failure from a combination of strong ground motion and reduced rock mass strength (Korup, 2004; Tatard and Grasso, 2013; Gallen et al., 2014; Gorum and Carranza, 2015; Fan et al., 2019; Bloom et al., 2021). Amplification of seismic shaking due to topographic effects and impedance contrasts can additionally enhance landslide susceptibility (Ashford et al., 1997; Martino et al., 2006). Furthermore, in some cases, tectonic structures may influence slope failure mechanism and volume as they provide an additional discontinuity set that can affect slope kinematics (Brideau et al., 2009; Stead and Wolter, 2015; Massey et al., 2018b). Large landslide complexes, in particular, have previously been observed in slopes affected by tectonic deformation (e.g., Amann, 2006; Penna et al., 2017; Agliardi et al., 2019; Glueer et al., 2019). Progressive strength degradation due to coseismic perturbations (i.e., repeated brittle rock damage due to seismic shaking) has been recognized as an important factor in slope destabilization (Wolter et al., 2015; Gischig et al., 2016; Agliardi et al., 2019; Massey et al., 2022; Brain et al., 2021); however, direct observations of coseismic failures within large landslide complexes are rare, despite their potential value in improving our understanding of coseismic and interseismic processes affecting slope deformation.

In this study, we present a detailed site characterization of the Half Moon Bay landslide complex, part of which failed as a rock avalanche during the 2016 Mw 7.8 Kaikōura earthquake in New Zealand’s South Island. We combined findings from geomorphological analysis, field mapping, measurements of coseismic (2016) and total (including pre-2016) landslide displacements, geophysical surveys, and a 60-m-deep borehole in the incipient portion of the landslide to develop an understanding of the structure and mechanism of this landslide. Drawing on this unique data set of geotechnical ground investigations, we built an integrated model of landscape and landslide evolution that elucidates the factors controlling the development of the landslide complex and, ultimately, the movement mechanisms leading to parts of the landslide mobilizing into a rock avalanche.

2.1. Regional Tectonics and Geology

The Half Moon Bay landslide is located in the northeastern South Island of New Zealand, ~800 m from the primary surface trace of the Hope fault (Fig. 1). The Hope fault is the southernmost tectonic structure of the Marlborough fault system, which consists of predominantly dextral strike-slip faults that link the Alpine fault transform plate boundary to the south with the west-dipping Hikurangi subduction zone to the north (Fig. 1; Litchfield et al., 2014). Net slip rates in the Marlborough fault system vary between 0.1 and 25 mm/yr on individual faults (Litchfield et al., 2014), while along the Hope fault itself, slip rates vary among different fault segments. Slip rates are as high as 23 ± 2 mm/yr on the Conway segment (Langridge et al., 2003) and decrease to ~2–5 mm/yr on the Mt. Fyffe and Seaward segments as slip is transferred onto the Jordan thrust and other faults in the Seaward Kaikōura Range (Fig. 1; Van Dissen and Yeats, 1991; Wallace et al., 2007; Collett et al., 2019). To the south of the Hope fault, relative plate motion is accommodated in a series of lower-slip-rate (as compared to the Marlborough fault system) strike-slip, reverse, and oblique faults and folds of the North Canterbury tectonic domain (Litchfield et al., 2014).

The geology of the study site is composed of graywacke of the Torlesse Supergroup in the Pahau terrane (Fig. 1C), which is typically characterized by closely jointed sandstone and argillite (Rattenbury et al., 2006). To the west of the Half Moon Bay site, Paleogene and Neogene sedimentary rocks are exposed in the Puhipuhi syncline, the fold axis of which is oriented NE-SW (Rattenbury et al., 2006), while Torlesse basement rocks are deformed in an anticlinal fold between the Puhipuhi syncline to the west and the Papatea fault to the east (Fig. 1C; Van Dissen and Yeats, 1991).

2.2. 2016 Mw 7.8 Kaikōura Earthquake

The 2016 Mw 7.8 Kaikōura earthquake initiated on the Humps fault in the North Canterbury domain (Fig. 1) and propagated to the north, causing surface rupture of more than 20 faults of the North Canterbury domain and the Marlborough fault system (Litchfield et al., 2018). During the earthquake, the “Papatea block” (i.e., the area bounded by the Papatea fault, Hope fault, Upper Kowhai fault, and Jordan thrust; Fig. 1C), in which the study site is situated, was uplifted up by 8 m and displaced to the south by 4–5 m (Hamling et al., 2017). Fault displacement along the Papatea fault (up to 11.8 m) is considered to have been an anelastic response to neighboring fault movement and a localized space problem (Diederichs et al., 2019; Langridge et al., 2022). Rupture along the Hope fault was only observed locally following the 2016 event, with fault displacement in the Seaward segment of the Hope fault reaching a maximum of ~0.3 m, while a maximum fault offset of ~7.5 m was measured along the Jordan thrust (Litchfield et al., 2018). During the 2016 Kaikōura earthquake, recorded peak ground accelerations (PGAs) exceeded 1 g in vertical and/or horizo ntal directions at the strong ground motion stations in Waiau (WTMC, ~75 km from site), Kekerengu (KEKS, ~35 km from site), and Ward (WDFS, ~45 km from site) (Fig. 1; Bradley et al., 2017; Kaiser et al., 2017). Overall, ~30,000 landslides were triggered during the event (Dellow et al., 2017; Massey et al., 2018b, 2020a). Landslide sizes varied from small rockfalls to large rockslides and avalanches (>1 × 106 m3; Massey et al., 2020a). Most landslides triggered by the event (~80%) were classified as soil or rock avalanches that tended to be relatively shallow, while some deeper translational and rotational slides were also observed (Massey et al., 2020b).

2.3. Site Overview

In this study, we focused on the slope comprising the Half Moon Bay landslide complex (Fig. 2). The landslide complex is located on a steep, SE-facing hillslope with ~550 m relief directly above State Highway 1 (SH1) and the Main North Line railway, which is a nationally important transport corridor. To the west of the Half Moon Bay hillslope, the Rakautara Stream is deeply incised, forming the steep western flank of the landslide (Figs. 2 and 3). To the north, the NW-SE–trending ridgeline terminates against a NW-SE–oriented saddle, while ridgelines and small river valleys make up the topography toward the east of the Half Moon Bay hillslope (Figs. 2 and 3). In the southeast, the hillslope ends in an uplifted shore platform, on which the transport corridor is situated (Figs. 2 and 3).

The 2016 Kaikōura earthquake triggered numerous landslides and rockfalls in the Half Moon Bay area, including the rock avalanche in the Half Moon Bay hillslope (Figs. 2 and 3). The extents of the source areas and debris trails of the 2016 coseismic landslides were mapped by Massey et al. (2018b, 2020a); see Figure 3. Slope geomorphology in historical imagery (Retrolens, 2023) shows preexisting slope deformation prior to the 2016 earthquake (Fig. 2). A 1942 aerial photo (the earliest publicly available historical image; Retrolens, 2023) shows evidence of a relict landslide headscarp and deposit as well as slope deformation in the upper portion of the slope (Fig. 2). Little change was observed in the slope between 1942 and March 2016 (Fig. 2). Historical earthquakes that predate aerial imagery are known to have caused strong to severe shaking (i.e., modified Mercalli intensity [MMI] > 6) at the location of the study site, e.g., the 1848 Marlborough earthquake (Mw 7.4) and the 1855 Wairarapa earthquake (Mw 8.2) (Geonet, 2023).

Following the Kaikōura earthquake in November 2016, debris was regularly remobilized by rainfall. Earthworks, visible in November 2019 aerial imagery, were conducted to reprofile the lower portion of the slope above the transport corridor and build debris-flow deflection bunds (Fig. 2).

3.1. Field and Remote Mapping

Morphological mapping was carried out according to the Geological Society of London guidelines (Griffiths, 1982). First, slope morphometry (i.e., breaks and changes in slope) was manually mapped based on a 1 m shaded relief model of the 2016 postevent digital elevation model (DEM) derived from light detection and ranging (LiDAR) data (Massey et al., 2018a). The morphology map, combined with a map of materials observed at the surface, provided the basis for interpreting the morphogenesis of features in the study area. A morphogenetic map (i.e., interpretation of how features formed) was then produced based on interpretation of the morphometric map and field observations.

Engineering geological logging of rock outcrops on and around the Half Moon Bay landslide complex was carried out to characterize the rock mass structure and quality (Fig. 3). Joint orientation, roughness, persistence, and spacing were described using International Society for Rock Mechanics and Rock Engineering (ISRM, 1978) and New Zealand Geotechnical Society (NZGS, 2005) guidelines. The geological strength index (GSI) was estimated according to Hoek and Brown (1997). Weathering grades were described following Brideau et al. (2020) and NZGS (2005).

Desktop mapping of lineaments on the shore platform between Te Ana Pōuri and Ōhau Point was carried out in a geographic information system (GIS) at 1:5000 scale (Fig. 3) due to restricted access to the eastern shore platform because of the Ōhau Point seal colony. We also carried out regional lineament mapping in the Torlesse basement rocks between the Puhipuhi syncline and the Papatea fault (Fig. 1C) to identify the orientation of bedrock structures over a larger scale (1:50,000). All lineament mapping was performed manually using 20 cm orthophotographs and 1 m LiDAR-derived shaded relief models (described in Massey et al., 2018a). While mapping of lineaments from remote-sensing imagery provides information on the strike direction of dominant trends in the bedrock, lineaments cannot be directly assigned to a specific discontinuity type (i.e., bedding, joints, or faults). We used lineament strike directions to compare joint set measurements from outcrop observations to dominant regional bedrock structures and set the results in context with structural geology and tectonics.

3.2. Landslide Displacement Analysis

Landslide displacement analysis was carried out to identify magnitudes and patterns of slope deformation at the site. In the following, we distinguish between “coseismic landslide displacement,” which refers to permanent slope displacement during the 2016 Kaikōura earthquake, and “total landslide displacement,” which encompasses any permanent slope displacement in the deformation history of the slope, including displacement from the 2016 event. In case of any preexisting slope deformation prior to the 2016 event, the total displacement would, therefore, be larger than the coseismic displacement.

3.2.1. Coseismic Landslide Displacement Analysis

The analysis of coseismic displacement followed the methodology described in Singeisen et al. (2022). Pre- and postevent (i.e., 2015 and 2017) digital surface models (DSMs) derived from digital stereopair aerial photographs (Massey et al., 2020b) were used to derive horizontal coseismic displacements using the digital image correlation (DIC) algorithm by Bickel et al. (2018). Technical specifications are provided in Table S1 of the Supplemental Material.1 The vertical component of the coseismic displacement vectors was obtained by sampling pre- and postearthquake DSM elevations at the locations of pre- and postearthquake displacement vector data points (Singeisen et al., 2022). Using the same method, we also assessed postseismic landslide displacement (i.e., DIC between the 2017 and 2019 DSMs).

Areas of high ground disturbance (i.e., ground cracking or debris) can cause noise in the displacement field. This noise was filtered with a spatial mask based on the orientation of displacement vectors and ground disturbance visible in imagery (e.g., the vectors facing upslope or those showing extreme variations in orientation between pixels were filtered out). The horizontal accuracy of the coseismic displacement measurements is assumed to be ~1 m based on the 1 m resolution of the underlying imagery. Reflectors must move at least one pixel for a change to be identified. For three-dimensional (3-D) displacements, vertical uncertainty from the LiDAR DEM (±20 cm) propagated into our estimates.

3.2.2. Estimate of Total Landslide Displacement

To place the coseismic landslide displacement from the 2016 earthquake into the context of the total amount of slope deformation (i.e., relative to the gravitationally undeformed state of the slope), we estimated total slope displacement based on the morphological mapping. We used the postevent 2016 LiDAR-derived DEM (1 m resolution) to estimate the total displacement by measuring the horizontal and vertical separation of the ground surface across morphological features such as scarps and counterscarps along two profiles (Figs. 4 and 5). Schematic cross sections of horizontal and vertical separation measurements for scarps and counterscarps are provided in section 3 of File S1 (see footnote 1). Cumulative total displacements (displacement component along the axis of the profile) were then compared to the magnitude of coseismic landslide vectors that were projected onto the same profile. The total landslide displacement values are estimates with relatively large uncertainties due to (1) the difficulty of accurately measuring subtle morphological features based on a 1 m DEM and (2) the assumption that topographic separation across features, such as scarps and counterscarps, is evidence of past slope displacement. However, such estimates were useful to compare coseismic landslide displacements in 2016 to the long-term pattern of slope deformation.

3.3. Geophysical Surveys

Electrical resistivity tomography (ERT) and seismic refraction tomography (SRT) surveys were carried out in the incipient portion of the landslide complex (location in Fig. 3). ERT surveys have been successfully applied in the past to image fracture density, weathering, and saturation conditions in a slope and have previously been used to reconstruct landslide depth and geometry (Heincke et al., 2010; Perrone et al., 2014). Intact rock is assumed to be resistive, with fracturing causing an increase in bulk porosity. In the absence of water and clay, a porosity increase would lead to an increase in resistivity. However, if water or clay is associated with the fracture zones, they can appear as low-resistivity zones in surveys due to electrical conduction through the pore spaces (Archie, 1942). In this study, an ERT survey was carried out using an Iris Syscal resistivity meter with a 2.5 m electrode spacing along a 240 m profile. Data were collected using Wenner-α and dipole-dipole electrode geometries. The locations of points along each survey line were recorded using a Trimble Geo7x global navigation satellite system (GNSS) receiver, and elevation data were extracted from the 2016 LiDAR DEM (Massey et al., 2018a). We displayed the recorded ERT data as apparent resistivity values in pseudosections for each of the two geometries. The pseudosections show the observed data in a relative position along the profile and at an approximate depth. In this format, it is easy to identify measurements that have anomalously high or low apparent resistivities compared to their neighbors. Anomalous readings were identified either by a large uncertainty in the individual measurement (random noise) or a rapid variation from its surrounding measurements that is not physically realistic for a diffusive phenomenon such as electrical current flow (bias). At shallow depths of investigation, more rapid changes in resistivity are allowed, but at larger depths in the pseudosections, the volume-averaging nature of the ERT measurements means only small variations both laterally and vertically should be present. The anomalous readings were removed from the data set, and it was then processed through a nonlinear, laterally constrained inversion process, which incorporated the topography to produce a depth section that was consistent with the observations. The specific software package used in this study was Res2Dinv (Loke and Dahlin, 2002).

Two overlapping SRT survey lines with 24 geophones at 4 m spacing were carried out at Half Moon Bay because previous studies have demonstrated the utility of SRT surveys in mapping discontinuities and identifying high-fracture-density zones in large landslides (Heincke et al., 2006, 2010; Cody et al., 2019). As fracture zones are characterized by reduced bulk density, these zones can be identified as low-seismic-velocity zones in a seismic survey. A 20 lb (9 kg) sledgehammer was used as a seismic source with a minimum stacking of three shots to increase the signal-to-noise ratio. The geometry of the geophone array and the shot locations are provided in Figure S10. The locations of the geophone array and all source positions were recorded using a Trimble Geo7x GNSS receiver. The seismic refraction data were processed using Reflex W 7.0, and first arrivals were picked manually (Sandmeier, 2012). Tomographic results were obtained through two-dimensional (2-D) inversion of the traveltimes for the first arrivals in Reflex W. Ray-path tracing was performed to understand the ray coverage within the model cells. Cells with low coverage result in poor resolution and were therefore masked.

3.4. Borehole

A vertical borehole with a depth of 59.6 m was drilled on the ridge crest, ~170 m behind the headscarp of the 2016 rock avalanche and ~10 m downslope of a downhill-facing scarp in a location of relatively flat terrain (<10° slope angle; Fig. S2E). The original target depth of ≥100 m could not be reached due to borehole wall instability and continuous loss of drill fluid, which prevented drilling at greater depths. The borehole was drilled using triple tube wireline rotary coring techniques, with a PQ (85 mm) diameter core recovered between 0 and 43.7 m depth and an HQ (63.5 mm) diameter core recovered between 43.7 m and 59.6 m depth. Detailed core logging was carried out off-site following the NZGS guidelines (NZGS, 2005). Optical televiewer and mechanical caliper borehole logging was carried out between 21.15 m and 46.2 m depth in the borehole; the instability of the borehole prevented further geophysical borehole logging. No groundwater was encountered during drilling.

The results from the multimethod approach used to characterize the ground conditions at the site were integrated into a three-dimensional ground model that is provided in the Supplemental Material (see footnote 1). This model can be used to help visualize the results described in sections 4.1–4.4.

4.1. Morphological Maps of the Landslide Complex

The morphological maps (Fig. 4) cover the extent of the main slope deformation observed, and the term “Half Moon Bay hillslope/slope” is used herein to refer to the entire hill mapped in Figure 4. We therefore use the term slope deformation to refer to any ground deformation observed on the hillslope. The terms upslope and downslope are used with reference to the southeast-facing slope aspect in which the coseismic rock avalanche is located. Field photographs additional to those presented in Figure 3 can be found in Figure S2.

In the “flat” portion (<30°) of the Half Moon Bay hillslope at and near the slope crest (Fig. 4A; Fig. S2E), a thin layer of soil covers graywacke rock mass, while most of the steep slopes (>30°) are characterized by densely vegetated blocky graywacke at the surface. In the southeast-facing aspect, the surface material comprises rock outcrops with partial debris cover as well as debris (Fig. 4A; Fig. S2C). The large graywacke rock mass outcrop mapped in the steep southeast-facing aspect (Fig. 4A) is delineated by an upper convex break in slope and a lower concave break in slope and represents the evacuated portion of the rock avalanche headscarp (Fig. 4B; Fig. S2C). Steep slope angles and an upper convex and lower concave break in slope were present in this part of the Half Moon Bay hillslope prior to the 2016 earthquake, the geomorphic expression of which indicates the occurrence of previous rock avalanches (Figs. 2 and 5A). Additionally, a debris lobe with an estimated volume of ~800,000 ± 200,000 m3 lies at the base of the slope, which consists of accumulated rock avalanche and debris-flow material that predates the 2016 coseismic landslides (Fig. 5A). This volume estimate was calculated by subtracting the predeposit surface from the current surface. The predeposit surface was estimated by interpolating elevation contours outside the landslide deposit.

Breaks in slope that are parallel to topographic contours are the predominant features of the southeast-facing aspect and parts of the upper flat terrain of the study area, which we interpreted as scarps (Fig. 3D), counterscarps (Fig. 3C), and grabens. Their relative movement appears to be associated with the dominant gravitational slope deformation mechanism (i.e., oriented downslope and related to the landslide mechanism as opposed to tectonic deformation). In the “flat” upper terrain of the Half Moon Bay hillslope, many of these breaks in slope are rounded in profile, whereas in the lower, steeper portions, sharp depressions and spurs are more common (Fig. 4A; Figs. S2E and S2F). The counterscarps and grabens are typically several meters wide, and scarp heights vary between several decimeters and ~1 m. Feature density analysis of the morphogenetic map (Fig. 4B) showed that the intensity of deformation is largest in the central portion of the slope (Fig. S1). Tensional ground cracks in the eastern portion of the Half Moon Bay hillslope as well as the main rock avalanche headscarp are also oriented ENE-WSW.

In the western portion of the Half Moon Bay hillslope, a sharp ridge-crossing depression and spur crosscut the breaks in slope that run parallel to topography (Fig. 4A; Fig. S2G). Due to its orientation, crosscutting nature, and larger scarp height (~1.5 m) as compared to the slope-parallel counterscarps, we interpret this feature as a tectonic fault (Fig. 4B). This fault has not previously been mapped, but its strike is similar to normal cross-faults mapped along the Seaward segment of the Hope fault by Simpson (1995) and along the Charwell segment of the Hope fault by Eusden et al. (2005). No indication of fault slip in the 2016 earthquake could be resolved. However, because this feature left-laterally displaces some of the slope-parallel counterscarps, it has likely been active during or after the formation of the counterscarps (Fig. 4; Fig. S2G).

Structures striking NNE-SSW (parallel to the lateral 2016 rock avalanche scarp) mainly represent recent deformational features such as ground cracks, lineaments, and scarps (Fig. 4B; Fig. S2H). In the lower portion of the southeast-facing aspect as well as in the west- and northeast-facing aspects, rounded convex and concave changes in slope are associated with local drainage lines and corresponding ridges, respectively (Fig. 4).

4.2. Landslide Displacements

4.2.1. Coseismic Landslide Displacements

The DIC and topographic sampling of pre- and postearthquake imagery revealed average coseismic landslide displacement of ~1.7 m horizontally and ~2.6 m three-dimensionally across the mapped area of coseismic slope deformation (i.e., in the incipient landslide above the headscarp of the coseismic rock avalanche). Coseismic displacements are lower in the upper portion of the slope where the terrain is relatively flat and increase toward the headscarp of the coseismic rock avalanche (Fig. 5). In the central-western portion of the slope (profile X-X′), displacement magnitudes reach up to ~13 ± 1 m, while in the eastern portion of the active slope (profile Y-Y′), maximum displacement magnitudes of ~3 ± 1 m were measured above the coseismic rock avalanche headscarp (Fig. 5).

Upslope ~230 m along profile X-X′, displacements could not be resolved as they were below the measurable 1 m limit of the DIC analysis. Some ground cracking was, however, observed near the concave break in slope following the earthquake, indicating that displacements <1 m did occur in this portion of the slope during the 2016 earthquake. We also note that the plunge of displacement vectors varies along the profile distance, with displacement vectors steeply plunging at ~60°–80° in the upper part of profile X-X′ and becoming gentler downslope (~30°–50°) (see Fig. S4).

No postseismic displacement greater than the 1 m measurable limit of the DIC analysis was identified in the slope between 2017 and 2019.

4.2.2. Estimate of Total Landslide Displacement

In profile X-X′, cumulative total landslide displacements of ~95 m were estimated above the headscarp of the 2016 rock avalanche (Fig. 5D). In profile Y-Y′, cumulative total landslide displacements reach a maximum of ~24 m (Fig. 5E). Based on the assumption that the measured topographic feature offsets are purely gravitational in origin, uncertainties in horizontal displacement at each measurement point are given by the DEM resolution of ±1 m. Vertical uncertainties are dependent on the topography within the horizontal uncertainty. Based on the number of measurement points, cumulative horizontal displacement uncertainty can be given as ±20 m and ±7 m in profiles X-X′ and Y-Y′, respectively. In both profiles, the pattern of displacement increase along the profile is strikingly similar between the coseismic and total landslide displacements, and in both profiles, the coseismic landslide displacement accounts for ~10%– 20% of the estimated total landslide displacement at the surface (Figs. 5D and 5E).

4.3. Rock Mass Structure and Properties

4.3.1. Geophysical Surveys

Results of the ERT survey were processed as outlined in the Geophysical Surveys section of the Methods (section 3.3). Pseudosections are available in section 4 of File S1 (Figs. S6–S9, see footnote 1). ERT survey results indicated variability in rock mass properties that allowed us to interpret bedrock structures and landslide geometry (Fig. 6). Given that ERT sections represent a bulk property that cannot be unequivocally assigned to a physical property, we describe the results in terms of resistivity zones (using reference depths) and discuss possible interpretations in conjunction with other observations. In general, low-resistivity values are consistent with higher fracture densities (assuming the fractures are moist or saturated with groundwater), zones of increased weathering, and/or zones of increased clay content, relative to intact rock masses (Archie, 1942; Perrone et al., 2014).

Between 0 m and ~220 m in the profile, high-resistivity values (>600 Ohm·m) were measured in an ~2–5-m-thick zone along the surface (Fig. 6D). Based on outcrop field observations, this zone likely represents unsaturated and dry, but highly fractured rock with soil. In the uppermost part of the profile (i.e., above the concave break in slope at ~40 m along profile), resistivity values range between ~400 and 1000 Ohm·m below the highly resistive zone and generally increase with depth. We interpret this high-resistivity zone as the (relatively) intact rock mass. At the concave break in slope (~40 m along profile), a low-resistivity zone extends to ~25 m depth. Downslope of the concave break in slope (~40–220 m along profile), resistivity values at depths between ~2 and 20 m are low (<400 Ohm·m). This low-resistivity zone increases in depth further downslope (>220 m along profile), with bulk resistivity decreasing in the downslope direction, and corresponds to a zone of low seismic velocities (Fig. 6C). The combination of low resistivity and low seismic velocity indicates that the zone is defined by high rock mass fracturing and, potentially, weathering. Due to the geometry of this low-resistivity and low-seismic-velocity zone pinching out in the upper part of the profile and increasing in depth downslope, we interpret its base as disturbed ground likely originating from coseismic sliding deformation within the landslide complex. Between ~40 and 220 m along profile, there is a zone of higher resistivity (up to 800 Ohm·m) below the low-resistivity zone (depth >25 m). The higher resistivity in this zone likely corresponds to less rock mass fracturing and can thus be interpreted to represent (relatively) intact rock mass.

Interestingly, several subvertical low-resistivity zones (<300 Ohm·m) crosscut the intact rock mass (at 40 m, 100 m, 200 m, and 240 m along profile; Fig. 6). Even though the dip angle of these low-resistivity zones cannot be directly interpreted from the ERT results, these low-resistivity zones likely represent subvertical structures in the rock mass. Consistent with the interpretations above, these low-resistivity zones are likely associated with near-vertical fracture zones. At the surface, the locations of these low-resistivity zones correspond to subtle topographic features (depressions and concave breaks in slope) that can be observed in the field (see the direct comparison of the locations of subvertical fracture zones with field and helicopter photos in Fig. S14).

The seismic refraction profile was located to the west of the ERT section and spanned the distance between ~120 m and 340 m along the ERT profile (Fig. 6). Examples of seismic traces and the first-arrival picks, as well as ray paths and modeled traveltimes, are available in section 4 of File S1 (Figs. S11–S13, see footnote 1). The stacking of shots provided a good signal-to-noise ratio, and first arrivals were consistently picked by hand as shown in the examples. For the tomography inversion, different initial models were tested and revealed similar results.

In the seismic refraction profile, P-wave velocities reach ~1500 m/s at a depth of ~25 m in the upper part of the profile (−18 m to 20 m along profile) and are as low as ~300 m/s to 400 m/s at depths shallower than ~5 m (Fig. 6C). The seismic velocities generally increase with depth. A slightly higher gradient in seismic velocities is located at a depth of ~10–15 m between −18 m and ~40 m along profile, and its depth subsequently increases until ~105 m along profile, where it reaches the limit of the survey’s depth of investigation as indicated by the modeled ray coverage (Fig. 6C). In conjunction with the ERT and borehole results, we interpret this higher-velocity gradient as the possible base of disturbed ground. Furthermore, the seismic refraction profile shows a relatively low P-wave velocity increase with depth between profile distances 100 and 192 m. This localized persistence of low P-wave velocities to depth corresponds with the decrease in bulk resistivity at depth in the same location on the ERT line. As mentioned above, the combination of lower P-wave velocities and lower resistivity provides strong evidence for an increase in rock mass fracturing.

4.3.2. Borehole and Downhole Geophysics

A borehole was drilled in the incipient portion of the landslide (borehole location latitude and longitude: 42.2449°S, 173.8145°E). The location of the borehole experienced ~1.6 m of coseismic landslide displacement (based on DIC; see section 4.2.1). Core recovery was very good to excellent, with an average recovery of 91.7% (Fig. 7). In the core, a thin soil cover of ~40 cm was overlying completely weathered to highly weathered Torlesse graywacke sandstone with rock mass weathering generally decreasing with depth. At ~36.5 m, rock mass weathering transitioned from completely to highly weathered to highly to moderately weathered (Fig. 7). Rock strength was highly variable, ranging from extremely weak to strong at depths shallower than ~35.2 m. Some moderately weathered intact sandstone fragments at greater depth were very strong (Fig. 7). While the lithology is predominantly fine- to coarse-grained graywacke sandstone, argillite beds were interbedded with sandstone at depths between ~43.0 and 59.6 m. Bedding dips subhorizontally with the steepest dip measured at ~20° in the core and between 8° and 27° in the downhole geophysics (Fig. S16). This is consistent with the subhorizontal NW-dipping bedding observed in other locations in and around the landslide (section 4.3.3; Fig. 8).

Defects in the borehole core were logged either individually or as defect zones, where joint spacing was very close to extremely close (NZGS, 2005). A detailed borehole log including a description of the defect zones is available in Figure S15. The borehole core was in general highly fractured, with the intensity of fracturing decreasing with depth (Fig. 7, columns showing defect zones and mean defect spacing). This is also reflected in the rock quality designation (RQD) (Fig. 7, column RQD), which is very poor to poor (<50%) between 0 m and 43.7 m depth and includes sections of fair to good (50%– 90%) between 43.7 m and 59.6 m. In the borehole core, most logged joints were open, while some sandstone fragments in the core also displayed incipient joints (Fig. 7). No healed (i.e., cemented or veined) joints were observed. Joint roughness was variable, with many joints displaying planar (rough or irregular) rough surfaces. Many of the joints showed black and brown iron-oxide coatings or were black or brown stained and weathered (Fig. S15). This degree of alteration suggests that many of the joints predated the 2016 earthquake. The structures identified through downhole geophysics between a depth of 21.15 m and 46.2 m generally reflected the same patterns in joint dip as the core measurements (Fig. 7) but showed a sampling bias in terms of joint dip (because vertical joints are more difficult to resolve in downhole geophysics). The mean defect spacing based on downhole geophysics is ~33 cm and therefore an order of magnitude higher than that measured from the core. This difference is likely due to the limited resolution of the downhole geophysics and the influence of extremely closely spaced jointing of defect zones logged in the core that were considered when calculating the mean defect spacing. The orientations of the structures identified in downhole geophysics are further analyzed in section 4.3.3.

Defect zones logged in the borehole core were varied in appearance and were classified into three main categories: zone recovered as gravel or cobbles, zone of very closely to extremely closely spaced joints, and weak zone/soil consisting of decomposed or crushed material and infilled zones. Each of these categories and their interpretation are described in Table 1. The presence of infilled zones as well as the high fracture density in the borehole highlight the dilated nature of the rock mass, which was also observed in outcrop (section 4.3.3). Decomposed and infilled zones displayed variable weathering grades, indicating progressive rock mass degradation.

4.3.3. Outcrop Observations and Structural Analysis

Graywacke rock mass is exposed in outcrops in and around the Half Moon Bay hillslope, at which we qualitatively described the rock mass and collected structural measurements (Fig. 8A). Here, we provide an overview of the rock mass characteristics and structures in outcrop. A detailed description is available in section 6 of File S1 (see footnote 1).

Bedding is oriented with a mean dip direction/dip of 321/15 in outcrops within the landslide and 245/20 based on downhole geophysics in the borehole HMB-01 (Table 2). In outcrops along the Rakautara Stream, bedding is steeper dipping with a mean orientation of 022/52 (Table 2). Cleavage was identified in outcrops along the Rakautara Stream as well as in HMB-01, with orientations varying distinctly (Table 2).

In outcrops within the landslide, two main joint sets (JS1 and JS2) were identified. Joint set 1 (331/71), oriented parallel to the strike of the Hope fault and antithetic to the slope, is dominant in terms of persistence (Table 2). Although no joint set was identified as dominant in the HMB-01 borehole, JS1 was identified in downhole geophysics (Table 2). Across the highly fractured outcrops and borehole, many random joint orientations were present. This random jointing and highly complex structural setting make the designation of joint sets and correlation of joint sets between measurement locations challenging. Dominant joint sets could be identified, however, by comparisons between locations, and they were qualitatively linked to tectonic structures and landslide mechanism as discussed in section 5.1.3.

On the shore platform, five joint sets were identified (Table 2; Fig. 8A). While the orientation of JS1 is not represented in the shore platform data, a joint set similar to JS2 was observed (Table 2; Fig. S17). Lineament mapping of the shore platform revealed a predominant NW-SE structural orientation coinciding with the strike direction of JS3 (Fig. 8A). To a lesser extent, lineaments oriented SW-NE were identified, which likely correspond to JS5 identified in the shore structural measurements (Fig. 8A). Along the Rakautara Stream, joint sets in orientations similar to the previously defined joint sets JS1, JS3, and JS4 were identified (Table 2).

4.4. Regional Morphotectonic Features

Regional lineament mapping was carried out in the Torlesse basement rocks between the Puhipuhi syncline and the Papatea fault (Fig. 8B). Evaluation of lineament mapping showed that, along the coast (~<4 km from the coast), the dominant lineament orientation strikes parallel to the Hope fault (SW-NE) (Fig. 8B). In the northeast of the lineament mapping area, lineaments oriented parallel to the orientation of the Papatea fault (NNW-SSE) become more common (Fig. 8B), and, to a lesser extent, in the western portion of the lineament mapping area, some lineaments, mainly expressed as ridgelines, strike parallel to the fold axis of the Puhipuhi syncline (NNE-SSW). Overall, the most dominant orientation of lineaments is parallel to the strike of the Hope fault (Fig. 8B).

5.1. Interpretation of the Landslide Failure Mechanism

5.1.1. Overview

We integrated the results from geomorphological mapping, coseismic and total landslide displacement analysis, geophysical surveys, borehole data, structural analysis, and regional morphotectonic mapping to interpret the slope deformation mechanisms. A 3-D model of the interpreted mechanism is available in the ground model in the Supplemental Material (see footnote 1). Based on the observations detailed in sections 5.1.2–5.1.5, we interpret the failure mechanism to comprise a combination of slope deformation through flexural toppling facilitated by the subvertical tectonic structures as well as joint-step-path sliding along preexisting, closely spaced discontinuities within the rock mass (Fig. 9; Singeisen et al., 2022). While slope deformation is thought to occur through both sliding and flexural toppling, the toppling mechanism likely accounts for most of the slope deformation at depth, whereas sliding along joint-step-path failure surfaces becomes predominant at shallower depths as rock mass dilation increases toward the slope’s free face. As such, the coseismic rock avalanche (Fig. 4B, headscarp) is likely to have resulted from sliding along multiple joint-step-path failure surfaces within the heavily deformed rock mass (Fig. 9). This interpretation agrees with and expands on the work by Clarke and Burbank (2010), who presented conceptual models of landsliding in regions with dominant tectonic and geomorphic fracturing, respectively. At Half Moon Bay, both models are superimposed. Tectonic fracturing is affecting the depth of slope deformation, and intense geomorphic fracturing at the surface is leading to shallow avalanching.

The evolution of the landslide complex results from (1) regional tectonic structures in the fault damage zone that increase susceptibility to coseismic slope deformation and (2) progressive rock mass deformation in oversteepened coastal slopes that can locally cause joint-step-path sliding, and, ultimately, lead to rock avalanching.

5.1.2. Geomorphological and Displacement-Based Evidence

Gravitational slope deformation above the coseismic avalanching headscarp at Half Moon Bay is geomorphologically expressed through scarps and counterscarps. Geomorphologic markers of deformation are subdued in northernmost and lateral portions of the slope (Fig. 4B, slope deformation subdued) and become more distinct in the central portion of the slope (Fig. 4B, slope deformation + sliding and predominantly sliding). Scarps north of the relatively “flat” slope crest portion of the slope mark the upslope extent of an area with intense deformation that is expressed through a high density of counterscarps, scarps, and ground cracks (Fig. 4B, slope deformation + sliding and predominantly sliding). In this central portion of the slope, coseismic displacements >1 m were measured with magnitudes increasing and plunge angles generally decreasing toward the headscarp of the coseismic rock avalanche (Figs. 5 and 6B; Fig. S4). This observation indicates that these scarps likely formed due to sliding, while counterscarps are typically recognized as a surface expression of toppling or “ridge renting” (Beck, 1968; Hungr et al., 2014), a gravitational faulting mechanism that can occur along ridgelines as a result of earthquake shaking (Beck, 1968; Korup, 2006). Based on the increase of scarp frequency toward the rock avalanche headscarp (Fig. 4B) and the decrease of displacement vector plunge, sliding appears to become dominant toward the rock avalanche headscarp, as indicated in Figure 9. Geomorphological features observed in the pre-2016 imagery (Figs. 2 and 5) also indicate that slope deformation prior to the 2016 earthquake was affected by similar failure mechanisms. The presence of a relict rock avalanche headscarp and deposit as well as scarp features in the “flat” incipient portion of the landslide suggest that the 2016 earthquake only represents the most recent evolution in a longer history of slope deformation at the site. The discrepancy of landslide displacement magnitudes between coseismic measurements (maximum of 13 m in 2016 event) and estimates of total landslide displacement (maximum of 96 m from total topographic scarp and counterscarp offsets including 2016) supports this interpretation.

5.1.3. Structural Evidence

The orientation of bedding was identified in the borehole as well as at several outcrops in and around the Half Moon Bay landslide complex (Fig. 8). The shallow dip angle of bedding into the slope at Half Moon Bay (14°–23° toward the NW in outcrops on the landslide and 5°–20° in the borehole), however, does little to explain landslide kinematics. The most dominant joint set orientation identified in outcrops within the incipient portion of the landslide (JS1) is oriented parallel to the strike direction of the most dominant regional lineaments mapped (Fig. 8). With an average dip angle of ~70° toward the NW (i.e., dipping into the slope), JS1 can also be associated with the subvertical fracture zones identified in the ERT section and can be correlated with the counterscarps in the slope. In the downhole geophysics, the joint set identified as JS1 is shallower dipping (~43°). This could be a result of either sampling bias toward structures of gentler dip angles or increased toppling rotation in this central part of the slope (i.e., at the borehole location) as compared to the structural measurement outcrop locations.

Because the strike of JS1 coincides with the strike direction of the Hope fault and is oriented antithetic to the slope, we infer that these structures likely originated as tectonic structures within the Hope fault damage zone (Fig. 9). Based on their orientation, either these structures could represent shear foliation, or they could have formed as synthetic strike-slip and normal faults within an internal fault wedge, similar to structures described by Eusden et al. (2005) along the Conway segment of the Hope fault. This interpretation is supported for the Half Moon Bay site by the observation that JS1 is absent in structural measurements on the shore platform (i.e., indicating that JS1 represents structures that developed in the uplifted hanging wall of the Hope fault). The orientation of this discontinuity set kinematically enables flexural toppling to occur in the slope. The prominence of counterscarps in the slope suggests that these features have been important in the development of the landslide complex prior to the earthquake in 2016 (Fig. 9). This interpretation is supported by dynamic numerical modeling conducted by Singeisen (2023). Modeling results showed that these defects facilitate flexural toppling and are required to replicate (1) the pattern of coseismic landslide displacements observed following the 2016 earthquake and (2) the type and depth of slope deformation observed in geotechnical and field investigations.

Joint sets JS2 and JS4 are oriented parallel to the strike direction of the Papatea fault and the fold axis of the Puhipuhi syncline, respectively (Fig. 8). The N-S–striking lineaments and ground cracks mapped in the slope morphology (Fig. 4) are potentially related to JS2 and JS4. These structures accommodate recent slope deformation and are expressed at the surface as minor scarps and ground cracks that crosscut and offset the slope-parallel counterscarps. The western lateral scarp of the 2016 rock avalanche is likely related to JS2 and JS4.

The orientation of JS3, which was identified as the dominant orientation of lineaments in the shore platform (Fig. 8), was not identified in outcrops within the landslide and only uncommonly in outcrops along the Rakautara Stream. Similarly oriented structures in the shore platform of the Kaikōura Peninsula (~22 km from the site) are predominantly associated with antithetic left-lateral strike-slip faults (Campbell et al., 2005). As such, the mapped fault within the landslide (Fig. 4) is oriented parallel to JS3 and is likely a conjugate Riedel structure associated with distributed NE-SW–oriented dextral shear. The observed left-lateral offset of counterscarps along the fault trace is consistent with this interpretation.

Joint set JS5 dips steeply toward the SE and therefore likely accommodates slope free-face displacement (i.e., gravitational extension), represented by the SW-NE–oriented tension cracks above the 2016 rock avalanche headscarp. Joint set JS6 and cleavage orientations could not directly be related to slope geomorphology.

Apart from the specific joint sets described herein, graywacke outcrops are characterized by close jointing, including random joint orientations not captured in the structural analysis. Sliding can therefore occur along centimeter- to decimeter-scale discontinuities in the rock mass. For other coseismic landslides in Torlesse graywacke rock mass triggered by the Kaikōura earthquake, Singeisen et al. (2022) identified that the geometry of failures cannot be explained by the identified joint sets of meter-scale persistence alone but is consistent with sliding along centimeter- to decimeter-scale discontinuities in the highly jointed rock mass. While flexural toppling can result in sliding (Adhikary et al., 1997; Hungr et al., 2014; Zheng et al., 2021), brittle failure of intact rock is not required at Half Moon Bay due to the closely fractured nature of the Torlesse graywacke rock mass. Instead, ridge renting and flexural toppling in this case are likely to contribute to slope destabilization through the opening of joints and fractures, which in turn leads to decreasing rock mass interlocking and allows for failure to occur more easily along centimeter- to decimeter-scale discontinuities (i.e., what we refer herein as joint-step-path failure; Singeisen et al., 2022).

5.1.4. Evidence from Geotechnical Ground Investigations

The borehole shows an intensely fractured and weathered rock mass containing many broken zones, some of which are extremely weak (i.e., strength of a soil [NZGS, 2005]; Fig. 7; Fig. S15). Mean defect spacing is low throughout the entire borehole (<5 cm), and particularly low (<2 cm) between ~15 and 35 m depth. This intense fracturing and the presence of weak zones with no systematic orientations in the borehole support the interpretation of a joint-step-path failure mechanism. In the closely jointed rock mass, no fracture propagation of intact rock is necessary for the development of a failure surface. Instead, internal distributed rock mass deformation can occur, and several failure planes might develop along multiple degrees of kinematic freedom. The absence of a single identifiable shear zone further supports the hypothesis that sliding likely occurs along multiple joint-step-paths within the rock mass. At a depth of ~36.5 m, rock weathering in the core transitions from highly to moderately weathered, and rock strength increases from extremely weak to moderately strong to predominantly strong to very strong with a few localized weak zones (NZGS, 2005). The increase of mean defect spacing and decrease of defect zone frequency from ~35 m down suggest that rock mass deformation decreases with depth. We therefore interpret that the base of the main sliding deformation at the borehole is located at ~37 m depth, even though sliding along localized weak zones may still occur at greater depths. Since localized weak zones are present to the full penetration depth of the borehole, the depth of slope deformation is likely deeper than 60 m. Furthermore, the presence of numerous fragmentation zones (shaking-induced or landslide-induced infilled and decomposed zones) with variable weathering grades provides further evidence of progressive (i.e., long-term) degradation of the rock mass.

The results of the geophysical surveys showed a low-resistivity and low-seismic-velocity zone that reaches down to the survey penetration depth of ~35 m (ERT survey) and 25 m (seismic survey), respectively (Fig. 6). Furthermore, subvertical low-resistivity zones were observed dipping into the slope. These physical characteristics are related to high rock mass fracturing (Heincke et al., 2010; Perrone et al., 2014) and increased chemical weathering (St. Clair et al., 2015). The structure of the low-resistivity and low-seismic-velocity zone does not correspond well with weathering models (St. Clair et al., 2015; Riebe et al., 2017). The zone is not homogeneously dept h-dependent, but it increases in depth downslope. We interpret the zones of low resistivity and seismic velocity as (1) subvertical fracture zones associated with tectonic structures along which ridge renting and flexural toppling can be accommodated and (2) zones of high rock mass damage and weathering in the sliding portion of the landslide complex, with increased fracturing toward the coseismic avalanche headscarp (Fig. 9). This interpretation correlates with the observations of a highly fractured and weathered rock mass in the borehole and outcrops that show increasing dilation toward the 2016 coseismic rock avalanche headscarp.

5.1.5. Slope Deformation Depth

Based on interpretations from the borehole and geophysical surveys as well as the elevation of lateral geomorphic markers, we can constrain the depths of the main sliding (i.e., flexural toppling + joint-step-path sliding) and deeper slope deformation (i.e., flexural toppling and ridge renting) mechanisms in the Half Moon Bay landslide complex (Fig. 9). In the upper part of the Half Moon Bay hillslope (i.e., north of the scarps delineating the main sliding deformation), the depth of the landslide complex is interpreted from the geophysical surveys (low-resistivity and low-seismic-velocity zone) to be ~20 ± 5 m (Fig. 6). The geophysical surveys and borehole indicate that the depth of the main sliding deformation increases toward the headscarp of the rock avalanche. While the depth of the main sliding deformation at the borehole location could be constrained to ~37 ± 1 m (Fig. 7), we estimate, based on the interpolation of lateral geomorphic markers, that the overall depth of slope deformation at the location of the borehole reaches at least ~120 m (assuming a planar failure surface based on the elevation of lateral scarps). According to the conceptual model (Fig. 9), we estimate the depth of slope deformation at the borehole location to be at ~150 ± 20 m.

Based on the 3-D ground model in the Supplemental Material (constrained by the considerations in sections 5.1.1–5.1.4), the volume of the deformed rock mass (flexural toppling + joint-step-path sliding) is estimated at ~60 ± 20 × 106 m3, while the volume of the main sliding deformation is estimated at ~10 ± 3 × 106 m3. Uncertainties are based on limited depth constraints as well as the depth uncertainties stated above.

5.2. Landslide Failure Complex Evolution

Given the location of the Half Moon Bay landslide complex proximal to the Hope fault and in an uplifted steep coastal slope, it is crucial to consider the tectonic history of the site in order to gain an in-depth understanding of the evolution of this landslide complex.

Exhumation initiated ca. 25–20 Ma in the southern Marlborough fault system and became widespread in the Seaward Kaikōura Range by ca. 5 Ma (Collett et al., 2019). Exhumation rates from thermochronology and 10Be have been measured between 0.3 and 0.5 mm/yr (Collett et al., 2019; Wilkinson et al., 2022), while Pleistocene uplift rates of ~1.1 mm/yr have been measured along the Kaikōura coast (Ota et al., 1996; Nicol et al., 2022). With a hillslope elevation of ~550 m and uplift rate of 1.1 mm/yr, the required time of uplift in the coastal hanging wall of the Hope fault at Half Moon Bay can be derived to a first order as ~500,000 yr at a minimum (wrongfully assuming no erosion) and ca. 1 Ma assuming a denudation rate of 0.5 mm/yr. Dextral movement on the Hope fault initiated ca. 1 Ma (Wood et al., 1994) and would have caused a total dextral offset of ~2–5 km along the Seaward segment of the Hope fault given a slip rate of 2–5 mm/yr (Van Dissen and Yeats, 1991; Wallace et al., 2007; Collett et al., 2019). The record of earthquakes on the Seaward segment of the Hope fault, however, is sparse, and it is thought that many ruptures transfer slip from the Hope fault to the Jordan thrust at their junction south of Half Moon Bay (Van Dissen and Yeats, 1991; Howell and Clark, 2022).

Regarding the overall seismicity at the site, the National Seismic Hazard Model indicates a 10% in 50 yr probability of exceedance for ground motions larger than ~0.7 g (Gerstenberger et al., 2022). This roughly corresponds to the ground motions at the site during the Kaikōura earthquake, which resulted in substantial slope deformation and partial failure as a rock avalanche. Based on this exceedance probability, an average number of 1000 seismic events >0.7 g could be expected at the site within a 500,000 yr time window, assuming that the seismic hazard has not changed over time. This illustrative calculation demonstrates that during the time of uplift, the Half Moon Bay hillslope experienced repeated seismic shaking. It is very likely that previous earthquakes also caused coseismic displacements, fracturing, and rock mass damage, contributing to increasing cumulative rock mass damage over time prior to the 2016 Kaikōura earthquake at Half Moon Bay. Rock mass damage in the fault damage zone and the locally deformed rock mass can also cause amplification of seismic shaking, further enhancing these effects (Martino et al., 2006). While healing effects are discussed in the literature (Marc et al., 2015), the rock type at the Half Moon Bay site (graywacke sandstone) typically behaves in a brittle manner. Laboratory studies have shown that failure preferentially occurs along preexisting joints, indicating that cumulative rock damage plays an important role in slope destabilization (Brideau et al., 2020) and in the postearthquake performance of the slope (Massey et al., 2022). The amount of rock mass damage along with the absence of healed (i.e., veined or cemented) joints observed in the borehole furthermore suggest that healing processes at this site do not play a significant role over the time periods between strong ground-shaking events (approximately every few hundred years at this site). Furthermore, past coastal erosion has likely steepened the Half Moon Bay site and contributed to landslide development (Bloom et al., 2023).

Other factors, such as vegetation and hydrological processes, could also affect slope stability. In the case of the Half Moon Bay landslide, the depth of slope deformation is too large for the failure process to be significantly affected by vegetation. However, circulation of meteoric water in the highly fractured rock mass is likely to enhance chemical weathering within the landslide, causing rock strength degradation. Short-term increases in pore-water pressure due to rainstorms could also affect slope stability and could cause slope deformation or trigger debris avalanches in the already dilated rock mass. This being said, despite the occurrence of Ex-Tropical Cyclones Debbie and Cook in 2017, no significant slope displacements (>1 m) were observed postearthquake. Based on postseismic observations, however, the effect of rainstorms is particularly significant in terms of secondary debris mobilization. Groundwater effects are unlikely to considerably affect the deeper-seated slope deformation because there is little upslope catchment feeding the slope. Currently, fracturing of the deformed rock mass is so intense that little to no groundwater is retained in the upper slope (water could not be retained in borehole HMB-01 as a result of high porosity and short water retention time).

Based on these considerations, the rapid rate of uplift, faulting, and the high frequency of strong seismic shaking likely played the most important roles in the evolution of this landslide complex. Considering the relative contribution of 2016 coseismic displacement (~10%– 20%) to the total estimated slope displacements, the landslide failure complex could have evolved in just the last ~5000 yr, if we assume that the relative contribution of slope displacement remained constant over time. However, it is likely that single-event contributions to overall slope displacements increased over time due to cumulative rock mass damage and changes in degrees of kinematic freedom.

Repeated seismic shaking likely allowed joint-step-path sliding to occur along centimeter- to decimeter-scale defects and weak zones, while enhanced weathering is likely to have affected the highly fractured rock mass during interseismic periods, causing further rock strength degradation (Brideau et al., 2020). The presence of defect zones with varying weathering grades in the borehole illustrates the progressive degradation of the rock mass. Amplification of ground motion in the damaged and weathered rock mass during seismic events (i.e., amplification caused by impedance contrasts in the fault damage zone and preexisting rock mass damage; Martino et al., 2006) likely additionally led to slope deformation and sliding along defects and weak zones. Dynamic numerical modeling by Singeisen (2023) showed topographic amplification at the slope crest at Half Moon Bay, with amplification increasing further when considering material strength effects and discontinuities. Based on the intense weathering and fracturing of the rock mass observed at Half Moon Bay (Fig. 7; Fig. S15), we infer that this positive feedback mechanism contributed to the development of joint-step-path sliding, and partial rock avalanching, in the hillslope during the 2016 Kaikōura earthquake. These same mechanisms are likely to continue into the future.

The structures facilitating slope deformation through flexural toppling (JS1) are likely linked to the Hope fault due to their parallel orientation and steep dip angle. Similar structures on the Charwell segment of the Hope fault have been interpreted as internal fault wedge structures that develop at the range front in response to transpressional deformation (Eusden et al., 2005). In case of the Charwell segment, normal movement on these structures due to gravitational collapse has resulted in the formation of counterscarps in the slope. In the glacially oversteepened, high-relief Torlesse graywacke terrain of the Southern Alps of New Zealand, counterscarps are often associated with ridge renting. Even though ridge renting is not commonly associated with fault damage zones or regional tectonic structures in particular, this process could have occurred and contributed to slope deformation in the oversteepened coastal Half Moon Bay hillslope. The counterscarps observed in the Half Moon Bay hillslope therefore likely resulted from normal-fault displacement along JS1 structures and ridge renting, as well as later-stage flexural toppling.

In the Half Moon Bay hillslope, the JS1 counterscarps are intersected and offset by a NW-SE–striking, ridge-crossing fault that is oriented parallel to JS3. Left-lateral antithetic faults of similar orientation are present in the Kaikōura shore platform (~22 km from the site; Campbell et al., 2005). Furthermore, Simpson (1995) identified similar NW-SE–striking normal cross-faults (approximately perpendicular to the main strand of the Hope fault) in the Seaward segment of the Hope fault, and Eusden et al. (2005) interpreted normal faults in a similar orientation relative to the Hope fault along the Charwell segment as late-stage secondary normal faults in the extensive fault wedge collapse. In the Half Moon Bay hillslope, the JS1 counterscarps as well as the NW-SE–striking, ridge-crossing fault are furthermore offset by NNE-SSW–oriented lineaments and ground cracks. Simpson (1995) also identified N-S– and NE-SW–striking cross-faults with normal components of movement, downthrown to the SE, between the Hapuku River mouth and Half Moon Bay (Fig. 1) and associated them with transtension in the eastern Seaward segment of the Hope fault, resulting in surficial gravitational collapse. We therefore interpret the geomorphology of the hillslope to be a result of tectonic and gravitational interactions (Fig. 10). The striking similarities between the total and coseismic displacement patterns (Fig. 5) suggest that the same structures repeatedly accommodate displacement and that the mechanism of slope displacement has remained similar over longer timescales (~5000 yr).

5.3. Implications for Landslide Susceptibility and Hazard

This case study demonstrates that rock mass damage and tectonic structures in fault damage zones, as well as seismicity, not only play a critical role in coseismic slope stability (Bloom et al., 2021), but they can govern the evolution of large slope instabilities over long timescales (Fig. 10). Several case studies of large rock slope instabilities around the globe have identified similar failure mechanisms of combined toppling and sliding in rock masses with a structural fabric steeply dipping into the slope (Amann, 2006; Bois et al., 2008; Chemenda et al., 2009; Brideau et al., 2011; Bouissou et al., 2012; Glueer et al., 2019), but here we can link the structural fabric to active faulting and landscape evolution. Furthermore, slopes in fault damage zones may be particularly predisposed to large landslides through progressive failure as rock mass damage accumulates over multiple seismic cycles (Gischig et al., 2016; Agliardi et al., 2019). In this context, the Half Moon Bay landslide complex provides a rare example where coseismic displacement and failure were observed in a recent earthquake and can be combined with a unique data set from geotechnical ground investigations and the tectonic site history to build an integrated model of landslide evolution (Fig. 10).

The landslide evolution model for Half Moon Bay can be extended to other failures around range-bounding strike-slip faults in New Zealand’s Torlesse graywacke, assuming a similar degree of rock mass deformation and similar structural fabric. Predisposition to deep-seated failure as well as rock mass strength effects in steep slopes affected by high seismicity may also, in part, explain the higher landslide source area density (i.e., percent of area covered by landslide source) identified around surface-rupturing faults during the Kaikōura earthquake (Massey et al., 2020a; Bloom et al., 2022).

Based on the detailed site characterization of the Half Moon Bay landslide complex and the past evolution of the slope, debris and rock avalanche failure of the currently incipient portion of the slope should be expected in future earthquakes and potentially during rainfall events. The future evolution of the Half Moon Bay landslide therefore presents a hazard to the transport corridor (State Highway 1 and the Main North Line railway) located on the shore platform at the base of the slope. The hazard posed to the infrastructure corridor is dependent on the response of the rock slope to seismic shaking as well as the volume of material mobilized in a single future event, which affects runout distance. The evaluation of potentially hazardous rock avalanche scenarios at Half Moon Bay lies beyond the scope of this study, but the results of the detailed site characterization presented herein can be used to inform hazard assessment. As such, in-depth understanding of the landslide failure mechanism, which could lead to potentially hazardous rock avalanching, as well as characterization of the extent and depth structure of the landslide complex can provide important constraints for hazard assessments.

Based on geomorphological mapping, coseismic and total landslide displacement analysis, geophysical surveys (ERT and SRT), a 60-m-deep borehole, and structural analysis, we developed a conceptual ground model of the Half Moon Bay landslide complex located ~1 km from the surface trace of the Hope fault in the Marlborough fault system of New Zealand’s South Island. Combining this detailed site characterization with the geomorphic and tectonic context of the study site, we interpreted the slope deformation and landslide failure mechanism to consist of combined flexural toppling along subvertical tectonic structures in the fault damage zone of the Hope fault and multiple joint-step-path sliding along centimeter- to decimeter-scale defects and highly weathered weak zones in the closely jointed and heavily damaged rock mass.

Our results suggest that a combination of high uplift rates, local slope oversteepening in proximity to the coast, faulting, seismicity, and weathering of the damaged rock mass led to the development of this area of slope deformation and rock avalanching. This study illustrates the importance of the interconnectedness among seismic, tectonic, coastal, and gravitational processes in tectonically active regions like New Zealand. While tectonic structures facilitated the development of a deep-seated landslide failure complex through flexural toppling, rock mass damage through a positive feedback cycle of high seismicity, weathering, rock strength degradation, and ground motion amplification resulted in sliding through a joint-step-path mechanism (favored by steep slope angles) and ultimately caused coseismic rock avalanches to occur. As such, the future evolution of this and other large landslide complexes should be given special attention in hazard and risk analyses of probable future earthquake scenarios.

1Supplemental Material. File S1: Additional information on methodology and results. File S2: Geomorphological and structural data. File S3: Ground model provides a 3-D model of the landslide including the two geophysical profiles and interpretation of landslide failure mechanism. The model can be viewed using the free desktop application Leapfrog viewer, which can be downloaded from https://www.seequent.com/products-solutions/leapfrog-viewer/. Please visit https://doi.org/10.1130/GEOS.S.24595647 to access the supplemental material, and contact editing@geosociety.org with any questions.
Science Editor: Andrea Hampel
Associate Editor: Andrew V. Zuza

This work was funded by the New Zealand government as part of the Ministry of Business, Innovation and Employment Endeavor project “Earthquake Induced Landscape Dynamics” and supported by the New Zealand Earthquake Commission. We would like to thank James Chapman and team (CW Drill), and Daniel Stevenson (South Pacific Helicopters) for their great efforts to realize the drilling project in this challenging location. We would also like to thank Oliver Gibson (Resource Development Consultants Limited) for his work on the downhole geophysics, and we would like to acknowledge the great help with geophysical data acquisition by Thomas Brakenrig (GNS Science), Caleb Gasston, Martin Brook, Julie Rowland, Adil Hameed, and Jivyde Despojo (University of Auckland). Thanks go to Valentin Bickel for his great advice and help regarding the derivation of coseismic landslide displacements. We would also like to thank Bill Medwedeff and an anonymous reviewer for their constructive comments. Last, but not least, this project would not have been possible without the support by the Department of Conservation, New Zealand.

Gold Open Access: This paper is published under the terms of the CC-BY-NC license.