Generating a seismogenic source zone model for the Pyrenees: A GIS-assisted triclustering approach

Generating a seismogenic source zone model for the Pyrenees: A GIS- assisted triclustering approach José L. Amaro-Mellado a,b, Laura Melgar-García c, Cristina Rubio-Escudero d,David Gutiérrez- Avilés c,∗ a Department of Graphic Engineering, University of Seville, ES-41092, Seville, Spain b National Geographic Institute of Spain, Andalusia Division, ES-41013 Seville, Spain c Data Science & Big Data Lab, Pablo de Olavide University, ES-41013 Seville, Spain d Department of Computer Languages and Systems, University of Seville, ES-41012, Seville, Spain Keywords: Seismic sources GIS Data science Triclustering A B S T R A C T Seismogenic source zone models, including the delineation and the characterization, still have a role to play in seismic hazard calculations, particularly in regions with moderate or low to moderate seismicity. Seismic source zones establish areas with common tectonic and seismic characteristics, described by a unique magnitude–frequency distribution. Their definition can be addressed from different views. Traditionally, the source zones have been geographically outlined from seismotectonic, geological structures, and earthquake catalogs. Geographic information systems (GIS) can be of great help in their definition, as they deal rigorously and less ambiguously with the available geographical data. Moreover, novel computer science approaches are now being employed in their definition. The Pyrenees mountain range – in southwest Europe – is located in a region characterized by low to moderate seismicity. In this study, a method based purely on seismic catalogs, managed with a GIS and a triclustering algorithm, were used to delineate seismogenic zones in the Pyrenees. Based on an updated, reviewed, declustered, extensive, and homogeneous earthquake catalog (including detailed information about each event such as date and time, hypocentral location, and size), a triclustering algorithm has been applied to generate the seismogenic zones. The method seeks seismicity patterns in a quasi-objective manner following an initial assessment as to the best suited seismic parameters. The eight zones identified as part of this study are represented on maps to be analyzed, being the zone covered by the Arudy–Arette region to Bagnères de Bigorre as the one with the highest seismic hazard potential. 1. Introduction Seismogenic source zones are necessary for applying the most widely used seismic hazard calculation method proposed by Cornell (1968) and McGuire (1976). Originally, this method contemplated the zones as an alternative to model the probabilistic space of distance-to- source relationships with the area under analysis, particularly when associations between seismicity and faults capable of triggering an earthquake are not apparent (García-Mayordomo,2015). In addition, the method conceived a fault as a linear element, but it has evolved to the use of fault as zones (even as 3D planes). This second option is much more widely used today. Currently, there is a tendency to make greater use of faults, but through a characterization collected in a database. These records include geological and seismic data, such as the slip rate, the fault length, the depth, recurrence periods, or maximum expected magnitude (García-Mayordomo et al.,2012a). ∗Corresponding author. The main assumptions of the Cornell-McGuire method regarding the occurrence of seismicity are the following (García-Mayordomo,2015): 1. Earthquakes are independent random events (Poissonian). 2. The probability of earthquakes occurrence is the same throughout the entire area. 3. The size of earthquakes is related to their frequency through the Gutenberg–Ritcher Law (Gutenberg and Richter,1944) of each source zone, usually limited to a maximum magnitude value. 4. The seismic activity rate of each zone is constant over time. These cannot be strictly fulfilled in practice since the mechanism of accumulation and release of energy in earthquakes is a memory process with long-range dependence (Barani et al.,2018), and the rate of seismic activity is variable over time. However, in seismic engineering, E-mail addresses: [email protected] (J.L. Amaro-Mellado), [email protected] (L. Melgar-García), [email protected] (C. Rubio-Escudero), [email protected] (D. Gutiérrez-Avilés). 2 seismogenic source zones are sufficiently adequate to estimate the probability according to different levels of ground shaking (García- Mayordomo, 2015). In this case, a memory process refers to the fact that shocks do not happen independently at any time. By contrast, they remember the occurrence (time and magnitude) of the last earthquake (Corral, 2006). The main strength of using zones is that they are efficient because their associated parameters are somehow easily obtained, and their calculations are fast. Thus, they have been used for over 50 years so far. Seismic source delimitation is a crucial task when a Probabilistic Seismic Hazard Analysis (hereinafter, PSHA) is conducted (Morales- Esteban et al., 2014). Once they have been geometrically outlined, their activity is characterized based on the recorded seismicity. Both stages define the seismic potential (geometric shape and size, maximum expected magnitude, seismic activity, and frequency distribution in size). Subsequently, the equations for predicting strong movements should be considered, and, finally, the annual probability of exceeding a certain magnitude earthquake should be obtained. Despite the importance of the seismogenic zonation to conduct a seismic hazard analysis, the criteria used for its definition are often not uniform. It is a very complex process and depends strongly on the analyst group’s criteria (Amaro-Mellado et al., 2017). There are two main groups within these criteria: one based exclusively on the distribution of seismicity and the other on geological domains or structures. The criterion, or its weight, depends on the background of the analyst who defines the zones. On the one hand, much weight is placed on the characterization of sources, it will mainly be based on seismicity. This can lead to inconsistent limits from a geological point of view but valid for the calculation of hazards. On the other hand, if the analyst relies on the location of the geological structure, two opposing situations arise. In the first, areas are based exclusively on the limits of large geological units, even if they have no relation to current tectonic activity. In the second, they only consider areas with active faults, even though the recorded seismicity does not allow for statistically significant characterization. Computational robust methods are becoming increasingly popular, as they can be applied in different geographical and geological contexts (Morales-Esteban et al., 2014; Martínez-Álvarez et al., 2015; Scitovski, 2018). Besides, when handling and representing georeferenced spatial information such as the distribution of epicenters, the use of a geographic information system (hereinafter, GIS) is powerful. Its ability to manage, represent, and generate geographical information is a useful tool for integrating different geographic data kinds. This work aims at proposing a seismic source model for the Pyrenees range in southwest Europe. To this end, the use of a GIS and a tricluster algorithm, called TriGen (Gutiérrez-Avilés et al., 2014), has been applied to derive seismicity patterns from a seismic catalog. Later, the source models have been drawn. Finally, seismic parameters have been obtained from these zones, and they have been analyzed. The range of seismicity in the Pyrenees is estimated as low to moderate and diffuse, where earthquakes with a magnitude larger than or equal to 5.0 (M5+) are infrequent. Besides, in these areas, obtaining a correlation between epicenters and faults is difficult (Drouet et al., 2020). Therefore, seismic zonings usage to drive a PSHA is adequate (Martínez-Álvarez et al., 2015; Amaro-Mellado and Tien Bui, 2020). 2. Related work In this section, research made in both domains presented in this study is described: firstly, the studies in seismogenic zonings and, secondly, the triclustering approach. 2.1. Previous seismic zonations for the Pyreenees The Pyrenees are located in a border-area between France, Spain, and Andorra, so each of these countries have studied their seismicity. They are also included in some European seismicity projects as SHARE (Woessner et al.,2011;Stucchi et al.,2013), or SIGMA (Pecker et al.,2017). In the Spanish seismic sources, the Pyrenees are an important region due to their activity, and they are included in different seismic zonations. For example, that conducted by Martín (1984), who established 27 zones for the whole Iberian Peninsula and adjacent area; Mezcua et al. (2011) addressed a PSHA as a combination from the one defined by Molina (1998) and the CODE one (from NCSE-2002 Ministerio de Fomento (Gobierno de España),2002, which is very similar to Martín, 1984). More recently, in the frame of the Updated Map of PSHA for Spain, both García-Mayordomo et al. (2012b) and Bernal (2011) defined a set of seismic sources. Later, Morales-Esteban et al. (2014) Morales-Esteban et al. undertook a zonation for the whole Peninsula from Mahalanobis distances, based only on an earthquake catalog (including information on the shock such as date and time, hypocentral location, or size) from 1978 to 2012. In Martínez-Álvarez et al. (2015), Martínez-Álvarez et al. used the triclustering methodology to obtain seismic sources for the whole Iberian Peninsula. French institutions have also defined different seismogenic zonations. In Metropolitan France, the first PSHA map was deployed in 2002 (Martin et al.,2002) and gave the foundations for the French zoning (Drouet et al.,2020). Moreover, in 2004, Marin depicted PSHA maps for France (Marin et al.,2004). Later, Baize et al. (2013) proposed a new seismotectonic zoning for Metropolitan France based on geological, seismological, and tectonic data. Recently, Drouet et al. (2020) constructed a seismotectonic model (GEOTER) for Metropolitan France, considering geological, structural, neotectonic, geophysical, and seismological data. Besides, some seismic zonations have been conducted specifically for the Pyrenees. In Njike-Kassala et al. (1992), the authors established nine zones to calculate the b-value of the Gutenberg and Richter relation (Gutenberg and Richter,1944). Later, Secanell et al. (2008) defined a seismotectonic zonation as a result of the ISARD project, held by France and Spain. Finally, Rigo et al. (2015) proposed eight zones to estimate stress tensors locally. Finally, the Institut d’Estudis Andorrans holds an internet site (https: //www.iea.ad/sismoweb, last accessed on November 2020) with information on Andorran seismicity. 2.2. Triclustering Regarding the second area of study of this research, the state of the art of triclustering techniques is presented. The traditional clustering approach is a data mining technique that works by grouping objects according to a predefined similarity in a one-dimensional space. When dealing with subspace clustering, the problem can be addressed with biclustering (Pontes et al.,2015) or triclustering (Zhao and Zaki, 2005) if the context data is structured over two or three dimensions, respectively. In the last decade, the tricluster algorithm has evolved, and many different algorithms have been proposed. For example, some triclustering algorithms are based on iterative searches, distribution parameter identification, biclustering algorithms, pattern mining procedures, or evolutionary multiobjective optimization (Henriques and Madeira, 2018). In particular, in Liu et al. (2008) and Bhar et al. (2015) some triclustering algorithms based on multiobjective techniques can be found. Furthermore, classic bio-inspired meta-heuristics as genetic algorithms (Holland,1992) or newer, as the COVID-19 propagation model (Martínez-Álvarez et al.,2020) emerge as a proper method to optimize multiobjective functions. In this sense, the TriGen algorithm presented in Gutiérrez-Avilés et al. (2014) is based on a population 3 of individuals and some genetic operators. Different fitness functions are also proposed in order to assess each research problem most accurately. TriGen is the triclustering algorithm used in this study and is further explained in 4.3. Besides developing triclustering methodologies, another open line of investigation is to measure the quality of the triclustering algorithms’ solutions. In this sense, we can find correlation-based measures inspired by Pearson and Filon (1898), and Spearman (Spearman, 1910) correlation indexes. In Gutiérrez-Avilés and Rubio-Escudero (2014) the authors propose a three-dimensional version of the mean squared residue, a classical biclustering evaluation measure (Cheng and Church, 2000). An evaluation measure based on the last squared lines depicted by the discovered patterns in the triclusters was presented in Gutierrez-Aviles and Rubio-Escudero (2014). Finally, in Gutiérrez-Avilés and Rubio-Escudero (2015), the authors present a new approach to measure the tricluster quality based on the graphical representation of it. Triclustering has many applications: biological as Li and Tuck (2009) that combines gene regulator information with expression data; medical as Melgar-García et al. (2020) that applies triclustering in the streaming environment to high-content screening images or the analysis of social networks as Gnatyshak et al. (2012). The domain of seismic data has not yet been deeply analyzed with triclustering techniques. To the best of our knowledge, only the work in Martínez- Álvarez et al. (2015) applies triclustering techniques, particularly, the TriGen algorithm, to data of the Iberian Peninsula, in order to discover possible seismogenic zones. 3. Seismicity and geological settings The Pyrenees mountain range, constituted from the collision between Iberian and Eurasian plates, was formed in Alpine times (Vissers and Meijer,2012); however, current deformation rates are relatively low. Spanning over 450 km in E-W direction and 150 km in N-S (Rigo et al.,2018), the Pyrenees have nearly cylindrical symmetry (Njike- Kassala et al.,1992). It is bordered to the north by the North Pyrenean Fault, which matches a 10–15 km vertical Moho offset, where the crust is thicker on the Iberian plate (Gallart et al.,1981). The North Pyrenean Fault extends along the whole chain in E–W direction, and it is supposed to be the boundary between the two plates (Njike- Kassala et al.,1992), and currently, between the well-known North Pyrenean Zone and the Axial Zone. Upon Vissers and Meijer (2012), the Pyrenees are characterized by, from north to south: first, a northdirected thrust-belt including Mesozoic and Tertiary sediments of the Aquitanian Basin; second, the North Pyrenean Zone, formed mainly by Paleozoic basement and Mesozoic sediments; third, the Axial Zone consisting of Paleozoic rocks; and, finally, the southern Pyrenees where the southward thrust Mesozoic and Tertiary sediments prevail, and they part of the Ebro foreland basin system. Regarding seismicity, it is very complex (Souriau et al.,2014), and it can be considered as low to moderate (Amaro-Mellado and Tien Bui,2020). In the western part, the seismicity is generally concentrated along the North Pyrenean Zone. In contrast, it is more diffuse in the eastern part and not particularly associated with the known faults (Njike-Kassala et al.,1992). Besides, the identification and parametrization of seismogenic structures have to face some challenges. The most important ones are the lack of surface faulting related to shocks in the last centuries, the low deformation rates, or the fact of being a political-border region (partially solved by the projects above mentioned) (Lacan and Ortuño,2012). Over time, as can be deduced from earthquake catalogs, some shocks have produced severe damages, mainly in the historical period, even with estimated magnitudes between 6 and 7 or with MSK intensities between VIII and IX (Lacan and Ortuño,2012). Among others, the 1373 earthquake (𝐼𝑜=𝑉 𝐼𝐼𝐼 −𝐼𝑋) by Maladetta Massif, in the Central Pyrenees; those happened during 1427 Amer and 1428 Queralbs, in Catalonia (up to 𝐼𝑜=𝐼𝑋), in the eastern Pyrenees; the 1660 Bigorre earthquake (𝐼𝑜 = 𝐼𝑋), near Lourdes, and the most recent (1750) in Juncalas (𝐼𝑜 = 𝑉 𝐼𝐼𝐼), close to the Bigorre one. In the instrumental period, some relevant shocks (magnitude equal to or larger than 5) have occurred, such as 1967 Arette (M5.3), in the western Pyrenees, and 1996 St-Paul Fenouillet (M5.0), in the Agly Massif, in the east (Amaro-Mellado and Tien Bui, 2020). Some specific studies have been driven to study the seismicity of the whole Pyrenees Njike-Kassala et al. (1992), Secanell et al. (2007), Rigo et al. (2018), Amaro-Mellado and Tien Bui (2020). Besides, the Arudy– Arette region has been analyzed in-depth by Gallart et al. (1980) and Sylvander et al. (2008). Fig. 1 shows a general view of the major events that have shocked the Pyrenees from the raw data obtained from the earthquake catalog published by the Instituto Geográfico Nacional (IGN)–National Geographic Institute of Spain (NGIS) (Insitituto Geográfico Nacional, 2020). The map only represents the events with a magnitude larger than or equal to 3.5, or macroseismic intensity, 𝐼𝑜, larger than or equal to ‘‘V’’ in macroseismic scale. If both data appear, the magnitude one prevails. 4. Methods In this section, the methodology driven to generate the seismogenic sources for the Pyrenees are described. Thus, the seismic dataset used in this study is presented; later, the data processing undertaken in this work is explained; and finally, a summary of the TriGen algorithm is set out. 4.1. Generation of the earthquake catalog The dataset employed in this study comes from the NGIS earthquake catalog (Insitituto Geográfico Nacional, 2020). In-depth research on the NGIS catalog is presented in González (2017). In particular, that compiled in Amaro-Mellado and Tien Bui (2020) has been chosen in this work. In this catalog, starting from the NIGS earthquake database, different steps were taken to generate an updated, reviewed, homogeneous, declustered, and extensive catalog, including events with a magnitude larger than or equal to 2.0; from 1373 to 2019; and the geographical extension was limited by 2.5◦ W and 3.5◦ E meridians, and 41◦ N and 44◦ N parallels. The steps driven to compile the catalog were the following: 1. Data from other sources, such as specific studies or other catalogs, were considered to re-evaluate the magnitude of some shocks, particularly the major ones. 2. To conduct a seismic analysis, the size of the events must be homogeneous. In this sense, the vast majority of the researches use 𝑀𝑤(Hanks and Kanamori,1979), due to its direct relationship with the released energy in the rupture by scalar seismic moments. Besides, it does not get saturated for the largest events (Das et al.,2019). Over time (from 1373 to 2019), different kinds of sizes have been assigned to the earthquakes (macroseismic intensities and several magnitude types). Herein, to homogenize the size, 𝑀𝑤was assigned to all the events, according to Cabañas et al. (2015). Besides, after estimating a lower magnitude of completeness (1.8), the cut-off magnitude was set as 2.0. 3. As stated in Section 1, the earthquakes must be independent events to fulfill the Poisson distribution. To this end, foreshocks, aftershocks, and seismic swarms must be eliminated. In seismicity, this process is called declustering. In that work (Amaro- Mellado and Tien Bui,2020), it was based on the definition of temporal and spatial windows following the method proposed by Gardner and Knopoff (1974). 4 Fig. 1. Seismicity in the Pyrenees from the raw NGIS earthquake catalog. Events instrumentally recorded (Circles); events with only macroseismic intensity (triangles). Above referred earthquakes have been numbered: 1—Maladetta(1373); 2—Amer(1427); 3—Queralbs(1428); 4—Bigorre(1660); 5—Juncalas(1750); 6—Arette(1967); 7—St-Paul Fenouillet(1996). 4. Although it is not a specific task in a catalog generation, the year of completeness related to different levels of magnitude (𝑀𝑐—year of completeness) were defined: M2—2013; M3— 1978; M4—1943; and, M5—1810. 5. The final catalog, after removing the events deeper than 65 km, since they are not important for the seismic hazard of the region (IGN-UPM-WorkingGroup,2013), consists of 7706 earthquakes. The present work aims at defining a source zone model, so such a low cut-off value as 2.0 should not be established. Therefore, the cut-off magnitude has been set as 2.5, as pointed out in other works (Njike- Kassala et al.,1992;Talbi et al.,2013). Besides, it has been considered in a recent work for the region (Drouet et al.,2020). This value is judged as a good trade-off between the seismic hazard and the low to moderate Pyrenees seismicity. The final catalog of this work consists of 3500 earthquakes, with 𝑀𝑤larger than or equal to 2.5, and it is shown in Fig. 2. 4.2. Dataset construction The seismic data, analyzed in the previous section, must be converted into a 3−𝐷dataset (or cube of data) in order to be interpreted by the TriGen algorithm. Herein, several new attributes must be generated to transform a 2𝐷dataset into a 3𝐷one. Firstly, using a GIS, all data included in the catalog have been sorted into 30 ×20 cells, each representing an area of approximately 16.7 km ×16.4km. Afterward, each cell is characterized as a set of features. These features, originally defined by literature (Scitovski and Scitovski,2013;Reyes and Cárdenas,2010), have been re-adapted to the Pyrenees seismicity and are: 1. 𝑀𝑚𝑎𝑥: maximum earthquake magnitude recorded in the cell. 2. 𝐷: mean depth of all earthquakes’ epicenters recorded in the cell. 3. 𝑀2.9: number of earthquakes occurring with a magnitude larger than or equal to 2.9 in the cell considered. 4. 𝑀3.3: number of earthquakes occurring with a magnitude larger than or equal to 3.3 in the cell considered. 5. 𝑀3.7: number of earthquakes occurring with a magnitude larger than or equal to 3.7 in the cell considered. 6. 𝑀4.1: number of earthquakes occurring with a magnitude larger than or equal to 4.1 in the cell considered. 7. 𝑁: total number of earthquakes occurring. As can be seen in Fig. 3, the result of applying this data transformation is a data cube (or 3𝐷dataset) where every cell is defined by three coordinates 𝑥,𝑦and 𝑓. The 𝑥coordinate is in [1,30] representing the relative latitude, the 𝑦coordinate is in [1,20] representing the relative longitude, and, the 𝑓coordinate is in {𝑀𝑚𝑎𝑥, 𝐷, 𝑀2.9, 𝑀3.3, 𝑀3.7, 𝑀4.1} and represents the particular feature. Therefore, as an example, the cell 𝐶1,2,𝑀3.7represents the number of earthquakes occurring with a magnitude larger or equal to 3.7 in the cell with relative latitude 1 and relative longitude 2. 4.3. The TriGen algorithm The triclustering algorithm used in this research is TriGen (Gutiérrez-Avilés et al.,2014), which stands for Triclustering Genetic based algorithm. As it can be appreciated in Fig. 4(a), the TriGen algorithm receives the 3D input dataset as well as the control parameters to yield a set of tricluster solutions. To do this, TriGen performs an evolutionary heuristic search with a population of individuals that evolves using genetic operators during a specific number of generations with the main goal of optimizing an evaluation function. Considering that an individual is a potential solution tricluster, in this case, each individual represents a zone of the Pyrenees. Focusing on the algorithm, TriGen repeats for each tricluster solution (𝑁control parameter) the following process. Firstly, it creates an initial population that is evaluated based on a fitness function that has to be optimized. Afterward, four genetics operators are applied to the initial population during a number of generations (𝐺control parameter) in order to select the best individual that is going to be the solution tricluster. The four genetics operators that are repeated are selection, crossover, mutation, and evaluation. Each operator has specific characteristics as were deeply analyzed in Gutiérrez-Avilés et al. (2014), Gutiérrez-Avilés and Rubio-Escudero (2015): 5 Fig. 2. Catalog of the work. The size of all the events is given as moment magnitude (𝑀𝑤). Fig. 3. TriGen input dataset. •The initial population is created differently to search for the first solution tricluster and all the other 𝑁-1 solution triclusters searches. The first one is produced with a random subset of all three data coordinates. For the next solution triclusters, some individuals are also randomly generated, and the other ones are created considering not to overlap the previously founded solutions. •The evaluation operator is very useful in promoting the best individuals to continue in the next generation and selecting the final solution tricluster. TriGen offers three different fitness function: the mean square residue measure (𝑀𝑆𝑅) (Gutiérrez-Avilés and Rubio-Escudero,2014), the least square lines (𝐿𝑆𝐿) (Gutierrez- Aviles and Rubio-Escudero,2014), and the multi slope measure (𝑀𝑆𝐿) (Gutiérrez-Avilés and Rubio-Escudero,2015). For this particular case, the 𝑀𝑆𝐿 fitness function is selected. In general, the 𝑀𝑆𝐿 function is based on the resemblances of the slope angles formed by the expression values of the tricluster coordinates. •The selection operator is the first operator computed in the loop of 𝐺generations. The tournament selection is implemented, and so the selected individual is promoted to the next generation. 6 Fig. 4. TriGen algorithm overview. •The crossover operator applied is the one point crossing technique. In particular, two selected individuals (by a probability of crossover parameter 𝑝𝑐) mix their coordinates randomly and generate two new individuals (their children). •The mutation operator also considers a probability of mutation parameter 𝑝𝑚that if it is satisfied, applies one of the six possible actions: add or remove one of the three coordinates of data. The goal of the mutation operator is to guarantee variability for the next generations. The algorithm ends when all the best solution triclusters are found. The produced tricluster solution set represents a partition of the 3𝐷 input dataset, whereas each tricluster is a subset of 𝑥coordinates, 𝑦 coordinates, and 𝑓coordinates. Therefore, as it is represented in the example in Fig. 4(b), a tricluster 𝑇 𝑟𝑁of the solution set is a dataset subset composed by the 𝑥coordinates from 7to 11, the 𝑦coordinates from 27 to 34 and, the {2,3,4,5} features, that is, the 𝑀2.9,𝑀3.3,𝑀3.7 and 𝑀4.1features. To yield the zonation discovered by the algorithm, for each tricluster, its 𝑥and 𝑦coordinates will be an area of the complete analyzed grid. 5. Results and discussion In this section, the experimental workflow carried out in this research is explained, as well as an analysis of the obtained results. For this purpose, the experimental setup of the TriGen algorithm, followed to produce the zonation, is described in Section 5.1. Finally, the discovered zone and a critical discussion are addressed in Section 5.2. 5.1. TriGen experimental setup The TriGen algorithm has been executed four times. As there is a component of randomness in each execution, different executions can lead to different solutions. From each tricluster solution set of each TriGen run, the non-overlapping zones are selected to yield the final zonation presented in the next Section 5.2. As can be seen in Table 1, for each run (from 𝑅1to 𝑅4), a particular parameter configuration has been used as the input dataset described in Section 4.2, in order to obtain tricluster solutions in the widest range possible. For all runs, TriGen has been set up to find ten solutions (𝑁) to guide the algorithm to cover the 𝑋− −𝑌space. The number of generations (𝐺) and the number of individuals in the first population (𝐼) has been {20,30}, and {15,20} respectively. This is considered a proper 7 Fig. 5. Zonation proposal (eight zones). The size of all the events is given as moment magnitude (𝑀𝑤). Source zone #1: Atlantic Pyrenees; #2 Arudy–Arette–Bagnères de Bigorre; #3: Central Pyrenees; #4: Andorra; #5: Easternmost Pyrenees; #6: Eastern Catalonia; #7: Central Iberian System; #8: Northeast Iberian System. Table 1 TriGen parameters setup. Parameter 𝑅1𝑅2𝑅3𝑅4 𝑁10 10 10 10 𝐺20 30 20 30 𝐼15 20 20 15 𝑝𝑐0.2 0.2 0.8 0.8 𝑝𝑚0.5 0.5 0.1 0.1 balance between intensification and diversification of the evolutionary process. This fact has been taken into account when the crossover and mutation parameters were set. Firstly, 𝑝𝑐varies in {0.2,0.8} to obtain high intensified and low diversified solutions in the 𝑝𝑐= 0.2runs and low intensified and high diversified solutions in the 𝑝𝑐= 0.8runs. In a similar way, with 𝑝𝑚varying in {0.1,0.5}, to get solutions with high and low variability. 5.2. Zonation obtained and its representation With the solutions obtained by the application TriGen, following the procedure described in Section 5.1, the zones with less than 25 events have been ruled out, according to Bender (1983) and Skordas and Kulhánek (1992). The use of a GIS (QGIS, https://www.qgis.org, last accessed November 2020) allows a rigorous representation, and the data can be handled for further analysis. The resulting zonation, as well as the epicenters distribution, is shown in Fig. 5. 5.2.1. Seismic parameters of the obtained zones Once the zones have been delineated according to solutions provided by TriGen, the next step is to seismically characterize each zone, in which some seismic parameters must be uniform, namely, the 𝑏-value, or the annual rate, as well as the maximum magnitude. Therefore, they must be estimated for the proposed zones. Thus, the Gutenberg–Ritcher law (Gutenberg and Richter,1944, 1954) holds: 𝑙𝑜𝑔10𝑁(𝑀) = 𝑎−𝑏𝑀 (1) where the b-value estimates the relationship between small and large earthquakes, and it is related to the physics of the source. The lower its value is, the more energy can be accumulated. Nis the number of events with a magnitude larger than or equal to a cut-off magnitude M, and arefers to the seismic productivity. The Maximum-Likelihood-Estimate proposed by Bender (1983) and Aki (1965) for binned data (𝛥) has been considered, in terms of 𝑏-value. 𝑏=1 ln (10)(  𝑀−𝑀𝑚𝑖𝑛 −𝛥 2) (2) where  𝑀is the average magnitude; 𝑀𝑚𝑖𝑛 is the magnitude of completeness, which means all the events of magnitude larger than or equal to this value have been recorded from the reference date; In this work, 𝑀𝑚𝑖𝑛 has been set as 3.0, and the corresponding year of completeness is 1978, as stated in Section 4.1, and the 𝛥 value for this work is 0.1. Besides, it is interesting to give uncertainty parameters in geology (Bárdossy and Fodor,2001). In this case, the b-value uncertainty (Sigma b) has been calculated by the approach suggested by Kijko and Smit (2012). Another relevant parameter is the annual rate referred to a unit area (𝑘𝑚2). It means the number of events exceeding a threshold per year and related to 𝑘𝑚2. 𝜆(𝑀𝑚𝑖𝑛) = 𝑁 𝑡(3) where 𝑡is the time-lapse. In this research, from the beginning of 1978 until the end of 2019, i.e., 42 years. When this value is normalized with the surface, by 𝑘𝑚2, and multiplied by 10,000, it is called 𝐴𝑅3, in this paper. 𝐴𝑅3 = 10,000 𝜆(𝑀𝑚𝑖𝑛) 𝑘𝑚2(4) Finally, the maximum magnitude 𝑀𝑚𝑎𝑥 is calculated for each zone. The resulting values can be found in Table 2. 8 Table 2 Obtained zones. Seismic parameters. Zone Name Area (km2)𝑏-value 𝐴𝑅3𝑀𝑚𝑎𝑥 𝐷𝑒𝑝𝑡ℎ𝑚𝑒𝑎𝑛 Sigma b 1 Atlantic Pyrenees 8157 1.01 5.43 5.2 8.4 0.07 2 Arudy– Arette– Bagnères de Bigorre 8157 0.99 14.54 6.7 6.0 0.04 3 Central Pyrenees 5478 1.20 3.30 6.2 5.7 0.14 4 Andorra 3287 1.32 3.19 5.0 5.7 0.20 5 Easternmost Pyrenees 4388 1.17 3.91 6.5 5.8 0.14 6 Eastern Catalonia 6931 1.34 1.89 5.8 6.3 0.18 7 Central Iberian System 6709 1.13 0.50 5.2 8.6 0.30 8 Northeast Iberian System 6630 1.40 0.61 5.7 8.0 0.34 Fig. 6. AR3 vs b-value. 5.3. Discussion After calculating some of the most relevant seismic parameters, the analysis of the results is driven. First, in order to check that the zones are seismically different from one another, a graphical representation of the AR3 and b-value is presented in Fig. 6. In addition, to support this discrimination, different color maps have been deployed. Herein, Fig. 7 shows the b-value variation; 𝑀𝑚𝑎𝑥; in Fig. 8 the 𝐴𝑅3is depicted; Fig. 9 represents the 𝑀𝑚𝑎𝑥; and finally, the mean depth is illustrated in Fig. 10. Second, the analysis of each zone is undertaken. Zone 1 is located in the Atlantic Pyrenees and presents one of the lowest b-values and the highest annual rate. However, the maximum magnitude is one of the smallest of the proposed zones (5.2), so the seismic hazard can be considered low to moderate in the Pyrenees region. The most seismically hazardous zone is undoubtedly the number 2, which runs from the Arudy–Arette region to the Bagnères de Bigorre environment. This is due to the combination of the minimum b-value, and both maximum annual rate (the highest in the Iberian Peninsula) and 𝑀𝑚𝑎𝑥. Therefore, this zone deserves special attention when a PSHA is driven in the Pyrenees. Zone 3 is situated in the Central Pyrenees, in the Axial Zone. This area, which includes the Maladetta Massif, has been shocked by some of the strongest events in the Pyrenees (M6.2). Besides, it presents medium values in both annual rate and b-value. Herein, the seismic hazard can be defined as moderate; it might be even said moderate to high. The smallest zone, numbered 4, is located in the Andorra environment, and not significant earthquakes have been taken place in it. Its b-value is high, and the annual rate is medium. Like Zone 3, the 5 one, in the easternmost Pyrenees, is characterized by a medium b-value and annual rate but has suffered the second most energetic earthquake in the whole belt. Therefore, the seismic hazard can be estimated as moderate to high. In zone 6, fully contained in eastern Catalonia (Spain), although the 𝑀𝑚𝑎𝑥 is not relatively high, four events with a magnitude exceeding 5.5 have occurred. Besides, it shows a high b-value, and the annual rate is low to moderate so that the seismicity can be rated as moderate. The epicenter distribution may indicate that this zone could be extended along the Catalonian Coastal Range to the southwest or consider a new one. The following zones, 7 and 8, could not be closely related to the Pyrenees, but in its environment, in the Iberian System. Besides, both zones are in the limit to be seen as a seismic zone, individually. Both zones have the lowest annual rate by far, and the maximum magnitude is not notably high (5.7 and 5.2). Regarding the b-value, zone 7 (Central) shows an average value, and zone 8 (Northeast) the highest in the region. Given the few events available for b-value calculations, these are not entirely meaningful. The epicenter distribution points out that both zones could correspond to the same seismic zone if a proper rotation were to be carried out. Although it is a debated topic, the highest seismic potential of source zone #2 (Arudy–Arette-Bagenères de Bigorre) can be related to the presence of several major faults. For example, the 25–30 km-long Herrère right lateral fault (Arudy) or the Lourdes reverse faults (60 km in three segments). Besides, the relevant seismicity of source zone #3 is related to the North Maladetta normal fault (30 km long) or Coronas normal fault (10 km) (Lacan and Ortuño,2012). It is also remarkable that in zone 5 (Easternmost Pyrenees), the more destructive earthquake could be caused by the Vallfogona or Ribes–Camprodon thrusts or the Amer normal fault. The latter fault generated the major 1428 Queralbs earthquake (in zone 6) (Perea,2009). Regarding the method limitations, one of them would be that it only produces ‘‘rectangular zones’’ along parallels and meridians arches, as the original TriGen algorithm works that way. Thus, to overcome this issue, the results could be improved using rotated polygons, which could have been useful to extend zone 6, or join zones 7 and 8. Another method limitation is related to the fact that conceptually, an active fault should not be cut by a source zone boundary (Hamdache and Peláez, 2019), so it could be considered in further research. 6. Conclusions In this research, a seismic zonation has been proposed for the Pyrenees mountain range. The zoning is purely based on the application of the TriGen algorithm to an earthquake catalog and a grid division. Herein, a GIS has been used to manage and represent both the input and the output data. The method aims to be as objective as possible, seeking patterns for the best considered seismic indicators and minimizing human interaction. The triclustering methodology has been the primary procedure to yield the zonation and applied to a pre-processed 3𝐷dataset using the TriGen algorithm. It is based on the genetic algorithm model. Thus, from an initial population of potential solutions, by applying genetic operators (crossovers, mutations), it evolves to optimal solutions based on the 𝑀𝑆𝐿 fitness function. The starting point is an updated, reliable, homogeneous, extensive, and with independent events earthquake catalog, coupled with a geographical grid defined by the GIS. This grid cell size has been set out as a nearly square shape in 16.7 km ×16.4km, a trade-off between the coverage and the low to moderate Pyrenees seismicity. Then, as a previous step to run the evolutionary algorithm, the suggested best indicators have been adapted to the Pyrenean seismicity. Finally, the algorithm 9 Fig. 7. Color map showing the b-value for the obtained source zones. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Fig. 8. Color map showing the 𝐴𝑅3for the obtained source zones. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) has generated the eight seismogenic zones, and their geometries have been integrated into a GIS. The following step is characterizing the zones seismically. To this end, the most relevant seismic parameters have been calculated for each source zone by a GIS to further analysis. This analysis reveals that the source zone covered by the Arudy–Arette–Bagnères de Bigorre region is the most relevant seismically by far, due to its minimum bvalue, the maximum annual rate, and the highest 𝑀𝑚𝑎𝑥. Therefore, this zone deserves special attention when a PSHA is driven in the Pyrenees. The seismic hazard is also remarkable in both the Axial Zone, that includes the Maladetta Massif, and in the easternmost Pyrenees, along the Spain–France border, including M6+ events over time. In the Catalonian Coastal Range, the zone rated as medium seismicity could be extended through this range to the southwest. In the Atlantic Pyrenees and the Andorra environment, the seismic activity is less pronounced. Finally, in the southwest of the studied region (even