Monogenetic volcanism in subduction settings: comparative statistical study of the Sierra Chichinautzin and Los Tuxtlas Volcanic Fields in Mexico

Monogenetic fields present significant diversity, yet this diversity has not been fully quantified, and its origin remains elusive. We studied two large subduction-related fields in Mexico, the Sierra Chichinautzin and Los Tuxtlas, that have distinct crustal stress regime and structures, magma compositions, vent types, and climatic conditions. Using recently available 5-m resolution topographical data, we located all the eruptive centers, studied their spatial distribution and analyzed scoria cone shapes in detail, calculating morphometric parameters for the best preserved. We then applied a set of statistical tools to analyze and compare the patterns of vent distribution, vent alignment, and diversity in cone shapes in these two fields. We observe that, despite their distinct setting, the two fields are similar in terms of vent distribution and cone morphology, which shows that this type of data cannot be used alone to infer the tectonic, magmatic, and climatic context of monogenetic fields. It also confirms previous results that the diversity in cone shapes (slope, height-to-diameter ratio) reflects processes that are common to all cones (e.g., ballistic emplacement followed by scoria avalanching on slopes), and hence do not vary significantly (at field-scale) with external parameters. Differences in the crustal stress regime had no apparent impact on vent distribution as the dikes followed active faults, irrespective of their motion. Climatic differences did not affect the shape variety of the studied cones probably because of their young ages (< 50,000 years old) and their location in a vegetated environment. The fields nevertheless differ in size and vent density, as well as scoria cone shape complexity and volume, which can be attributed to differences in the geometry of the magma source for its impact on the closeness of the dikes feeding the activity. Differences in the relative proportion of small cones in both fields are likely due to factors impacting eruptive style such as magma-water interaction, magma composition, and/or fissure lengths.


Introduction
. In turn, those processes vary with the tectonic setting that determines the conditions of magma generation, ascent and eruption, and, to a lesser degree, with the climatic conditions that control the availability of surface water and rates of degradation of the edifices (e.g., Wood 1980b;Fornaciai et al. 2012;. All these factors have multiple effects and complex interactions, some of which are highlighted here. The geometry of the magma source has probably a major role on the overall size and shape of the fields as well as the occurrence of vent clusters, as it determines where magmas are generated at depth (Connor 1990;Valentine and Perry 2006;Le Corvec et al. 2013). Nevertheless, the crustal stress regime, thickness and fracture network, as well as magma composition, control if and where dikes are likely to reach the surface, acting on vent alignment and spacing (e.g., Haug and Strecker 1995;Conway et al. 1997;Parsons et al. 2006;Gaffney et al. 2007;Cassidy and Locke 2010;Valentine and Hirano 2010;Mazzarini et al. 2010;Maccaferri et al. 2014). The size of scoria cones, the most common monogenetic edifice, is determined by magma batch size, dike geometry or fissure length, and the explosivity index that control, respectively, how much magma is erupted, how many edifices form during a single event, and what is the relative proportion of pyroclastic vs effusive products (e.g., Wood 1980a; Bemis et al. 2011;Dóniz-Páez et al. 2012). The magnitude of the explosive activity, influenced by magma composition, volatile content, and magma-water interaction, affects the trajectory, size, and degree of welding of the pyroclasts, which primarily control scoria cone shapes, in addition to substrate topography, cone spacing, and faultand lava-induced breaching (e.g., Wood 1980a;Head and Wilson 1989;Tibaldi 1995;Bemis and Ferencz 2017). The climate, in particular rainfall and temperature, also influences scoria cone morphology due to the redistribution of tephra by erosional processes (e.g., Wood 1980b;Hooper and Sheridan 1998;Dóniz et al. 2011). The high number of factors involved suppose a significant diversity in monogenetic fields.
Comprehensive studies of monogenetic fields in distinct settings are necessary to capture the entire diversity of monogenetic activity and understand the related mechanisms of magma generation, crustal storage and eruption, all of which have great relevance for hazard assessment and mitigation. A deeper knowledge of the elements that control the distribution of vents has allowed to improve forecasts on the location of future eruptions (e.g., Bebbington and Cronin 2011;Kereszturi et al. 2017;Connor et al. 2018;Sieron et al. 2021). In particular, qualitative and quantitative studies on vent alignment and clustering have informed about the role of faults on vent location (e.g., Conway et al. 1997;Germa et al. 2013;Le Corvec et al. 2013). The morphological characterization of the vents has provided valuable information about the nature of the most probable future events (e.g., Siebe et al. 2004;Dóniz et al. 2008;Becerra-Ramírez et al. 2022), which is being integrated into hazard maps (Marrero et al. 2019).
The central and southern Mexico region is a particularly interesting case of study as monogenetic fields are large and numerous in this area and display significant diversity with respect to their size, chemical composition, prevalent volcano type, and association or not with a central volcano (e.g., Guilbaud et al. 2009;Chédeville et al. 2020;Sieron et al. 2021). To assess more quantitatively the variety that exists in this region and identify the controlling factors, we chose to focus on two relatively well-studied fields that contain preserved scoria cones of similar age (Late Pleistocene to Holocene, mainly) but are located in widely separated sectors of the volcanic arc: The Sierra Chichinautzin volcanic Field in the central-eastern sector of the Trans-Mexican Volcanic Belt and the Los Tuxtlas Volcanic Field situated to the SE of the belt (Fig. 1). Although there have been debates, most data support that both fields are related to the subduction of the Cocos plate, albeit with differences in the geometry of the slab and in the crustal stress field and associated structures (e.g., Pardo and Suárez 1995;Pérez-Campos et al. 2008;Ferrari et al. 2012;Calò 2021).
This study aims to test whether these two fields which have different settings have different spatial vent distribution and cone morphologies or not. Accordingly, using recently available high-resolution topographical data, we located all the eruptive centers in these two fields, studied their spatial distribution and analyzed scoria cone shapes in detail, calculating morphometric parameters for the best preserved. We then applied statistical tools to analyze the patterns of vent distribution (density and alignment) and cone morphometry and to test whether the fields are statistically different in these aspects. The morphometric data has been partly presented elsewhere (Vivó Vázquez 2017;Sieron et al. 2021) but was entirely revised here using similar methods for both fields, in order to make meaningful comparisons. The statistical tools we use, specifically principal components analysis (PCA), and two techniques for comparing the covariance matrixes of morphometric parameters from both volcanic fields testing for similarity, have been rarely applied in volcanology and we briefly discuss here their usefulness. PCA is often used in morphological studies in biology for example (Arnold et al. 2008), and we hope that this study will contribute to their more ample application in this research field. In the following, we first summarize the differences in setting of the two fields. We then use the previously mentioned statistical techniques to assess in a semi-quantitative manner the differences and similarities of two fields in terms of vent distribution and scoria cone morphometry, and discuss how the elements of their setting may be involved in causing such variations. Based on this, we can then conclude on the factors that have significant impacts on the morphological variability of monogenetic fields.

Volcanic, tectonic, and climatic context
In this section, we introduce the key background characteristics of the two fields that are likely to affect their vent distribution and the morphology of the scoria cones that occur in them. The contrast between the two fields is highlighted in a summary table (Table 1).

Sierra Chichinautzin Volcanic Field (SCVF)
The SCVF forms an E-W elongated field constructed on an elevated plateau (the Sierra Chichinautzin) that stretches between two stratovolcanoes, the Nevado de Toluca to the west and the Popocatepetl to the east ( Fig. 2A). It is crossed roughly in the middle by an approximately N-S range of deeply eroded domes (Sierra de Las Cruces, Fig. 2a)). Márquez et al. (1999) identify a minimum of 221 individual vent structures in the SCVF. Volcanoes consist mostly of scoria cones with lava flows, with fewer medium-sized shields crowned by a cone and thick lava flows and domes (e.g., Bloomfield 1975;Martin del Pozzo 1982;Siebe et al. 2004Siebe et al. , 2005Agustín-Flores et al. 2011;Lorenzo-Merino et al. 2018). Noteworthy, there is no phreatomagmatic construct within the field, except for one tuff cone located in the Mexico basin to the north. The ages of the volcanoes range from Holocene, with Xitle volcano being the youngest dated volcano at ca. 1670 y BP (Siebe 2000) to 1.2 Ma for a deeply eroded cone at the foot of Nevado de Toluca volcano (Arce et al. 2013). In general, older edifices are located on the borders of the volcanic field (Arce et al. 2013) and the central part of the field is mostly Late Pleistocene to Holocene in age (Jaimes-Viera et al. 2018). Magma compositions range from basaltic to dacitic, with andesitic being the most common and abundant (Wallace and Carmichael 1999;Schaaf et al. 2005). Structurally, the SCVF can be divided into two parts, separated by the Sierra de las Cruces (Fig. 2a)). The western part of the field is transected by NNW-SSE, NE-SW, and E-W faults attributed to the Taxco-Queretaro, San Antonio, and Tenango fault systems, respectively (Fig. 2a), Garcıa-Palomo et al. 2000Norini et al. 2006Norini et al. , 2008Bellotti et al. 2006). The Tenango system is active since the Late Pleistocene and consists in transtensive left-lateral faults, with an average extension rate of 0.3-0.5 mm/year in the Holocene . By comparison, faults are not expressed at the surface to the east of Sierra de las Cruces. Morphological and geophysical data suggest that this area is principally affected by subparallel E-W normal faults (e.g., La Pera, Xochimilco) defining a horst-and-graben structure (De Cserna et al. 1988;Siebe et al. 2004;García-Palomo et al. 2008;Campos-Enríquez et al. 2015;Arce et al. 2019, Fig. 2a)), but the location and direction of dip of these structures are poorly defined.
Previous studies of the distribution of volcanoes in the SCVF include an analysis of the alignment of volcanic centers that shows that these display a broad range of orientations whose main directions coincide with those of tectonic lineaments in the area (Márquez et al. 1999). Mazzarini et al. (2010) present vent intensity maps showing a main E-W trend and a weak N-S trend for the eastern portion. They estimate an average vent density of 0.09 vent/km 2 and relate the self-similar clustering of scoria cones within the field to the crustal thickness in this area. Doing so, they assume that the magmas ascend directly from a partially melted area located just below the crust and had limited crustal storage, which is questionable based on the petrological data (Schaaf et al. 2005).

