Out of the 2652 Holocene and Pleistocene volcanoes in the Global Volcanism Program (2023) Database, we identified 370 edifices as volcanic islands. In this study, we investigated 47 of these in three groups.

Group A consisted of 26 volcanic islands in 4 archipelagos (Aeolian Islands, Canary Islands, Hawaii, and South Sandwich Islands) that were almost completely covered by high-resolution bathymetry data ( Table 1 ). The only unmapped part in the archipelagos of Group A was the coastal white ribbon, which has proven to be a major challenge to investigate (Weymer et al., 2022). The islands were chosen for their high data coverage, as well as for representing a variety of volcanic island types and sizes. The ship-based high-resolution bathymetry for the Aeolian Islands, the Canaries, and Hawaii were downloaded from the online platform of the Global Multi-Resolution Topography (GMRT) (Ryan et al., 2009), at the highest available resolution, which was between 120m and 240m grid cell size. Besides this global compilation, several regional platforms such as EMODnet (Europe) or JAMSTEC (Japan) as well as bigger more general platforms such as PANGAEA provide bathymetric data. We chose to use only the high-resolution data (masked) and interpolated the coastal white ribbon for these datasets. The data from the South Sandwich Islands was provided by Fretwell (2015). This is a data compilation (200m resolution) of various sources, including GEBCO and, below 60°South, the International Bathymetric Chart of the Southern Ocean (IBCSO), where direct measurements were not available.

In addition, we investigated 18 volcanic islands from the Marianas and the Galapagos Islands, which form group B. The islands of these archipelagos were only partly covered with high-resolution bathymetry. To study the effect straightforward interpolation has on the geomorphological parameters we investigated, compared to using GEBCO in data gaps (unmasked), we used both versions (unmasked and interpolated) grids of the Marianas. These, along with the unmasked bathymetry of the Galapagos Islands, were also downloaded from GMRT.

The GEBCO grid is a combination of direct measurements and indirect measurements such as predicted bathymetry from satellite altimetry (Weatherall et al., 2015; Tozer et al., 2019). For the volcanic islands of Group A and B, a major part of the grid covering the islands is based on direct measurements. To isolate the difference between indirect measurements and ship-based high-resolution bathymetry, three volcanic islands from Vanuatu (Group C) were included. The high-resolution bathymetry was provided by Daniel (1993) for these islands, however, the GEBCO grid in the same area is based on indirect measurements only.

Several automatic algorithms for the delineation of volcanic edifice outlines have been published with the intent to decrease the time and subjectivity a manual delineation is subject to, such as the Matlab-based program MBOA (Bohnenstiehl et al., 2012) or the ENVI IDL-based NetVolc (Euillades et al., 2013). Nevertheless, manual delineation is often less time-consuming and less sensitive to artefacts in the DEM. Thus, morphometric studies of volcanic edifices are often performed with manual delineations. In an attempt to decrease the subjectivity of this critical step, we modified a semi-automated approach by Van Wees et al. (2021), where the 3° slope contour line is used as a starting delineation and is manually edited where necessary. To make the delineation more robust and comparable across different volcanic edifices and less sensitive to artefacts, we resampled all datasets in this study to 1000m grid cell size for the delineation. For the parameter calculation, the unaltered DEM is used. For isolated islands on a flat seafloor and good data coverage, the initial delineation did not need to be modified considerably. For the remaining islands, the three main challenges were a) islands whose edifices share a plateau and need to be delineated separately, b) undulating terrain at the volcanic base due to spreading processes or sediment waves that result in a complex wavy pattern of contour lines and c) other complex terrain at the surrounding seafloor that results in a puzzle of shorter contour lines, rather than a single closed one. Several derivatives of the DEM (hillshade, slope, aspect etc.) were used as visual guides for the manual editing. For the separation of islands, catchment delineations proved the most useful. Of the 47 volcanic islands of this study, 8 were left with minimal manual edits, for 18 the only edits were the separation of edifices, and 21 require > 50% of the outline to be edited manually.

