M2.4/3.3 Submission of developed models for publication

Noname manuscript No. (will be inserted by the editor) Comparison of Geostatistical and Machine Learning Methods for1 Spatial Analysis of Natural Resources Data2 Maria Konstantina Germanou ·Andreas Pavlides ·3 Emmanouil A. Varouchakis4 5 Received: date / Accepted: date6 Abstract This research work employs advanced geostatistical and machine learning methods to analyze7 and model spatial data, with a focus on zinc concentration measurements. The primary objectives are8 to evaluate and compare these methods in terms of generating accurate spatial predictions, quantifying9 uncertainty, and identifying critical spatial patterns. The suggested approach used Self-Organizing Maps10 (SOM) to augment the applicability of Ordinary Kriging (OK) to a considerably large dataset. The11 classical OK approach and the Gaussian Process Regression (GPR) method were applied for compari-12 son purposes, as they are widely used methods. All three methods satisfy both fields explored herein,13 geostatistics and machine learning.14 In preliminary analysis, a variety of kernels were tested, including the novel Harmonic Covariance15 Estimator (HCE). The exponential kernel was selected for the comparisons among the three methods.16 GPR is a flexible Bayesian approach that is exceptionally efficient in capturing complex spatial patterns17 and providing robust uncertainty estimates. While both GPR and OK aim to achieve accurate spatial18 predictions, SOMs enhance kriging by classifying data and facilitating predictions based on the best-19 matching neuron’s data, effectively adapting kriging to localized spatial patterns and improving the20 interpretability of spatial dependencies.21 This integration of methodologies demonstrates the combination of stochastic geostatistical methods22 with machine learning for an improved understanding and prediction of spatial phenomena. The results23 highlight that this modeling approach enhances the management and evaluation of natural resources,24 specifically in the case of a zinciferous ore deposit. More specifically, the research findings indicated that25 the conventional selection of neighborhoods in OK tends to favor certain validation measures. In contrast,26 the SOM-guided selection of the neighborhood enhances the predictive-observational correlation.27 1 Introduction28 The integration of geostatistical methods with machine learning techniques has emerged as a trans-29 formative approach for analyzing and predicting complex geospatial data in various fields, including30 mining, environmental monitoring, and resource management. These datasets often exhibit intricate31 spatial patterns and dependencies, which are essential for accurate analysis but challenging to model.32 Traditional methods such as Ordinary Kriging (OK) are widely used for spatial prediction, since they33 effectively utilize the spatial correlations in observed data to make unbiased predictions through co-34 variance or variogram modeling. However, while effective, OK and its variants suffer from significant35 limitations regarding non-stationary processes that are common in many complex datasets. These have36 been extensively discussed and analyzed by (Cressie 1993) and (Goovaerts 1997).37 Rasmussen and Williams (2006) introduced Gaussian Process Regression (GPR) as a flexible prob-38 abilistic framework that captures nonlinear spatial patterns and provides robust uncertainty quantifica-39 tion, thereby addressing several challenges faced by traditional methods. Further emphasizing its utility,40 (Hengl et al. 2004) demonstrated the application of GPR in predicting soil properties such as organic41 carbon, texture classes, and nutrient distribution. Similarly, (Liu and Weisberg 2005) discussed the flexi-42 bility of GPR in modeling both groundwater quality and ocean current variability, thereby demonstrating43 its ability to adapt across diverse geospatial contexts.44 School of Mineral Resources Engineering, Technical University of Crete, Chania 73100, Crete, Greece E-mail: apa[email protected] 2 Maria Konstantina Germanou et al. An expansion of spatial analysis was proposed by (Kohonen 2001), who developed Self-Organizing45 Maps (SOM), an unsupervised learning algorithm. SOM transforms high-dimensional data into clustered46 and interpretable spatial structures. (Bação et al. 2005) demonstrated how SOM could serve as a47 substitute for k-Means clustering in the analysis of soil variability. Similarly, (Vesanto and Alhoniemi48 2000) illustrated the use of SOM to cluster high-dimensional environmental data for identifying the49 sources of groundwater contamination. More generally, (Liu and Weisberg 2011) provided an extensive50 review of the applications of SOM in meteorology and oceanography, demonstrating its wide applicability51 in environmental sciences.52 The potential of hybrid approaches that synergistically combine the strengths of GPR, SOM, and53 traditional geostatistical methods, such as OK, has been tried and tested, demonstrating their efficiency54 in enhancing geospatial data analysis. (Hengl et al. 2004) presented regression-kriging frameworks to im-55 prove predictions of soil variables. Further extensions in advanced geostatistics, integrated with machine56 learning techniques, have been made in hydrology and water resource management, where SOM-based57 clustering has been instrumental in assessing pollution patterns, as demonstrated by (Liu and Weisberg58 2005) and (Vesanto and Alhoniemi 2000).59 OK remains the most widely established stochastic geostatistical method, preferred by the geoscien-60 tific community (Cressie 1993; Chiles and Delfiner 2012) since its original formulation (Krige 1951;61 Matheron 1965). GPR can be regarded as a supervised machine learning counterpart to kriging, grounded62 in similar probabilistic principles. Given that SOM fall within the field of unsupervised learning, they63 were contrasted with GPR to compare their efficiency on spatial modeling. (Varouchakis et al. 2023) pro-64 posed a hybrid Geostat-SOM framework, which integrates SOM with OK to classify and then interpolate65 groundwater levels across complex hydrogeological systems.66 The work examines the combined application of OK, GPR, and SOM to develop an integrated frame-67 work for spatial data analysis. From the various Kernels investigated in preliminary analysis (including68 the Harmonic Covariance Estimator), the Exponential Kernel performed better and is used in this work.69 The current research introduces an important improvement in the predictiveness and interpretation of70 spatial patterns by involving the complementarities of these methodologies. It examines OK and GPR71 for efficient spatial interpolation and uncertainty quantification, while SOM is used to guide kriging72 estimations by enhancing data classification. This integrated approach overcomes the intrinsic shortcom-73 ings of conventional geostatistical methods and provides a flexible methodology for geospatial analysis74 supporting sustainable resource management. The framework proposed in this research can be extended75 to different datasets so that follow-up research can augment geostatistical methodologies with machine76 learning models.77 2 Methodology78 2.1 Ordinary Kriging79 Kriging refers to a group of stochastic interpolation methods that are by construction linear, unbiased,80 and minimum variance estimators (Goovaerts 1997; Varouchakis et al. 2018).81 The value of the random field X(u)at an unmeasured location u∈Gis estimated by a linear82 combination of the measurements at n(u)nearby points s1(u), . . . , sn(u), where si(u)∈ {s1, . . . , sN}is a83 neighbor of ufor all i= 1, . . . , n(u). A map of the spatial distribution of the field is obtained by repeating84 the estimation process at every node of a prediction grid. Such maps can be accompanied by estimates of85 the prediction variance at each point. The variance can adequately represent the prediction uncertainty86 if the data probability distribution is Gaussian. If the data follow a skewed probability distribution, the87 Kriging-based uncertainty estimates are not reliable (Pavlides et al. 2022).88 To reflect the spatial dependence of uwith each point in n(u), the variogram is used. The variogram89 function γX(r)for stationary random fields is connected to the covariance function cX(r)as shown in90 Eq. (1).91 γX(r) = cX(0) −cX(r)(1) where ris the distance between two points. Thus, the variogram function can be used interchangeably92 with the covariance function if the random field is considered a stationary field.93 The variogram function is defined as follows.94 γX(r) = 1 2Var [X(s)−X(s+r)] (2) Comparison of Geostatistical and Machine Learning Methods for Spatial Analysis of Natural Resources Data 3 where Var is the variance operator, i.e., Var(X) = E[X2(s)] −E2[X(s)] (Hristopoulos 2020).95 OK assumes that the random field is stationary and its mean, mX, is unknown, but constant within96 the kriging neighborhood. The OK estimator is obtained by means of the following equation (Krige 1951;97 Cressie 1990):98 ˆ X(u) = n(u) X i=1 λi(u) (Xi)(3) In Eq. (3), λi(u)are the Kriging weights at the target point u. The Kriging weights are calculated99 by solving the following linear system.100 C XΛ=Cu(4) where C Xis the covariance (or variogram) matrix of the fluctuation random field X0(·)at the data101 locations, i.e., [C X]i,j =c(si−sj)for i, j = 1, . . . , n(u),Λ= (λ1(u), . . . , λn(u))>is the vector of the102 Kriging weights, and Cu=c(u−s1), . . . , c(u−sn(u))>is the covariance (or variogram) vector between103 the unknown point uand each of the n(u)neighbor points that contribute to the estimate (Chiles and104 Delfiner 2012). The Kriging estimate ˆ X(u)at uis obtained by replacing in Eq. (3) the Xiwith the105 respective data values, mXwith the estimate of the mean, and the Kriging weights with the solution of106 Eq. (4) for Λ.107 2.2 Gaussian Process Regression108 GPR is a nonparametric, flexible regression that captures complex patterns in spatial data by modeling109 correlations between data points using kernel functions. AGaussian Process (GP) is a collection of110 random variables, any finite subset of which follows a joint multivariate Gaussian distribution. Formally,111 a GP {f(x)}x∈X , indexed by an input space X, is characterized by its mean function m(x)and covariance112 function c(x,x0), and can be written as:113 f(x)∼ GP(m(x), c(x,x0)) ,(5) where:114 –m(x) = E[f(x)] is the mean function, which gives the expected value of the process at input x,115 –c(x,x0) = E[(f(x)−m(x))(f(x0)−m(x0))] is the covariance function, which defines the relationship116 (similarity) between function values at different input points xand x0.117 These two components fully specify a Gaussian Process—the mean function and the covariance func-118 tion. GPR provides a flexible, non-parametric framework for modeling functions, as they do not assume119 a specific form of the function but instead defines a distribution over all possible functions that could120 explain the observed data. GPR assumes that the data points are drawn from a multivariate normal dis-121 tribution and models the spatial structure of the data through the covariance (kernel) function mentioned122 above.123 Due to its analytical tractability, the GP is commonly used to model real-valued functions, particularly124 in cases where uncertainty needs to be quantified. In practice, one can think of a GP as a distribution125 over functions, where any collection of function values has a multivariate Gaussian distribution.126 2.2.1 The Covariance Function (Kernel)127 The covariance function, also known as the kernel function, is a fundamental component of GPR. It128 defines the covariance between function values at different input points, controlling key properties of129 the resulting function such as smoothness, periodicity, and behavior over various scales. Formally, the130 covariance function c(x,x0)takes two input points xand x0and returns a scalar that represents the131 covariance between the corresponding function values, f(x)and f(x0).132 A covariance (kernel) function c(x,x0)is a non-negative, real-valued, integrable function that repre-133 sents a mapping Rd×Rd→R. Let xand x0be two points in the RdEuclidean space, and let u=kx−x0k134 represent the Euclidean distance between them (Hristopoulos 2020). A covariance function satisfies the135 following properties:136 1. Non-negativity: c(u)≥0for all u.137 4 Maria Konstantina Germanou et al. 2. Symmetry: c(u) = c(−u)for all u.138 3. Normalization: R∞ −∞ du c(u) = 1.139 4. Finite second-order moment: The integral R∞ −∞ du u2c(u)exists.140 5. Mode at the origin: c(u)takes its maximum value at u= 0.141 6. Continuity: c(u)is a continuous function of u.142 7. Scaling: If c(u)is a covariance function and h > 0is a non-negative number, then the function143 ch(u) = 1 hcu his also a covariance function.144 Models of Covariance Functions:145 Commonly used models of covariance functions are listed below. These functions depend on the normal-146 ized distance u=kx−x0k h, where h > 0represents the kernel bandwidth.147 For radial kernels, the covariance functions depend only on the magnitude u=kx−x0k, not the148 direction of the Euclidean distance between two points.149 1. Uniform:150 c(u) = 1 2I|u|≤1(u).(6) 2. Spherical:151 c(u) = 4 31−1.5u+ 0.5u3I|u|≤1(u).(7) 3. Gaussian:152 c(u) = 1 √2πexp(−u2).(8) 4. Exponential:153 c(u) = 1 2exp(−|u|).(9) In the above, I|u|≤1(u)denotes the indicator function, which satisfies:154 I|u|≤1(u) = (1if |u| ≤ 1, 0if |u|>1.(10) The first two covariance functions are compactly supported, meaning they are zero outside the interval155 |u| ≤ 1, while the last two have unbounded support.156 2.2.2 Harmonic Covariance Estimator (HCE)157 A novel covariance function, the Harmonic Covariance Estimator (HCE) has been developed and detailed158 in (Varouchakis et al. 2025). HCE is based on the Spatiotemporal covariance functions from spatial159 physics and specifically from the Linearly Damped Harmonic Oscillator (LDHO) (Hristopulos 2024).160 The equation of HCE is shown in eq. (11) below:161 c(u) = η2 0Γd+1 2β πd+1 2(u2+β2)d+1 2 ,(11) where dis the number of dimensions, η2 0is a scaling factor, and βis a parameter that controls the rate162 of the decay of the kernel with distance. HCE has been used in various natural resources applications,163 utilizing Euclidean and non-Euclidean distances.164 2.2.3 Maximum Likelihood Estimation (MLE)165 Maximum Likelihood Estimation (MLE) is a commonly used method for selecting the parameters θ∈Θ166 of a parameterized covariance kernel cθ(x,x0)(Rasmussen and Williams 2006). Under the Gaussian167 Process model GP(m, cθ), the probability density function of the observed data ygiven θis expressed168 as:169 p(y|θ) = 1 pdet(2πC)exp −1 2y>C−1y,(12) Where:170 –y: The observed target values,171 Comparison of Geostatistical and Machine Learning Methods for Spatial Analysis of Natural Resources Data 5 –C: The covariance matrix computed using the kernel cθ(x,x0)for the observed inputs x,172 –θ: Parameters of the kernel (e.g., length scale, variance, nugget).173 Maximizing this likelihood function is equivalent to minimizing the negative log-likelihood function.174 The log-likelihood function is given by:175 log p(y|θ) = −1 2y>C−1y−1 2log det C−n 2log(2π),(13) Where nis the number of observed data points.176 For computational convenience, we often minimize the modified log-likelihood function, which177 removes constant terms:178 l(θ|y) = y>C−1y+ log det C.(14) 2.2.4 Optimization of Kernel Parameters179 The parameters θare selected by solving:180 θML ∈arg min θ∈Θl(θ|y),(15) Where:181 –The first term y>C−1ymeasures the goodness-of-fit of the model to the data,182 –The second term log det Cpenalizes model complexity to prevent overfitting.183 By minimizing this objective function, MLE balances the trade-off between model fit and complexity,184 ensuring that the selected kernel parameters θgeneralize well to unseen data. This process is critical185 for accurately modeling the covariance structure c(x,x0)between observed data xand unobserved test186 points x0.187 2.2.5 Constructing the Gaussian Process188 To model a set of data T={(xi, yi)}n i=1, we can employ a GP. This provides a flexible, nonlinear189 alternative to simple linear regression models, such as Y=β0+β1X+, which may not fully capture190 the complexity of the data (Rasmussen and Williams 2006). Instead of assuming a linear relationship,191 we hypothesize a nonlinear function fthat maps inputs to outputs via the model:192 yi=f(xi) + i,(16) where irepresents noise or observation error, typically assumed to be independent and identically193 distributed (i.i.d.) according to:194 i∼ N(0, σ2 n),(17) with σ2 nbeing the noise variance.195 In GPs, the function values f(x)and f(x0)for any two distinct known inputs xand x0are treated196 as random variables that follow a joint Gaussian distribution. This assumption allows us to model the197 correlation between function values at different input locations.198 The joint distribution of the function values at known inputs xand x0can be written as:199 f(x) f(x0)∼ N m(x) m(x0),c(x,x)c(x,x0) c(x0,x)c(x0,x0),(18) where m(x)and m(x0)are the mean function values at the known inputs xand x0, and c(x,x0)is the200 covariance between the function values at xand x0.201 A GP is formally defined as a collection of random variables, any finite subset of which follows a202 joint multivariate Gaussian distribution. Thus, for a finite set of known inputs {x(1),x(2), . . . , x(n)}, the203 corresponding function values {f(x(1)),f(x(2)), . . . , f(x(n))}are jointly distributed as:204 f=     f(x(1)) f(x(2)) . . . f(x(n))     ∼ N (m,C),(19) 6 Maria Konstantina Germanou et al. where mis the mean vector and Cis the covariance matrix derived from the covariance function c(x,x0),205 also known as the kernel. This covariance function governs the relationships between the function values206 at different known inputs.207 If the target function f(x)is estimated as a GP, then one needs to check if the consistency requirement208 of the estimator is fulfilled (Papoulis 1991). One can easily do this by the marginalization property of209 multivariate Gaussian distributions. The marginalization property tells us that if210 (f1,f2)∼ N(m,C),(20) Then we also have211 f1∼ N(m1,C11),f2∼ N(m2,C22),(21) where C11 and C22 are sub-matrices of C.212 2.2.6 The Prior in GPR213 In GPR, the prior encodes assumptions about the underlying function before any data is observed214 (Rasmussen and Williams 2006). Two key components specify the prior:215 –Amean function,m(x), typically indicating no prior knowledge about the function’s value before216 data is observed. Often, it is assumed to be zero, i.e., m(x) = 0.217 –Acovariance function or kernel, c(x,x0), which defines the correlation between function values at218 different input points.219 The function we wish to model is assumed to be drawn from a Gaussian process:220 f(x)∼ GP(m(x) = 0, c(x,x0)),(22) Where f(x)represents the function, m(x) = 0 is the mean function (often assumed to be zero), and221 c(x,x0)is the covariance function, which controls the structure and smoothness of the function.222 The GP prior describes a distribution over possible functions, before any observations are made.223 Specifically, the prior assumes that any finite set of function values follows a joint Gaussian distri-224 bution, meaning that the values of f(x1), f(x2), . . . , f(xN)for different input points x1,x2, . . . , xNare225 jointly distributed according to a multivariate Gaussian.226 In summary, the GP prior encapsulates our beliefs about the range of possible functions the model227 can represent, with the covariance function dictating how function values are related across the input228 space.229 2.2.7 Constructing the GPR230 GPR assumes that the function we wish to model, denoted f(·), is drawn from a Gaussian process231 (Murphy 2012). Given a set of training inputs X={x1,x2, . . . , xN}and corresponding observed function232 values f={f(x1), f(x2), . . . , f(xN)}, the goal is to predict the function value f(x0)at a new test input233 x0.234 To achieve this, we model the joint distribution of the function values at the training points and235 the test point. This distribution is Gaussian and described by a mean function m(·)and a covariance236 function c(·,·)(also known as the kernel function).237 The joint distribution of the function values f(x)at the training points {x1,x2, . . . , xN}and the test238 point x0can be expressed as:239      f(x1) . . . f(xN) f(x0)     ∼ N           m(x1) . . . m(xN) m(x0)      ,     c(x1,x1)··· c(x1,xN)c(x1,x0) . . ..... . .. . . c(xN,x1)··· c(xN,xN)c(xN,x0) c(x0,x1)··· c(x0,xN)c(x0,x0)           ,(23) Where:240 –m(xi)is the mean of the Gaussian process at the point xi,241 –c(xi,xj)is the covariance between the function values f(xi)and f(xj), and it is determined by the242 chosen kernel function,243 –N(m,C)represents the multivariate normal distribution with mean vector mand covariance matrix244 C.245 Comparison of Geostatistical and Machine Learning Methods for Spatial Analysis of Natural Resources Data 7 2.2.8 Posterior Distribution246 To predict the function value f(x0)at the test point x0, we compute the posterior distribution by con-247 ditioning the joint distribution of the training data and test point on the observed values at the training248 points (Duvenaud 2014). The posterior distribution of f(x0), given the observed training data f(X), is249 a Gaussian distribution with the following mean and variance:250 251 Posterior Mean:252 E[f(x0)|f(X)] = m(x0) + C(x0, X)C(X, X)−1(f(X)−m(X)),(24) where:253 –m(x0)is the prior mean of the function value at the test point x0,254 –C(x0, X)is the vector of covariances between the test point x0and the training points {x1, . . . , xN},255 i.e., C(x0, X) = [c(x0,x1), . . . , c(x0,xN)],256 –C(X, X)is the covariance matrix of the function values at the training points, with (i, j)-th entry257 given by c(xi,xj),258 –f(X) = [f(x1), . . . , f(xN)]Tis the vector of observed function values at the training points,259 –m(X) = [m(x1), . . . , m(xN)]Tis the mean vector of the function values at the training points.260 The posterior mean E[f(x0)|f(X)] represents the best estimate of the function value at the test point261 x0based on the training data. It is a linear combination of the observed values f(X), where the weights262 depend on the kernel function and the relative distances between the test point and the training points.263 264 Posterior Variance:265 Var[f(x0)|f(X)] = c(x0,x0)−C(x0, X)C(X, X)−1C(X, x0),(25) where:266 –c(x0,x0)is the prior variance (i.e., the variance without observing any data) of the function value at267 the test point x0,268 –C(x0, X)C(X, X)−1C(X, x0)is the reduction in uncertainty due to the information provided by the269 observed data.270 The posterior distribution incorporates information from observed data points (red dots), resulting271 in low variance and high confidence near these points. Further from the training points, the model’s272 uncertainty increases, as indicated by the higher variance. Thus, the posterior variance reflects the273 model’s confidence, decreasing near known data and remaining high in less-sampled regions.274 2.2.9 Prediction275 The posterior distribution enables predictions at new test points x0based on training data {x1, . . . , xN}.276 The posterior mean E[f(x0)|f(X)] serves as the point estimate, while the posterior variance Var[f(x0)|f(X)]277 reflects prediction uncertainty.278 With a Gaussian process and a covariance function like the exponential kernel of Eq. (9) or Eq. (11),279 predictions are smooth and continuous. This dual output of prediction and uncertainty makes GPR ideal280 for applications like natural resources reserves estimation, where high uncertainty indicates areas that281 may need further data.282 2.3 Self-Organizing Map (SOM) Algorithm283 The Self-Organizing Map (SOM), developed by Teuvo Kohonen, is an unsupervised learning algorithm284 that projects high-dimensional data onto a lower-dimensional grid while maintaining the topological285 properties of the input space (Villmann and Bauer 1998).286 8 Maria Konstantina Germanou et al. 2.3.1 SOM Architecture287 The SOM is composed of one input layer and one feature map layer in a two-dimensional space, where288 every neuron in the feature map is connected to all the input dimensions via weight vectors (Kohonen and289 Honkela 2011; Kohonen 2005). The SOM projects high-dimensional input data onto a lower-dimensional290 map while maintaining topological relationships. In training, for every input vector, the Best Matching291 Unit (BMU) is determined by using a similarity measure, such as Euclidean distance, and the weights of292 the BMU and its neighbors are iteratively updated. Over time, the SOM self-organizes so that similar293 input vectors are mapped to neighboring neurons, thereby preserving the topology and density of the294 input space.295 Fig. 1: Self-Organizing Map (SOM) architecture (modified after (Karimi et al. 2023)). Fig. 1 shows how the SOM maps the high-dimensional input space into an orderly two-dimensional296 map: Input vectors are mapped to a two-dimensional feature map. The BMU and its neighboring neurons297 are updated iteratively, enabling the SOM to approximate the input space topology. That makes SOMs298 applicable to clustering and reduction of dimensionality problems, as well as the visualization of complex299 datasets containing nonlinear relationships, which appear more effective than PCA due to its linearity.300 2.3.2 SOM Algorithm301 The SOM algorithm iteratively organizes the input data onto a grid. Below is presented a step-by-step302 layout for the implementation of the algorithm:303 –Step 1: Each neuron’s weight vector wjis initialized randomly over the input space.304 –Step 2: For each iteration k, a random input vector x(k)is selected from the dataset. Here, k305 represents the current iteration number in the training process.306 –Step 3: The BMU is identified by finding the neuron whose weight vector is closest to x(k)in terms307 of Euclidean distance:308 wBMU(k) = min jkx(k)−wj(k)k.(26) –Step 4: The weights of the BMU and its neighbors are updated to move closer to x(k):309 wj(k+ 1) = wj(k) + θ(j, BMU, k)·α(k)·(x(k)−wj(k)),(27) Where:310 –α(k)is the learning rate, decreasing over iterations,311 –θ(j, BMU, k)is the neighborhood function, which decreases with distance from the BMU and with312 k.313 –Step 5: Repeat the steps for all iterations, gradually reducing α(k)and the neighborhood radius314 σ(k)over time until the map stabilizes.315 Comparison of Geostatistical and Machine Learning Methods for Spatial Analysis of Natural Resources Data 9 where krepresents the iteration index of the training process. Each iteration corresponds to one step316 in adapting the SOM to the input data, and kprogresses from 1 to the total number of iterations.317 Once training is complete, the SOM provides an organized map that reflects the structure of the input318 space, grouping similar data and preserving topological relationships. This enables effective clustering319 and visualization of high-dimensional data (Haykin 2009).320 2.3.3 Optimal Map Size321 For proper clustering and correct representation of the data set, different SOM grid configurations should322 be considered. When choosing a suitable map size, the aim is to achieve optimal preservation of spa-323 tial relationships in the dataset by minimizing corresponding errors. Therefore, any map size has been324 considered in terms of two main performance metrics:325 –Quantization Error (QE): Quantization Error measures the average Euclidean distance between326 each data point and its nearest neuron, also known as the Best Matching Unit (BMU). A lower327 QE indicates that the neurons within the grid effectively capture and represent the data’s intrinsic328 characteristics, improving clustering quality.329 –Topographic Error (TE): Topographic Error measures the extent to which the SOM preserves the330 topological structure of the dataset. It is defined as the proportion of data points for which the first331 and second closest neurons (BMUs) are not adjacent. A lower TE suggests that the map maintains332 the spatial continuity of the data, which is crucial for accurate spatial pattern recognition.333 2.3.4 OK Estimation with SOM Neurons334 One of the goals of this research is to present how SOM can be combined with kriging methods to335 improve the neighborhood selection, specifically with OK. The methodology for OK estimation using336 SOM neurons involves:337 1. Data Grouping with SOM: Apply SOM to cluster data points into neurons based on similarity,338 organizing spatial data into a grid for local pattern capture.339 2. Selecting Nearest Neurons: For each prediction location (Xp, Yp), identify the nearest neurons.340 Use data from these neurons to inform the kriging estimation, ensuring local influence is maintained.341 3. OK with Localized Data: Apply OK using data from selected neurons. Model spatial correlation342 with a variogram to weight data points according to their distance from (Xp, Yp), yielding localized343 predictions.344 4. Generating the Spatial Prediction Map: Create a continuous prediction map by applying kriging345 across all target locations.346 Sparse regions in the dataset can result in empty neurons, leading to unreliable or incomplete predic-347 tions. Aggregating data from the knearest neurons (herein k= 5) mitigates this issue by incorporating348 information from neighboring regions, ensuring robust predictions even in areas with fewer data points.349 3 Case Study and preliminary steps350 The data of this case study are drawn from the Geoscience Regulation Office (GSRO) of the Department351 of the Environment, Climate and Communications of the Republic of Ireland. GSRO was established in352 January 2021 and provides a large list of historical mining data (Geoscience Regulation Office Information353 Hub 2021).354 The dataset used in this case study consists of values from 633 exploratory drill holes, representing355 Zinc concentration with a range from 30.0 to 6988.3 ppm, a median of 189.0 ppm, and a mean of 158.5356 ppm. The coordinates of the original data are in the Irish Transverse Mercator (ITM) system, expressed357 in kilometers, and have been converted to the same unit for ease of use. The drill locations are presented358 in Fig. 2. The data in question are in http://gofile.me/75XIE/uBX1jglCy.The analysis is conducted359 with the logarithm of the concentration with units of log of ppm. Thus, the range of the variable values360 is from 3.40 to 8.85, with a median of 5.24 and a mean of 5.07, as shown in Table 1.361 16 Maria Konstantina Germanou et al. Fig. 8: SOM Sample Hits Plot showing the number of times each neuron was selected as the Best Matching Unit (BMU) during training. of similar data points, and finds sparse regions. This visualization ensures that the SOM effectively509 represents the patterns in the dataset, capturing both dense clusters and less populated areas.510 SOM Input Planes511 The Self-Organizing Map (SOM) input planes provide a visual representation of how each input feature512 contributes to the neurons in the SOM grid.513 Fig. 9: SOM Input Planes showing the mapping of different input variables to the SOM neurons. Each plane corresponds to an individual input variable. The third input plane corresponds to the zinc concentration. Each input plane corresponds to a specific input variable and illustrates how the values of that variable514 are distributed across the neurons. The color of each hexagonal cell represents the weight associated with515 a particular input feature, with darker colors indicating higher weights. This visualization helps identify516 the influence of each input feature on the SOM.517 Comparison of Geostatistical and Machine Learning Methods for Spatial Analysis of Natural Resources Data 17 In Figure 9, the input planes represent three features. The first input plane shows the influence of the518 Xspatial coordinate, where darker regions indicate neurons associated with higher X-values. Similarly,519 the second input plane corresponds to the Yspatial coordinate, and darker regions highlight neurons520 influenced by higher Y-values. The third input plane represents zinc concentrations, with darker regions521 revealing areas of higher zinc levels. These input planes collectively illustrate how the SOM captures522 spatial and attribute-based variability within the dataset.523 By comparing these input planes, the relationships between features can be observed. For those524 regions with similar patterns across planes, such as overlapping darker areas, it means that there are525 correlated features that influence the same neurons. On the other hand, regions of dissimilar patterns526 indicate that different features contribute independently to the structure of the SOM. Such comparisons527 give insight into how spatial and concentration patterns interact within the dataset. The SOM input528 planes serve as a validation tool to understand the model’s behavior. They show the features’ relevance529 to the SOM’s structure and allow the cluster identification, where some neurons respond vigorously530 to specific feature values. Moreover, the comparison between input planes enables the identification of531 correlations/contrasts between variables, making it useful to understand which relationships underlie the532 data.533 SOM Weight Positions534 The SOM Weight Positions plot in Figure 10 visualizes the final locations of the weight vectors associated535 with each neuron after training.536 Fig. 10: SOM Weight Positions plot showing the final locations of the neurons’ weight vectors in the input space. The green points represent the input data, blue circles represent the SOM neurons (weight vectors), and red lines indicate the connections between neighboring neurons. These weight vectors represent positions in the input space that the SOM neurons adjust to fit the537 distribution of the input data during training.538 The proximity of blue circles to green points indicates how well the SOM represents the input data.539 Neurons close to dense clusters of green points reflect accurate mapping of data-rich regions, while540 sparse regions with fewer neurons show that the SOM has adapted to areas with less data. The red lines541 connecting neighboring neurons illustrate the SOM’s topological structure, with shorter lines representing542 similar regions and longer lines highlighting transitions between clusters or sparse areas.543 The arrangement of weight vectors in the plot demonstrates how well the SOM has adapted to544 the data’s structure. Regions with dense clusters of green points typically have neurons (blue circles)545 positioned closer together, indicating successful mapping of similar data points. The red lines between546 18 Maria Konstantina Germanou et al. neighboring neurons ensure that the topological relationships in the input space are preserved during547 the SOM training process.548 4.3.2 OK Prediction with SOM Neuron Grouping549 The OK Estimation with SOM Neuron Grouping (OK-SOM), as shown in Figure 11, illustrates550 the application of OK for predicting zinc concentration across the study area by utilizing SOM to group551 the data into neurons.552 Fig. 11: OK-SOM showing the spatial distribution of zinc concentration (ppm) across the study area. Values in ppm. This method enables spatial prediction of zinc concentration in ppm after back-transformation from553 the log scale. That is, the estimation is made in log scale and then transformed into ppm as shown554 in the figure. This plot shows the kriging estimation, where regions with higher zinc concentration are555 notably in the central part of the area, indicating possible hot spots of mineralization that merit further556 investigation. These areas reflect significant zinc accumulations, indicating favorable geology for mineral557 deposits. The areas on the periphery yield lower concentrations, indicating a slight chance of significant558 mineralization.559 With this approach, OK interpolates the zinc concentrations in unsampled locations based on the560 spatial autocorrelation of the measured data points, providing continuous predictions for the whole area561 under study, even for areas with no direct measurements. Moreover, in this plot, the grid structure is562 influenced by the fact that the SOM neuron grouping divides the region into distinct clusters or neurons,563 each of which applies kriging. This grouping ensures that the local characteristics of each neuron are564 considered in the estimation, enhancing the accuracy and reliability of the predicted zinc concentrations565 throughout the region.566 The plot provides a continuous spatial prediction, which enables the identification of areas with567 elevated or reduced zinc concentration. Additionally, it helps in identifying the specific area of interest568 for further exploration, especially the areas with high estimated zinc concentrations that may require569 closer investigation. OK-SOM indicates the spatial variability in zinc concentration, showing possible570 local anomalies and broader regional trends.571 Table 5 presents the validation measures derived from LOOCV with OK-SOM. Compared to the572 range of values (30-7000 ppm), RMSE and MAE are acceptably low. Correlation ρ= 91% shows a robust573 Comparison of Geostatistical and Machine Learning Methods for Spatial Analysis of Natural Resources Data 19 Table 5: Performance of OK-SOM. RMSE: Root Mean Square Error, MAE: Mean Absolute Error, ME: Mean Error, ρ: Pearson’s correlation coefficient. Metric OK-SOM RMSE (ppm) 490.93 MAE (ppm) 153.26 ME (ppm) -115.46 ρ0.91 correlation between the observed and predicted values. ME of −115.46 shows a slightly negative bias, at574 -1.65% of the range.575 4.4 Comparisons576 The engaged methods for the mathematical modeling of the spatial dependencies are compared to provide577 a better overview of the results. Table 6 summarizes the LOOCV metrics for all the methods. Table 6: Validation measures. RMSE: Root Mean Square Error, MAE: Mean Absolute Error, ME: Mean Error, ρ: Pearson’s correlation coefficient. Metric OK GPR OK-SOM RMSE (ppm) 485.79 532.14 490.93 MAE (ppm) 186.97 190.46 153.26 ME (ppm) -71.9 -95.88 -115.46 ρ0.76 0.73 0.91 578 The neighborhood selection is different for each method. GPR takes into account the entire data for579 each prediction point. OK uses a neighborhood of 300 m, based on the kernel bandwidth (correlation580 range). OK-SOM takes as neighbors all the data points in the cluster. All methods produced RMSE of581 the same magnitude, with OK achieving slightly better RMSE than the rest. The correlation coefficient582 is satisfactory for OK and GPR, with OK showing marginally better results. OK-SOM outperforms583 them all with a strong correlation of ρ= 91%.The significant improvement of the last-mentioned metric584 highlights the combination of a stochastic approach with a classification method.585 ME presents an interesting low value for the OK method. Since a non-linear transformation is per-586 formed on the data, a bias error is introduced in the back transformation (Goovaerts 1997). This can be587 attributed to the benefits of the local OK neighborhood selection, which also explains the lower RMSE588 of OK. However, this leads to a contradiction, as the BMU classification system of OK-SOM specifically589 assists OK in this regard. So it is a profound insight that the neighborhood selection for the OK pre-590 diction supports the ME and thus RMSE. In contrast, the OK-SOM neighborhood selection supports591 the correlation of the prediction with the observation. This outcome should be further investigated in592 different random fields.593 4.4.1 Spatial Difference Visualization: OK vs. GPR594 Figure 12 illustrates the spatial distribution of differences in zinc concentration predictions (in ppm)595 between the models generated by OK and GPR. The color bar illustrates that differences range from596 approximately zero ppm in the blue areas to values close to 1300 ppm in the dark red areas. Warmer597 colors indicate regions with greater differences in predictions. White values indicate that OK and GPR598 values were effectively the same (within 0.1 ppm difference). In most of the mines, the differences are599 below 30 ppm, but there are substantial differences in the Zinc-rich areas of the mine. Visual comparison600 of Figs. 3, 4 shows that the smaller neighborhood used in OK results in less smoothing than the GPR.601 The differences in predicted zinc concentrations that are observed provide some critical insights into602 the performance of each model. GPR gives considerably higher predictions than OK in the central regions603 of the study area, where dark red values are present. This pattern indicates that GPR effectively captures604 localized variations in zinc concentration, underlining the ability of the model to respond to small-scale605 spatial trends. Contrarily, peripheral blue areas exhibit a minimal difference, at which OK’s prediction is606 close to GPR’s. This is a consequence of the kriging smoothing effect, compared with the more adaptive607 GPR approach.608 20 Maria Konstantina Germanou et al. Fig. 12: Absolute values of Spatial differences in zinc concentration predictions (ppm) between OK and GPR. Values in ppm. 4.4.2 Scatter Plot Comparison OK vs. GPR609 Figure 13 compares OK and GPR predictions with the observed zinc concentrations (ppm). The scatter610 plot compares the predictive performance of OK and GPR against the observed zinc concentrations (in611 ppm). The y=xblack dashed line represents perfect prediction and thus serves as a benchmark. OK612 performs quite well for lower to moderate concentrations (<1000 ppm); its predictions are close to613 the observed values, which reflects the ability of OK to produce stable and smooth estimates. However,614 OK gives a slight underestimation at higher concentrations (>3000 ppm). On the other hand, GPR is615 more variable for all concentrations and better reflects localized spatial trends and finer heterogeneity.616 However, GPR gives a greater underestimation of extreme concentrations (>3000 ppm), as reflected by617 wider deviations, reflecting a problem in generalizing sparse or highly variable data. Overall, the scatter618 plot highlights the complementary strengths of the methods: OK excels in capturing broader spatial619 trends with consistency, while GPR captures detailed variability but may introduce higher uncertainty620 for extreme values.621 4.4.3 Spatial Difference Visualization: OK vs. OK-SOM622 Figure 14 presents the spatial distribution of absolute prediction differences in zinc concentrations (ppm)623 between OK and OK-SOM. The color scale ranges from minor differences (near zero ppm) to substantial624 discrepancies (up to 1300 ppm). The visualization underscores the predictive tendencies of the two625 methods. OK provides smoother, continuous predictions suited for broader regional analyses, whereas626 OK-SOM highlights localized anomalies and spatial clusters.627 4.4.4 Scatter Plot Comparison OK-SOM vs. OK628 Figure 15 depicts a scatter plot of OK and OK-SOM against observed zinc concentrations. The predictions629 by OK reflect relatively stable and consistent performance for all concentration ranges. There is a slight630 underestimation (negative bias) in the OK-SOM predictions (shown in Tables 5, 6), especially where the631 observed concentrations are greater than 1000 ppm.632 Comparison of Geostatistical and Machine Learning Methods for Spatial Analysis of Natural Resources Data 21 Fig. 13: Scatter plot comparing original zinc concentrations (ppm) with predictions from OK and GPR. Blue points: GPR predictions. Red points: OK predictions. Back dashed line: perfect prediction. Fig. 14: Absolute values of Spatial differences in zinc concentration predictions (ppm) between OK and OK-SOM. Values in ppm. In contrast, the OK-SOM predictions show larger variability, indicating a more detailed mapping633 of the data. SOM performs well at lower to moderate concentrations (< 3000 ppm), showing closer634 agreement with the observed values and resolving finer spatial trends. This variability reflects the OK-635 SOM’s ability to account for localized spatial heterogeneity. However, at higher concentrations (> 4000636 ppm), the OK-SOM predictions are showing greater deviation from the perfect prediction line, which637 explains the slightly elevated RMSE shown in Table 6 compared to OK, as RMSE is affected more by638 outliers.639 22 Maria Konstantina Germanou et al. Fig. 15: Scatter plot comparing original zinc concentrations (ppm) with predictions from OK-SOM. Red points: OK predictions. Green points: OK-SOM predictions. Black dashed line: perfect predictions. Overall, the results indicate the complementary strengths of the two methods. While OK gives640 smoother and more generalized predictions suitable for broader spatial scales, OK-SOM is superior in641 capturing local variability and fine-scale patterns in the data. This duality underlines the added value642 of combining these approaches, with OK-SOM providing detailed predictions in data-rich areas and OK643 ensuring robust estimates across larger regions.644 4.5 Uncertainty quantification645 An essential aspect of natural resources analysis is uncertainty quantification. Both GPR and OK can646 produce uncertainty maps, reflecting the variance of the predictions. However, both methods require647 Gaussianity for the uncertainty to have a physical meaning (Pavlides et al. 2015).648 The logarithmic transform moves the data closer to the normal distribution, as can be seen from649 Table 1. The first step of the back-transform is650 ˆ X(u) = exp( ˆ X∗(u)),(29) Where ˆ X∗(u)is the estimation of the logarithm of the concentration in the prediction point u. In651 the case of the presence of bias, a bias correction factor Bcan be used in eq. (29) as652 ˆ X(u) = Bexp( ˆ X∗(u)).(30) (Goovaerts 1997; Hristopoulos 2020)653 However, while back-transforming the predictions to the original scale works well seamlessly, the same654 cannot be said for the Error Variance, as it is additive to the log-transformed value, not the original655 value. Thus, if εis the estimation error, so that the true log‐value is656 X(u) = ˆ X∗(u) + ε(u) rather than ˆ X∗(u), then under the back‐transform,657 expX∗(u) + ε(u)= expX∗(u)expε(u)6= expX∗(u)+ expε(u). Hence, it is not possible to map the original scale uncertainty by simply exponentiating the log‐scale stan-658 dard deviation. Conditional Simulations can be used to create maps of uncertainty, but this was beyond659 the scope of this research. For the simulations to effectively capture uncertainty, Gaussian anamorphosis660 should be conducted first (Pavlides et al. 2022).661 While maps of uncertainty are not included, to give a quantification of the uncertainty of the esti-662 mation of each method, relative absolute errors (RAE) have been calculated, as shown in Eq. (31).663 Comparison of Geostatistical and Machine Learning Methods for Spatial Analysis of Natural Resources Data 23 RAEi=ˆ X(ui)−X(ui) X(ui), i = 1, . . . , n. (31) E20 and E80 represent the relative absolute error at specific percentiles, shown in Table 7. In geostatistics,664 such percentile relative errors are used to assess the spread and magnitude of errors through quantiles:665 E20 as the 20th percentile (the relative absolute error below which 20% of errors lie), E80 as the 80th666 percentile (the relative error below which 80% of errors lie).667 Table 7: Percentile Relative Errors (E20 and E80) for OK, GPR, and SOM‐Grouped Kriging Metric OK GPR OK-SOM Neurons E20 (Relative Error at 20th Percentile) 31.59% 32.54% 45.50% E80 (Relative Error at 80th Percentile) 84.04% 85.31% 92.89% Both OK and GPR exhibit similar performance in terms of percentile errors. GPR shows an E20 of668 32.54% and an E80 of 85.31%, while OK achieves an E20 of 31.59% and an E80 of 84.04%. This indicates669 that both methods perform comparably in predicting values within these specified percentile ranges.670 Considering the results of Table 6, it would be unexpected that OK-SOM would show higher relative671 errors. However, as Eq. (31) shows, Relative errors are based on the initial value X(ui). Thus, if X(ui)672 is low, a slight deviation in ˆ X(ui)could lead to a higher relative error. As can be seen from Fig. 14,673 there are differences in the range of 15 to 50 ppm in areas of low Zinc concentration (30-100 ppm). Thus,674 while RMSE is only slightly affected, the relative errors are higher. Furthermore, while OK-SOM has675 a better correlation coefficient, it has somewhat higher Mean Error, which would also indicate higher676 relative error.677 5 Conclusions678 This work presents an overview of research on the integration of geostatistical methods with machine679 learning techniques to enhance spatial data modeling of natural resource reserves, focusing on OK, GPR,680 and SOM. Several kernels were tested to determine the data spatial dependence, including the novel HCE.681 In preliminary analysis, the performance of HCE was found to be comparable to the exponential kernel;682 however, the latter was selected due to its slightly superior cross-validation performance. The main goal683 of enhancing prediction accuracy and quantifying spatial uncertainty was successfully achieved.684 Gaussian Process Regression effectively captured the non-linear and complex spatial patterns, yielding685 satisfactory validation errors. As all methods employed are stochastic-based, GPR produced uncertainty686 bounds comparable to those of Ordinary Kriging (OK), and slightly narrower than those obtained with687 SOM-OK. OK also provided satisfactory estimates. Moreover, applying a 300 m search neighborhood in688 OK helped reduce excessive smoothing in the prediction map.689 Self-Organizing Maps introduced a novel clustering approach for spatial data, based on the similarity690 between neurons. Combined with kriging, SOM increased the predictability, reflected in a higher corre-691 lation coefficient and improved error metrics. Fine-scale anomalies and spatial clustering were captured692 more effectively with the presented approach, allowing for improved localized accuracy and a more nu-693 anced representation of spatial variability. An important result from the comparison between OK-SOM694 and OK, is that OK gives significantly less bias, (ME=-71.9 ppm compared to ME=-115.5 ppm for OK-695 SOM). At the same time, OK-SOM primarily enhances the ρ, giving significantly better values (ρ= 91%)696 compared to the other methods that give ρ= 73% (GPR) and ρ= 76% (OK).697 In conclusion, the methods have different strengths: OK and GPR provide broad and stable predic-698 tions, and Self-Organizing Maps enhance localized insight and clustering to deal with complex spatial699 patterns. These methods together set up a robust framework for addressing intricate spatial datasets,700 balancing accuracy, interpretability, and the representation of spatial variability.701 Future work could focus on further uncertainty quantification through conditional simulations, more702 efficient methods for the Gaussian anamorphosis of the data, as well as exploring more kernels, like the703 Advanced Harmonic Covariance Estimator.704 24 Maria Konstantina Germanou et al. Acknowledgement705 The research project is implemented in the framework of H.F.R.I call ”Basic research Financing (Hori-706 zontal support of all Sciences)” under the National Recovery and Resilience Plan ”Greece 2.0” funded707 by the European Union – NextGenerationEU (H.F.R.I. Project Number: 16537).708 Conflict of Interest709 The authors declare no conflicts of interest.710 Comparison of Geostatistical and Machine Learning Methods for Spatial Analysis of Natural Resources Data 25 References711 Bação F, Lobo V, Painho M (2005) Self-organizing maps as substitutes for k-means clustering. Computers712 & Geosciences 31(5):545–553713 Chiles J P, Delfiner P (2012) Geostatistics: Modeling Spatial Uncertainty. Wiley, 2 edition714 Cressie N (1990) The origins of kriging. Mathematical Geology 22(3):239–252, ISSN 0882–8121715 Cressie N (1993) Spatial Statistics. New York: John Wiley and Sons716 Duvenaud D (2014) Automatic model construction with Gaussian processes. Ph.d. dissertation, Univer-717 sity of Cambridge718 Geoscience Regulation Office Information Hub (2021) Mineral exploration and mining web applications.719 Technical report, Department of the Environment, Climate and Communications, accessed: 2025-03-1720 Goovaerts P (1997) Geostatistics for Natural Resources Evaluation. Oxford University Press721 Haykin S (2009) Neural Networks and Learning Machines. Prentice Hall, 3 edition722 Hengl T, Heuvelink G B M, Stein A (2004) A generic framework for spatial prediction of soil variables723 based on regression-kriging. Geoderma 120(1-2):75–93724 Hristopoulos D T (2020) Random Fields for Spatial Data Modeling: A Primer for Scientists and Engi-725 neers. Advances in Geographic Information Science, Springer Netherlands726 Hristopulos D T (2024) Non-separable covariance kernels for spatiotemporal gaussian processes based727 on a hybrid spectral method and the harmonic oscillator. IEEE Transactions on Information Theory728 70(2):1268–1283729 Karimi E, Haghighi F, Sheykhfard A, Azmoodeh M, Shaaban K (2023) Self-organized neural network730 method to identify crash hotspots. Future Transportation 3(1):286–295731 Kohonen T (2001) Self-Organizing Maps. Springer, 3 edition732 Kohonen T (2005) Intro to som. http://www.cis.hut.fi/projects/somtoolbox/, retrieved 2025-06-18733 Kohonen T, Honkela T (2011) Kohonen network. Scholarpedia 2(1):1568, bibcode: 2007SchpJ...2.1568K734 Krige D G (1951) A statistical approach to some basic mine valuation problems on the witwatersrand.735 Journal of the Chemical, Metallurgical and Mineral Society of South Africa 52:119–139736 Liu Y, Weisberg R H (2005) Patterns of ocean current variability on the west florida shelf using the737 self-organizing map. Journal of Geophysical Research: Oceans 110(C6)738 Liu Y, Weisberg R H (2011) A review of self-organizing map applications in meteorology and oceanog-739 raphy. In Self-Organizing Maps: Applications and Novel Algorithm Design, 253–272740 Matheron G (1965) Les variables régionalisées et leur estimation: une application de la théorie de fonctions741 aléatoires aux sciences de la nature. Paris: Masson et Cie, École Nationale Supérieure des Mines de742 Paris743 Murphy K P (2012) Machine Learning: A Probabilistic Perspective. The MIT Press744 Papoulis A (1991) Probability, Random Variables, and Stochastic Processes. New York: McGraw-Hill, 3745 edition746 Pavlides A, Agou V D, Hristopulos D T (2022) Non-parametric kernel-based estimation and simulation747 of precipitation amount. Journal of Hydrology 612:127988748 Pavlides A, Hristopulos D T, Roumpos C, Agioutantis Z (2015) Spatial modeling of lignite energy reserves749 for exploitation planning and quality control. Energy 93:1906–1917750 Rasmussen C E, Williams C K I (2006) Gaussian Processes for Machine Learning. MIT Press751 Varouchakis E A, Corzo G A, Karatzas G P, Kotsopoulou A (2018) Spatio-temporal analysis of annual752 rainfall in Crete, Greece. Acta Geophysica 66(3):319–328753 Varouchakis E A, Koltsidopoulou M D, Pavlides A (2025) Designing robust covariance models for geo-754 statistical applications. Stochastic Environmental Research and Risk Assessment755 Varouchakis E A, Solomatine D, Corzo Perez G A, Jomaa S, Karatzas G P (2023) Combination of756 geostatistics and self‐organizing maps for the spatial analysis of groundwater level variations in complex757 hydrogeological systems. Stochastic Environmental Research and Risk Assessment 37(8):3009–3020758 Vesanto J, Alhoniemi E (2000) Clustering of the self-organizing map. IEEE Transactions on Neural759 Networks 11(3):586–600760 Villmann T, Bauer H U (1998) Applications of the growing self-organizing map. Neurocomputing761 21(1):91–100, ISSN 0925-2312762