Landslides in the Upper Submarine Slopes of Volcanic Islands: The Central Azores

Small landslides in the upper submarine slopes of volcanic islands present potential hazards locally because of their high frequency. We examine evidence for landsliding in high‐resolution bathymetric data from Faial, Pico, São Jorge, and Terceira islands of the Azores. Because the rugged morphology of the upper slopes makes landslides difficult to interpret, we develop two classification schemes for the 1,227 identified slope valleys. One scheme addresses how recognizable the valleys were as originating from landslides (whether scarps are prominent or indefinite), whereas the other scheme addresses valley types (whether apparently produced by single or multiple failures). Size distributions are used to assess the relative occurrence of large versus small landslides. Thirteen landslides are predicted to have generated tsunami heights at source of >1 m and one with height of >7 m. Some slopes have gradients far above 30°, the angle of repose of incohesive clastic sediment, so the seabed in those areas is strengthened perhaps by carbonate cementation, by seismic shaking or by the presence of coherent lava or lava talus. Using all types of slope valleys, Faial and Pico have smaller affected volumes per unit slope area than those of São Jorge and Terceira. These differences could be associated with varied seismic activity, with more frequent earthquakes beneath Faial and Pico preventing the build‐up of sediments on their slopes. Submarine landslide statistics are therefore potentially useful for assessing long‐term earthquake hazards of volcanic islands in seismically active environments such as the Azores.

and Hawaiian Islands (Garcia et al., 2006). Although tsunamis they generate are smaller, those tsunamis may nevertheless be significant geohazards because of their greater frequency. According to historical reports analyzed by Andrade et al. (2006), at least 23 tsunamis have affected the Azores since the human settlement of the islands and one with wave height >10 m. Some of them were not explained by earthquakes so they may have originated from submarine landslides. One motivation for this research was therefore to investigate the sizes and geometries of upper-slope submarine landslides among the islands, in order to assess whether the landslides were large enough to generate tsunamis.
Earthquakes also pose significant geohazards on volcanic ocean islands and can include events of magnitude M > 7.0 (Crossen & Endo, 1982;Furumoto & Kovach, 1979;Madeira et al., 2015). Unfortunately, the instrumental record is not sufficient to characterize seismic hazards of small areas (Scholz, 2002). As the central and eastern Azores islands lie in a region of distributed deformation of the Eurasia-Nubia plate boundary (Lourenço et al., 1998; Figure 1), they are seismically active. The question of adequacy of the instrumental record is, therefore, important locally. In particular, onshore mapping and sonar data has revealed potentially active faults of 10-20 km in length (Casalbore et al., 2015;Lourenço et al., 1998;Madeira et al., 2015;Mitchell et al., 2012), whereas the north slope of São Jorge, which appears fault controlled (Walker, 1999), is 60 km long and the walls of the Terceira Rift ( Figure 1a) reach 100 km in length. Madeira et al. (2015) estimated that earthquakes up to M = 6.8 could occur based on the lengths of active faults mapped on the islands, but this is likely an underestimate as the submarine faults are longer. From data compiled by Scholz (1994), a strike-slip rupture of 110 km length (the approximate maximum of the faults outlined later) could produce an earthquake with seismic moment ∼5 × 10 20 N m (i.e., moment magnitude Mw ∼ 7.4). Such an earthquake would be similar in size to events that have occurred in the Azores historically (e.g., Mw ∼ 7.2; Gaspar et al., 2015). If such a large event were to occur in the Azores, its effects would be exacerbated by the poor preparedness of legacy buildings (Maio et al., 2017). Unfortunately, assessing the frequency of large earthquakes is difficult because paleoseismological methods are not easily applied to volcanic sequences and many of the faults lie underwater. Alternative assessment methods are therefore needed.
High-resolution bathymetric data have revealed evidence of widespread submarine mass wasting along the upper submarine slopes of islands in the Azores (Mitchell et al., 2008Quartau et al., 2012Quartau et al., , 2014Quartau et al., , 2015Ricchi et al., 2020). Unfortunately, those data such as shown in Figures 1b-1e present a difficulty of interpretation. Whereas landslides in continental slopes often occur over decollements parallel to the seabed and create landslide scars that are clearly defined (e.g., ten Brink et al., 2006), this is not the case in the Azores data, which instead show a rugged morphology likely due to past volcanic history, transfer of slope sediments and multiple small-scale landsliding. We address this problem by adopting the general term "slope valley" for topographic depressions without connotations of origin. Within a database of such identified features, we classify them first according to the degree to which they adhere to the classical features of landslides (Varnes, 1978). A second classification is then used to separate them according to whether they appear to form single landslide scars, amphitheaters created by multiple landslide events (e.g., Pratson & Coakley, 1996) or potential retrogressive failures. We then interrogate this database using the different classes depending on the objective being addressed. The results indicate a significant variation in submarine slope landsliding and hence hazard among the islands.

