Machine Learning in Discriminating Active Volcanoes of the Hellenic Volcanic Arc

Abstract

Identifying the provenance of volcanic rocks can be essential for improving geological maps in the field of geology and providing a tool for the geochemical fingerprinting of ancient artifacts like millstones and anchors in the field of geoarchaeology. This study examines a new approach to this problem by using machine learning algorithms (MLAs). In order to discriminate the four active volcanic regions of the Hellenic Volcanic Arc (HVA) in Southern Greece, MLAs were trained with geochemical data of major elements, acquired from the GEOROC database, of the volcanic rocks of the Hellenic Volcanic Arc (HVA). Ten MLAs were trained with six variations of the same dataset of volcanic rock samples originating from the HVA. The experiments revealed that the Extreme Gradient Boost model achieved the best performance, reaching 93.07% accuracy. The model developed in the framework of this research was used to implement a cloud-based application which is publicly accessible at This application can be used to predict the provenance of a volcanic rock sample, within the area of the HVA, based on its geochemical composition, easily obtained by using the X-ray fluorescence (XRF) technique.

1. Introduction

The Hellenic Volcanic Arc (HVA) in Southern Greece consists of four volcanic regions (Figure 1), namely the Methana region, which includes Methana, Aigina, and Poros in the west; the Melos region, consisting of Melos, Antimelos, and Kimolos in the southwest; the Santorini region, consisting of Santorini, Kolumbus, and Christianoi in the southeast; and the Nisyros region consisting of Nisyros, Giali, and Pergousa in the east. All active volcanoes of the HVA were formed by a converging tectonic setting, producing Island Arc Basalts (IABs) belonging to the calc-alkaline and high-K calc-alkaline series. The first activity took place during the Pliocene (5 m.y.) and continues to the present day [1,2,3]. The extent of the HVA is well-mapped on the surface [1], but little is known about its extent under the seafloor. A large part of the HVA is buried under the sediments of the Southern Aegean Sea. Future deep-sea drilling for mineral prospecting will provide rock samples to produce more precise geological borders for the HVA. In order for these rock samples to be used for mapping, they have to be correlated to one of the four active volcanic regions. Furthermore, volcanic tocks from the HVA have been used to produce artifacts like millstones and anchors from prehistoric times. These artifacts can be found scattered all around the Mediterranean. One of the questions about these artifacts in archaeology is the source of the original material. By answering that question, ancient trade routes can be revealed.

One obvious approach used to correlate volcanic rocks and fingerprint ancient artifacts is to use the geochemistry of these rocks to identify their provenance. By studying the binary K₂O versus SiO₂ diagram [4] of the HVA compiled from 1500 data samples (Figure 2) [2], major overlaps of the data points can be seen. It is obvious that no discrimination between the volcanic regions of the HVA can be established based on this diagram.

The diagram in Figure 2 is based on binary and ternary diagrams, known as discrimination diagrams, which were introduced in the 1970s to discriminate geotectonic settings based on the trace element [5] and major element [6] geochemical data of basalts and sandstones, respectively. The main geotectonic settings are the boundaries of either converging or diverging lithospheric plates where volcanic activity is observed, resulting in the formation of basalts, the most common rock formed by magma pouring out from volcanoes originating from the Earth’s mantle. Basalts can also be formed at hot spots far away from tectonic boundaries, which are known as Ocean Island Basalts (OIB). Given the small amount of data available and the limited areas under examination at the time, discrimination diagrams performed well and this is why they are still used and are continuously being evolved [7,8,9]. Since 1970, however, the amount of rock data acquired from all over the world has rapidly increased. Huge datasets containing information on the chemical analysis of rock samples have become openly available from databases such as GEOROCK [10] and PetDB [11]. These datasets have been used to prove that conventional discrimination diagrams are inaccurate in many cases [12,13] because of the overlap of the sample points, as is the case in the diagram in Figure 2.

Machine learning (ML) has been used in many fields of science to study the hidden patterns in large datasets [14]. After 2000, attempts have been made to use this technology to discriminate tectonic settings. The first successful attempt to discriminate tectonic settings using ML was made using Decision Trees (DT) in 2006 [15]. Ten years later, a Support Vector Machine (SVM) was tested on a dataset using the geochemical and isotopic data of basalts to discriminate tectonic settings [16]. One year later, SVM, Random Forest (RF), and Sparse Multinomial Regression (SMR) were compared using 20 elements and five isotopic ratios of basaltic rocks. The accuracy achieved was over 83% [17]. In 2019, chemical data derived from olivine, a major mineral found in basalts, were used to compare the accuracy of Logistic Regression (LR), Naïve Bayes (NB), RF, and Multi-Layer Perceptron (MLP) for discriminating tectonic settings [18]. The same group introduced a new method named SONFIS incorporating fuzzy logic [19]. The latest published work used the Logistic Regression Analysis (LRA) algorithm with the geochemical data of monazite, a mineral found in basalt rock [20]. In 2017, SVM was used to discriminate volcanic centers in Italy [21].