Using the software QGIS and Python, we created an algorithm that is able to extract the morphometric parameters described in Table 2 and in the following section. Though the MorVolc algorithm created by Grosse et al. (2012) was unavailable to us at the time of processing, we did use a modified version of their list of geomorphometric parameters. As we were primarily interested in the shoreline-crossing and submarine characteristics of the volcanic islands, we focus on the basic size and shape and omit any information about the summit region that are important in a terrestrial environment. In addition to the parameters we extracted for the entire edifice, we looked at the subaerial and submarine parts separately. The inputs were the interpolated DEM and the edifice boundary.

The previous steps were taken separately for each dataset mentioned in Table 1 , resulting in a vector file for each volcanic island and DEM. The results were summarized in *.csv files and further processed in Python. We then calculated the deviation (D) of the results of the parameter extraction using GEBCO for both steps (G) from those of the high-resolution data (B) in percent:

D is negative when GEBCO underestimates values compared to the high-resolution bathymetry, and positive when overestimating. Additionally, we produced a table of parameters using the GEBCO DEM with the delineation derived from the high-resolution bathymetry (C). This enabled us to calculate the deviations of the results produced by the difference of the edifice boundary and the DEM used in parameter extraction separately, using the following equations:

In addition to the comparison of the parameters, we calculated the Jaccard Index (aka Tanimoto Index), a measure of the similarity of sample sets in statistics, or in this case, polygons. This will yield information about how different the edifice boundary is, depending on which DEM was used for the delineation. The index is the ratio of Intersection and Union of the polygons enclosed by the edifice boundaries resulting from the delineations with both DEMs, respectively:

Where G and C are the polygons enclosed by the edifice boundaries produced with GEBCO and high resolution bathymetry, respectively. This results in an index ranging from 0 (no overlap) to 1 (perfect overlap).

Figure 3 shows an overview of the parameters calculated from high-resolution bathymetry data for all three data groups investigated. In the figure, the x-axis follows a non-linear sorting of the volcanic islands descending by total height H_tot for all subplots.

Our study covers a great range of sizes of volcanic islands from Genovesa (Galapagos), with 78m the smallest island by elevation above sea level to the great Hawaiian Islands which reach more than 4000m a.s.l. ( Figure 3A ). There is a trend of volcanic islands with higher elevation above sea level also reaching into lower depths (-623m, Floreana, Galapagos; -5864m Maui edifice group, Hawaii), resulting in total heights from 1138m (Pinzon, Galapagos) to 9878m (Hawaii main Islands, Hawaii) from the seafloor to the summit. The distribution of the area is shown in Figure 3B . The area above sea level mostly ranges from 0.5km² (Ka’ula, Hawaii) to 2468km² (Maui edifice group, Hawaii), with Hawaii Island being the exception with 10500km². The area below sea level A_bsl contributes between 79km² (Pinzon, Galapagos) and 25780km² (Hawaii Island, Hawaii), resulting in total areas between 104.8km² (Pinzon, Galapagos) and 36289km² (Hawaii Island, Hawaii).

In Figure 3C the distribution of the volume of the islands is shown. V_asl ranges from 0.2km³ (Ka’ula, Hawaii) to 12000km³ (Hawaii Island, Hawaii) and V_bsl from 25km³ (Pinzon, Galapagos) to 69000km³ (Hawaii Island). This results in total edifice volumes (V_tot ) from 29km³ (Pinzon, Galapagos) up to 81000km³ (Hawaii Island, Hawaii). Additionally, there are two bends visible in Figure 3 , separating the islands into three groups by volume; those larger than Kauai (three islands from the Hawaiian archipelago), those larger than Agrihan (islands from Hawaii and the Canaries), and the remaining islands. The radius parameters are visualized in Figure 3D ). The average radius of the islands calculated from the area varies from 5776m (Pinzon) to 107477m (Hawaii Island). Its values lie between those of R_min (3932m - Pinzon, 66572m Hawaii) and R_max (7957m - Alicudi - Aeolean; 163006m - Hawaii Island). Additionally, we calculated R_A,asl producing values between 415m (Ka’ula) and 57838m (Hawaii) and R_bsl ranging from 2912m to 59844m (Maui edifice group). Some of the smaller islands show significantly larger R_max than expected from their height, producing peaks in the graph in Figure 3D .

