Identifying Probable Submarine Hydrothermal Spots in North Santorini Caldera Using the Generalized Moments Method

Abstract

The presence of active hydrothermal vent fields near residential areas and their possible link to volcanic activity poses a potential hazard to the environment, society, and the economy. By capitalizing on Autonomous Underwater Vehicle sampling methodologies and applying the Generalized Moments Method model for geological and physical processes in these environments, we shed light on the underlying dynamics shaping the physicochemical characteristics of the vents. In this study, we focus on the Northern Caldera of Santorini and, more specifically, on the recorded CTD data (Conductivity, Temperature, Depth). The data sets were collected in 2017 in Santorini using an Autonomous Underwater Vehicle during the GEOMAR POS510 mission. Our research shows that the active vent field within the caldera probably follows a multifractal behavior and exhibits a weak memory effect. Depth Profiles and Time Series show similar behavior among conductivity and temperature. The variance and moments of both parameters underline the existence of two different mechanisms governing the behavior of the vent field. Finally, the structure function shows that changes in the time series are described by a Cauchy–Lorentz distribution.

1. Introduction

The active Hellenic Volcanic Arc (HVA; e.g., [1]), also known as the South Aegean Volcanic Arc (e.g., [2]) in the Eastern Mediterranean is the result of the subduction of the Eastern Mediterranean lithosphere beneath the active Hellenic margin of the European plate [1,3,4,5]. It extends from the Saronic Gulf in the west to the Kos-Nisyros Complex in the east [6]. The HVA is characterized by active volcanism both onshore and offshore [7], with notable volcanic centers such as the Kolumbo volcanic chain at the Christina-Santorini-Kolumbo volcanic field [8] and the Kos-Nisyros volcanic Complex (Figure 1) [9,10].

Santorini Volcano is located in the central part of the Christiana-Santorini-Kolumbo volcanic field [8]. It consists of three islands surrounding a flooded 10 × 7 km diameter caldera. The volcanic complex of Santorini includes the islands of Thera, Therasia, and Aspronisi, arranged in an oval shape around the caldera. Thera, the largest of the islands, has a semicircular shape with its concave side towards the west, where Therasia and Aspronisi lie (Figure 2). Within the caldera are two smaller islands, Palea Kameni and Nea Kameni (the Kameni Islands) [11].

The caldera walls rise to 300 m above sea level, while the maximum depth of the caldera seafloor is about 390 m. There are three distinct flat-floored morpho-sedimentary basins [12], as shown in Figure 3. The Northern Basin is the largest and deepest of the three (389 m). On the other hand, the Western Basin, the smallest, has a depth of 325 m, while the Southern Basin, south of the Kameni Islands, forms the southern part of the Santorini caldera and has the shallowest sea bottom at 297 m deep. The caldera is connected to the sea via three breaches: one in the NW and two in the SW (Figure 3).

Figure 2. Aerial photograph of Santorini caldera from the north [13].

Figure 2. Aerial photograph of Santorini caldera from the north [13].

In the recent volcanic history of Santorini, spanning approximately 0.36 million years [15], there have been twelve significant eruptions, numerous minor eruptions, and at least four caldera collapses [11,13]. The most significant explosive eruption, known as the Late Bronze Age (LBA) or the Minoan eruption, occurred around 3600 years ago [16] and is regarded as one of the most significant eruptions in the Holocene [15,17,18].

Following the LBA, volcanic activity in Santorini continued but with diminishing intensity, characterized by several less intense eruptions. Over the past 600,000 years, the region has witnessed at least seven significant explosive eruptions [11,17,19,20]. To this day, Santorini remains the most active volcanic field in the HVA [21], with the most recent eruption occurring early in 1950.

So far, no high-temperature hydrothermal vents have been found in the caldera floor during Remotely Operated Vehicle (ROV) explorations [7,22]. However, a low-temperature hydrothermal vent field in the northeastern part of the northern basin covers an area of over 200 × 300 m² [7]. Additionally, CO₂-rich fluid pools with slightly higher temperatures than the surrounding waters have been discovered at the base of the caldera wall and at water depths 200–250 m, connected by narrow channels downslope [23]. Monitoring changes in temperature and the chemistry of these areas can provide valuable indicators of volcanic unrest alongside traditional methods such as seismicity and fumarolic activity.

