3.1 GENERAL CONSIDERATIONS

The main problem underlying the identification of “natural” clusters, is that in a general sense there is not a unique definition for what should be considered as “natural”. In the context of the application of methods of cluster analysis to the identification of natural groups of volcanic vents (or of magmatic centers), the non-uniqueness of a “natural” division relies on the various parameters that may be of geological significance. For example, from the point of view of a petrologist, it may seem natural to define the groups of volcanoes in terms of the mineral composition of the erupted products. This criterion alone might provide important information concerning the source of the magma feeding the regional volcanic activity. Nevertheless, this criterion, alone, does not provide information that can serve us to assess the rates of eruption characteristic of that region. To this end, we need to rely on a different source of information, as for example isotopic dates or stratigraphic position of the erupted products (^{Martí et al., 2018}; ^{Németh and Palmer, 2019}). Of course, if attention is focused exclusively on the age-related information, we would be able to define groups based on age of activity. Certainly, each of the petrologically defined groups and each of the age-defined groups may seem very natural. Nevertheless, there is no guarantee that we will be able to establish a neat correlation between the petrologically-defined and age-defined groups. Consequently, the “natural” groupings become less natural when both sources of information are considered simultaneously. If it is further considered that petrologic and age-related information are not the only pieces of information that are important on the study of volcanism, the departure from a “natural” group becomes more marked. For example, a third source of information that might be useful to consider is the size of the volcanic edifices. As it is the case with each of the petrologic and age-related criteria, the size-related classification should be enough to identify a set of “natural” groups if considered in isolation, but the number of those groups, and the volcanoes that belong to each of them, may not have a simple correlation with the petrologic and age-related partitions. Evidently, if other criteria, such as morphology of the edifices, volume of erupted products, proportion of pyroclastic/lavas erupted, etc. are to be considered at some time during our analysis, the case for a “natural” definition of a cluster of volcanic vents becomes ever weaker.

Even if attention is focused on only spatial location of eruptive centers, it might be difficult to identify a “natural” partition. For example, let us think of a polygenetic volcano that has experienced eruptive episodes from three radial rift zones (Figure 1). The question here is what should be considered to be more natural: a) one cluster that includes all the vents and suggests a single, almost concentric distribution, b) two clusters located on top of the two rift zones that have the larger number of eruptive vents, c) three clusters located each on top of a rift zone, d) four clusters, three each on top of a rift zone and a fourth group surrounding the vents near the central conduit, e) a larger number of clusters that delineate the location of the three rift zones and the central conduit, as well as also serving to identify departures from a simple radial pattern in one of the rift zones, f) a structure that reveals the same components as e) plus a few clusters in the outskirts of the main edifice, etc. As it turns out, all of those spatial structures can be documented in one specific case (Mauna Kea volcano, Hawai’i). Thus, each of the mentioned alternatives is physically possible, and has a relevant geologic interpretation. Therefore, it is unclear which of the alternative grouping schemes should be considered more “natural” than the others.

Figure 1 Diagrams showing different clustering options that can be achieved when analyzing the vent distribution of cinder cones on top of Mauna Kea Volcano (^{Cañón-Tapia, 2022}). Each grouping scheme is found using a Gaussian Kernel with a smoothing factor as indicated by the value of h on top of each diagram. A similar range of h values is obtained when using different automatization algorithms. The location of eruptive centers on top and around Mauna Kea were obtained from GoogleEarth, and saved as a MATLAB file. Background image produced by using a section of the SRTM 1arc_v3 geotif downloaded from the USGS EarthExplorer server. A total of 242 vents were used on the analyses. 

3.2 SOME TECHNICAL DIFFICULTIES

In addition to the problems found when attempting to define what should be the natural spatial structure of vent distribution at a given place, there are other issues related with the methods of analysis that need to be examined. Historically, most works documenting the spatial distribution of volcanic vents published prior to 2020 adopted an approach that privileged the statistical interpretation of results, assuming that the method selected for the spatial analysis had no influence on the outcome (^{Bleacher et al., 2009}; ^{Kiyosugi et al., 2012}; ^{Kiyosugi et al., 2009}; ^{Magill et al., 2005}; ^{Martin et al., 2003}; ^{Weller et al., 2006}). Although indeed, not all of those works explicitly aimed to identify clusters of vents, the groupings used on each of those works were justified by the use of a unique method of analysis, often based on an “optimal” value of the smoothing factor in a kernel function. This is not the place to examine whether the conclusions reached on each of those studies is truly geologically significant or not. Nevertheless, the example of Mauna Kea volcano described at the end of the previous section should suffice to make us realize several issues: 1) even when it might be tempting to focus in only one solution, such procedure not necessarily captures all the relevant characteristic of the volcanic activity on the region of study, 2) the interpreted outcome may not be the only situation of geological significance, and 3) other interesting relationships not revealed by the “optimal” statistical approach might be waiting for discovery.