Los Tuxtlas Volcanic Field (LTVF)
The LTVF is located in the Tuxtlas area, near the coast of the Gulf of Mexico, about 200 km south of the easternmost limit of the TMVB (Fig. 1). A general description of the field can be found in previous works (Nelson and González-Caver 1992;Sieron et al. 2014Sieron et al. , 2021. The LTVF covers an area of ca. 2800 km 2 and contains at least 368 clearly identifiable monogenetic vents, mainly scoria cones, but also 67 phreatomagmatic constructs (maars and tuff rings/ cones) that mostly have Pleistocene to Holocene age based on radiometric dating and morphology (Reinhard 1991;Nelson and González-Caver 1992;Martínez and Milan 1992;Sieron et al. 2014). Of specific interest for this study, this field covers the flank of San Martín Tuxtla volcano, which is a broad shield volcano that last erupted in 1793 (Fig. 2b), Moziño 1870;Espíndola et al. 2010;Carrasco-Núñez et al. 2010;Sieron et al. 2021), and can be hence considered contemporaneous to the monogenetic activity. The climate is tropical, with high precipitations (1272-4201 mm, up to 7000 mm at higher elevations: Gutiérrez-García and Ricker 2011) that promote a wide range of mass wasting processes and may be related with the importance of phreatomagmatic activity (Nelson and González-Caver 1992;Santley 2007;Sieron et al. 2014Sieron et al. , 2021. The volcanic products are mostly primitive, alkaline and sub-alkaline basalts (Nelson and González-Caver 1992;Nelson et al. 1995;Ferrari et al. 2005).
The structural framework of the LTVF is relatively well constrained from Andreani et al. (2008) and Ferrari et al. (2005), both of which include results from previous studies (e.g., Jennette et al. 2003). Despite the dense vegetation cover, Andreani et al. (2008) identifies numerous faults across the study area ( Fig. 2b)), using PEMEX seismic profiles, digital elevation models, satellite imagery, and microstructural analysis at some sites. Major structures are two parallel regional NW-SE left-lateral, subvertical strike-slip faults called the Sontecomapan and Catemaco faults (Fig. 2b)), that may join at depth to form a large transpressive flower structure called the Veracruz fault (Andreani et al. 2008). Secondary E-W trending faults have been interpreted as synthetic Riedel shears (Andreani et al. 2008). Such transpressive deformation started during the Miocene with en-echelon folds which then converted into faults (Andreani et al. 2008). Seismic studies indicate that the fault system is still active (e.g., Suárez 2000;Franco et al. 2013). Active and extinct polygenetic edifices, together with the surrounding monogenetic volcanoes, form alignments that follow the Veracruz fault (Nelson et al. 1995;Fig. 2 Vent distribution and fault pattern in the studied areas. a) The Sierra Chichinautzin Volcanic Field and b) Los Tuxtlas Volcanic Field. Red symbols indicate the vents mapped in this study (see details in text). The DEMs on this figure have been created by using INEGI 20 m elevation contour lines. Faults for the SCVF are from Garcıa-Palomo et al. (2000), Norini et al. (2006Norini et al. ( , 2008; Bellotti et al. (2006), for the western part (Toluca region) and Arce et al. (2019) for the central and eastern parts. Faults in the LTVF are from Andreani et al (2008) ◂ Andreani et al. 2008;Sieron et al. 2021), suggesting a faultcontrolled magma emplacement that will be discussed here in more details.

Inventory of vent location and density analysis
An inventory was made of the centers of eruptive vents in the two studied fields. Those were determined using hillshade images constructed from a digital elevation model (DEM), which was created from LIDAR data (5 m pixel size) published by Mexico's National Institute for Statistics and Geography (INEGI). We also used information from the literature, satellite imagery (Google Earth), and fieldwork (see also Sieron et al. 2021, andVivó Vázquez 2017). Locally, public pictures shared on the Google Earth platform as well as the street view feature proved to be useful to verify the nature of the topographic features. For both fields, all identifiable craters were considered as individual vents, noting that several closely spaced vents often corresponded to a single volcano or eruptive event. Nevertheless, for the LTVF, only clearly distinguishable craters (separated from each other) were treated as separate ones; when one edifice had several small, coalesced craters only one of them was counted (examples will be presented in the result section). This results in an underestimation of the number of vents relative to the SCVF, which will be considered in the discussion of the results. Note that the databases used in this study are corrected and completed versions of the ones presented by Sieron et al. (2021) and Vivó Vázquez (2017).
In order to compare the spatial distribution of the volcanic fields, we first normalized the UTM coordinates with the min-max transformation: This normalization transformed the data to a scale in [0,1], preserving the relative position between the volcanoes. Applying min-max normalization to UTM longitude and UTM latitude, the volcano locations fell within the unit square (Fig. 3).
Next, their kernel density was estimated using a 2D Gaussian kernel density estimation (KDE) and probability contours were calculated, following Connor et al. (2019). A KDE plot is a method for visualizing the distribution of observations in a dataset (in this case x-y coordinates of volcanic vents), analogous to a histogram. KDE represents the data using a continuous probability density curve in one or more dimensions. As the volcanic vents are defined (1) data − min{data} max{data} − min{data} by 2 coordinates, in this case the KDE a bivariate density function. The probability of a volcanic center occurring in a certain area is estimated (mathematically, the values correspond to a "volume" underneath an "x-y-blanket," being an integration).
We then calculated the smoothing bandwidth matrix H with the SAMSE pilot estimator, implemented in the functions Hpi and Hpi.diag of the ks R package (Duong 2020). This way, we estimated the spatial density of the volcanoes using the kde function, also implemented on the ks package.
When calculating the bandwidth matrices H for both fields, we obtained the following smoothing matrices: The determinants of these matrices named det(H) are 753 km and 21 km for SCVF and LTVF respectively. Furthermore, we used the "point density" tool of the Arc-Map software suite to determine the number of vents per unit area, which is named vent spatial intensity in the following. We used unit cells of 1 × 1 km (1 km 2 ), 5 × 5 km (25 km 2 ), and 10 × 10 (100 km 2 ) in order to evaluate clustering.

Vent alignments
We used the Two-point azimuth MATLAB GUI (Thomson and Lang 2016) to identify vent alignments in both fields. This GUI allows us to analyze the orientation of the vents by using two different methods: the Lutz algorithm (Lutz 1986) that considers all possible pairs of vents (N(N-1)/2) and the Cebriá algorithm (Cebriá et al. 2011 pairs of points that are located at less than a threshold distance d m to each other. The threshold distance is empirically defined as d m ≤ x − ∕3 (Cebriá et al. 2011), where x is the mean value of all distances between vents and is the standard deviation. For both methods, we report the results obtained with the Monte Carlo (MC) option available in the MATLAB GUI. The MC option allows for reducing the influence of the shape of the polygon containing the vents and identifying significant alignment angles (Thomson and Lang 2016 and references therein). Under the MC approach, an angle range is considered significant if the total normalized counts for the corresponding histogram bin is higher than a threshold given by the mean of all MC runs plus two standard deviations (Lutz 1986).

Morphological analysis and morphometric parameter correlation
The morphology of all scoria cones identified in the vent location inventory was studied in detail on the LIDARconstructed images. In this context, we followed Bishop (2009) to classify cones as simple (single edifices), compound (formed of two or more cones which coalesced or overlapped), and complex (associated with a different edifice type, such as a maar). The morphometric parameters of a majority of those cones (see below) were determined "manually" using the procedures that are detailed in Vivó Vázquez (2017) for the SCVF and Zarazúa-Carbajal and De la Cruz-Reyna (2021) for the LTVF. The method and tools used are summarized in Table 2. All measurements were initially made using tools available in ArcGis software. However, for about 50% of the LTVF cones, the slope and volume values that were derived from measurements in ArcGis (zonal statistics tool) were anomalous, probably due to the existence of errors in the Crater depth was not considered in the formula because of high uncertainties in vertical measurements Slope (S) S = 2 * Hco∕(Wco − Wcr) Corresponds to average external slope Flatness (f) f = Wco∕Wcr Height-to-diameter ratio (Ar) Ar = Hco∕Wco DEM. Sieron et al. (2021) attribute those to the high vegetation cover in a major part of this field, which resulted in the incorrect correction of this factor during the processing of the Lidar dataset (causing a "rough" artificial surface). For this reason, for both fields, these morphometric parameters were calculated using a simple formula of a truncated cone, deriving the volume from cone width and crater width values, and the average external slope from cone height and cone width (Table 2). Additionally, the dimensionless parameters height-to-diameter ratio (Ar) and flatness (f) (e.g., Bemis et al. 2011) were obtained for each cone (see Table 2).
The morphometric dataset does not include all cones, but only those for which all the parameters could be estimated (a requirement for the statistical studies), excluding the edifices that are strongly eroded, extensively quarried, rafted, or deeply buried by younger lavas. A critical element in determining whether the cone was sufficiently preserved to be included in our analysis was the occurrence of a clear crater. It follows that our method mainly considers the youngest edifices, which is taken into account in the discussion. For the SCVF, Lidar data was only available for the part of the field that is east of the Sierra de Las Cruces. In addition, the compound cones of the SCVF (Vivó Vázquez 2017, see also results section) were not considered in the morphometric analysis.
To identify differences between the fields, we first assessed for each field separately the correlation among the different morphometric parameters, by calculating the covariance and correlation matrices. The correlation matrix corresponds to the normalized covariance matrix. The correlation matrix was then used to perform several statistical analyses which are described in the following section.

Statistical analysis
First, we applied a Kolmogorov-Smirnov test to evaluate the equivalence of the two datasets corresponding to the two fields. We tested all the morphometric parameters and all values are calculated using the Matlab's two-sample Kolmogorov-Smirnov test function (kstest2). The null hypothesis describes both datasets to be equivalent at a significance level of 0.05. Hence, h = 0 means no rejection of the null hypothesis and h = 1 would mean the rejection of the k-stest statistic ks2stat. P represents the probability of observing a test statistic as extreme as, or more extreme than, the observed value under the null hypothesis. Values of p > 0.05 indicate that the two sampled populations compared could share the same population distribution function.
In the following, five of the morphometric parameters (excluding f and Ar as they show linear relationship with their dependent variables) ( Table 2) were transformed into a logarithmic scale. We did so, because this scale is common in morphometric studies since the classical works of Jolicoeur and Mosiman (1963) and Jolicoeur (1960), and this way the differences in scale and units do not affect the analyses. Cone width (Wco), crater width (Wcr), and cone height (Hco) are measured directly, whereas slope and volume are a product of formulas (slope and truncated cone) of the other variables/parameters (Table 2). Nevertheless, the relation between slope and Hco and Wco is inversely proportional and the relation between volume and Wcr, Wco, and Hco is squared. Hence, the validity of the analysis is not compromised, which would be the case of a lineal relation, such as the one that exists for f and Ar.
First, we conducted principal components analysis (PCA) separately for each covariance matrix of the volcanic fields. Then, we compared the two fields using common principal components analysis (CPCA) (Flury 1984). PCA is one of the oldest and most widely applied multivariate data analysis tools for identifying and studying patterns in multivariable data sets. It reduces the dimensionality of a multivariate dataset while minimizing the loss of information. (e.g., Abdi and Williams 2010;Jolliffe and Cadima 2016). CPCA is a generalization of PCA to multiple groups of observations. Given a centered n × p data matrix data resulting from observing p variables in n objects, PCA produces a new set of p variables, called principal components, by orthogonally projecting the observations along the eigenvectors of the sample covariance matrix of the data. These new variables have maximal variances and are uncorrelated among them. For presentations on PCA, see Jackson (1991) and Jolliffe (2002). Recent reviews of PCA can be found in Jolliffe and Cadima (2016) and Björklund (2019).
We also use biplots for data visualization. A biplot is a joint display of the units and the variables of a multivariate data matrix (Gabriel 1971). A PCA biplot is a scatterplot of the sample units on the plane generated by the first two principal components (PC1 and PC2), which explain most of the variation in the data (Jolliffe and Cadima 2016) and are hence used for data visualization. The variables are represented on this same plane (that is the reason for the prefix "bi") with arrows. The lengths of these arrows are approximately proportional to the standard deviations of the variables, and the cosines of the angles between the arrows approximate the correlations of the variables. Thus, the more parallel two arrows, the higher the correlation between the variables, while orthogonal arrows indicate uncorrelated variables. See Gower and Hand (1996) for a detailed presentation of biplots.
Let Σ 1 and Σ 2 be the population covariance matrices of the parameters in log scale of LTVF and SCVF, respectively. Flury (1984), see also Flury (1988), developed the common principal components model as where is the orthogonal matrix whose columns are the eigenvectors of Σ 1 and Σ 2 ; Λ 1 and Λ 2 are diagonal matrices with the eigenvalues of Σ 1 and Σ 2 in their main diagonals, respectively. So, the common principal components (CPC) model establishes that the covariances matrices Σ 1 and Σ 2 share the same eigenstructure. The null CPC hypothesis is H CPC ∶ Σ 1 = BΛ 1 B � ,Σ 2 = BΛ 2 B � and is tested against the alternative that the covariance matrices are arbitrarily different H a ∶ Σ 1 ≠ Σ 2 . In simpler words, two systems with the same principal components are interpreted as having equivalent covariance matrices. In the particular case of morphometric parameters of scoria cones, the equivalence of the principal components of two different volcanic fields would imply that the morphometric parameters vary in the same way among each other (our null hypothesis). Flury (1984) developed maximum likelihood estimation of the covariances matrices under multivariate normality. Let, j = 1, 2, be the maximum likelihood estimators of Σ 1 and Σ 2 under H CPC . The estimators of Σ 1 and Σ 2 under H a are the sample covariances matrices, S 1 and S 2 , of the volcanic fields. The log-likelihood ratio test statistics for testing H CPC against H a is (3): where n 1 = 110 and n 2 = 100 are the sample sizes. Under H CPC and multivariate normality, the test statistics 2 follows asymptotically a chi-square distribution with df = (k − 1)p(p − 1)∕2 degrees of freedom, where k is the number of populations and p is the number of variables. In our case, k = 2 and p = 5 , so the degrees of freedom are 10. The asymptotic p-value is, therefore, P( 2 df ≥ x obs ) where 2 df denotes a chi-square random variable with df degrees of freedom and x obs is the observed value of the test statistic 2 (likelihood ratio).
Since we do not want solely to rely on the multivariate normality assumption, we also tested H CPC using a permutation procedure. We used 2 as a test statistic for this procedure since large values of 2 represent departures from H CPC . Permutation tests are powerful counterparts of parametric tests. Asymptotically they are, in general, equivalent to parametric tests (e.g., Berry et al. 2014). Furthermore, permutation tests are appropriate for nonrandom samples or when analyzing entire populations, see Berry et al. (2014). Also, permutations tests are essentially nonparametric in nature, because they do not assume an underlying distribution of the data. A sufficient condition for a permutation test is the interchangeability of the sample units. In our case, this assumption establishes that the two populations of scoria cones, constrained within the LTVF and SCVF, are equal in terms of H CPC . Thus, the name of the volcanic field attached to the volcanoes is just an uninformative label. We can pool the 210 volcanoes of both fields together and randomly resample without replacement 110 volcanoes and label them where I denotes the indicator function. A value of 1 is added to the numerator and the denominator to take account for the observed value x obs . The algorithm of our permutation is available as supplementary material.
We tested H CPC with the function cpc.test of the R package CPC of Pepler (2019). The R code developed is available upon request to the authors. We used the R function princomp for PCA (R Core Team 2013).

Spatial distribution of volcanic centers
KDE-derived probability contours are drawn on a map of both fields in Fig. 4 and spatial vent intensity is shown on Fig. 5. The LTVF covers a smaller area (1/5) compared to SCVF. The LTVF has a greater spatial intensity of vents (Fig. 5). The maximum vent intensity per 25 km 2 is 25 in the LTVF, compared to 10 in the SCVF. This difference is significant, considering that this value is under-estimated in the LTVF because of the lack of consideration of closely spaced craters in complex cones (see method section). Interestingly, when reducing the unit area to 1 × 1 km 2 , both fields present a maximum intensity of 5, but those areas with higher intensity (dense spots) are more abundant in the LTVF. This means that in both fields spatial clustering exists (up to 5 vents per 1 km 2 ), but in the SCVF they are less abundant and single (separated) vents predominate over the whole area.
In both fields, the KDE probability contours display a prominent elongation in the direction of the main topographic high, which is the Sierra Chichinautzin for the SCVF and the Tuxtlas massif for the LTVF. This elongation is more marked for the LTVF that displays ellipticalshaped contours than for the SCVF whose contours are more rectangular shaped. For the SCVF, the 5 to 25% probability contours correspond to an area with a high concentration of young volcanoes, many of these belonging to an alignment composed of several partly imbricated vents or series of craters (from E to W: Guespalapa, Chichinautzin, Cerro del Agua, Pelagatos). For the LTVF, the highest probability contours correspond to the central area comprised between the two main faults (Sontecomapan and Catemaco) of the Veracruz fault system, where young volcanoes occur, but also volumetrically larger cones, from which long lava flows were emitted (Sieron et al. 2021).
In both fields, the probability contours of 5 to 25% form regular ellipses, while the 50 to 95% probability contours show convoluted shapes. Those include clusters of vents that are located outside of the main topographic high zone. For the SCVF, these clusters are made of aligned vents with similar morphology and inferred ages (e.g., Sierra Santa Catarina, Joquicingo, Tenancingo, Tetillas). Those range from 1 Ma old (Tenancingo, Arce et al. 2013) to recent (Santa Catarina, Jaimes-Viera et al. 2018). Within the LTVF, vent probability diminishes with distance from the NW-SE central zone comprised within the two major faults. The highest concentration of phreatomagmatic vents (mainly maars) lies in the SE part (Sieron et al. 2021).
When comparing the bandwidth matrices of both fields, it becomes clear that, the determinant of the SCVF matrix (753 km) is much greater than the one of LTVF (21 km). Hence, the LTVF is much more densely clustered than the SCVF. This also means that in the SCVF volcanic vents are grouped in different clusters over the whole area, while in the LTVF volcanoes seem to belong to a single cluster. ranging from 95 to 5%. Spatial density can be interpreted as a probability function (Connor et al. 2019), meaning that within the 95% contour there is a 95% chance of a volcano to form in the future This clustering in the SCVF most probably originates from considering the Santa Catarina vents that are located in the basin, at some distance from the main mountain range (Fig. 2).
Results from the Two-point azimuth analyses are displayed in Table 3. For SCVF, the average distance between vents is x ~32 km, the standard deviation is ~16 km, and the threshold distance for the Cebriá et al. (2011) method is d m ~5 km; whereas for LTVF x ~16 km, ~ 9 km, and d m ~ 2 km. In the LTVF, vents show a frequent alignment of NW-SE (N145°), which corresponds to the direction of the main Veracruz fault (Fig. 2). There is a subordinate alignment direction of almost N-S (N85°), detected by the Lutz method. In the SCVF, the predominant alignment is E-W, with a subsidiary contribution of ENE-WSW (N25°). The main alignment trends correspond clearly to the local and regional fault systems in the respective study areas (Fig. 2).

Scoria cone morphology: general descriptions
Scoria cones in the SCVF are morphologically diverse. They consist of simple cones (78%) and compound cones (22%). Simple cones have clearly defined outlines and a unique crater (Fig. 6a)). Compound cones range from single cones that have two or more craters (Fig. 6b)) to coalesced edifices and craters, generally aligned, that formed during the same eruptive vent (Fig. 6c)). The cones typically have gullies on their external slopes, which range from thin, shallow, and closely spaced in young cones (Fig. 6a)) to deep and widely separated in older cones (Fig. 6d)). However, some eroded cones with partly filled craters and smooth outer slopes do not display regular gullies (Fig. 6e)). Some cones are not closed constructs but present breaches or notches, causing a horseshoe shape. Lava outflows are often associated to the open part of the crater, which in some cases is probably related to rafting and in others to contemporaneous effusive and explosive activity as described by Bemis and Ferencz (2017). Older cones are often buried partially to nearly completely by lavas from younger vents, especially when they occur downslope of other vents (Fig. 6f)). Finally, cones are mostly circular and elliptical in plan view.
In comparison with the SCVF, the LTVF contains a much lower percentage of simple cones (about 36%) compared to compound ones (about 64%). There is some uncertainty however in these estimates because of the very high density of edifices, which results in a wide variety of atypical shapes (see below) and in the large amount of overlapping between cones that may be contemporary or very distinct in age. Simple cones occur as round (ring-) shaped edifices (Fig. 6g) or they may be notched, with lava flows emerging from the notch in places (Fig. 6h). Compound edifices range from simple cones with several coalesced craters to cones with complex shapes (elongated in one direction, or angular, Page 13 of 24 14 coalesced edifices) that mostly enclose rows of craters and may include notched cones (Fig. 6i, j).
Additionally, in the LTVF, the scoria cones may be associated with maars (e.g., Fig. 6k and m), which form complex structures (e.g., the Nixtamalapan-Las Animas complex: Fig. 6l). Scoria cones sometimes form exactly at the center of a maar structure (Fig. 6k) and, in this case, may be very small in size (Fig. 6m). Also, a maar may cut a previously formed scoria cone (Fig. 6l). The proximity in age of the maar and the cone was in cases confirmed by field observations (directly overlapping deposits of neighboring volcanoes, with no soil formation in between). As mentioned in Sieron et al. (2021), highly vegetated lower-altitude parts of the field present low-resolution topographic data, which obscure the morphology of the cones (Fig. 6n). In addition, large scoria cones clusters often contain remnants of scoria cones that are of undefined shape and have been excluded (Fig. 6o).