In this research, the aforementioned methodology used for discriminating tectonic settings was applied in order to evaluate its suitability for distinguishing volcanic regions belonging to one tectonic setting. This was accomplished by comparing the performance of 10 MLAs using six variations of the initial dataset retrieved from the GEOROCK databases. Only major elements were used, as usable trace element data were only available for 10% of the samples of the dataset. Finally a cloud-based application (CBA) was implemented using the best performing model with the ability to correlate a volcanic rock sample from the HVA with one of the four volcanic regions of the South Aegean, based on the geochemical composition of 10 major elements (SiO_2, TiO2, Al₂O_3, Fe₂O₃T, CaO, MgO, MnO, K₂O, Na₂O, and P₂O₅). This application is publicly accessible at www.geology.gr (accessed on 6 September 2021).

2. Materials and Methods

2.1. Machine Learning Algorithms

The goal of ML is to construct a model for making predictions [22]. There are two essential parts to constructing a model: (1) the Machine Learning Algorithm (MLA) [23] and (2) the data. A wide variety of MLAs have been developed since the 1960s, and their improvement has never ceased. The choice of the best-performing algorithm is determined heuristically for a given problem. The best learning algorithm is considered to be the one that has low bias and accurately relates the input to the output values, achieving low variance, meaning that it can successfully predict the same result using different datasets. The latter is also called generalization. Ten MLAs, both linear and non-linear, were used in this study (Figure 3).

Of the linear algorithms, Logistic Regression (LR) [24] and Linear Discriminant Analysis (LDA) [25] were used. The non-linear algorithms tested were k-Nearest Neighbor (kNN) [26], Gaussian Naive Bayes (GNB) [27], and Support Vector Machine (SVM) [28]. Furthermore, Multi-Layer Perceptron (MLP) [29,30] and Extreme Learning Machine (ELM) [31,32] belonging to the Neural Network Group were evaluated. Finally, decision trees were tested, namely Classification and Regression Trees (CART) [33], Random Forest (RF) [34,35,36], and Extreme Gradient Boost (XGB) [37,38].

2.1.1. Logistic Regression

The LR algorithm [24] is an extension of the linear regression algorithm [39], where the goal is to find the ‘line of best fit’ in order to map the relationship between the dependent variable and the independent variables. The LR algorithm predicts whether something belongs to a class or not, i.e., if it is true or false. LR is mainly used for binary classification but can be extended for multiclass classification. In our case, Multinomial Logistic Regression was used.

2.1.2. Linear Discriminant Analysis

LDA is a supervised algorithm which reduces dimensionality. LDA is a two-class classifier but can be generalized for multi-class problems. The linear discriminants are the axes, which maximize the distances of multiple classes [25].

2.1.3. k-Nearest Neighbor

In the kNN classifier [26], the distances of the sample points are calculated from a random central point, and each data point is grouped according to this distance. It is regarded as one of the simplest MLAs. The k parameter indicates the number of points that will be taken into account. Larger values create high bias and should be avoided. Another parameter to be selected is the way the distance is computed. Some options are: (1) Euclidean distance, (2) Manhattan distance, and (3) Minkowski distance.

2.1.4. Gaussian Naïve Bayes

The GNB is a probabilistic non-parametric classifier and, as the name implies, is based on Bayes’ theorem [27].

2.1.5. Support Vector Machine

SVM is a supervised classification technique and was initially developed for binary classification problems. The main idea is that input vectors are mapped in a non-linear way to a very high-dimension feature space. A linear decision surface is constructed in this feature space. When a dataset consisting of three features is trained by an SVM, each point is plotted in the three-dimensional space. The algorithm’s goal is to classify the points by applying three planes. In the case of n features, the points are projected in an n-dimensional space, which cannot be visualized. The algorithm’s goal is to classify the objects by dividing the n-dimensional space with the use of what is known as hyperplanes. Predictions are achieved by projecting new points in these divisions [28].

2.1.6. Multi-Layer Perceptron

An MLP is a type of feedforward Artificial Neural Network (ANN). ANNs are based on the assumption that the way the human brain operates can be mimicked artificially by creating artificial neurons named perceptrons [30]. This is achieved by providing input and output data to the algorithm, which tries, by changing the weights of its neurons, to create a function that successfully maps the relationship between the input and the output data. In the case of the MLP, at least three layers are required: the input layer, the hidden layer, and the output layer. Each node is a neuron that uses a non-linear activation function. In this way, the MLP can classify separable non-linear data [29].