Consequently, it is of interest to undertake a closer examination of the technical aspects concerning the form in which the “optimal” spatial structure that minimizes “error” has been commonly selected. To facilitate their presentation, the technical aspects concerning the selection of an “optimal” spatial structure will be separated below in two broad groups: those concerning the selection of a method of analysis to isolate the unique spatial structure that will be further analyzed and those concerning the selection of an optimal parameter to be used with a single method. Among the second group, I will focus on the form in which the selection of an optimal smoothing factor has been decided in a kernel method.

4.1 CONVERGENCE OF DIFFERENT METHODOLOGIES

Although it might be considered that there are fundamental differences that distinguish the many methods available for the study of spatial distributions (alignment detection, cluster identification, spatial statistics, etc.), all of these methods share the objective of providing a summary of the distribution that should facilitate its association with a probable causal agent. Thus, despite their inherent technical differences, it might be possible to use two or more of those methods to detect the underlying pattern that might point to the presumed causal agent. Evidently, combination of methods that belong to a similar category (for example different algorithms used to define clusters) is easier than the combination of methods that belong to different categories (for example one method of cluster analysis with one designed to detect alignments, or one method that yields results in the form of numerical statistics with another method that provides results of a more graphical nature). Nevertheless, if attention is focused on the common aspects of both methods, a fruitful combination can be devised.

An example of the form in which different methods belonging to the same category can be used to constrain the possible number of outcomes deemed of relevance, is the definition of the so-called core-groups. A core-group was defined by ^{Cañón-Tapia (2020)} as the common output of more than two clustering methods applied to obtain partitions from a collection of vents whose geographic coordinates are known. As shown on that paper, the common outcome of several clustering algorithms (agglomerative with different linkage definitions, k-means, k-medoids, gaussian-model) yields more significant results in a geologic context than either of those algorithms used alone. Consequently, the location of the groups that might have more relevance in a geologic context can be constrained more or less efficiently by adopting such approach. Unfortunately, a condition required for the use of coregroups is that the number of groups needs to be specified at the beginning of the analysis, and in most situations of volcanic interest, the number of relevant groups is hard to be known with certainty at the beginning. Furthermore, as exemplified by the case of Mauna Kea volcano illustrated above, such number might be different depending on the scale at which the analysis is made.

4.1.1 CONSTRAINING THE NUMBER OF CLUSTERS

To constrain the number of groups that should be considered relevant, it is possible to resort to one method whose results are relatively easy to visualize: the gaussian kernel estimator. As a general rule, a large smoothing factor in the gaussian kernel should yield a unimodal distribution (marked by a single set of concentric contour lines, as in Figure 1f), and a small smoothing factor should yield a large number of groupings (marked by many different sets of concentric contour lines, and some contours that are not concentric in relation to others, as in Figure 1a). The rate of increase of the number of modes in the resulting distribution is not lineal, nor monotonic, but commonly four or five attempts should suffice to constrain the values of the smoothing factor that are likely to bracket the number of clusters that seems convenient to consider with more detail. An additional advantage of using this method to constrain the number of clusters of interest is that it also provides some indication of the location of those clusters, location that can also be compared at a later stage with that obtained with the use of the core-groups, as shown in Figure 2. Agreement between the kernel and core-groups methods therefore increases the probability of identification of the groups with more geological relevance.

Figure 2 Comparison of the groupings identified using the Core-Groups method (CG) and the zones of higher vent density identified with the Gaussian kernel (PDF), when analyzing the vent distribution in the island of Jeju, North Corea (^{Cañón-Tapia, 2021b}). Some of the core groups coincide with the kernel clusters, but some do not. This result suggests the possibility of a sub-clustered structure. Vent location obtained by examination of Google Earth images, combined with the location of vents indicated in the maps by ^{Song et al. (2018)}, ^{Sohn and Park (2005)}, and the Digital Elevation Model STRM 1 arc-second downloaded from https://earthexplorer.usgs.gov. A total of 304 edifices were included in the analysis. All data were processes using MATLAB codes. 