Morphometric analysis
The scoria cones selected for advanced morphometric analysis amount to 180 for the LTVF and 100 for the SCVF, which account for more than 50% of the total number of scoria cones in both fields. The range and mean value of cone´s morphometric parameters are similar for the two fields (Table 4). They are in average ~ 100 m in height, have a 200-m-wide crater, a 500-600-m-wide base, a volume of ca. 0.04 km 3 and an average external slope of 25°. Those values agree with estimates of the average size of cones worldwide (Hargitai and Kereszturi 2015). The average slope value is within the range observed at historic cones (Table 1 in Bemis et al. 2011).
The morphometric dataset is reported in Table 4 and plotted in Fig. 7 as histograms. It can be observed that, although the ranges of the cone dimensions are similar in both  . a), b), c), d), e), f) SCVF cones. The white line deliminates the cone base that was considered for deriving the morphometric parameters; g), h), i), j), k), l), m), n), o) LTVF cones. All images: North is up and the scalebar (white rectangle) is 200 m long. a) Example of a large simple cone with a unique deep crater and smooth external slopes. Xitle volcano, Holocene. Note a smaller segundary crater (Xicontle) that was formed to the west during the same eruptive event and was counted as a vent in the density analysis; b) example of cones with multiple craters. The scars on the southern flanks of the cones are quarries. Cuatepel and Aholo volcanoes, Late Pleistocene age; c) example of a complex cone made of imbricated or coalesced cones and craters formed during a single, fissure-fed event. Tlalóc volcano, Holocene age; d) older cone with deep gullies on external slopes and partly filled crater. Huipilo volcano, estimated Mid to Late Pleistocene age; e) simple cone with smooth and low-angle outer slopes, partly-filled crater and single, deep gully to the SE. Tulmeac volcano; f) example of cones standing on the slope of a shield, which were partly buried by young lava flows. g) Simple cone (young). Perfect ring-shaped, symmetrical single scoria cone; h) Notched cone with coalesced craters and small lava flow (also exhibits traces of mining-small quarry-on the E flank); i) scoria cones with multiple coalesced craters, and cone with wide crater, interpreted as having an higher phreatomagmatic component, and small notched cone; j) older scoria cone with coalesced craters and deep erosion rills (left) and not clearly recognizable older structures (center and right); k) compound scoria cone formed in a maar and truncated cone to the left. l) Contemporary 2 maar-2 scoria cone complex; m) very small scoria cones (red dot) within a maar crater formed after an initial phreatomagmatic eruptive phase; n) scoria cone cluster with heavy vegetation in the center area of the field; poor Lidar quality. Example of a case where it is difficult to workout the cones' shape and degree of complexity; o) example of a high density area near the center of the fault system with several shapeless, older scoria cones that were excluded from the analysis ◂ Table 4 Range of morphometric parameters in the studied areas (parameter abbreviations as in Table 2 (Table 4), all of them higher than the horizontal resolution (5 m) and vertical accuracy (1-2 m) of the DEMs used to calculate them. Standard deviations are high, an indication of the shape diversity. Slope and flatness mean values are also smaller in the LTVF (Fig. 7), but the difference in slope is too small (2°) to be significant (Table 4). This suggests that cones in the LTVF have smaller craters relative to their basal diameter than those in the SCVF. However, the heightto-diameter ratio (Ar) distribution is the same in both volcanic fields (Fig. 7). Additionally, most of the cones in the LTVF have values lower than the mean value as can be seen in Table 4 (i.e., histograms are skewed to the left). Flatness seems to have a bimodal distribution in both fields, which is more accentuated in the SCVF (Fig. 7f), Table 4). The first (lowest) modal value in the flatness histogram is equivalent in both volcanic fields, whereas the second (highest) modal value is greater in SCVF (Fig. 7f). Figure 6h illustrates the difference in the average size proportions (flatness and Ar) of a cone in the LTVF compared to the SCVF. It shows that LTVF cones have smaller craters to base ratios in average, compared to those of the SCVF. Figure 7i and j illustrates the differences in the same parameters but for the modal values, considering two modes for the flatness, which suggest the existence of two groups of cones. These 2 groups exist in both fields, and represent cones with similar heightto-diameter ratios but strong differences in crater diameter; this is much more evident in the SCVF.