Despite its importance for risk assessment and mitigation, monitoring the activity of submarine volcanoes is impeded by the remoteness and extreme conditions of underwater volcanoes [24]. An efficient method to study those characteristics is to employ sensors aboard Remotely Operated underwater Vehicles (ROV) or Autonomous Underwater Vehicles (AUV) able to record CTD data (Conductivity, Temperature, Depth) in the water column [25,26]. Such data were collected for the first time in the northern basin of the caldera of Santorini in 2017 [27]. One of the mission’s objectives was to observe the previously unexplored hydrothermal fields in the caldera using AUV. In the present work, results from the entire investigation are reported for the first time regarding high-frequency CTD measurements using sensors mounted on an AUV.

The CTD data have been sorted and analyzed using the Generalized Moments Method (GMM), a relatively recently developed mathematical method with various applications in finance time series [28,29,30] and environmental cases [29,31,32]. Our group has successfully used the same modeling method previously to identify the underlying stochastic processes driving the activity of other volcanoes, such as the Kolumbo volcano [33] and Nisyros volcano [34,35], followed by others, e.g., in the case of the shallow underwater volcano El Hierro [36,37]. Previous studies yielded exciting aspects of the underlying dynamics of hydrothermal vent fields during rest and unrest periods.

2. Materials and Methods

2.1. The Mission

In March 2017, the cruise of POS510 embarked on Catania, Italy, and headed to Heraklion, Greece (500 NM transit) to investigate the initiation of arc rifting and associated back-arc hydrothermal activity [27]. The GEOMAR-led cruise included the following objectives: (i) an evaluation of the NE-trending Christiana-Santorini-Kolumbo line, (ii) the investigation of the thermal and structural history of the Anydros Basin and its link to hydrothermal activity, and (iii) a time series and lateral facies analysis of Pleistocene to Holocene tephra.

Using integrated geophysical and geological approaches, the mission aimed to answer critical questions related to the arc rifts’ thermal and structural evolution. The cruise was held by the Research Vessel (R/V) POSEIDON and lasted 25 days, 19 of which were on-site operations. In comparison, the remaining 6 days were necessary for mobilization and demobilization, personnel transfers, and transits.

The AUV Abyss from GEOMAR has an operation capacity of up to 6000 m in water depth. The system comprises the AUV itself, a control and workshop container, and a mobile Launch and Recovery System with a deployment frame that was installed at the afterdeck of R/V POSEIDON [27]. According to the cruise report, the AUV missions were planned based on the previously obtained swath data [7,38,39].

Typically, the AUV Abyss (Figure 4) uses a Long Baseline system to support its navigation. Due to the low depth of the area and the close coastlines, dives were made possible without transponders. During the dives, the pre-calibrated CTD sensor continuously acquired data at a rate of 4 Hz. During the POS510 cruise, fifteen missions were completed by the AUV Abyss. The turbidity and the CTD sensors ran simultaneously and served as secondary sensors. The CTD data sensor was a Seabird SBE49 FastCAT, S/N 4948793-0168, and its last calibration before the mission was on 28 January 2016. The CTD sensor is capable of recording data up to the fourth decimal digit, providing ample precision to capture even subtle variations. All data have a time stamp and are related to the vehicle’s position [27]. Apart from conductivity and temperature, the sensors also recorded salinity and sound velocity.

2.2. The Generalized Moments Method

The GMM has been previously applied to study the underlying mechanism driving the CTD times series in two cases: the first recorded over a very active submarine hydrothermal vent field (Kolumbo, Santorini) [33] and the other recorded over a less active area (Nisyros caldera) [34,35]. Details of the novel application of the GMM are presented in those references, where the reported results have also established the validity of the novel approach. In the following few paragraphs, we briefly overview the method.

Hydrothermal activity may be detected by measurements of the physical and chemical properties of the fluid obtained from CTD depth profiles and time series. Such a sequence of events (measured property versus time) can be assumed as a manifestation of a stochastic process. In our case, such a process may be the complex hydrothermal vent activity, which shapes the physical and chemical parameters of the vent field. A first estimate of the behavior of a stochastic process, x(t), can be made in terms of its variance. For discrete data time series, the variance reads:

$$W\left(\Delta ,T\right)=\frac{1}{M}{\sum }_{i=1}^M\frac{1}{\Delta -T}{\sum }_{n=1}^{T-\Delta }{\left(x_i(\Delta +n)-x_i\left(n\right)\right)}^2-{(\frac{1}{M}{\sum }_{i=1}^M\frac{1}{\Delta -T}{\sum }_{n=1}^{T-\Delta }\left(x_i(\Delta +n\right)-x_i\left(n\right)))}^2,$$

(1)

where $M$ is the number of realizations, and for $M=1$, the first summation is omitted. The total length of the trajectory, $T=N\times \tau$, is expressed in terms of the elapsed time between two consecutive measurements, τ, or minimum lag time, which is the reciprocal of the sampling rate, $f$. The lag time, $\Delta$, takes the role of the time, and its values are in the range $\tau <\Delta <\frac{T}{10}$.

The variance in either continuous or discrete space is expressed in terms of the first two moments of the probability density function. This suffices when we deal with random processes drawing steps from a Gaussian distribution. Nevertheless, if extreme events contribute significantly to the probability density function, which a Gaussian distribution can no longer describe, more moments of the probability density function are necessary to draw conclusions.

By fitting the graphs of the variance with the equation W = bΔ^γ, we gain the value of the exponent γ, which classifies the type of process. More precisely, if $\gamma =1$, then we have normal or Brownian motion; if $0<\gamma <1$, the behavior is subnormal or subdiffusional; and lastly, if $1<\gamma <2$ and supernormal or superdiffusional.

More insights into the mechanisms governing a stochastic process can be extracted by using sophisticated methods appropriate for time series analysis [29,40]. Among them, the GMM, which has been successfully applied in numerous fields, is robust and works well even for short time series [31,41].

In brief, the GMM works as follows:

• Construct a time series with different lag times.
• Estimate the statistical moments. See Equation (3).
• Show the moments scale according to Equation (4).

In more detail, firstly, we construct a time series with different lag times Δ.

$$y_n\left(\Delta \right)=\sqrt{{{x\left(n+\Delta \right)}^2-x\left(n\right)}^2}.$$

(2)

Then, the moments of $y_n\left(\Delta \right)$

$$\rho \left(q,\Delta \right)=\frac{1}{T-\Delta }{\sum }_{n=1}^{{\rm T}-\Delta }{\left(y_n\left(\Delta \right)\right)}^q,$$

(3)

can be obtained. If $q$ stands for the order of the moment and $\Delta$ is the lag time, then according to this method, the moments, $\rho \left(q,\Delta \right)$, are expected to scale with $\Delta$ as

$$\rho \left(q,\Delta \right)\approx {\Delta }^{z\left(q\right)},$$

(4)

where $z\left(q\right)$ is the structure function, which identifies the stochastic mechanism or mechanisms that drive the random process. This equation considers moments lower and higher than two. Lower moments are responsible for the core of the probability density function. Higher moments contribute to the wings of the probability density function. Furthermore, the form, $z\left(q\right)$, gives information about the types of motion [42].

Stochastic processes can be of varying degrees of complexity. When the structure function takes the simple form $z\left(q\right)=Hq$, there is a direct relation between the exponent $\gamma$ and the cofactor of the structure function, namely $H=2\gamma$, where H is the Hurst exponent. If $z\left(q\right)$ is a linear function of q or is linear in different portions (bilinearity), then the mechanisms that likely govern the stochastic sequence are additively decomposable, e.g., fractional Brownian motion, fractional Gaussian noise, and Lévy flights and/or walks. The linear dependency of $z\left(q\right)$ on q describes monofractal processes. Instead, $z\left(q\right)$ has a convex shape as a function of the order of the moment, $q$, then multiplicative (multifractal) processes drive the stochastic sequence [43], which means that $H$ is a function of $q$. For some q values, the exponents $H\left(q\right)=z\left(q\right)/q$ are associated with unique features. For $q=1$, the value $H\left(1\right)$ describes the scaling behavior of the absolute value of the increments. It is identified as the Hurst exponent if the process is monofractal. For $q=2$, the value of $H\left(2\right)$ is associated with the scaling exponent of the autocorrelation function and is related to the power spectrum [44]. In the context of the behavior of the two parameters in this study, multifractal properties of the conductivity field indicate that the volcano continuously provides material (ions) to the marine environment. At the same time, multifractality in temperature signifies the continuous flux of energy into the environment. Among multifractals, universal multifractals are likely to be ubiquitous [31,45,46], and the structure function reads [31]:

$$z\left(q\right)=Hq-\frac{C}{a-1}\left(q^a-q\right),$$

(5)

where $H$ defines the scaling of the mean-field of $\Phi$, and for $H=0$, the field is conservative, $C$, where $C>0$ measures the mean homogeneity of the field, and $\alpha$, where $0\le \alpha <2$ is the Lévy index, that is, the α–stable Lévy distribution from which the process draws its changes. The value of α shows how fast the inhomogeneity increases with the order of the moments. For $\alpha =1$, $z\left(q\right)=Hq-Cqlog\left(q\right)$, and the distribution draws changes from a Cauchy–Lorentz distribution. For $\alpha =2$, Equation (5) reads $z\left(q\right)=Hq-C(q^2-q)$, and the distribution draws changes from a lognormal distribution.

In the present study, the GMM is performed on the temperature and conductivity time series recorded over the seafloor of the northern basin of the caldera of Santorini. All measurements considered in the analysis presented below have been carried out near the seabed directly above the vent sources, which, despite being weak, are ostensibly active and present anomalies in the CTD depth profiles. The results from applying the method are reported in the next section.

3. Results

3.1. Depth Profiles

Santorini was surveyed on three separate days: 11 March 2017, Area A; 25 March 2017, Area B; and 26 March 2017, Area C (Figure 5). Although a few vents have been discovered previously in Area A (Sigurdsson et al., 2006), we focused on the entire basin, wherever the GMM was practicable. More specifically, the dataset was split into tracks, six of which the GMM is considered directly applicable (see maps in Figure 6). Two tracks came from Area A (Tracks A1, A2), three from Area B (Tracks B1, B2, B3), and one from Area C (Track C1); see also

Table 1. Tracks used for the GMM in the case of Santorini.

Tracks	Time	Total Time (in Mins)
Track A1	24:50–25:20	30
Track A2	26:00–26:45	45
Track B1	19:40–20:30	50
Track B2	20:50–21:45	55
Track B3	22:00–23:00	60
Track C1	20:40–21:10	30

. Tracks analyzed were those with an AUV interlude of 30 min or more to ensure the reliability of the results. Each of the tracks chosen for the GMM satisfies different prerequisite parameters. Only Track B3 was an hour long, which is why it will be used here to present the results. The areas where the AUV passed over an underwater relief were excluded from the analysis.

Images from the seafloor revealed a lack of significant fluid outflow from seabed edifices, bubble formation, and black smoke, compared to previous studies in Santorini’s neighboring submarine volcano, Kolumbo, where intense CTD anomalies were found [25]. Nevertheless, the mathematical model was applied to data recorded in separate locations within a period of half an hour or more. At first glance, anomalies were detected at a depth of 342–346 m in almost all vertical profiles (depth profiles) for both conductivity and temperature. Salinity and Sound Velocity, on the other hand, showed a straight line in all cases. Maps of the separate tracks are depicted in Figure 6, while all depth profiles appear in Figure 7, Figure 8, Figure 9 and Figure 10. On 27 March 2017, while the AUV was going through Area B, there was a track where the vehicle covered the shortest distance near the seafloor for an hour exactly (Track B3; 22:00–23:00; Figure 6). This track will be used in the next section to present the results of the GMM.

3.2. The Generalized Moments Method Results

CTD time series were constructed from the raw data for each track shown in Table 1. A particular focus was given to the temperature and conductivity time series, as they are physical properties directly measured by the sensors and can provide helpful information on the thermodynamic state of the system. Moreover, these parameters provide information about the underlying mechanisms governing the active volcanoes.

In the case of Track B3, which is used here as an analysis example, both conductivity and temperature variances exhibit two discrete regimes characterized by different slopes, with a turning point at Δ = 5 s. This value for the turning point is practically the same as the turning point found in the case of Kolumbo [33] (Δ = 5.1 s), where intense vent activity was observed in 2010. A turning point indicates that the behavior of the parameters is different through time, in the sense that there is not one single cause that drives the behavior of the parameters but perhaps multiple. In these data, two causes are evident, with the existence of two regimes, respectively, and two different behaviors of the same element.