2.1.7. Extreme Learning Machine

ELMs are essentially neural networks with only one hidden layer. According to the original authors [32], this algorithm performs faster than MLPs with good generalization. In an ELM, the hidden nodes’ learning parameters, such as the input weight and biases, are randomly assigned and, according to the author, no tuning is needed. Furthermore, the output weights are determined by a simple inverse operation [31]. The only parameter tuned is the number of neurons in the hidden layer.

2.1.8. Classification and Regression Trees

When CART is used, classification of the samples is achieved by the binary partitioning of a feature [33]. It is not only the ease and robustness of construction that make CART a very attractive classifier. The big advantage is the possibility of interpreting the classification, which is either difficult or not possible at all with the aforementioned classifiers. Another advantage is the ability of classification trees to handle missing data. The disadvantage, however, is the overfitting problem. The tree created works well with the data used to create it but does not perform that well with new data. This is also called low bias–high variance. No parameters are set in this algorithm.

2.1.9. Random Forest

To deal with the overfitting problem of CART, the ensemble classifier RF was developed [35,36]. The main idea is to use many trees instead of one. The trees are built by using randomly selected subsets of the original data. The only constraint is that all the vectors must have the same dimensions. This is called bootstrapping. The results obtained from each tree are then put to a vote. The aggregate is the final result. This technique is known as bagging (from bootstrap aggregation) [34]. The split at each node is made by randomly selecting a subset of features. As the number of trees grows, better generalization is achieved. The parameters that can be set are: (1) the number of trees, (2) the number of features used to make the split at each node, and (3) the function used to measure the effectiveness of the split.

2.1.10. Extreme Gradient Boost

XGB is the latest classifier in the evolution of tree-based classifiers. It introduced tree gradient boosting [38]. Furthermore, XGB can handle sparse data [37]. According to the authors, XGB runs 10 times faster than other solutions.

2.2. Data Preparation

A collection of volcanic rock geochemical data, covering the active volcanic regions (4) of the HVA and most non-active pre-Pliocene volcanic regions (32) of Greece, was retrieved from the GEOROC database [10]. The initial dataset provided consists of 6229 entries containing 181 features. The first 41 features are descriptive and cover general information about the rock samples. Most notable is the year of publication of the rock sample data and the sampling location. Initially, data reduction was performed using the Python libraries Pandas [40] and Numpy [41] in nine sequential steps (Figure 4). In the first step, all rock sample entries published before 1985 [12] were deleted because analytical methods—the way the chemical composition was analyzed—prior to this year are considered unreliable. This resulted in a new data matrix with the dimensions of 5506 × 181. In the second step, 39 of the first 41 features were deleted. The ‘LOCATION’ feature was preserved, as it was essential for creating the classes in the eighth step. The ‘CITATION’ feature was used to select the validation data in the ninth step. In this way, the data matrix was reduced to 5506 × 142. The first feature represents the location and the second represents the scientific paper where the sample’s chemical composition was published, while the other 140 features contain the geochemical data of the rock samples. The chemical composition is expressed was the mass fraction (wt%) of the major elements, trace elements in parts per million (ppm) or parts per billion (ppb), and isotopes as ratios. In the third step, locations outside the research area were deleted, resulting in a 4153 × 142 data matrix.

Moreover, in the initial dataset, many NAN fields were present, which could be computed by the MLAs; therefore, in the fourth step, these were replaced with zeroes and, consequently, no data matrix dimension change occurred.

The chemical compositions provided by GEOROC were determined using the XRF method. Depending on the methodology used, the iron oxides’ mass fraction (wt%) was provided either in FeO, FeOT, Fe₂O₃, or Fe₂O₃T. Thus, in the fifth step, to calculate the mass fraction uniformly, these four features were unified into one (Fe₂O₃T) through the following procedure: all oxides measured as FeOT were transformed to Fe₂O₃T by first multiplying the wt% of FeOT’s chemical composition by 1.1111. The chemical composition of FeO was multiplied by 1.1111 to obtain its chemical composition in Fe₂O₃, then, this value was added to the initial value of Fe₂O_3, resulting in the Fe₂O₃T. In this way, the table’s new dimensions were 4153 × 139.

In the sixth step, features where zeroes represented more than 65.29% (this value corresponds to the oxide P₂O₅) were deleted. This was necessary, as features beyond this threshold had zero values higher than 90%, making them useless for the training process of a ML model. This produced a 4153 × 12 data matrix. In the seventh stage, 199 entries had to be deleted as all 10 of their geochemical data fields were equal to zero. The emerging new data matrix’s dimensions were 3954 × 12. In the eighth step, the classes were prepared based on the four volcanic groups of the HVA [1,2]. Methana along with Aigina and Paros were labeled as Class 1. Melos together with Antimelos and Kimolos were labeled as Class 2. Santorini with Christianoi and Kolumbo were labeled as Class 3, and finally, Nisyros with Giali, and Pergoussa were labeled as Class 4.