Figure 3E shows the distribution of several slope parameters. Unlike the parameters previously mentioned in this section, none of the slope parameters show any trend correlating with the height, the volcanic islands on the x-axis are sorted by. The average slope above sea level Sl_asl,mn is ranging from 4.3° (Santa Cruz, Galapagos) to 13.1° (Wolt, Galapagos). The average across all the islands of 9.6° is 2° higher than that of the average submarine slope Sl_bsl,mn (7.6°), which shows values between 2.7° (Santa Cruz) and 10.9° (Wolt). For the average slope of the entire edifice, Sl_mn ranges between 4.5° (Santa Cruz) and 13.5° (Aoba, Vanuatu). The Height over Radius ratio for the entire edifice H/R, on the other hand, shows higher values, both on average (9.9°) and in a number of high values of up to 17.5° (Filicudi, Aeolian Islands). The most shallow values are in a similar range (4.4°, San Christobal). H/R_asl and H/R_bsl (data not shown) have a similar relationship as the corresponding slope parameters. H/R_asl is ranging from 2.0° (Genovesa, Galapagos) to 28.3° (Alicudi). However, there is one island (Wolt, Galapagos) with an even higher value of 57.7°. We consider this to be an outlier with questionable accuracy. On average the ratio on land (11.6°) is 1.2° higher than the values in the submarine part (10.4°). It ranges from 5.5° (O’ahu, Hawaii) to 17.6° (Salina, Aeolian Islands).

In Figure 3F the dimensionless parameters are summarized. ΔM ranges from 0.07 (Southern Thule; South Sandwich Islands) to 0.84 (Genovesa, Galapagos) with an average of 0.25. The ellipticity index ei describes the elongation of an edifice. The majority of islands show values between 1.003 (Lipari, Aeolian), which is close to a perfectly circular edifice boundary, and 2.72 (Maui group, Hawaii). There are, however, outliers on both sides with Genovesa showing an ellipticity of 3.69 on the one hand and values below 1 for three islands, including Santa Cruz at 0.73. Values under 1 are errors since no edifice can be more circular than a circle. There is no visible trend in Figure 3F regarding ei, but unusually high values usually correlate with high values of R_max in Figure 3D . Finally, ii is ranging from 1.05 (Alicudi, Aeolian Islands) to 2.56 (Hawaii Island) with an average of 1.63. In Figure 3F , a trend of larger islands having higher values of ii is visible. In addition to the aforementioned parameters, we calculate the Jaccard Index, which is a measure of similarity between the polygons enclosed by the edifice boundaries produced with GEBCO and the high-resolution bathymetry, respectively. It is ranging from 0.79 to 0.97 with a median of 0.91.

Figure 4 shows the deviation of the values calculated using GEBCO from those using high-resolution data (D) as well as the influence of the delineation on the deviation (D_Delineation ) and the influence of the DEM (D_DEM ).

The total volume V_tot and area A_tot as well as volume V_bsl and area A_bsl below sea level show a median total deviation D of less than ±1% and quartiles lower than ±5%. The influence on the deviation by the delineation D_Delineation shows similar values, while D_DEM quartiles are below ±1% for these parameters. The total height (H_tot ) and the depth below sea level (H_bsl ) show a slight systematic underestimation of the parameters by GEBCO, mostly visible at the shifted quartiles, rather than the median values (-1.7% and 0.7%, respectively). D_Delineation is symmetric around zero, with quartiles lower than ±1% for both of these parameters, while D_DEM follows the same trend as D, but with smaller quartiles. Where the outliers of the aforementioned parameters lie well within the ±30%, visualized in Figure 4A , the three parameters calculated above sea level (H_asl , A_asl , V_asl ) show significantly higher values of up to ±120% deviation. These outliers are not shown in Figure 4 for reasons of better visualization. While the median deviation of A_asl is close to zero (-0.4%), the height (H_asl ) is slightly underestimated (-4.1%) and the volume (V_asl ) slightly overestimated (+2.3%) by the values calculated with GEBCO. Here, D_DEM is identical to D, while D_Delineation is zero for all three parameters.