The distinct regimes in the plots are treated separately by fitting the variances of conductivity and temperature with a function of the form W = bΔ^γ. Figure 11 shows the variance for conductivity and temperature for Track B3. In the first regime for conductivity, for Δ < 5 s, a value of γ = 0.282 ± 0.013 is found, while in the second regime, γ = 1.115 ± 0.003. For temperature, for Δ < 1 s, we obtain a value of γ = 0.173 ± 0.010 in the first regime, whereas in the second γ = 0.931 ± 0.005.

Table 2. γ exponents extracted from the two regimes of the variance for conductivity and temperature (Santorini).

Tracks	1st Regime		2nd Regime
	Conductivity	Temperature	Conductivity	Temperature
Track A1	0.264 ± 0.021	0.136 ± 0.010	0.635 ± 0.002	0.335 ± 0.003
Track A2	0.207 ± 0.014	0.141 ± 0.010	0.508 ± 0.004	0.337 ± 0.007
Track B1	0.342 ± 0.020	0.210 ± 0.013	1.219 ± 0.005	1.278 ± 0.005
Track B2	0.339 ± 0.019	0.196 ± 0.013	1.238 ± 0.003	1.216 ± 0.002
Track B3	0.282 ± 0.013	0.173 ± 0.010	1.115 ± 0.003	0.931 ± 0.005
Track C1	0.524 ± 0.028	0.346 ± 0.023	1.113 ± 0.002	0.972 ± 0.005

Note: The values in red indicate superdiffusive behavior, while those in black show subdiffusion.

compares these values of γ with those obtained from the rest of the tracks. Most of the exponents appear to have a value of $0<\gamma <1$, which means that the behavior is subdiffusive. The values marked with a red color, however, show a superdiffusive state, $1<\gamma <2$. More specifically, in the case of conductivity, we have a subdiffusive state in about half the tracks until Δ = 5 s. At the same time, in the second regime, the behavior is superdiffusive (Tracks B1, B2, B3, C1). In Area A, the behaviors appear to be subdiffusive before and after the turning point (Tracks A1, A2). On the other hand, temperature seems to have the same behavior as conductivity in the first regime (subdiffusive state until Δ = 5 s). In the second regime, however, the behavior is superdiffusive in Tracks B1 and B2, while most tracks reveal a subdiffusive behavior before and after the turning point (Tracks A1, A2, B3, C1).

As for the time series, the behavior found in the graphs is consistent in both conductivity and temperature throughout all tracks, meaning that spikes appear at the same time in both parameters. However, the intensity of the spikes may differ. Figure 12 shows the time series for Track B3, while the time series for the rest of the tracks appear in Appendix B.

The moments ρ(q,Δ) in relation to the lag time Δ for Track B3 are presented in Figure 13. A total of 10 different values were used for q for Equation (3); these values are q = 0.25, 0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0. We first notice that conductivity and temperature seem to have different behaviors. Conductivity has a slow curve and an evident turning point in all graphs. The turning point is most prominent with higher values of q (q ≥ 2), whereas with lower, the plots show an almost flat line (parallel to the lag time). This is also true for temperature. The turning point for both parameters appears to be the same in the moments as is in the variance plots at Δ = 5 s. Fitting the moments for Δ > 5 s and observing the difference between those fits give us a clearer picture of how the graphs change as the order of the moments increases.

The structure function (Figure 14) has a convex shape and best fits Equation (5) for a = 1, corresponding to a Cauchy–Lorentz distribution. The overall field is non-conservative (H > 0), while C has a very small value close to 0, indicating that the mean field is slightly inhomogeneous.

Table 3. H and C parameters for conductivity and temperature as derived from the structure function.

	Conductivity		Temperature
a = 1	H	C	H	C
Track A1	0.200 ± 0.011	0.118 ± 0.021	0.219 ± 0.003	0.048 ± 0.005
Track A2	0.170 ± 0.002	0.079 ± 0.004	0.349 ± 0.002	0.135 ± 0.004
Track B1	0.538 ± 0.005	0.047 ± 0.010	0.211 ± 0.004	0.111 ± 0.009
Track B2	0.559 ± 0.023	0.286 ± 0.046	0.210 ± 0.006	0.065 ± 0.011
Track B3	0.480 ± 0.007	0.132 ± 0.014	0.417 ± 0.004	0.109 ± 0.007
Track C1	0.438 ± 0.005	0.129 ± 0.011	0.418 ± 0.006	0.021 ± 0.012