Statistical analyses
The Kolmogorov-Smirnov test was used to test the equivalence of the datasets of the two volcanic field in terms of each of the morphometric variables (Table 5). The test shows no-equivalence for all datasets, but Ar as mentioned above, which means that the datasets of the other morphometric parameters do not belong to the same distribution. The correlations and covariances matrices of each volcanic field, shown in Fig. 8, are noticeably similar. In both, LTVF and SCVF, the strongest positive correlations exist between volume and cone basal diameter (SCVF: 0.97, LTVF: 0.99), followed by volume and cone height (SCVF: 0.92, LTVF: 0.93), and cone height and cone basal diameter (SCVF: 0.81; LTVF: 0.87) (Fig. 8). Crater width and volume are well correlated in the LTVF (0.80) and moderately well correlated in the SCVF (0.66) (Fig. 8). Slope and cone height are moderately well correlated, whereas slope is poorly correlated with all the other parameters (Fig. 8).
The results of a PCA on each of the covariance matrices shown in Fig. 8 are displayed in Table 6. Even though it is possible to numerically identify similarities and differences between the PCA of each covariance matrix, those are small. When using a biplot representation of the same data (Fig. 9), a similar covariation structure is immediately visually identified. The variance and covariance patterns of the parameters of scoria cones, represented by vectors drawn as lines in the biplots, have very similar behavior: same orientation, same length, and a similar angle between them. The volume vector is sub-horizontal and sub-parallel to the first component (PC1 in Fig. 8), while the slope vector is the most vertical and is sub-parallel to the second component (PC2 in Fig. 8). The volume vector is also the one that most contributes to the PC1 (greater length along the x-axis), whereas for PC2, the crater width has the greater contribution (greater length along the y-axis). This means that the variance of the morphometric parameters of both fields, which is a measure of the dispersion, is  The similarity of the biplots of the LTVF and SCVF implies that the morphometric parameter-relations vary in the same way, and thus the processes responsible for these variations could be similar. To confirm this hypothesis, we conducted further formal testing using the CPC model.
In Fig. 10, we compare the two null distributions of the test statistic for H CPC , the blue line is the chi-square density with 10 degrees of freedom and the permutation distribution is the density histogram, based on 10,000 resamples. The observed value of the test statistic is x * = 32.31 . Hence, the asymptotic p-value is P 2 df =10 ≥ 32.31 = 0.000356 while the permutation p-value is as follows: with I as the indicator function. The red line is located on the value of the observed test statistic, so we can appreciate the difference between both p-values to the right of the red line, which correspond to the area under the curve or histogram considering the asymptotic test and permutation test respectively. The asymptotic p-value strongly indicates that the CPC hypothesis has to be rejected (blue curve almost at 0 to the right of the red line). This decision (5) p * = 1 + ∑ 10,000 i=1 I(x * i ≥ 32.31) 10, 000 + 1 = 1 + 2, 319 10, 001 = 0.232  . 9 Biplot showing the correlation between individual morphometric parameters of each field. a) The correlation matrixes of both fields are plotted in distinct colors. b) Same as a) but, in addition, the vectors calculated combining the database of both fields are reported (black line), to show their similarity stands in contrast with the evidence provided by the permutation p-value where the decision would be not to reject the CPC hypothesis (histogram representing the permutation test is still present to the right of the red line). Why is this happening? First of all, we have to keep in mind that the chi-square distribution is an asymptotic approximation to the true, but unknown, null distribution of the test statistic 2 (likelihood ratio). The validity (usefulness) of this approximation also depends on the normality assumption.
In simpler words, it depends on how well the multivariate normal distribution describes the distributional behavior of the cone's parameters. Therefore, we do not discard a potential poor approximation by the chi-square distribution to the true null distribution of the test statistic. This poor approximation, together with the non-normality of the data, may therefore be affecting the exact p-value. On the other hand, the permutation test does not require the normality of the data. Most importantly, it is an exact test (when all the possible permutations are resampled) in the sense that the p-value is exact. Thus, for the CPC hypothesis, we will trust the permutation test over Flury's asymptotic test. We show that the asymptotic test is not adequate for treating the matter of comparison of morphometric parameters of two volcanic fields as it supposes an underlying normal distribution of the data, while on the contrary, the permutation test is adequate. Table 7 summarizes the common principal components output. These results were obtained using the R package multigroup from Eslami et al. (2020). In principle, the maximum likelihood estimates of the covariance matrices should be preferred over the sample covariance's matrices ( Table 7). The common eigenvectors and their respective specific eigenvalues are provided too. Since we now require only one set of eigenvectors, rather than two sets, we can claim that a more parsimonious description of the covariances is being provided. Furthermore, as long as the CPC model holds for the covariances of LTVF and SCVF, it is reasonable to think that the covariances of the morphometric characteristics of the scoria cones in both volcanic fields are governed by similar underlying processes.