In the final ninth step, 40 entries were randomly handpicked from the dataset to compile the validation dataset. This was done with the following constrains: 10 samples from each class, with each sample from a different publication and year. Furthermore, the ‘CITATION’ feature was removed, as it was no longer needed. This resulted in a 3914 × 11 matrix dimension for the final dataset.

Following data reduction, six datasets (All datasets are available upon request) were created. The initial dataset was labeled Dataset 1 and consisted of 3914 entries and 11 features. The first 10 features were the geochemical data of the major elements and the last was the class. A snapshot of Dataset 1, the initial dataset, is shown in

Table 1. Snapshot of Dataset 1.

Feat.#1	Feat.#2	Feat.#3	Feat.#4	Feat.#5	Feat.#6	Feat.#7	Feat.#8	Feat.#9	Feat.#10	Class
SiO2	TiO2	Al2O3	Fe2O3T	CaO	MgO	MnO	K2O	Na2O	P2O5
67.72	0.32	15.7	2.75	3.63	1.23	0.07	3.36	3.16	0.11	1
67.5	0.36	15.15	3.30	4.85	1.69	0.08	3.69	2.55	0.17	1
60.55	0.52	17.26	5.61	6.26	2.49	0.11	2.25	2.98	0.14	1
71.58	0.25	12.72	2.07	1.49	0.5	0.06	4.19	3.01	0.03	2
69.13	0.38	14.06	2.79	2.23	0.67	0.09	3.93	3.47	0.07	2
71.41	0.38	14.05	2.56	2.73	0.77	0.04	3.32	4.02	0.08	2
67.89	0.39	14.11	2.84	2.43	0.64	0.11	3.70	3.73	0.08	2
54.28	0.84	16.97	8.11	8.52	5.56	0.15	1.57	3.07	0.12	3
52.21	0.84	17.79	8.70	10.81	6.17	0.16	0.60	2.71	0.10	3
52.3	0.85	17.67	8.52	10.73	6.13	0.16	0.60	2.78	0.11	3
55.31	0.92	19.19	7.18	8.55	4.11	0.12	1.09	3.61	0.20	4
55.52	0.92	18.99	7.15	8.51	4.30	0.12	1.09	3.42	0.22	4
54.91	0.92	19.51	7.17	8.61	4.34	0.12	1.08	3.61	0.19	4

. These 3914 data samples were used to plot SiO₂ vs. K₂O [42] (Figure 5) and SiO₂ vs. Na₂O + K₂O [43] (Figure 6) to check if the four volcanic regions could be discriminated by the use of this bigger dataset. The diagrams show that this was not the case, as major overlapping of the sample points was present. In Dataset 2, the zeroes were replaced by the mean value. In Dataset 3, the zeroes were replaced by the median. Datasets 4–6 were created from Datasets 1–3 correspondingly by deleting the outliers. The outliers were removed by using the stats class of the Scipy library and the Numpy library. The number of entries of each class in every dataset is shown in

Table 2. Classes and entries of Datasets 1–6.

	Methana (Class 1)	Melos (Class 2)	Santorini (Class 3)	Nysiros (Class 4)	Total
Dataset 1	274	115	1022	2503	3914
Dataset 2	274	115	1022	2503	3914
Dataset 3	274	115	1022	2503	3914
Dataset 4	135	104	908	2275	3422
Dataset 5	132	104	906	2263	3405
Dataset 6	130	104	905	2263	3402

. The basic statistic data of Datasets 1–6 are presented in

Table 3. Basic statistics of Dataset 1.

	Feat.#1	Feat.#2	Feat.#3	Feat.#4	Feat.#5	Feat.#6	Feat.#7	Feat.#8	Feat.#9	Feat.#10
	SiO2	TiO2	Al2O3	Fe2O3T	CaO	MgO	MnO	K2O	Na2O	P2O5
mean	50.62	1.77	15.18	11.33	7.01	8.18	0.18	1.05	2.96	0.23
median	53.43	0.33	15.52	6.44	6.09	1.33	0.12	0.43	3.12	0
min	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
max	81.58	50.47	68.35	100.42	53.85	52.47	1.81	9.62	9.92	44.94
std	17.47	6.34	10.79	18.41	6.48	13.44	0.26	1.40	2.25	2.70
var	305.37	40.24	116.33	339.1	42.04	180.73	0.07	1.96	5.08	7.27

,