shows the values of H and C for all six tracks. These results show that volcanic activity inside each caldera has nearly the same multifractal behavior everywhere despite the difference in diffusivity.

4. Discussion

The exploration of the Santorini caldera spanned a duration of three days. The survey focused on the entire northern basin, and it included six tracks where the GMM was applicable (Tracks A1, A2, B1, B2, B3, C1). These tracks were considered the most reliable compared to the rest of the dataset since the vehicle remained very close to the seafloor for more than 30 min.

The variance, where both conductivity and temperature exhibit two distinct regimes, with a turning point at Δ = 5 s in all tracks, suggests that multiple underlying mechanisms govern the parameters’ behavior over time. In most cases, the behavior of conductivity in the first regime is subdiffusive, with a value of 0 < γ < 1. However, in the second regime, Tracks B1, B2, B3, and C1 show a superdiffusive state with a value of 1 < γ < 2. This suggests that the behavior of conductivity is not consistent across all tracks in the second regime and that there are variations in the way that the parameters are affected by the underlying mechanisms. The moments also appear to have two different regimes, with a turning point at Δ = 5 s. Interestingly, the turning point is more prominent in temperature with higher values of q, while lower orders of the moments show a flat line. This seems to suggest that the behavior is more random and less predictable at shorter time scales, and it would be interesting to see the behavior of temperature in a more extended time series.

Both parameters show dependence on different lag times. A unique fit for each moment and all lag times captures the dynamics exclusively for large times. This can be interpreted by assuming that the respective dynamics for each regime are different. Data in the second regime have been fitted with Equation (4) for every moment. The convex shape of the structure functions implies that each time series corresponds to multifractal processes.

For conductivity data in the first regime, which covers data points with a maximum lag time of 5 s, the consistent turning points indicate that ions are produced by volcanic activity at the same time. The subdiffusive behavior of Area A’s first and second regimes underlines that outgoing products share common properties. However, this is not the case with the tracks from Areas B and C, where the second regime appears to behave superdiffusive. This difference in behavior between the different areas studied is in complete agreement with the depth profiles, where tracks from Areas B and C show a more scattered profile in the two parameters (conductivity and temperature).

The subdiffusive behavior underlines the presence of a magnetic field, which results from the motion of the outgoing ions from the crust to the water through the vents. This makes the conductivity mean field non-conservative. On the other hand, superdiffusion exhibits persistent behavior that can be explained by a temporal increase in conductivity, likely followed by another increasing step. In the second regime, the mean inhomogeneity of the conductivity field is almost unchanged, the mean field remains non-conservative, and a Cauchy–Lorentz distribution describes its changes. As was somewhat expected for a wide caldera, no trapping effects are observed, suggesting that memory effects are much weaker in the formation of repeated patterns than in a very active field.

Temperature expresses energy dissipation, the process of converting a certain form of energy into thermal energy that cannot be fully converted back into the original form of energy. In our results, the variance in the different temperature regimes indicates a continuous energy flux into the seawater, resulting directly in a temperature rise since no value goes to 0—a continuous flow of energy points to a system that is not in thermodynamic equilibrium, which agrees with the slight anomalies observed in the CTD profiles.

The vent field in the Santorini caldera exhibits multifractal behavior, just like the previous fields studied, Kolumbo and Avyssos. Submarine volcano Kolumbo is located 7 km NE of Santorini [7,22,38] and is considered one of the most dangerous submarine volcanoes in the Mediterranean [47]. On the other hand, Avyssos is a submarine caldera at the easternmost part of the HVA, N of Nisyros volcano [10,48,49].

Comparing the results presented here with results from our previous studies, we find that Santorini has few similarities and differences with both Kolumbo and Nisyros. The comparison here is made between data from the POS510 mission (present work), the Avyssos caldera [34,35], and the data from the Kolumbo area during the 2010 unrest period and during the following year when the activity had subsided (rest) [33].