Discussion
Monogenetic fields generally present a great diversity in morphology whose origin is poorly known (e.g., Connor and Conway 2000;Valentine and Gregg 2018;Le Corvec et al. 2013;Németh and Kereszturi 2015). We studied two large subduction-related fields in Mexico that have distinct tectonic, compositional, and climatic characteristics (Table 1). In the following, we first briefly summarize the similarities and differences of these fields in terms of vent distribution and cone shapes and then discuss the factors that can explain those.

Similarities and differences
The two fields have numerous similarities. They both present elongated outlines and density probability distributions (Fig. 4). The main direction of elongation coincides in both cases in the orientation of active crustal faults, irrespective of their type (normal for SCVF vs. strike-slip for LTVF). Another similarity is the tendency of vent clustering in the fields, expressed by large spatial variations in vent intensity (Fig. 5). The main directions of alignment of vents are parallel to the main faults in both fields (Table 3; Figs. 2 and 4). Both fields display, as well, a wide variety of scoria cone shapes that vary from simple to compound and commonly show crater rows, notched cones and cone overlapping (Fig. 6). Cone shape parameters present similar ranges (Fig. 7) and strikingly similar correlations (Fig. 8) in both fields, and the volume and crater width are the parameters that account for most of the shape variety in the two fields (Fig. 9). Cones have similar height-to-diameter ratios and slopes in both fields (Fig. 7), despite external factors which vary strongly (precipitation values, vegetation cover). The age range of the studied scoria cones in both fields is similar.
The two fields present, nevertheless, some significant differences. The LTVF is smaller in area than the SCVF, even though the Tuxtlas massif itself is as large as the SCVF, but the active part comprises approximately an area half as large. The LTVF has however a much higher number of vents, and hence a significantly higher vent intensity, evidenced in both kernel and vent intensity exercises (Figs. 4 and 5). This is associated to more common vent clustering, creating a higher complexity in cone shapes (compound cones). This clustering in the LTVF seems to be, in many cases, related to the formation of several edifices during a single event. In addition, there is a higher percentage of small cones in the LTVF than in the SCVF (histograms skewed to low size values), which causes lower averages in cone sizes. An important difference of the LTVF is the common occurrence of phreatomagmatic constructs that are often associated with scoria cones of the same age. Additionally, kernel density estimations show that in the SCVF, there are outliers to the general vent distribution, which is caused by the occurrence of vents outside of the main mountain range.