Table 4. Basic statistics of Dataset 2.

	Feat.#1	Feat.#2	Feat.#3	Feat.#4	Feat.#5	Feat.#6	Feat.#7	Feat.#8	Feat.#9	Feat.#10
	SiO2	TiO2	Al2O3	Fe2O3T	CaO	MgO	MnO	K2O	Na2O	P2O5
mean	52.33	2.39	16.71	11.57	7.38	9.52	0.23	1.30	3.49	0.38
median	53.43	1.35	15.52	6.80	7.01	4.53	0.18	1.05	3.12	0.23
min	0.01	0.00	0.00	0.02	0.00	0.00	0.00	0.00	0.00	0.00
max	81.58	50.47	68.35	100.42	53.85	52.47	1.81	9.62	9.92	44.94
std	14.70	6.22	9.53	18.34	6.27	12.96	0.24	1.28	1.79	2.69
var	216.00	38.75	90.8	336.39	39.32	168.01	0.06	1.63	3.22	7.21

,

Table 5. Basic statistics of Dataset 3.

	Feat.#1	Feat.#2	Feat.#3	Feat.#4	Feat.#5	Feat.#6	Feat.#7	Feat.#8	Feat.#9	Feat.#10
	SiO2	TiO2	Al2O3	Fe2O3T	CaO	MgO	MnO	K2O	Na2O	P2O5
mean	52.43	1.89	16.74	11.46	7.33	8.40	0.22	1.15	3.52	0.34
median	53.43	0.33	15.52	6.44	6.10	1.33	0.12	0.43	3.12	0.17
min	0.01	0.00	0.00	0.02	0.00	0.00	0.00	0.00	0.00	0.00
max	81.58	50.47	68.35	100.42	53.85	52.47	1.81	9.62	9.92	44.94
std	14.69	6.31	9.52	18.36	6.28	13.32	0.24	1.33	1.79	2.69
var	215.94	39.85	90.70	336.92	39.39	177.41	0.06	1.78	3.19	7.22

,

Table 6. Basic statistics of Dataset 4.

	Feat.#1	Feat.#2	Feat.#3	Feat.#4	Feat.#5	Feat.#6	Feat.#7	Feat.#8	Feat.#9	Feat.#10
	SiO2	TiO2	Al2O3	Fe2O3T	CaO	MgO	MnO	K2O	Na2O	P2O5
mean	54.58	0.6	16.93	6.83	7.75	6.6	0.13	1.11	3.36	0.06
median	54.72	0.34	16.5	4.98	7.21	1.20	0.10	0.55	3.45	0.00
min	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
max	81.58	7.04	36.91	46.88	25.44	48.48	0.96	5.18	8.25	2.68
std	12.16	0.84	9.87	7.94	5.75	11.24	0.17	1.24	2.1	0.11
var	147.98	0.7	97.4	63.08	33.04	126.31	0.03	1.53	4.43	0.01

,

Table 7. Basic statistics of Dataset 5.

	Feat.#1	Feat.#2	Feat.#3	Feat.#4	Feat.#5	Feat.#6	Feat.#7	Feat.#8	Feat.#9	Feat.#10
	SiO2	TiO2	Al2O3	Fe2O3T	CaO	MgO	MnO	K2O	Na2O	P2O5
mean	55.44	1.23	18.00	6.97	7.95	7.98	0.19	1.28	3.67	0.20
median	54.79	1.09	16.53	5.53	7.27	4.33	0.18	1.05	3.46	0.23
min	30.42	0.00	0.00	0.02	0.08	0.00	0.00	0.00	0.00	0.00
max	81.58	7.04	36.91	46.88	25.44	48.4	0.94	4.81	8.25	2.68
std	10.02	0.82	8.77	7.62	5.61	10.73	0.14	1.14	1.79	0.07
var	100.41	0.66	76.92	58.02	31.48	115.05	0.02	1.31	3.19	0.01

and

Table 8. Basic statistics of Dataset 6.

	Feat.#1	Feat.#2	Feat.#3	Feat.#4	Feat.#5	Feat.#6	Feat.#7	Feat.#8	Feat.#9	Feat.#10
	SiO2	TiO2	Al2O3	Fe2O3T	CaO	MgO	MnO	K2O	Na2O	P2O5
mean	55.49	0.72	18.03	6.84	7.94	6.71	0.17	1.18	3.68	0.06
median	54.8	0.34	16.54	5.52	7.28	1.33	0.12	0.56	3.47	0.00
min	30.42	0.00	0.00	0.02	0.08	0.00	0.00	0.00	0.00	0.00
max	81.58	7.04	36.91	46.88	25.44	48.35	0.94	4.81	8.25	2.68
std	10.00	0.77	8.76	7.58	5.61	10.98	0.15	1.18	1.78	0.11
var	100.05	0.59	76.78	57.51	31.52	120.51	0.02	1.4.0	3.17	0.01