Geological Setting
It is useful to know aspects of Azorean geology in order to assess how frequently volcanic materials are fed to the coasts during eruptions and earthquakes, both of which may affect landslide frequency and geomorphology. The central Azores islands of Faial, Pico, São Jorge, Graciosa, and Terceira lie east of the Mid-Atlantic Ridge (Figure 1). Rift zones, faults on land and submarine ridges commonly lie perpendicular or oblique to Nubia-Eurasia plate motion suggesting that they were produced by trans-tensional deformation, continuing through the Quaternary to the present-day (Lourenço et al., 1998;Madeira & Brum da Silveira, 2003). Recent 40 Ar/ 39 Ar and K/Ar dating results suggest that Terceira is generally young, with the earliest-formed edifice in the east younger than 0.4 Ma and the most recent edifice in the west ∼50 ka (Calvert et al., 2006;Hildenbrand et al., 2014). São Jorge Island can be subdivided into old and young volcanic complexes (Figure 1). The oldest complex in the southeast is 1. to young volcanic complex in the northwest is younger than 0.75 Ma based on 40 Ar/ 39 Ar and K/Ar dating (Hildenbrand et al., 2008;Marques et al., 2018;Pinto Ribeiro et al., 2010). The northwest of São Jorge also has had historical eruptions in 1580 and 1808 (Madeira, 1998). Although K/Ar dates suggest the maximum age for Faial Island is 0.85 Ma (Hildenbrand et al., 2012), most parts of the island are widely blanketed by a 10-440 ka volcanic complex (Pacheco, 2002;Pimentel et al., 2015). The youngest part of Faial in the west (<10 ka) also experienced the historical Capelinhos eruption in 1957/58 (Machado et al., 1962). K/Ar dates suggest an oldest age of 0.19 Ma for Pico Island (Costa et al., 2015), which has had three historical eruptions since 1562. Lava flows of Pico volcano in the west of the island originated from the summit area and some from flank sources. The flows include both pahoehoe and 'a'a lava types. In contrast, 'a'a lava flows are more common in eastern Pico ( Figure 1; Scarth & Tanguy, 2001).
Neotectonic structures (Figure 2) can be explained by regional trans-tension (Lourenço et al., 1998). Active strike-slip faults on São Jorge and Pico islands are revealed by sheared cinder cones (Madeira & Brum da Silveira, 2003). Some Holocene lava flows have been displaced by faults , which is evidence of pre-historic seismic activity. From 1522 to 1964, at least 10 destructive earthquakes occurred in the central Azores ( Figure 2) responsible for hundreds of deaths (Gaspar et al., 2015). Since 1980, seismometers have been installed on several Azores islands (CIVISA, 1998), greatly improving earthquake detection threshold and location accuracy (Carvalho et al., 2001). At least 30 earthquakes of Mw > 4.5 were recorded ( Figure 2).  Ryan et al. (2009). Island elevations above sea-level are shown in gray with black outlines. Red circles locate the sediment cores analyzed by Chang et al. (2021). Inset shows the regional tectonic setting of the area in the diffuse boundary between the Eurasia (EU) and Nubia ( (Gaspar et al., 2015, their Table 4-1). Red and black lines represent active and inferred faults on land, respectively . F, P, SJ, and T are Faial, Pico, São Jorge, and Terceira, respectively. The deposits forming the upper submarine slopes have not been directly observed or sampled, but their origins can be speculated on. Quartau et al. (2012Quartau et al. ( , 2015 showed that the shelves of the Azores islands have recently been accumulating volcaniclastic particles from subaerial and coastal erosion and carbonate particles from in situ biogenic production. Such sediments likely accumulated since the Last Glacial Maximum and are generally only 10 ms thick, although a seismic reflection image in Marques et al. (2018, their Figure  8) reveals an unusually thick 250 ms (∼225 m) of sediment in the outer shelf of northeast São Jorge. During sea-level lowerings, such sediments are likely to be remobilized by coastal and subaerial erosion, and some of those sediments redeposited in the upper submarine slopes. Sonar images of the shelves of Pico and Faial reveal Holocene dendritic lava flows, lobes and deltas (Mitchell et al., 2008;Quartau et al., 2012Quartau et al., , 2015. During sea-level low-stands, seaward extending shorelines probably allowed such lava flows to reach the outer shelves and hence feed the uppermost slopes. Sedimentary particles released by erosion by streams and waves also likely fed the upper submarine slopes of the islands at such times. However, during sea-level high-stands, sediment can also be produced and exported from shelves due to greater biogenic production (Droxler & Schlager, 1985). It is presently unclear exactly how the flux of sediment to the island slopes has varied over time.
Offshore sediment gravity cores have been collected around these islands, which contain turbidites likely originating from slope landslides such as those we describe here. Analyses of four of those cores located in Figure 1 (Chang et al., 2021) has revealed that those turbidites comprise mainly silt to sand grade volcaniclastic and minor bioclastic particles. Glass shard analysis and other data suggest that more than half of the turbidites originated as primary volcaniclastic events (i.e., pyroclastic flows). Much of the uppermost slopes of the islands may therefore also include such volcaniclastic deposits.

Research Materials
High-resolution bathymetric data were collected on RV Arquipélago in 2003 and RV l'Atalante in 2011, providing the coverage outlined in Figure 1. The Arquipélago survey was carried out with a portable multibeam sonar (Mitchell et al., 2008) and covered the submarine slopes of Faial, Pico and São Jorge islands. These data suffer from imperfect rigidity of the transducer mounting, which led to an artificial high-frequency cross-track corrugation of the data, which needs to be borne in mind in interpretation (e.g., in the class C example in Figure 2). The 2011 survey was carried out during project Features of Azores and Italian Volcanic Islands (FAIVI). Bathymetry data were collected around Terceira Island, with a survey boat deployed into coastal areas. Those combined data have a resolution of ∼1 m in the first 100 m water depth increasing to 50 m at 2,000 m depths (Casalbore et al., 2015;Chiocci et al., 2013;Quartau et al., 2014).
The earthquake event data were acquired from the International Seismological Centre (ISC) un-reviewed catalog for the years 1964-2019. Figure 2 shows the 30 largest earthquakes with Mw > 4.5. Historical earthquakes from 1522 until 1964 are from Gaspar et al. (2015) and the locations and lengths of active and inferred faults on land are from Madeira et al. (2015).

Landslide Identification and Calculations of Area and Volume
An idealized landslide comprises both an embayment or scar created by removed material and an emplaced mound of deposits, sometimes also with a chute in between (Clare et al., 2018;Hampton et al., 1996). Slope valleys that are potential landslide scars were identified by inspecting multiple perspective views of the shaded bathymetric data and seafloor gradient maps, also using tightly spaced contour lines as an indicator of steep gradients. Given the varied data resolutions and poor imaging at deep water depths where the deposits lie, mapping complete landslides was considered to be unrealistic. Therefore, our mapping focused on the heads of landslides only. Even with this focus, subjective decisions were necessary because some headwall escarpments are variably reworked or obscured by later sedimentation. Moreover, headwalls and boundaries between adjacent valleys can become obscured by subsequent slope failures and some escarpments appeared likely formed from multiple events.
The slope valleys were classified in two ways ( Figure 3; see detailed mapping results in Supporting Information S1). The first classification addresses the varied level of recognition of these features as landslides and hence conveys our confidence in their landslide origins. The second classification addresses whether their morphologies suggest they formed by single slope failure events or by multiple events (either subsequent different events or retrogressive failures occurring over finite periods). For the recognition classification, five classes (A-E) were chosen as follows ( Figure 3a). Class A: headwalls and sidewalls can be clearly identified by sharp gradient changes and commonly with significant depth changes. Class B: headwalls and sidewalls mostly can be identified by sharp gradient changes and commonly with significant depth changes, though parts of them are less well defined. Class C: headwalls can be mostly identified from gradient changes but sidewalls are obscured or missing. Class D: headwalls and sidewalls are delineated based on modest gradient and depth changes (valleys belonging to this class commonly still have a depression). Class E: headwalls and sidewalls are delineated based on mild gradient and depth changes.
Examples illustrating our second classification scheme are shown in Figure 3b. The slope valleys were subdivided into three types (1-3) based on their morphometric features. Type 1 valleys were interpreted as produced by multiple failure events. In this cases, multiple-branching valleys were identified and each of them appears important. Later erosive features can also overprint the main rupture surface. Type 2 valleys were interpreted as each produced by a single major failure event. Each valley contains a prominent headscarp. Type 3 valleys lie above the headscarps of type 1 or 2 valleys, so we interpreted them as the result of retrogressive failures. They have gentler gradients and form depressions of the sediment surface. Although those depressions could be artifacts, such as due to sediment aggradation around them or other factors, landsliding has occurred in similar shallow gradients and layered carbonate-rich deposits (Schwab et al., 1988), so we suggest these were potentially also earthquake-triggered failures. Type 3 valleys typically occur in the outermost shelves of the islands, which are fully mapped in the multibeam data.
The areas and volumes of valleys were measured following the methods of ten Brink et al. (2006). We first manually digitized the boundary of each valley encompassing the potentially collapsed region to acquire valley area. We then created a smooth upper surface by interpolating over the perimeter of each digitized polygon. The present-day bathymetry (lower surface) was then subtracted from the interpolated surface (upper surface) to acquire valley volume. As the investigated submarine slope area varies amongst the islands and hence also affects the abundances of landslide valleys, we normalized landslide valley areas and volumes by dividing them by the submarine slope areas.
These measurements have the following potential uncertainties.
(1) Later sediments deposited in slope valleys bias volume estimates.
(2) Further erosion may be possible from sediment movements during or shortly following each slope failure. (3) Oblique artifacts due to the non-rigid mounting of the 2003 data can influence the topography of headscarps. Uncertainty (3) is minimized, as most slope valleys were located in the upper slopes. We mitigate against uncertainty (1) when assessing tsunami risks by only studying the more distinct landslide scars (classes A-C and type 2). That uncertainty may affect our inter-island comparisons of landslide volume, although we also check variations using scars with prominent scarps. Uncertainty (2) can potentially also affect the landslide volumes (e.g., Sun et al., 2018). However, the variations in landslide volumes can be verified if similar variations also occur in the landslide areas, which are less affected by either depositional or erosional modifications. Furthermore, our results below suggest that the material has significant cohesion and is thus more resistant to erosions by small sedimentary flows in the uppermost island slopes.

Landslide Size-Frequency Analysis
The complementary cumulative distribution function (CCDF) has been used to represent the probability that a landslide greater than a particular size will occur and is helpful for summarizing the distribution of landslide sizes as it emphases the larger features (e.g., Malamud et al., 2004). CCDFs of area and volume were generated from slope valleys best representing single-event landslide scars (classes A-C and type 2). Their correspondence with models (goodness-of-fit) was tested using Kolmogorov-Smirnov test (K-S test) and the p-value suggested by Clauset et al. (2009). The K-S test is a nonparametric test of the equality of continuous or discontinuous variations. It is used to quantify an interval between the empirical cumulative distribution function of the sample and the reference distribution. The p-value can vary from 0 to 1. If it is higher than 0.1, the observed data fit the tested distribution, whereas a p-value equal to or less than 0.1 suggests that the data are unlikely to follow the distribution. These tests have been applied in other submarine landslide studies (Casas et al., 2016;Geist & ten Brink, 2019).

Estimates of Peak Horizontal Accelerations During Earthquakes
Earthquake-induced ground motions have been widely suggested to trigger landslides (Fine et al., 2005;Gràcia et al., 2003;Heezen & Ewing, 1952;Tappin et al., 2001). Earthquakes with Ms > 6.5, which can create a local peak horizontal acceleration (PHA) >0.5 g, occur every 70 years on average in the Azores (Nunes & Ribeiro, 2001;. To compare with the landslide data, we estimated PHA as follows. (1) Based on an empirical study of ground motion recordings in western North America, Boore et al. (1997) suggested the following dependence of PHA (Y in fraction of g) on Mw and closest horizontal distance from source (r jb in km).
and V A are determined by nonlinear regression (Joyner & Boore, 1994), which in this study were derived from Boore et al. (1997, their Table 8) corresponding with period 0.75 s as suggested by ten Brink et al. (2009). The 200 m/s value used here for Vs is a representative average for 0-30 m depth in marine sediments from data compiled by Hamilton (1976). Equation 4 was derived from larger magnitude earthquakes than we study here and represents effects of attenuation and scattering over continental crust, so the derived PHA could be biased, though we use it here as it is based on a large amount of strong-motion data and lack of suitable alternatives.

Sediment Cohesion Estimates
Geotechnical data are not available to assess the stability of these slope sediments. However, remarkably steep deposits in headwalls and, in places, chutes of landslides (Figures 1b-1e) around the volcanic islands imply that the slope sediments are cohesive and that in turn has implications for how we interpret them. In contrast, if the slopes were instead purely incohesive sediment, they would be expected to have gradients no greater than the angle of repose of ∼30°. We estimated cohesion by inverting the conventional equations used in static limit-equilibrium analysis using the seabed gradients found. Such analysis is still widely applied due to its simplicity and ease of usage but tends to be conservative (Jibson, 2011), implying that cohesion estimates based on it will be minima. According to that analysis, failure occurs when the factor of safety (FS) (the ratio of resisting to driving stresses on the potential failure surface) is lower than 1. We use the simplified equations for failure on a surface parallel to the dipping seabed (Figure 4), assuming the structure extends infinitely. FS is then simply given by: where c = sediment cohesion, H = average thickness of failure, γ' = unit weight of submarine deposit in water, ρ s = sediment wet bulk density and ρ w = water density, g = gravitational acceleration, θ f = the angle of failure, and Φ = sediment friction angle. Using parameters given later, we estimated c as the value satisfying FS = 1. As earthquake shaking was not considered but some of these steep sediments may have survived such shaking, this analysis provides a minimum value for c.

Estimates of Tsunami Height
Coastal runup heights of tsunamis are difficult to reconstruct without sophisticated numerical modeling.
There is also the uncertainty over whether the valleys were created by single or multiple landslides (Hunt et al., 2011(Hunt et al., , 2013. Nevertheless, Watts et al. (2003) derived analytical expressions for the initial waves above the landslides that are potentially useful for applying to the Azorean valleys likely produced by single events. Those expressions were derived from laboratory experiments simulating landslides that are large along-slope compared with the tsunami wavelength. McAdoo and Watts (2004) updated the expressions to allow for landslides that are narrower along-slope than the tsunami wavelength. In their model, the sliding body was idealized as elliptical with specific gravity 1.85. With further corrections by De Lange and Moon (2004), this wave amplitude at source for translational slides is:   (Figure 4) T = head scarp height, w = along-slope landslide width, λ = tsunami wavelength, θ = average slope angle, b = initial length of the slide measured downslope, and d = depth of landslide initial center of mass. Sediment density has not been measured here. If, based on the offshore cores, our suspicion that the sediment comprises detrital particles is correct, a wet bulk density of ∼2.0 g/cm 3 is appropriate (Hamilton & Bachman, 1982) so the specific gravity of 1.85 should not be far in error. Spreading during wave propagation toward the coasts can reduce wave amplitudes, whereas shoaling typically increases wave amplitudes. Local bathymetry can refract waves, leading to either convergence or divergence. Depending on direction, the incoming waves can also interfere with wave fronts emanating from the landslide deposit area (Geist et al., 2009). Local wave impacts are therefore likely to vary greatly around coasts and differ from these values of A. Nevertheless, Casalbore et al. (2011) found that tsunami runup heights in the Tyrrhenian Sea roughly compared with those predicted by Equation 6. Here, the estimated values of A were used mainly to provide a rough sense of relative scale, in a similar way to how earthquake magnitudes are used to assess relative size.

Seismic Ground Shaking
The map of PHA ( Figure 5) calculated from the ISC data ( Figure 2) suggests that the submarine slopes of the four islands studied here have all experienced >0.05 g and locally >0.4 g since 1964. Three areas of high acceleration occur around the slopes of Faial and São Jorge associated with the largest earthquakes in the ISC catalog. However, if our suspicion is correct that the landslides have formed over the Holocene, all of these slopes may have experienced as much as 0.5 g PHA over that timescale for at least two reasons. First, the historical (1522-1964) record includes 10 large events (Figure 2), many of which occurred on Terceira where seismicity since 1964 has been modest ( Figure 2). This suggests a discrepancy with the instrumental data. Although magnitudes of such old events are difficult to estimate accurately, they are likely to have been M > 6.5 based on magnitude assessments of historical records elsewhere (Nunes & Ribeiro, 2001;. Second, active faults are widespread amongst the islands. Madeira et al. (2015) estimated that >20-km-long faults on the islands could produce events of M = 6.5 if they ruptured along their whole lengths. However, mapping the full extents of active faults is difficult on land, where they can be masked by subsequent deposits . Identifying active faults between the islands and continuations of faults from land underwater is another challenge. Deep-tow sonar images (Mitchell et al., 2018) show many fault escarpments that lack erosional features seen in older inactive faults (Tucholke et al., 1997), and the FAIVI data also show many "fresh" submarine fault scarps around Terceira (Casalbore et al., 2015;Chiocci et al., 2013;Quartau et al., 2014). We thus suspect that the submarine areas host active faults, as also implied by the broad distribution of seismicity. Many of these faults have greater lengths than those mapped on land by Madeira et al. (2015) so they likely generate earthquakes larger than the M = 6.5 that Madeira et al. (2015) suggested as a maximum. Consequently, the pattern of seismic accelerations over the long-term remains uncertain.

Submarine Landslide Inventory and Volume-Frequency Distribution
We mapped 1,227 slope valleys in total. The counts of the different classes and types are shown in Figures 6a and 6b. Volumes range from 10 2 to 10 8 m 3 with sample mean 7.9 × 10 5 m 3 . The volume distribution for all valleys (Figure 6c) has a logarithmic bell shape. The volume distributions for individual islands (Figure 6d) differ from each other. The modes for Pico and Faial (10 4 -10 5 m 3 ) are smaller than those for São Jorge and T SJ P F Terceira (10 5 -10 6 m 3 ), but the distributions all have similar bell shapes. Figure 6d also reveals that São Jorge and Terceira have higher percentages (20%-25%) of the larger valleys (V > 10 6 m 3 ) than Faial and Pico islands (5%-10%). The largest 10 valleys (all types) all occur around São Jorge and Terceira. São Jorge and Terceira also have the largest valley cumulative area and volume (Figures 7a and 7c), an observation that persists when the data are normalized for slope area (Figures 7b and 7d).
The relationship between area (A l ) and volume (V) for valleys most likely to be single landslide scars (type 2 and classes A-C) is shown in Figure 8. Regressing the logarithmic variables suggests  1.3365 0.234 l E V A (R 2 = 0.924). Mean thicknesses were also obtained by dividing individual landslide scar volume by area (also for type 2 and classes A-C valleys). The mean thicknesses are 9.2 ± 0.5, 7.3 ± 0.5, 12.6 ± 0.8, and 10.6 ± 1.3 m for Faial, Pico, São Jorge, and Terceira, respectively (uncertainties of means are 1σ). CCDFs of area and volume were also generated from these valleys (type 2 and classes A-C) as shown in Figure 9. Both CCDFs follow log-normal distributions (p = 0.98 and 1.0).

Assessment of Cohesions of Steep Deposits
Steep gradients in some headwalls and sidewalls exceed 60° locally and generally reach 30°-40° (Figures 1b  and 1c). Using these slope gradient values for the failure gradient angle, we inverted Equation 5 for cohesion c. For H, we used the average thickness (10.6 m) of the scars associated with single landslides (type 2 and classes A-C) found by dividing their volumes by their areas. A bulk density of 1.0 g/cm 3 was used for (ρ s -ρ w ) representing sand in water (Hamilton & Bachman, 1982;Hamilton et al., 1956). We used 25°-40° for internal friction angle Φ typical of volcaniclastic sands and calcareous ooze (Boldini et al., 2009;Lee et al., 1994). Repeating the calculation for gradients of 30°-60°, we estimated c between 9 and 33 kPa ( Figure 10). However, these c-values are minima because the infinite-slope calculation is conservative. Furthermore, if some of these steep sediment deposits have survived shaking by earthquakes, such accelerations have not been The upper-slope sediments lie in water that is saturated with respect to aragonite and calcite (Wisshak et al., 2010(Wisshak et al., , 2015, so cementation may be involved in stabilizing the sediments. Carbonates are precipitated from seawater on the edges of clasts (Gabitov et al., 2019), binding the sediments (e.g., Tucker et al., 2020).
These cohesion values ( Figure 10) noticeably overlap with laboratory results of Nafisi et al. (2020), who formed calcium carbonate cements in sand by microbial precipitation and tested the results under 10-100 kPa effective stresses. They showed that cohesion was 9-12 kPa after moderate cementation (1-4 wt% CaCO 3 ) and 56-65 kPa after high cementation (4-10 wt% CaCO 3 ). Our 9-33 kPa estimates of sediment cohesion therefore correspond with moderate to high cementation.
These cohesion values are also compatible with test results on some volcaniclastic deposits (e.g., lava-breccia ∼30 kPa) found in volcanic islands as here but lower than the ∼600 kPa of pyroclastic deposits (Di Traglia et al., 2018). Volcaniclastic deposits lying at up to 40° were also found off the Sciara del Fuoco slope of Stromboli . Some of these cohesions could therefore be compatible locally with lava deltas (Mitchell et al., 2008) underlying parts of these slopes. The cohesions we estimated are logically smaller than found in intact rocks (c > 25,000 kPa; Bieniawski, 1975;Marinos & Hoek, 2000).

Assessment of Tsunami Wave Heights at Source
Tsunami heights A immediately above landslides were derived for 83 single-failure landslides (type 2) with V > 10 6 m 3 using Equation 6 (detailed morphometric data in Supporting Information S1). The CCDF of wave height in Figure 11a also follows a log-normal model (p = 0.43). We define an initial wave height at source >1 m as potentially hazardous as such a wave could be amplified to several meters when inundating coasts  Area (m 2 ) l (e.g., Tinti et al., 2005). At least 18 landslides were found capable of generating tsunamis with A > 1.0 m, and one with A > 7 m ( Figure 11). These landslides were located around São Jorge, Terceira and Faial Islands.
Sophisticated modeling of landslide-triggered tsunami wave propagation is beyond the scope of this paper. Coastal topography, refraction and wave steepening with shallowing water all affect final runup height (tides are a secondary consideration here as they are only ∼1 m [Cruz & Silva, 2001]). However, a first-order estimate of tsunami height at the coast can be made using Green's law (e.g., Federici et al., 2006), which predicts runup heights relative to source heights. In our study area, the shoaling factor to 1 m water depth is between 3 and 5, suggesting that there is an underappreciated risk of tsunamis around volcanic islands such as the Azores.

Assessing Varied Submarine Landslide Features in the Central Azores
The above analysis has revealed lower cumulative volumes and areas of slope valleys around Faial and Pico islands than around São Jorge and Terceira islands (Figure 7). The density of landslides around São Jorge is comparable with those of Faial and Pico, so landslides around São Jorge are generally larger than those around the other two islands. Accounting for the low density but high affected slope volume of landslides in Terceira Island, the landslides there are large as well. As might be expected from their different average volumes, therefore, landslides are modestly thicker on average around São Jorge than around Faial and Pico. We explore these contrasts in this section.

Earthquake Triggering
ten Brink et al. (2016) suggested that slopes become more stable with decreasing sedimentation rate and increasing frequency of earthquakes. Such increased sediment shear strength arises due to seismic shaking, causing shear-induced compaction (e.g., Sawyer & DeVore, 2015;ten Brink et al., 2016) or sediment internal structure changes (Wu et al., 2021). This may help to explain the comparatively high 9-33 kPa estimated cohesion of these Azorean sediments besides carbonate cementation and prompts us to ask if earthquakes may explain other aspects of the landslides data also. To represent the regional occurrences of landslides, we prefer the  cumulative volumes or areas of all slope valleys (Figure 7) as they form the largest database, although the restricted groups (type 2 and classes A-C) also show a similar tendency for volumes and areas of valleys around São Jorge and Terceira to be larger than those around Faial and Pico islands.
The pattern of extreme shaking in Figure 5 is dominated by the three largest earthquakes occurring since 1964, with two near São Jorge and one in northeast Faial. The pattern does not correspond with the landslide volume and area differences between the islands. However, the lack of correspondence could arise if the ISC catalog is unrepresentative of the long period over which these submarine landslides likely developed (>ky). For example, patterns of seismicity can change location over decades due to changes in stress after the largest events (Stein, 1999). The historical record of earthquakes is also subject to uncertainties. The wide distribution of active faults (Figure 2) implies that there is a potential for large earthquakes that is more widespread than suggested by the simplistic reading of either instrumental or historical records.
Sediments accumulating on the upper island slopes originate from various sources; from coastal and subaerial erosion, biological production and volcanic eruptions. Their fluxes are hard to quantify and most likely vary around the islands. However, a simple explanation for the larger landslides around São Jorge and Terceira compared with Faial and Pico could involve a larger time interval between destructive earthquakes, leaving more time for deposits to accumulate. This in turn implies that there could be a higher frequency of larger earthquakes around Faial and Pico islands and thus a greater longer-term threat of earthquakes there. This is an important inference concerning the seismic hazards of the Azores islands. We however temper this inference by assessing other factors below that may have also affected the landslide volumes and areas.

Varied Lava Flow Types
Lava flows on east Pico island are mainly 'a'a whereas those of west Pico are more mixed pahoehoe and 'a'a types (Figure 1). Pahoehoe flows typically form stacks of thin sheets whereas 'a'a lava flows are commonly thicker (Macdonald, 1953). Although 'a'a flows can have friable clinkers above and below them, their interiors are commonly more massive and less porous or fractured (Macdonald, 1953). During sea-level low-stands, lava flows would have more easily reached the outer shelves, supplying the uppermost slopes directly. Although we anticipate there being mainly sediments underlying the uppermost slopes because turbidites in offshore cores contain mainly silt-sand grade volcaniclastic and bioclastic particles (Chang et al., 2021), different landslide statistics may occur if there were widespread lava underlying the upper slopes because of differences in geotechnical properties between the two lava flow types.
However, landslide volumes ( Figure 12) differ only slightly between west and east Pico. We suggest that, even if lavas are common, volcaniclastic deposits of 'a'a and pahoehoe are probably similarly susceptible to failure, rapidly disaggregating on the steep upper submarine slopes (cf. Sansone & Smith, 2006). Furthermore, varied effusion rates of lava flows and the pre-eruption seafloor morphology can also lead to uneven  loading of slopes by lava flows, which ultimately affect slope stability (Bosman et al., 2014;Casalbore et al., 2021), so factors other than merely flow type likely affect landslide occurrences.

Volcanic Ages of Adjacent Land and Shelf Widths
The ages of the volcanic sequences comprising each island might be expected to influence landsliding if more recent volcanic building led to steeper slopes of friable material that were, susceptible to failure. On the other hand, broad shelves typical of older volcanic islands (Menard, 1986) may host thick sediments and those sediments may be exported to the island submarine slopes by wave erosion and/or mass wasting (e.g., Dengler et al., 1984;Fornari et al., 1979). Based on high-resolution bathymetry data, Quartau et al. (2010Quartau et al. ( , 2014Quartau et al. ( , 2015 confirmed that most of the older volcanoes in the Azores have wider shelves than younger edifices, though shelf width can be affected by vertical tectonic motions, varied wave climate and other factors as well (Ramalho et al., 2013). As sediments are continually produced by biogenic production as well as subaerial erosion and wave abrasion of sea cliffs, we might expect the upper slopes of islands with wide shelves to be capable of accumulating thicker sediment deposits during sea-level high-stands. For instance, the easterly older shelf of São Jorge has a broad sedimented terrain where seismic data reveal an extreme ∼250 ms of layered sediment (Marques et al., 2018). However, shelf width is likely not the sole pre-conditioner for slope failures. For instance, Santa Maria, which is the oldest island and has the widest shelf of the Azores, has deposits on its shelf of only a few m thick (Ricchi et al., 2020). This is partly because it is a low relief, semi-arid island and hence contributes little sediment from stream erosion. In addition, wide shelves can attenuate waves crossing them, protecting cliffs from erosion (Ricchi et al., 2020).
Although São Jorge has a higher cumulative volume of slope valleys than either Faial or Pico islands, there are also differences between the volcanically old and young parts of São Jorge. The slope valleys in old southeast São Jorge have nearly twice the volume per unit slope area than those in the young northwest part of the island (Figure 12). This favors an explanation involving export of shelf sediment from the broad shelf in southeast São Jorge, leading to thick slope deposits susceptible to failure. Sediment production rates are difficult to anticipate, given that the sediments originate from varied sources such as biogenic particles created on the shelves, subaerial and coastal erosion, and pyroclastic fallout from eruptions (Quartau et al., 2012), though collecting sediment cores and other data around the islands could in principle help resolve this in the future. Nevertheless, the thick sediment imaged seismically beneath the shelf of São Jorge (Marques et al., 2018) suggests this could indeed partly explain its high cumulative volume of slope valleys. Despite this effect, the cumulative volumes of landslide valleys of west São Jorge (∼1 m 3 /m 2 ), where the island is volcanically young and the shelf is narrow, is still larger than the volumes for Faial and Pico, which are both <0.5 m 3 /m 2 , so the difference between the islands still needs an explanation, such as the difference in long-term seismicity.

Comparing the Azores Landslide Area-Volume Relationship With Those of Other Settings
The relationship between volume (V) and area (A l ) is generally written as    l E V A where α and β are constant parameters for each data set. Landslides that are perfectly self-similar or show no systematic variation in ratio between the vertical and horizontal dimensions should uniformly have β = 1.5 (Guzzetti et al., 2009;Klar et al., 2011). With accurate landslide dimensional measurements, Guzzetti et al. (2009) reported β = 1.449 for subaerial landslides in central Italy, which is close to β = 1.5. Table 1 shows α and β values for six studies of data from different submarine geological environments for comparison with the Azores results. ten Brink et al. (2006)   homogeneously deposited, though likely having layered physical properties. Failure along seabed-parallel layers tends to form tabular landslides. A Mediterranean result (Urgeles & Camerlenghi, 2013) has intermediate β = 1.251. The Mediterranean slopes have varied geological structure and sedimentary inputs, leading to diverse landslides types and thickness.
The Azores data have intermediate β = 1.337, though closer to a size-invariant ratio (β = 1.5). Landslides have β < 1.5, perhaps because of later infilling of valleys by sediments or resistance to deep failure (e.g., because of increasing sediment compaction or harder volcanic rocks at depth). Figure 13 shows CCDFs of landslide volume from seven different submarine settings. CCDFs are used to assess the maximum size in a region and investigate the relative occurrences of features of differing size, though are affected by completeness and other factors. Substantial sediments can accumulate on open passive continental margins as the slopes are broad and homogenous, with only occasional earthquakes. However, once landslides occur, movements can transport sediment far from their origins so landslides can be large, reflected in their CCDFs (ten Brink et al., 2006). Though landslides in tectonically active areas are usually small (Urgeles & Camerlenghi, 2013), some giant submarine landslides still can be produced there (e.g., large slumps in Hawaii; Masson et al., 2002).  Small landslides (<10 5 m 3 ) tend to be identified more readily in higher-resolution data sets. Hence, at least 10 times more landslides were identified in bathymetry grids with 20-m cell sizes than in those with 100-m cell sizes ( Figure 13). Comparing data sets of similar grid cell sizes, some differences occur that could be explained by differing sediment physical properties. For instance, the smallest mapped scars in the Azores have volumes 5 times larger than the smallest scars in the southern Tyrrhenian Sea, despite comparable grid cell sizes (<20 m). One possible explanation could be different sediment cohesions. Limit-equilibrium slope stability modeling has suggested that incohesive sediments tend to form shallow, elongate landslides with no lower size limit, whereas cohesive sediments tend to develop landslides with a minimum size (Frattini & Crosta, 2013). The 9-33 kPa or larger cohesion of the Azores sediments may explain the larger minimum volumes there compared with those of the Tyrrhenian Sea. However, the variance of geological setting needs to be considered as well, as most of the small landslide scars in the Tyrrhenian Sea were found at the heads or flanks of submarine canyons. Confined spaces there may also have limited landslide development (Casas et al., 2016).

Lessons for Assessing Geohazards of Volcanic Islands
We have suggested that differences in the cumulative volumes and areas of submarine landslides may imply differences in long-term seismicity. Such differences could potentially help to overcome the inadequacy of the instrumental seismic record (Scholz, 2002). Hence, submarine landslide mapping can help in assessing earthquake hazards as well as tsunami hazards. This analysis could easily be repeated around other volcanic islands with suitable bathymetric data, although effort is still needed to assess the other possible causes of varied landslide properties identified above, for example, fluxes of sediment from erosion and biogenic production.
Geotechnical data would help evaluating slope responses to earthquakes. Obtaining samples for such tests by vibrocoring would be challenging on the steep (30°) slopes of islands, but more indirect assessments of sediment rigidity may be possible, for example, velocities from shallow seismic refraction.
More efforts could be put into assessing landslide frequency. Long-term records can be obtained by dating landslide-origin turbidites (e.g., Hunt et al., 2013). Short-term records could be obtained by repeating multibeam sonar surveys (e.g., Casalbore et al., 2012Casalbore et al., , 2020Kelner et al., 2016;Soule et al., 2021), with timings of events obtained from tide gauge records and acoustic recordings (Caplan-Auerback et al., 2001) to help identify possible formative events such as earthquakes.

Conclusions
We have mapped 1,227 submarine slope valleys in the central Azores volcanic islands. To overcome difficulties of interpretation arising from the rugged morphology of the upper submarine slopes, valleys were first categorized based on the levels of recognition as landslides and whether they appear to have formed by single, multiple or retrogressive failures events.
Considering slope valleys of all types and classes, São Jorge and Terceira islands have greater valley volumes and areas per unit slope area, compared with Faial and Pico. We highlight this observation as suggesting that Faial and Pico potentially have greater earthquake hazard. In this interpretation, frequent earthquakes prevent the build-up of unstable sediment deposits on slopes, leading to mostly smaller landslides around Faial and Pico. This suggestion is tempered by an observed greater valley volume in easterly São Jorge where thick sediments are also observed on its shelf. Such sediments are likely exported to slopes during sea-level lowstands, which suggests that greater sediment accumulation has also affected landslide volumes there. Nevertheless, westerly São Jorge has a narrower shelf and its valley volume is still greater than Faial and Pico islands. Thus, differences in factors other than sediment input such as seismicity are still needed to explain the contrast in slope valleys statistics.
Cumulative area and volume distributions for valleys most likely to be single landslides fit log-normal probability density functions, as found elsewhere. Based on an analytical formula, at least 13 of the valleys most likely to have been produced by single landslides would have generated tsunamis with heights >1 m at source. One may have produced a wave of ∼7 m. Those heights also follow a log-normal probability density function.
Static slope stability analyses suggest that some steep slopes have cohesive strengths of at least 9-33 kPa, much larger than 0.08-2 kPa expected of typical superficial incohesive seafloor sediment. This could be explained by moderate to high carbonate cementation, by earthquakes shaking and/or perhaps by the presence of coherent lava or talus.
Overall the study suggests that mapping submarine landslides around volcanic islands in general could also be useful for investigating seismic and tsunami hazards that are not well characterized by other methods, as well as uncovering other aspects of their submarine slope sediments.

Data Availability Statement
The 2003 multibeam data are available at 100 m grid resolution from the Marine Geoscience Data Portal (http://www.marine-geo.org/index.php). The multibeam data for Terceira are not publicly available due to the sensitive nature of the area and environmental sustainability concerns, but can be made available to researchers with appropriate credentials. Background bathymetric data are from the GeoMapApp (GMRT; www.geomapapp.org, Ryan et al., 2009). The seismic catalog is from International Seismological Centre (2019) (https://doi.org/10.31905/EJ3B5LV6; Di Giacomo et al., 2014). Yu-Chun Chang thanks the editor Peter van der Beek and reviewers Daniel Casalbore and Uri ten Brink for their constructive comments, which led to significant improvements of this manuscript. He also thanks the Taiwanese government for providing his PhD scholarship and the International Association of Sedimentologists for a travel grant to the annual British Sedimentological Research Group meeting, where some of these ideas were explored. The European Communities 7th Framework Programme under EUROFLEETS funded the bathymetric data collection around Terceira used in this study. Iasmina Cristina Popa is thanked for a preliminary attempt at the peak ground horizontal analysis. Thor Hansteen, Christoph Beier, Adriano Pimentel, Ting-Wei Wu, Sung-Ping Chang, and Zhongwei Zhao are thanked for informative discussions. Figures were prepared with the GMT software system (Wessel & Smith, 1991).