Tectonic regime and pre-existing faults
According to the KDE contours and intensity pattern of vents (Figs. 4 and 5), most vents are concentrated in the center of both fields and coincide with the orientation of active faults. The alignment of vents also coincides with such faults (Table 3), which could be studied further by a fully quantitative analysis. This suggests that there was a coincidence between the dominant direction of magma emplacement and the orientation of the active faults, which may be due to the capture of dikes by faults at shallow levels or to the independent rise of dikes according to the active stress regime (e.g., Connor and Conway 2000;Valentine and Gregg 2008). The structural framework is different in the two fields (Table 1), which suggest different dike emplacement mechanisms. In the SCVF, the vents are mostly distributed E-W (Table 3), which coincides with the active N-S extension and is coherent with preferred magma rise along the maximum horizontal stress direction (Nakamura 1977). Hence, it is possible that dikes formed their own pathways in this area, only following existing faults when they occurred in their path. Evidence for this could be the lack of fault expression in the most-active central part of the field (Fig. 2) which may be due to the accommodation of extensional forces by frequent dike injections (Valentine and Gregg 2088). This may be also caused by the burial of fault traces by younger volcanic products; however, we note that this does not occur in the LTVF, despite a higher vent spatial intensity there. The situation is different for the LTVF that occurs in a transpressive regime (Andreani et al. 2008). The formation of a monogenetic field in such context is somewhat unusual as most distributed fields occur in extensive settings (Valentine and Gregg 2008) but it has been reported elsewhere (Grosse et al. 2020;Uslular et al. 2021). In this field, fault expressions are abundant. In such context, it is more likely that the dikes used pre-existing fault structures rather than making their own pathways. The faults are believed to join at depth, to form a flower structure (Andreani et al. 2008). This geometry could create a flexure of the upper crust, which could have opened tension gashes that may have favored dike injections. Another simpler hypothesis is that sheared, fractured zones along the faults provide low-density areas that allow dike ascent (Sieron et al. 2021). Related mechanisms are described in Tibaldi (2005). Other hypotheses include the formation of small pull-part basins (Espindola et al. 2016) but this is unlikely provided as pull-apart basins usually form in strike-slipextensional settings (e.g., Farangitakis et al. 2021), which is not consistent with the compressive setting and the structure of the uppermost crust constrained by seismic data (Franco et al. 2013;Andreani et al. 2008;Sieron et al. 2021).