. The histograms of Datasets 1–3 with outliers are depicted in Figure 7. The histograms of Datasets 4–6 with the outliers removed are depicted in Figure 8. The histogram of Dataset 5 shows the best distribution of the data.

2.3. Simulations

The training of the models was implemented using Python 3.8 (Guido van Rossum, The Hague, The Netherlands) [44] joint with Scikit-learn (David Cournapeau, Paris, France) [45], the xGboost wrapper for Scikit-learn(https://xgboost.readthedocs.io/, accessed on 30 June 2021) and elm.ELM Kernel (Augusto Almeida, Brazil). The goal was to create the most accurate model, with the smallest standard deviation and the least computing time, by using one out of the 10 MLAs and one of the six datasets. Training of the model was conducted on a Windows 10 System with an Intel Core 2 Quad Q8400 CPU and 8 GB of RAM. Tenfold cross-validation was applied (i.e., the original sample was randomly partitioned into 10 equal size subsamples). The datasets were split randomly into 70% training data and 30% testing data. In this way, the experiments were repeated 10 times by using 10 different versions of the dataset. Finally, three stages of testing and three improvement procedures were applied (Figure 9).

The initial metrics were computed by using the six datasets on the 10 algorithms with their default parameters.

The threshold values for the first testing stage were an accuracy of >88% and a standard deviation of <1.5%. In the following initial improvement procedure, a grid search was applied to the MLAs of the models passing the first testing stage to optimize the selected model’s algorithms.

The second testing stage had a threshold accuracy of >92% and the standard deviation remained at <1.5%. The models outperforming all the others were then put into the second improvement procedure, where feature reduction was applied using the results of a heatmap.

In the third testing stage, the goal was to improve accuracy by keeping the standard deviation of <1.50%. Following this, a confusion matrix [46] was compiled for the best-performing model. This revealed that a third improvement procedure was necessary to augment the data. Namely, the Synthetic Minority Over-sampling Technique (SMOTE) [47] was applied to the best-performing model’s dataset. Finally, the best-performing model was tested using the testing data, i.e., the 40 entries left out of the training and validation.

3. Results

The results of the first testing stage showed accuracy ranging from 48.39% (2.24% standard deviation, 2.01 s) for the model using GNB with Dataset 2 to 92.75% (1.36% standard deviation, 11.48 s) for the model using XGB with Dataset 5 (

Table 9. Testing results for the first testing stage (LR: Logistic Regression, LDA: Linear Discriminant Analysis, SVM: Support Vector Machine, KNN: K Nearest Neighbor, GNB: Gaussian Naïve Bayes, MLP: Multi-Layer Perceptron, ELM: Extreme Learning Machine, CART: Classification and Regression Trees, RF: Random Forest, XGB: xGBoost, std: standard deviation).

	LR			LDA			SVM			KNN
Dataset	Accuracy	Std	Time	Accuracy	Std	Time	Accuracy	Std	Time	Accuracy	Std	Time
			(Sec)			(Sec)			(Sec)			(Sec)
1	69.82%	2.15%	2.66	66.84%	2.02%	2.34	69.75%	2.23%	5.05	88.02%	1.97%	3.49
2	68.03%	1.85%	2.7	64.74%	2.18%	2.2	70.20%	1.64%	5.32	87.77%	1.82%	3.66
3	68.64%	1.87%	2.68	65.98%	2.32%	2.14	70.17%	1.75%	5.43	88.57%	1.38%	4.14
4	67.67%	2.08%	2.52	65.99	1.47%	2.32	68.94%	0.85%	3.88	88.43%	1.91%	3.6
5	71.00%	1.74%	2.52	65.42%	2.54%	2.32	69.31%	1.12%	5.25	87.96%	1.57%	3.86
6	71.22%	1.60%	2.48	66.12%	1.60%	2.45	71.00%	1.06%	5.42	88.58%	1.08%	3.77
	GNB			MLP			ELM			CART
Dataset	Accuracy	Std	Time	Accuracy	Std	Time	Accuracy	Std	Time	Accuracy	Std	Time
			(Sec)			(Sec)			(Sec)			(Sec)
1	48.48%	2.89%	2.00	84.64%	1.80%	30.00	80.46%	1.30%	3.41	87.37%	2.18%	3.63
2	48.39%	2.24%	2.01	87.61%	1.74%	30.12	85.52%	1.50%	0.95	87.66%	1.48%	3.59
3	48.74%	2.45%	2.17	83.65%	1.86%	30.42	74.27%	3.45%	1.44	87.66%	1.62%	3.61
4	62.07%	1.60%	2.01	86.33%	0.98%	22.59	74.78%	5.22%	0.44	88.02%	1.89%	3.41
5	66.30%	3.32%	2.03	85.72%	2.03%	31.22	77.54%	1.78%	1.35	87.58%	2.30%	3.6
6	66.19%	1.60%	2.21	86.81%	1.77%	32.35	76.13%	1.92%	1.01	87.86%	2.12%	3.59
	RF						XGB
Dataset	Accuracy	Std				Time	Accuracy			Std		Time
						(Sec)						(Sec)
1	91.60%	1.60%				10.3	91.57%			1.34%		13.19
2	91.68%	1.76%				10.18	91.42%			1.60%		13.02
3	91.64%	1.24%				10.39	91.97%			1.55%		13.02
4	92.02%	1.54%				9.64	91.82%			1.99%		11.44
5	92.73%	1.46%				9.26	92.75%			1.36%		11.48
6	92.69%	1.48%				9.05	92.73%			1.68%		11.33

, Figure 10).

Six models met the requirements of the first testing stage (>88% accuracy and <1.5% standard deviation). These were Model 1: kNN combined with Dataset 6 (88.58% accuracy, 1.08% standard deviation, 3.77 s); Model 2: RF with Dataset 3 (91.64% accuracy, 1.24% standard deviation, 10.39 s); Model 3: RF with Dataset 5 (92.73% accuracy, 1.46% standard deviation, 9.64 s); Model 4: RF with Dataset 6 (92.69% accuracy, 1.48% standard deviation, 9.05 s); model 5: XGB combined with Dataset 1 (91.57% accuracy, 1.34% standard deviation, 13.19 s); and Model 6: XGB combined with Dataset 5 (92.75% accuracy, 1.36% standard deviation, 11.48 s). All MLAs were used with their standard parameters (

Table 10. Best performing models (first testing stage).

Model	MLA	Dataset	Accuracy	Std	Time (Sec)	Parameters
1	kNN	6	88.58%	1.08	3.52	n_neighbors = 5, weight = uniform, metric = Euclidean
2	RF	3	91.64%	1.24%	10.39	max depth = 3, n_estimators = 100
3	RF	5	92.73%	1.46%	9.64	max depth = 3, n_estimators = 100
4	RF	6	92.69%	1.48%	9.05	max depth = 3, n_estimators = 100
5	XGB	1	91.57%	1.34%	13.19	max depth = 3, n_estimators = 100, learning_rate = 0.1
6	XGB	5	92.75%	1.36%	11.48	max depth = 3, n_estimators = 100, learning_rate = 0.1

).

The performance of a MLA can be optimized by changing a set of parameters. Each MLA has its own parameter that can be fine-tuned. The six models were improved in the first improvement procedure by using the gridsearch() function [45], which allows bulk testing of models on a selected set of parameters. This procedure is known as hyperparameter tuning. In the second testing stage (accuracy >92%, standard deviation <1.5%), only three models reached these stricter criteria, namely Model 3: RF with Dataset 5 (92.79% accuracy, 1.29% standard deviation, 16.13 s); Model 4: RF with Dataset 6 (92.75% accuracy, 1.35% standard deviation, 15.91 s); and Model 5: XGB combined with Dataset 5 (92.80% accuracy, 1.48% standard deviation, 18.97 s). A side effect of the improvement in accuracy was a rise in computational time during the training (

Table 11. Performance improvement of the models after hyperparameter tuning (first improvement procedure).

Model	MLA	Dataset	Accuracy	Std	Time (Sec)	Parameters
1	kNN	6	91.26%	1.70%	3.21	n_neighbors = 2, weights = distance, metric = Manhattan
2	RF	3	91.85%	0.95%	30.16	max depth = 26, n_estimators = 200
3	RF	5	92.79%	1.29%	16.13	max depth = 17, n_estimators = 200
4	RF	6	92.75%	1.35%	15.91	max depth = 14, n_estimators = 200
5	XGB	1	91.93%	1.56%	23.03	max depth = 9, n_estimators = 106, learning_rate = 0.1
6	XGB	5	92.80%	1.48%	18.97	max depth = 8, n_estimators = 101, learning_rate = 0.1

).

In the second improvement procedure, the best models from the second testing stage were analyzed to accomplish feature reduction. A heatmap was used to rank the features of Datasets 5 and 6 and further improve Models 3, 4, and 6. The features of Dataset 5 were reordered as follows: Al₂O₃, TiO₂, CaO, P₂O₅, Na₂O, MnO, SiO₂, Fe₂O₃T, MgO, and K₂O. The features of Dataset 6 were reordered as follows: Al₂O₃, CaO, Na₂O, TiO_2, MnO, P₂O_5, SiO₂, Fe₂O₃T, K₂O, and MgO. After using all the rearranged features, the three models showed minor improvements (

Table 12. Performance improvement of the models after applying a heatmap (second testing stage, second improvement procedure).

Model	MLA	Dataset	Accuracy	Std	Time (Sec)	Parameters
3	RF	6	92.81%	1.3%	15.28	max depth = 17, n_estimators= 200
4	RF	5	92.90%	1.25%	14.43	max depth = 14, n_estimators= 200
6	XGB	5	92.94%	1.15%	13.85	max depth = 8, n_estimators= 101, learning_rate = 0.1

). Model 6 achieved the best results. Rearranging the features resulted in an improvement in computational time. By removing features based on the aforementioned ranking, no increase in accuracy was achieved.

In the third improvement procedure, a confusion matrix was created for Model 6 to evaluate the accuracy of each class (Figure 11a). The low prediction accuracy results of Melos can be explained if we take into account that an extremely imbalanced dataset was used: 132 entries for Methana, 104 for Melos, 906 for Santorini, and 2263 for Nisyros (Table 2).

For this reason, SMOTE was applied to augment Dataset5. One new dataset was created, labeled Dataset 5—balanced 900. The confusion matrix (Figure 11b) showed an improvement. The final model was devised by using the optimized XGB algorithm and Dataset 5—balanced 900 created by this imputation. Accuracy improved to 93.07% and the standard deviation was calculated as 1.21%.

This final model was tested with the 40 handpicked entries left out of the training and testing data. In this validation sample, 90% accuracy was achieved. The Methana class scored 100% and the Melos class scored 80%. One sample was misclassified as Santorini and one as Nisyros. The Santorini class scored 90%, with only one sample being misclassified as Melos. The Nisyros class scored 90%, one sample being misclassified as Santorini.

4. Discussion

The SiO₂ vs. K₂O diagram plotted in this study using 3914 data points showed that it is not possible to discriminate the tectonic regions by using this methodology. This is in accordance with the SiO₂ vs. K₂O diagrams of the HVA’s volcanic regions using 1500 data points in 2005 [1]. In order to classify the volcanic regions, the simulations showed that the best result was achieved by Model 6, which used the XGB classifier and a dataset where imputation was applied by replacing the zero values with the mean value and where the outliers were removed (Dataset 5). This final model performed at an accuracy of 92.75%. After hyperparameter tuning, feature rearrangement, and data augmentation had been applied, the accuracy reached 93.07%. This shows that this procedure can augment accuracy and should always be applied. It is notable to mention that Model 1, which used the kNN classifier with Dataset 6, performed at an accuracy of 91.26% with a very low computational cost of only 3.21 s. The final validation was conducted using the validation dataset consisting of 40 entries left out of the training and testing which resulted in an accuracy of 90%.

This study is in accordance with the positive results obtained from the Calabrian Volcanic Arc [21]. In this study, a similar dataset was used to discriminate the volcanic groups of Central Italy in the framework of tephrochronological studies using the SVM classifier. In this research, 98% accuracy was achieved.

The use of the model developed in the present study should not be used as a black box. Other means of examining the rock samples, such as mineralogy and the age of crystallization, must be taken into account to cross-validate the model’s results. The misclassification of four samples of the validation dataset shows that further improvement of the methodology described is necessary. This will be possible by retraining the selected classifier with new data which will emerge in the future.

5. Conclusions

The results of this research prove that the ML methodology used for discriminating tectonic environments can also be used to discriminate volcanic regions belonging to the same tectonic setting. In the case of the present study, four volcanic regions of the HVA producing IABs in a converging tectonic environment were successfully discriminated using the XGB classifier with an accuracy of 93.07% during training and 90% during validation.

The final model compiled was used for the implementation of a CBA [48] written in Python code and deployed in the flask framework. The CBA implemented in the framework of this study can be used as a valuable tool in the field of geology to classify volcanic rocks found in the Southern Aegean Sea in one of the four active volcanic regions of the HVA using 10 major elements. This will also be a useful tool for mapping of the HVA in the future. Furthermore, in geoarchaeology, the CBA can be used as a valuable tool for tracking the origin of artifacts found all over Greece, such as millstones and anchors made from volcanic rock dating from prehistory to antiquity, a procedure known as geochemical fingerprinting. The CBA is publicly accessible at www.geology.gr (accessed on 6 September 2021).

The authors’ future plans include applying this methodology for the two active volcanic arcs in the Mediterranean, namely the Hellenic Volcanic Arc in the Aegean Sea and the Calabrian Volcanic Arc in the Tyrrhenian Sea, as similar datasets are available in the GEOROC database. This would produce a tool for applying geochemical fingerprinting of ancient artifacts to a broader geographic region.