A clustering algorithm noted by its versatility, known as HDBSCAN (Hierarchical Density Based Clustering Algorithm with Noise), is similar in many respects to an exploratory analysis made with the gaussian kernel using different values of the smoothing factor. In both cases zones of high density of data are identified. The difference in the methods consists in the form in which the evolution of those zones of high density are handled and displayed. To obtain results with the HDBSCAN it is necessary to specify a minimum number of members to be considered as a relevant group, whereas such specification is not required on the kernel approach. In the HDBSCAN algorithm, the minimum membership is used to trim the hierarchical tree of clusters, and only those branches that have at least the required number of members are retained. The visualization of this process is thus made in a two-dimensional diagram in which clusters and subclusters are shown as branches on a dendrogram, but no clear reference is made in regards the spatial location of those branches. In contrast, in the kernel approach the branches are not trimmed, and the spatial location of each group is always visible. In any case, despite those differences the results of these two methods of spatial analysis (HDBSCAN and Gaussian kernel) tend to be similar to each other.

Yet a different method that can be useful to determine the number and location of groups of vents is to adopt a map representation of the modified two-vent azimuth method as described in more detail by ^{Cebriá et al. (2011)}. Although in this case the separation between vents that is deemed significant needs to be specified beforehand, the results of the groupings outlined by the zones with more lines joining adjacent vents provides an independent estimation of the location and number of possible clusters that is not based on density structure alone. Thus, both methods (density based and line-identification with length limitation) should yield somewhat similar results at small spatial scales, as illustrated in Figure 3.

Figure 3 Comparison of the results obtained when applying a two-point azimuth method with distance limitation and the Gaussian kernel during the study of vent distribution in the Pinacate area, Northwest Mexico (^{Cañón-Tapia and Jacobo-Bojórquez, 2023}). Note the zones of higher vent density marked by the larger concentration of line segments, and also the possible connection between clusters (indicating a probable sub-clustered structure) marked by a few line segments. 

4.2 GEOLOGIC CONSTRAINS

In the rare situation when precise information concerning composition and age of each vent is available, criteria for the identification of spatial patterns might be easy to devise. Alternatively, such abundant information might serve to test the relevance of the results of the various methods of spatial analysis used. Unfortunately, in most situations it is common to lack detailed information about the vents, and most commonly we have to deal with scenarios in which the only information available comes from a satellite image, or a list of coordinates indicating the position of the vents. Fortunately, even in those situations it is possible to include some constraints that might be useful for the interpretation of results. The most important constraint is the distance separating two vents that are likely to be connected at depth by the same conduit. Although the extension of some eruptive fissures can be very large, the separation between adjacent vents along that fissure is commonly less than 5 km (^{Churikova et al., 2015}; ^{Thordarson and Self, 1993}; ^{Tibaldi, 1995}). Thus, limiting the search of neighboring vents to this distance might be helpful to constrain the location of relevant clusters. In addition, knowing that the orientation of feeder dykes is always perpendicular to the orientation of the minimum horizontal stress (except from very local deviations), but that such local orientation might change systematically with the distance from a central conduit or magma reservoir at depth, might help to separate geologically significant interpretations from other possibilities less unlikely. Another aspect to be considered is that in general larger zones of magma storage are likely to be found deeper than shallow reservoirs, and that larger reservoirs are likely to feed more than one eruption than smaller ones. For that reason, a cluster-subcluster structure might become significant in some occasions. Finally, acknowledging that the shape of reservoirs of magma at depth might evolve with time, and sometimes younger eruptions may take place overlapping zones of older activity might prove to be an important consideration that helps to discriminate the relevant from the not-so relevant vent clusters. A summary of conceptual models is shown in Figure 4.

Figure 4 Summary of conceptual models illustrating the relation of volcano distribution at the surface and subvolcanic components (magma reservoirs, dikes) beneath the surface. a) system in which magma has the same probability to be extracted from every location within the zone of storage, travelling to the surface in conduits that might have a range of inclinations relative to the vertical. The distribution of volcanoes at the surface should be relatively uniform, and the lateral extension of the zone covered by volcanism reflects the lateral extension of the zone of magma storage (ZMS) with some accuracy. b) Similar to a, but with the depth to the main ZMS at a shallower depth, which results in a decreased lateral extension of the system at the surface. c) System with multiple, intermediate depth ZMS. Even if it is assumed that magma has the same probability to be extracted from every location within each of the multiple ZMS, the distribution at the surface is likely to experience some degree of clustering related to the presence of those zones. The lateral extension of the zone of volcanism in this case not longer reflects the dimensions and location of a unique ZMS, but still provides a good approximation to the integrated system defined by those zones. d) Similar to c), but showing the effect of a concentration of activity towards the center of the deepest ZMS , and the construction of a polygenetic volcano. E) Effects of a time gap in a system similar to the one displayed in c. 