Magma source geometry
The differences in size (area) and vent spatial intensity of the two studied field (Fig. 5) are probably due to a different geometry of the magma source (e.g., Connor and Conway 2000). In the SCVF, geophysical data and modeling indicate that the magmas rise directly from an area of partial melting in the mantle wedge generated by the dehydration of the plunging subducting slab (e.g., Ferrari et al. 2012). Hence, the field location and size seem to be controlled mostly by the slab geometry, rather than by crustal structures that only provide favorable pathways for the rising magmas. Nevertheless, the thick crust probably causes significant trapping of the dikes on route to the surface, as evidenced by the dominant intermediate magma compositions (Schaaf et al. 2005).

Fig. 10
Chi-square density with 10 degrees of freedom and permutation distribution based on 10,000 resamples. The red line is located on the observed value of the test statistic x * = 32.31. The blue line represents a parametric Chi-square density with 10 degrees of freedom, showing the importance of using a nonparametric test This context can explain the relatively large size of the field (extended magma source) and dike stalling in the crust may account for the relatively low vent spatial intensity, although this could also be caused by a relatively low magma supply rate. The LTVF lies in a very different context. The occurrence of the active San Martín Tuxtla polygenetic volcano within the monogenetic field suggests the existence of one or several large magma bodies underneath, which is supported by the gravity data (Espindola et al. 2016) and recent seismic data (Castellanos et al. 2018;Calò 2021). These probably act as magma reservoirs from which the dikes feeding the monogenetic vents radiate, following vertical strikeslip faults created by crustal stresses. The magma source is hence more localized and shallower than for the SCVF, which can account for the smaller field size, and dike injections are more closely spaced, explaining the higher vent spatial intensity and complex overlapping cone shapes (large number of compound cones and complex cone and maar structures). A thinner crust, in comparison to the SCVF, can also affect the clustering (Mazzarini et al. 2010) and facilitate magma ascent (lower lithostatic pressure). A higher magma supply rate is also a possible cause but there are no data to test this. Moreover, the ascending magma seems to be distributed to form multiple edifices along the ascending dikes, which probably reduces the final edifice size as the total volume of pyroclastic material produced by a single eruption is split into several edifices. Recent studies have shown that feeder dikes can follow complex paths, segment near the surface and within the edifices, to create complex volcanic edifices (e.g., Blaikie et al. 2014;Heimisson et al. 2015;Sigmundsson et al. 2015).