Figure 4A also contains the results of the parameters related to the radius of the island. All of these parameters show median deviations of less than ±3% and percentiles of less than ±9%. R_min is slightly overestimated and R_max underestimated when using GEBCO for our calculations, mostly visible in the shifted quartiles, rather than the median values (1%, -0.5% respectively). Consequently, the range of radius R_range has a slightly negative median deviation as well (-2.4%). For all three parameters, D_Delineation is identical to D and D_DEM is zero. The average radius R_A shows quartiles of ±1% or less for all deviations. For R_A,asl D and D_DEM are identical and show quartiles of -5.0% to 1.2% with a median of -0.2% while D_Delineation is zero. For R_A,bsl the median deviation and quartiles are within ±2.3% with a slightly positive shift. D_DEM shows a small upward and D_Delineation downward shift, with minimal overall deviations of less than ±1.6%. Three outliers ±30% of R_range are not shown in Figure 4A .

Figure 4B shows the results of slope-related parameters. Here, we find the highest systematic median deviations of our investigation with -10.3%, -17.4%, and 10.5% for the mean slope of the entire edifice (Sl_mn ), above sea level (Sl_asl,mn ) and below sea level (Sl_bsl,mn ), respectively. The quartiles for all three lie between -6.8% and -22%. These values are almost identical for D_DEM . The median and quartiles of D_Delineation , however, remain < ± 1% for all slope-related parameters. While the height-over-radius ratios HR, HR_(asl) and HR_(bsl) show significantly smaller median deviations, they are still all negative with -1.7%, -2.3%, and -1.8%, respectively. For HR and HR_(bsl) D_DEM shows similar trends with smaller overall deviations, while it is identical to D for HR_(asl) . D_Delineation on the other hand, is zero for HR_(asl) and shows quartiles of < ± 2% around a median of almost zero. Additionally, these parameters produce stronger outliers than the slope parameters. Two outliers of HR_(asl) were omitted for the creation of Figure 4B .

In Figure 4B the results of the shape parameters are shown. All three parameters show median deviations of less than 2%, and quartiles of < ± 10%. D_DEM is almost zero for both the ellipticity (ei) and irregularity index (ii), while D_Delineation is almost identical to D. They show a median deviation of -0.6% and -1.9% as well as quartiles of -7.5%, 7.1% and 5.1%, 4.2%, respectively. The parameter ei produced a single outlier of 133% that was omitted in Figure 4B ). ΔM, the distance between centroid and summit shows median deviations of ±0.3%, and quartiles of -8.9% and 6.5% as well as many outliers, several of which were higher than 50% and were thus not shown in Figure 4B . D_DEM and D_Delineation follow similar trends but with lower overall deviations.

In summary, 20 out of 24 parameters show < ± 2.5% median deviation, and quartiles < ± 10%. Only the slope parameters (Sl_mn, Sl_asl,mn, Sl_bsl,mn ) and H_asl have systematic deviations greater than that.