4.3 FUZZY CLUSTERING, OVERLAPPING SYSTEMS, OUTLIERS AND SUB-CLUSTER STRUCTURES

The complexity of volcanic phenomena may result in “noisy” data that are difficult to separate in different groups with sharp boundaries. Several methods of clustering analysis have been devised to deal with such situations, including the so-called fuzzy clustering techniques, as well as finite mixture models (^{Everitt et al., 2011}; ^{McLachlan and Peel, 2000}). As is the case with the non-fuzzy (aka crisp) partitioning methods, it is necessary to specify the number of expected clusters beforehand, and selection of a different number of clusters result in different cluster partitions. Nevertheless, unlike the crisp-partition counterparts, in fuzzy clustering the partition results are accompanied by information that denotes the grade of membership of each element to each of the identified clusters. In this form, one observation can be interpreted as belonging to two or more groups simultaneously, although with a different level of probability in each case. Ambiguity in the association of each observation with one specific group may be convenient to describe situations in which two or more volcanic systems overlapped spatially, being the result of eruptions fed from different magma reservoirs, or as the result of reactivation of activity in the same general region after a large period of repose after which stress orientations may have changed notably.

One advantage of the fuzzy methods is that they allow us to define a numerical threshold level dividing the clustered observations from those outside a given cluster. Nevertheless, deciding which threshold level should be applied in each particular instance is best left as a subjective decision that ideally should incorporate other sources of information (age, composition, etc.) even if such information is only scarcely distributed among all the observations under examination. Thus, deciding whether a threshold value should be applied to all the clusters in a given region, or it might take different values as a function of geographic position, are examples of decisions that need to be made in a case-by-case basis, and always under the guidance of the local geological features observed in the region of study.

One aspect of the fuzzy methods that may be very useful is the possibility to define somewhat qualitatively observations that can be outliers (isolated vents that do not belong to any particular cluster in the region). Interpretation of such types of vents is always complicated because it is uncertain whether they are extreme events within a dispersed cluster, they represent a first event that could lead to the formation of a different cluster in the region (therefore marking a change in the conditions leading to volcanic activity), or they should be interpreted as isolated eruptions that occurred by mere chance. From an analytical point of view, the outliers influence the location of other clusters in the region, and could lead to biases in the sense that they can be considered by some algorithms to be one group, therefore changing the boundaries and locations of the other groups in the region. Due to their isolated character, this type of observations have very small probabilities of being associated with the other clusters in the region, and therefore can be filtered out from the analysis based on such a probability. Nevertheless, their geologic interpretation or relevance is not always straightforward.

One aspect of the results of fuzzy clustering methods that is not user-friendly concerns the identification of a cluster-subcluster structure. Although all the observations of a subcluster might be characterized by a smaller probability of membership association than the vents in a different subcluster, those differences might be subtle and difficult to detect unless a more complex series of tests is completed. A simpler approach that might be adopted to assess the possibility of such an spatial structure was proposed by ^{Cañón-Tapia (2021b)}. The procedure is applied to groups of vents that have been identified by any other methods, but that for some reason might seem to be part of a subclustered structure. The assessment is made based on the calculation of the nearest neighbor distances of the observations. From these, a reference distance for each suspected group (called the within cluster distance) is defined by finding the largest nearest neighbor separation for the vents thought to be part of one group. This reference distance is then compared with another reference distance (called the between cluster distance) defined as the shortest nearest neighbor distance found between one vent in one suspected group and one vent in a different suspected group. Then, for each pair of suspected groups the following rule is applied: If the two within cluster distances are smaller than the corresponding between cluster distance, then the two suspected groups can be considered as independent from each other. On the contrary, if one of the within cluster distances is larger than the corresponding between cluster distance, then both groups should be considered to be subclusters of a larger cluster that includes at least those two. Displaying the between and within cluster distances in a triangular matrix therefore facilitates the identification of possible subclusters, and provides a quantitative basis that could be then used as a starting point to make more stringent tests within finite mixture models, if so desired. Nevertheless, the possibility of having subclusters, and their possible geologic implications can be assessed simply from the examination of the mentioned matrix.