Eruptive style
Despite these distinct magmatic and structural frameworks, scoria cones in the SCVF and LTVF are surprisingly similar, in terms of shape parameters (volume, cone height). It confirms previous results that the diversity in cone shapes reflect processes that are common to all cones (e.g., ballistic emplacement followed by scoria avalanching on slopes) and hence do not vary with external parameters (Bemis and Ferencz 2017). An important outcome is that those characteristics cannot be used on their own (without over independent information), to infer characteristics such as the tectonic regime, magma composition, and climatic conditions (Tibaldi 1995;Bemis et al. 2011). Nevertheless, despite these similarities, our study reveals that there is a higher proportion of small cones (Hco < 100 m, Wcr < 200 m, Wco < 500 m) in the LTVF, which decreases the average values of the cones and craters dimensions but does not seem to affect height-to-diameter ratio and slope ( Fig. 7; Table 2). We can exclude weathering processes as a cause because these should also affect the slope and aspect ratio (Hooper and Sheridan 1998). In addition, they seem to cause significant changes only over long timescales (Dohrenwend et al. 1986), especially in vegetated areas where cone slopes are rapidly stabilized after eruption (Bemis et al. 2011). Cone burial by later-erupted lava can significantly decrease the size of the cones (Favalli et al. 2009), which may be an important factor in the LTVF, due to the high vent intensity and pronounced slopes (due to its construction on a large shield). Yet, it is very unlikely to affect the crater width, which is much lower in average in the LTVF. Besides, it affects more importantly the height than the width of the cones, and hence also would change the height-to-diameter ratio (Wood 1980a,b;Favalli et al. 2009), which does not fit with our observations. It is then more likely that this difference in size reflects primary processes, which are those that take place during the eruptions. Based on the previous section (magma source geometry), we can hypothesize that the eruption of smaller cones in the LTVF is linked to the generation of smaller magma batches from a shallow magma reservoir,in comparison to a deeper source for the SCVF, where smaller magma batches are more likely to be trapped on route. Nevertheless, we acknowledge that lavas are a dominant component in the total erupted magma in both fields (e.g., Sieron et al. 2021;Siebe et al. 2004), which means that the smaller cones may only mean a lower general explosivity (more volume emitted as lavas).
Another hypothesis arises from the observation of the occurrence of tiny scoria cones in previous-formed maar craters in the LTVF (Fig. 6b). This is interpreted as resulting from a short-lived dry magmatic phase at the end of a mostly phreatomagmatic eruption (e.g., Kshirsagar et al. 2016;Ort and Carrasco-Núñez 2009). Such features are non-existent in the SCVF due to the lack of maars. The occurrence of initial phreatomagmatic activity may hence be a cause of the higher proportion of small cones in the LTVF, which would then reflect lower-magnitude, shorter-lived gas-driven explosive activity. Sieron et al. (2021) note however that maar features are restricted to a specific area of the field but a recent revision of this database suggests a much higher abundance (note that the maars were excluded from the morphological analysis herein presented). For example, phreatomagmatic activity either at the start or end of an eruption is often underestimated (see Lorenz et al. 2020). In the LTVF, approximately half of the area is strongly affected by phreatomagmatism and the other half only slightly. Hence, this type of activity is probably sufficiently frequent to affect average cone sizes.
Another hypothesis may be a link with the magnitude of the eruptions (explosivity), as hinted before. Sieron et al. (2021) note that smaller cones tend to be built of coarser material and hence formed by less-energetic eruptions. In comparison, they note that larger cones are made of finer grained pyroclastic material without bombs, which, together with their generally smaller crater size and large height and volume, points towards more effective fragmentation (violent Strombolian activity) and higher magma eruption rates (e.g., Wood 1980a;Pioli et al. 2008;. Although this should be confirmed by a more systematic study, this correlation suggests that the cone-forming activity in the LTVF may be less energetic, in general, than in the SCVF. Magma composition can impact eruption intensity by acting on magma viscosity, mass flux or ascent rate, and gas pressure built-up (e.g., Rutherford and Gardner 2000). Therefore, the more evolved composition of magmas feeding cones in the SCVF in comparison to the LTVF (Table 1) may be responsible for their larger sizes in average (higher eruption explosivity and mass flux) and possibly, a lower proportion of magma emitted as lavas, which should nevertheless be tested by further studies. The high volatile content of magmas in the SCVF (Cervantes and Wallace 2003;Roberge et al. 2015) may also have favored higher explosivity, but recent studies have shown that magmas in the LTVF are also relatively wet (Díaz-Bravo et al. 2022). As noted before, in the LTVF, cones are commonly arranged in long alignments of coalesced edifices that suggest high fissure lengths and a more distributed eruption of magma at the surface, which may also lead to smaller cone sizes.

Conclusions
The Sierra Chichinautin and Los Tuxtlas volcanic fields occur in distinct sectors of the Trans-Mexican Volcanic Belt, in regions with different tectonic regimes, magma compositions, volcano types, and climatic conditions. Our study reveals that, in terms of vent distribution and cone morphology, the two fields present significant similarities, which shows that this type of data cannot be used alone to infer the tectonic, magmatic, and climatic context of monogenetic fields. Volume and crater diameter are the parameters that express the overall morphological variety (principal components PC1 and PC2 in the statistical analysis) and could be used as a first-hand indicator for cone diversity in a volcanic field. When considering all the existing information, it is however possible to infer the respective role of different factors in the morphology and spatial distribution data. Differences in the geometry of the magma source probably caused the distinct size, vent density, and degree of cone complexity of the fields because it determined the closeness (spacing) of the dikes feeding the activity. Differences in the crustal stress regime (transtensional vs. transpressive) had no apparent impact on vent distribution as the dikes followed active faults, irrespective of their motion (strikeslip vs normal), indicating that, in this case, data on crustal stress cannot be provided by vent alignments. Climatic differences between the two fields (temperate vs. tropical) did not affect the shape variety of the studied cones probably because of their young and similar ages (< 50,000 y old) and their location in a vegetated environment. Therefore, systematic differences in the cone sizes (higher proportion of small cones in Los Tuxtlas) are more likely due first, to differences in erupted volume, as controlled by the magma source and storage depth (shallow reservoir in Los Tuxtlas vs. deep mantle source in Sierra Chichinautzin) and possibly the magma supply rates, but this latter cannot be tested because of the lack of estimates for both fields. Then, it can depend on eruptive style caused either by differences in initial magma/water interaction, magma composition (basaltic vs andesitic) and explosive/effusive behavior, and/or fissure lengths. This study shows how statistical studies on vent distribution and cone morphology in monogenetic fields allow to improve our knowledge on the factors that control the location and style of eruptive activity in these areas, contributing to improve hazard assessment.
Acknowledgements Gustavo Vivo Vasquez and Andrea De San Jorge Cárdenas are thanked for contributing to the inventory of vents in the SCVF and the preliminary analysis of cone shapes. Project UNAM-DGAPA-PAPIIT IN103421 supported the studies in the SCVF. We further acknowledge the contribution of the Universidad Veracruzana in terms of a full salary sabbatical grant for Katrin Sieron. The constructive comments of the editor and three anonymous reviewers were greatly appreciated and contributed highly to the quality of this paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.