The graph of the size parameters in Figures 3A–D follows the trend of the “law of high numbers” or “law of large events” (Verma et al., 2006) where large-scale events, or in this case large volcanic islands are less common than small volcanic islands. This is especially evident in Figure 3C , depicting the volume, where we observe the two bends in the graphs. This confirms that the chosen sample of volcanic islands is indeed representative. Figure 3E illustrates how the ratio of height over radius shows higher peaks and a higher average value than the equivalent slope parameter derived directly from the DEM. This is likely due to the resolution dependency of slope parameters (Carrera-Hernández, 2021), which is further discussed in the following subsections. This does not necessarily mean that the H/R is more accurate, we find high outliers, such as Wolt (Galapagos) with a subaerial H/R of 57.7°. Even though Wolt is known for extremely steep slopes (Mouginis-Mark et al., 1996), this proxy for the average slope exceeds the values commonly found for the steepest parts of the island. We therefore consider this to be an error. Figure 3F shows the dimensionless shape parameters. The irregularity index ii appears to be dependent on the size of the volcanic island. This is due to the enormous difference in the size of the edifices, being investigated with the same resolution. Larger islands such as Hawaii are represented in detail, compared to their size. Consequently, their outline will be more complex than a small island delineated with the same resolution. We also see that islands with high values of ei usually also show higher values of R_max and ΔM, which confirms the ability of these parameters to give information about the asymmetry of a volcanic island.

In the following section, we will repeatedly refer to parameters being underestimated or overestimated by using GEBCO. This is relative to the higher resolution ship-based bathymetry data, which in itself can be subject to artefacts or other misrepresentations of the bathymetry. It is, however, currently the most accurate data available and for the sake of this evaluation, we assume it to be the best approximation of reality. This subsection refers to the investigation of all three data groups (A, B, C), in subsection 4.3 and 4.4, we investigate Group B and Group C in more detail, respectively. Even though most of the parameters show only small deviations, there are general trends.

In our investigation, we included two archipelagos with partial high-resolution coverage: the Marianas and the Galapagos Islands. We tested the effect of using a nearest neighbour interpolation of the data gaps in the Mariana Islands compared to using the unmasked (i.e. underlain with GEBCO) version of the dataset provided by GMRT. We found random differences of about ±1%, so this does not seem to have a significant effect on the delineation or the resulting parameters. A major challenge with composite datasets, however, is the artefacts on the borders between high-resolution and background information. These cannot be avoided and are often directly visible in the GEBCO grid. They can create problems, both during the delineation, where the edifice boundary falsely follows the artefact, and during parameter extraction, where they can create distortions in the elevation and slope and thus create a bias in some parameters. Figure 5 shows a plot to compare the distribution of the parameters of the Marianas and Galapagos Islands to the other archipelagos. The Marianas do not show a significantly different distribution compared to the rest of the data, except perhaps for the slope parameters. Here, the islands of this archipelago tend to cluster slightly above the median. Since they are only partially covered with high-resolution bathymetry, the interpolated areas are more similar to GEBCO, resulting in overall lower deviations. The Galapagos islands, on the other hand, do show a slightly different picture, especially concerning the most extreme negative outliers. The Galapagos Islands were the most difficult to delineate, not only because of their partial coverage but also because they are situated in shallow water depths and edifices tend to overlap underwater. They contained one of the smaller islands of the study, Wolt, which is responsible for the negative outliers of about -100% in Figure 5 . This raises the question of whether there is a size threshold for volcanic islands to be investigated with GEBCO. We are unable to define a clear threshold because too many factors besides the size of the island (proximity to other islands, undulating terrain, artefacts, the morphology of the edifice) influence the success of the automatic part of the delineation.

In addition to the two archipelagos with partial high-resolution coverage, we included three examples from Vanuatu. The high-resolution data from Daniel (1993) for these three islands is not freely available and not included in GEBCO. Figure 6 shows a strip plot in which these three islands are highlighted in pink. Though a sample of three is admittedly small, there appears to be no significant systematic error except for the slope parameters, which show a clustering below the median. This is likely the effect of the indirect satellite measurements, which form the base of GEBCO here, already being smoother compared to those areas where high-resolution bathymetry is incorporated in GEBCO. This data is then additionally smoothed during the processing of GEBCO, increasing the difference in slope between the high-resolution data and the satellite-based GEBCO. We, therefore, conclude that even though our dataset mostly consists of volcanic islands covered by bathymetry created from direct measurements, our findings will be applicable to areas where this is not the case.