In the present work, the vent activity is relatively low, which is also the case in Avyssos [33,34]. In both these cases, this is confirmed by visual images captured by the cameras aboard the AUV and ROV. This contrasts with the findings from Kolumbo in 2010, where the activity was evident, and the images show ongoing plumes [38,50]. The turning point in the variance and the moments in Santorini is at Δ = 5 s, similar to that of Kolumbo in 2010, as opposed to Avyssos, where it is at Δ = 1 s. Interestingly, in Santorini, we have a general subdiffusive behavior in the first regime, which turns into a superdiffusive behavior at later lag times. On the other hand, the behavior in Avyssos starts as superdiffusive and then turns subdiffusive afterward.

The case of the submarine volcano Kolumbo in 2011 showed the existence of multiple time scales for conductivity while the temperature was stationary [33]. On the other hand, the results from Avyssos [26] showed that when the volcano is at rest, conductivity is constant, revealing a stationary process. Comparing these data to findings in Track B3 and the Avyssos caldera (results published in 2021; [27]), we see that the new results differ from all previous cases. In

Table 4. A comparison of the values of parameters H and C (Equation (4)) for data gathered during the 2010 unrest of Kolumbo, the 2011 rest period of Kolumbo, and the 2010 Nisyros data.

Conductivity
	Track B3	Track 7	Avyssos10	Kolumbo10	Kolumbo11
H	0.480	0.162	-	0.185	0.371
C	0.132	0.052	-	0.030	0.134
Temperature
	Track B3	Track 7	Avyssos10	Kolumbo10	Kolumbo11
H	0.417	0.137	0.601	0.051	-
C	0.109	0.041	0.070	0.189	-

, we present the values of H and C, side by side, from Track B3 in Santorini (present study), Track 7 in Avyssos [27], the data published in the 2017 Kolumbo study for the unrest period (Kolumbo10) and the rest period (Kolumbo11) [33], as well as the 2018 Nisyros study (Avyssos10) [26]. All data compared in the table were selected based on having similar periods, with similar distances crossed by the vehicles. In all cases, the vehicle was hovering above the seafloor, making the comparison easier. Note that all tracks presented in Table 4 had a period of 60 min, except for Avyssos10, which is 90 min. The missions in Kolumbo and Nisyros shared the same equipment (ROV) and survey methodology, while the missions in Santorini used different equipment (AUV).

Kolumbo is considered one of the most dangerous submarine volcanoes in the Mediterranean [47], and the results from Kolumbo10 confirm this. However, in Kolumbo11, the submarine volcano appears at rest [33]. Tracks B3 and 7 differ from the results in the other studies, making them less of a risk than Kolumbo10 but also more prominently active than Kolumbo11 and Avyssos10. During the resting period of Kolumbo [33], the values of temperature parameters are null, while Nisyros gave null conductivity values, which is not the case in Track B3 and Track 7.

Moreover, the structure function has a convex shape in all cases, though changes in Track B3 and Kolumbo10 are described by a Cauchy–Lorentz distribution. In contrast, a lognormal distribution describes changes in Track 7 and Avyssos10. The moments show a behavioral difference in temperature among the volcanoes, with the temperature in Avyssos remaining unchanged as the order of the moments rises. The difference in behavior could be attributed to the different magma chambers that lie beneath the volcanoes or the different geomorphologies. At the same time, the different equipment used in the mission should also be taken into consideration here. Even though this may seem somewhat speculative, the differences in the geomorphology of the studied chambers are known to have a direct effect on the vent activity, which is the main focus of this study. Even indirectly, the present results confirm the existing evidence from other types of measurements, establishing at the same time the GMM application as a useful approach for obtaining information where more expensive and complicated techniques are necessary. Finally, Kolumbo10 is the only case where actual particle trapping occurs, where outgoing products with large differences in lag times share common properties. This means that the combined effect of seawater and volcanic activity retains a memory. In all the other cases, the memory effect is weak.

5. Conclusions

In this study, we have demonstrated that the hydrothermal field of Santorini exhibits multifractal characteristics. These characteristics arise from two distinct underlying mechanisms. However, it is noteworthy that over time, the behavior of the hydrothermal field in the entire northern basin of the caldera becomes inconsistent. In certain areas, conductivity and temperature exhibit a subdiffusive behavior, while others display a superdiffusive behavior. Furthermore, Santorini’s behavior differs from that of its neighboring submarine volcano, Kolumbo, located NE of Santorini, and the submarine caldera, Avyssos, situated at the easternmost part of the HVA.