The Hydrothermal Vent Field at the Eastern Edge of the Hellenic Volcanic Arc: The Avyssos Caldera (Nisyros)

Abstract

Almost three-quarters of known volcanic activity on Earth occurs in underwater locations. The presence of active hydrothermal vent fields in such environments is a potential natural hazard for the environment, society, and economy. Despite its importance for risk assessment and risk mitigation, the monitoring of volcanic activity is impeded by the remoteness and the extreme conditions of many underwater volcanoes. The morphology and the activity of the submarine caldera, Avyssos, at the northern part of Nisyros volcano in the South Aegean Sea (Greece), were studied using a remotely operated underwater vehicle. The recorded time series of temperature and conductivity over the submarine volcano have been analyzed in terms of the Generalized Moments Method. This type of analysis can be used as an indicator for the state of activity of a submarine volcano. Here, we expand the work conducted for the first time in 2018. We present the findings of the geological exploration and the mathematical analysis, obtained from the data collected in October 2010. The temperature and conductivity time series show minor fluctuations in a rather stable environment. Based on these results, the impact of developing appropriate mechanisms and policies to avoid the associated natural hazard is expected to be important.

1. Introduction

The Hellenic Volcanic Arc (HVA) resulted from the northward subduction of the African plate underneath the active margin of the European plate [1,2,3,4]. The HVA extends from volcanic centers of the Saronic Gulf in the west (W) to the Kos–Nisyros Complex in the east (E) [5,6].

The HVA volcanoes have been especially active in the Late Pleistocene–Holocene and occur onshore [7,8,9] the peninsula of Methana and the islands of Poros, Milos, Santorini, Kos, and Nisyros [10], as well as offshore [6], with several submarine volcanoes detected in the Epidauros Basin in the W Saronikos Gulf [11,12,13] NE of Santorini, the submerged Kolumbo Volcano [14,15,16,17], and in the submarine area around Nisyros [18,19,20,21,22,23,24,25,26,27] (Figure 1).

Volcanism related to the modern South (S) Aegean volcanic arc in the Kos-Nisyros area began in Pliocene times [10,28]. Moreover, in Plio–Pleistocene times, lavas were emplaced on the centers of Pachia and Perigusa (a few km W of Nisyros) and in the lower succession of Nisyros [29,30,31]. In Quaternary times, major magmatic activity in the area led to the eruption of the Kos Plateau Tuff, which covers an area of more than 3000 km${}^2$. The Plateau is thought to have been produced about 160 ka from a volcanic source located between Kos and Nisyros, with its center close to the young volcanic island of Yali [30,31,32,33,34,35,36,37]. According to Pe-Piper et al. [38], the Kos Plateau Tuff was erupted from an original andesite stratocone, the remnants of which are present beneath Yali, Nisyros, and Pachia. The center of this eruption is not known with accuracy, but it is probably located in the submarine area, north of the Yali islet, a few kilometers NW of Nisyros [21,25].

The main onshore morphological features of the Kos-Nisyros-Tilos volcanic field are three zones of positive relief comprising Kos in the NW (843 m elevation in Mt. Dikeos), Nisyros (698 m in Mt. Prophitis Ilias) and surrounding islets in the middle, and Tilos (654 m, Prophitis Ilias) in the SE [23,24]. These zones divide the submarine area between Kos and Tilos into two basins with an average sea bottom depth of 600 m. Consequently, the topographic differences between the mountain ranges and the submarine basins are of the order of 1–1.5 km. This results in very steep slopes and a narrow zone of shallow depths (0–300 m) in the S coast of Kos Island. By contrast, the area of the islets around Nisyros is characterized by extended shallow water depths [21,26,39] (Figure 1).

Volcanism developed in the marine area around Nisyros with the formation of individual volcanic centers, each of them exhibiting a special geometry and evolutionary stage [19,21,23,24,26]:

(i)

The Nisyros Volcano, which is exclusively made of Quaternary volcanic rocks and has been characterized by a stratovolcanic cone, whose eruption created the Nisyros caldera. The largest present-day phreatic crater of Stephanos occurs within the Nisyros caldera along with other phreatic craters (such as Alexandros, Polyvotis, Phlegethon, and Achelous) that are still emitting fumaroles.

(ii)

The Yali volcano, which exhibits a partly submerged caldera. The two parts of Yali are dislocated by a postcaldera N–S fault. Western Yali is made of two successive pumice formations, whereas Eastern Yali is made of obsidian glass covered by a pumice formation (equivalent to the upper pumice of the Western Yali).

(iii)

The Strongyli volcanic cone, which starts at a −600 m depth (sea-bottom) and reaches +120 m at the top of the present day Strongyli Islet, where a 300 m-diameter volcanic crater is observed (pre-caldera stage).

(iv)

The Pergousa volcanic cone, which is made of stratovolcano type formations, with alternating lava flows and pumice layers (pre-caldera stage).

(v)

The volcanic domes of Pachia Islet and the submarine volcanic domes to the E of Kondeliousa Islet.

No volcanic activity is documented on the island for at least 25 ka after the formation of the volcanic domes of Prophitis Ilias; the only reported historical explosions are related to the formation of several phreatic craters inside the caldera. Violent earthquakes, gas detonations, steam blasts, and mudflows accompanied the most recent hydrothermal eruptions in 1871–1873 and 1887 AD (see the work in [10] and references therein). During this activity some people were injured and minor damages were caused to the infrastructures.

The crustal seismicity of Nisyros exhibits the basic properties of active volcanic regions in quiescent periods [40]. Nisyros’ seismicity has shown episodic unrest with local swarm-like earthquake sequences separated by months and years of quiescence. This is a characteristic which makes Nisyros comparable to other large calderas, such as Campi Flegrei (Naples, Italy), Rabaul (Papua, New Guinea), Yellowstone (WY, USA), and Long Valley (CA, USA) [41]. At the end of 1995, the Nisyros volcanic region entered an episodic unrest period with long-lasting, highly clustered earthquake activity, which seems to be a characteristic of the region rather than a precursor of a volcanic eruption [40].

From 1996 to 1998, the island experienced an episode of unrest that included a series of shallow earthquakes up to a M = 5.5 magnitude [40], a considerable temperature increase in the hydrothermal system [42], and ground deformation [43]. Seismic activity related to the Mandraki fault caused significant damage on the island along the W edge of the town of Mandraki (located in NW Nisyros) [26]. The fault extends N under the sea.

Despite its importance for risk assessment and mitigation, monitoring of the activity of submarine volcanoes is impeded by the remoteness and the extreme conditions of underwater volcanoes [44,45]. An efficient method to study those characteristics is to employ sensors aboard Remotely Operated underwater Vehicles (ROV), able to record CTD data (Conductivity, Temperature, Depth) in the water column (see, e.g., in [46]). Previous studies of the underwater features around Nisyros have focused on the geomorphology, using ROV to collect data, and have explored the volcanic field offshore, mainly in the Nisyros-Yali-Strongyli area [26].

Such data were collected for the first time in the underwater caldera Avyssos at the NE part of the volcanic islet of Strongyli, N of Nisyros volcano (Figure 2). The size of Avyssos caldera stretches for 3 km and 4 km along NW-SE and NE-SW directions, respectively, with a maximum depth of 680 m [23,24]. In addition, there is an underwater hill, about 1 km long at the center of Avyssos with a positive relief of about 60–70 m above a flat sea bottom. A dedicated mission conducted in 2010 focused on exploring the unknown volcano-tectonic characteristics of the Nisyros-Yali-Kos region and the associated volcanic field around Nisyros [26]. The main focuses of that mission included the observation of the previously unexplored hydrothermal field of Avyssos. In the present work, results from the full investigation are reported for the first time in terms of high-frequency CTD measurements using sensors mounted on a ROV (Figure 3).

A small portion of the full data set has been reported and analyzed by our group in [47] using a novel mathematical method, based on the Generalized Moments Method (GMM) [48]. The same modeling method had been successfully applied to identify the underlying stochastic processes driving the activity of the Kolumbo volcano in the Aegean sea [49]. This study yielded interesting aspects of the underlying dynamics of Kolumbo’s hydrothermal vent field during rest and unrest periods.

Due to the non-zero probability of eruption, and thus of its potentially devastating effects on the environment, society, and economy, the dynamic state of Avyssos is studied utilizing the data collected in October 2010, significantly extending our initial work in 2018, as reported in [47], to the full hydrothermal vent field in Avyssos.

2. Materials and Methods

2.1. Data

Recent multibeam mapping and seismic profiling have revealed several previously undetermined features near Nisyros and Yali, including a number of submarine volcanoes located NE of the Strongyli islet. The ROVs Argus and Hercules were deployed from E/V Nautilus in October 2010 to explore the slopes of the crater and the volcanic domes identified on the bathymetric map acquired from the multibeam data (Figure 2).

The CTD sensor was a calibrated SeaBird FastCAT 49Plus CTD probe mounted on ROV Hercules (Figure 3), which was tethered to ROV Argus receiving power from the ship and transferring data via a fiber-optic cable. During dives, the CTD sensor continuously acquired data at a rate of 16 Hz [50]. All data used in the mathematical analysis described in Section 2.2 were recorded at the deepest possible locations inside the Avyssos caldera. These CTD data sets from Nisyros/Avyssos are presented for the first time in literature.

The precise position of the vehicle was measured by highly calibrated pressure sensors on the vehicle. Argus is the tow sled that hung directly beneath the operating vessel and helped reducing ship motion effects from the hovering of Hercules over the seafloor. This vessel additionally provides working lights and “birds-eye” views for Hercules, which was used to gather the images from the seafloor during the mission (Figure 4, Figure 5 and Figure 6).

The main data collected were conductivity, temperature, and depth, as well as salinity and sound velocity. Vertical (depth) profiles and times series have been reconstructed from the raw data. In the following paragraphs, the focus of analysis is on temperature and conductivity, which are directly measured by the mounted sensors, as in the case of earlier applications [46,47,49].

2.2. The Generalized Moments Method

In recent works, a novel mathematical approach based on the GMM has been 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) [49], and the other recorded over a less active area (Nisyros caldera) [47]. Details of the method are presented in those references, where the reported results have also established the validity of the novel approach. In the next few paragraphs, a brief overview of the method is given.

Hydrothermal vent 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 (generally, but not necessarily) assumed as a manifestation of a stochastic process. A first estimate of the behavior of a stochastic process, $x\left(t\right)$, can be made in terms of its variance. For discrete data time series, by taking the time average, the variance reads

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

(1)

where M is the number of experiments repeated under the same conditions, to express the ensemble averaged over all records, whereas for $M=1$, the first summation is omitted. T is the length of the recorded trajectory, and $\Delta$ is the lag time, which is the elapsed time between any two measurements and plays the role of time for discrete data sets. $K_{\gamma }$ is a generalized coefficient expressed in proper units. It can serve as a tool for a first classification of the stochastic process based on the value of the exponent $\gamma$. A stochastic process is classified as sub–normal or antipersistent for $0<\gamma <1$, Brownian or normal under certain conditions for $\gamma =1$, super-normal or persistent for $1<\gamma <2$, ballistic for $\gamma =2$, and stationary for $\gamma =0$. Furthermore, $n=1,2,\dots ,T-\Delta$ and $\Delta =\tau ,2\tau ,\dots ,{\tau }_{max}\tau$ [51].

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

Considering a discrete time series in the form $x\left(n\right)$ with $n=1,2,\dots ,N$, where N is the total number of steps (measurements), and if the minimum lag time defined as the reciprocal of the sampling rate, f, or in other words the time elapsed between two consecutive measurements, is $\tau$, then the total length of the trajectory is $T=N\times \tau$. In brief, the “recipe” works as follows:

(1)

Construct time series with different lag times.

(2)

Estimate the statistical moments (Equation (2)).

(3)

Show that moments scale according to Equation (3).

In more detail, we create new time series, $y_n(\Delta )={(x{(n+\Delta )}^2-x{\left(n\right)}^2)}^{1/2}$, then we estimate the moments of the time series according to

$$\begin{array}{c}\hfill {\displaystyle {\rho }_{(q,\Delta )}=\frac{1}{T-\Delta }\sum _{n=1}^{T-\Delta }{\left(y_n(\Delta )\right)}^q}\end{array}$$

(2)

where even fractional values of the moment, q, are taken into account. Note that we use only positive values of the moments. The values of the moments in the range $0<q\le 2$ are responsible for the core of the probability density function (pdf), whereas higher moments, $q>2$, contribute to the tails of the pdf.

In general, the moments scale according to Equation (3), in which $z\left(q\right)$ is the structure function whose shape gives information on the stochastic mechanism(s) governing the motion.

$$\begin{array}{c}\hfill {\rho }_{(q,\Delta )}\approx {\Delta }^{z\left(q\right)}\end{array}$$

(3)

When the structure function takes the simple form $z\left(q\right)=Hq$, there is a direct relation between $\gamma$ and the co-factor of the structure function, namely, $\gamma =2H$, where H is the Hurst exponent [53]. 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. The linear dependency of $z\left(q\right)$ on q describes monofractal processes. Instead, if $z\left(q\right)$ has a convex shape as a function of the order of the moment q, then multiplicative processes drive the stochastic sequence. This also implies that H is a function of q. For some values of q, the exponents $H\left(q\right)=z\left(q\right)/q$ are associated with special features. For $q=1$, the value of $H\left(1\right)$ describes the scaling behavior of the absolute value of the increments and is identified with 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 [55]. Of all types of multifractals, universal multifractals are likely to be ubiquitous [56,57], and the structure function reads [54]

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

(4)

For $a=1$, Equation (4) takes the form $z\left(q\right)=Hq-Cqlog\left(q\right)$, and the distribution diverges from a Cauchy–Lorentz distribution. For $a=2$, Equation (4) reads $z\left(q\right)=Hq-C(q^2-q)$, and the distribution diverges from a lognormal distribution.

In the present study, the GMM is performed on the temperature and conductivity time series recorded over the hydrothermal vent field of Avyssos. As the ROV hovers over the extended hydrothermal vent field along various tracks, the time series present a variability in locations on the horizontal plane. 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 (temperature and conductivity). The ROV hovering motion (<1 m/s) perturbs the exact nature of the mathematical formalism, but it has been considered a safe approximation due to the widespread, overall homogeneous, vent activity along the studied tracks and the very slow horizontal translation. The results from the application of the method are reported in the next section.

3. Results

3.1. Mapping

CTD data were collected in Avyssos between October 10 and 12 2010 at a 16 Hz sampling rate. Raw data integrity was ensured first. Data sets representing measurements on deepest locations were selected, close to the seabed, where hydrothermal vent activity is possibly expected. The mission adapted to the geomorphology of the seafloor, a variance of 2–5 m, as per avoiding collision of the ROV with the submarine rocks and abnormal relief. In all cases, the ROV was kept at an almost fixed distance from the seafloor. Figure 7 marks the approximate area where the ROV collected those data sets, which were further grouped in tracks with the same maximum depth reached by the ROV within particular time slots and locations over the seabed (see

Table 1. List of the ROV tracks and their correspondence to data sampling time intervals.

Track Name	Surveys on 10 October 2010	Track Name	Surveys on 11 October 2010
Track 1	12:30–13:00	Track 7	01:00–02:00
Track 2	13:00–14:00	Track 8	02:00–03:00
Track 3	17:20–18:00	Track 9	03:00–04:00
Track 4	18:00–19:00
Track 5	19:00–19:48
Track 6	20:15–21:00

for the details of all tracks). The areas where the ROV passed over an underwater relief were excluded from the analysis.

The ROV Hercules explored the northern part of the caldera rim in the volcanic hill (SITE 1) from the south, and then moved to the other two volcanic hills (SITE 2, SITE 3), as shown in Figure 4, Figure 5 and Figure 6. High-resolution images from the seafloor display sediments and rocks without signs of apparent activity. Visual survey during the ROV movement in Avyssos (Figure 8) revealed a lack of significant fluid outflow from seabed edifices, bubble formation, black smoke, etc., whereas in the case of other active hydrothermal vent fields investigated by the same mission, such as the Kolumbo submarine volcano, intense CTD anomalies were found [46]. Nevertheless, despite being hardly detectable in direct optical observations, weak anomalies are evident in the CTD data, confirming an advantage of sensor-based measurements over “traditional” surveying methods.

Mathematical modeling was applied for data recorded in separate locations, within the time span of about an hour or less. At a first glance, anomalies were detected at a depth of 677–679 m (1–3 m above the seafloor) in almost all vertical profiles (depth profiles) for both conductivity (Figure 9) and temperature (Figure 10).

On 11 October 2010 between 1:00 a.m. and 2:00 a.m. (Track 7), the ROV crossed a short distance in the N part of the caldera, as shown in Figure 11. The conductivity and temperature vertical profiles obtained from these data show activity at approximately 679 m water depth, which is likely the deepest ever recorded activity in Avyssos to the best of our knowledge.

3.2. Time Series

CTD time series were reconstructed from the raw data for each track shown in Table 1. A special focus was given on the temperature and conductivity time series used to apply the GMM, as they are physical properties directly measured by the sensors and can provide useful information on the thermodynamic state of the system. In the case of Track 7, which is used here as a typical example of the analysis, both conductivity and temperature variances follow the same pattern in relation to the lag time ($\Delta$) in the time series, as shown in Figure 12. Both variances exhibit two discrete regimes characterized by different slopes, with a turning point at $\Delta \approx 1$ s. This value for the turning point is practically the same of the one reported in our earlier work where a single time series in a different location inside Avyssos was monitored [47]. However, it differs from the corresponding value found in the case of Kolumbo [49] ($\Delta =5.1$ s), where intense vent activity was observed.

Those distinct regimes in the plot are treated separately by fitting the variances of conductivity and temperature with a function of the form $W(\Delta )=b{\Delta }^{\gamma }$. In the first regime for conductivity, for $\Delta <1$ s, a value of $\gamma =0.954\pm 0.003$ is found, while in the second regime, $\gamma =0.211\pm 0.001$. For temperature, for $\Delta <1$ s, we obtain a value of $\gamma =1.172\pm 0.055$ in the first regime, whereas in the second $\gamma =0.171\pm 0.001$.

Similarly to Track 7, all other tracks appear to have two distinct regimes for temperature and conductivity with a clear turning point at $\Delta \approx 1$ s. The values of the exponents of $\Delta$ for each regime in the case of conductivity appear in

Table 2. Values of the exponent $\gamma$ for conductivity.

Track No	Regime 1	Regime 2
Track 1	$0.624\pm 0.005$	$0.290\pm 0.001$
Track 2	$1.028\pm 0.001$	$0.236\pm 0.001$
Track 3	$1.212\pm 0.001$	$0.069\pm 0.001$
Track 4	$1.364\pm 0.001$	$0.123\pm 0.001$
Track 5	$0.991\pm 0.001$	$0.100\pm 0.002$
Track 6	$1.363\pm 0.001$	$0.171\pm 0.001$
Track 7	$0.954\pm 0.003$	$0.211\pm 0.001$
Track 8	$1.250\pm 0.015$	$0.101\pm 0.001$
Track 9	$1.301\pm 0.001$	$0.292\pm 0.001$

, and in the case of temperature in

Table 3. Values of the exponent $\gamma$ for temperature.

Track No	Regime 1	Regime 2
Track 1	$1.349\pm 0.013$	$0.124\pm 0.001$
Track 2	$1.149\pm 0.069$	$0.190\pm 0.001$
Track 3	$1.471\pm 0.031$	$0.080\pm 0.001$
Track 4	$1.462\pm 0.039$	$0.110\pm 0.001$
Track 5	$1.321\pm 0.032$	$0.082\pm 0.002$
Track 6	$1.431\pm 0.041$	$0.154\pm 0.002$
Track 7	$1.172\pm 0.055$	$0.171\pm 0.001$
Track 8	$1.401\pm 0.082$	$0.092\pm 0.001$
Track 9	$1.264\pm 0.071$	$0.274\pm 0.001$

. In both cases, a variety in values of obtained exponents implies that the vent field fluid infuses and mixes with the surrounding water in different ways at different places. Despite the fact that the data cannot support it, the particular geomorphology of Avyssos (wide caldera), which allows water mobility in and out of the caldera, may be a major factor for the present observations. This is further supported by the CTD data from Kolumbo, which are affected by its unique morphology as discussed in [46].

Construct Time Series with Different Lag Times

Times series were constructed for both conductivity and temperature as shown in Figure 13. The similarity in behavior across these parameters is evident with spikes appearing at the exact same points in time. Moments for temperature and conductivity were found by using Equation (2). The results are illustrated in Figure 14. The values of q used for the estimation of the moments are 0.25, 0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, and 4.0. Only positive values of the moments were considered [58]. For both conductivity and temperature, all the moments showed the existence of two distinct regimes with a turning point at $\Delta =1$ s. For lag times greater than $\Delta =1$ s, all moments are almost parallel to the time axis in all tracks. This indicates that the parameters of the hydrothermal fluid equilibrate with the parameters of the surrounding water.

The moments of conductivity and temperature were fitted using Equation (3) for $\Delta >1$ s. The obtained exponents were used to find the form of the structure function, $z\left(q\right)$. The structure function poses a convex shape as a function of the order of q (see a typical example in Figure 15), which reveals that multiplicative processes driving the stochastic sequence. $z\left(q\right)$ was best fitted by Equation (4) for $a=2$, lognormal distribution. In

Table 4. Parameters H and C values obtained from Equation (4) for conductivity and temperature.

	Conductivity		Temperature
Track No	H	C	H	C
Track 1	$0.225\pm 0.023$	$0.098\pm 0.012$	$0.119\pm 0.014$	$0.042\pm 0.006$
Track 2	$0.169\pm 0.013$	$0.044\pm 0.006$	$0.165\pm 0.026$	$0.054\pm 0.011$
Track 3	$0.056\pm 0.005$	$0.016\pm 0.002$	$0.073\pm 0.009$	$0.022\pm 0.004$
Track 4	$0.105\pm 0.014$	$0.036\pm 0.006$	$0.099\pm 0.016$	$0.033\pm 0.007$
Track 5	$0.086\pm 0.009$	$0.023\pm 0.004$	$0.078\pm 0.012$	$0.019\pm 0.005$
Track 6	$0.118\pm 0.009$	$0.029\pm 0.004$	$0.115\pm 0.014$	$0.030\pm 0.006$
Track 7	$0.162\pm 0.004$	$0.052\pm 0.002$	$0.137\pm 0.012$	$0.041\pm 0.005$
Track 8	$0.077\pm 0.008$	$0.022\pm 0.004$	$0.084\pm 0.014$	$0.027\pm 0.006$
Track 9	$0.189\pm 0.007$	$0.041\pm 0.003$	$0.189\pm 0.013$	$0.046\pm 0.006$

we present the values of the parameters H and C obtained per track from Equation (4) for conductivity and temperature.

The GMM was also applied on a longer time series for a total of three hours during Tracks 7, 8, and 9, thus from 1:00 a.m. to 4:00 a.m. of 11 October 2010 (Figure 16). These three hours were chosen over the rest of the data sets because the ROV remained uninterrupted at the deepest location possible during this time interval. The activity demonstrated in the depth profiles in both conductivity and temperature for these total of three hours appears at about 678.5–679.5 m deep. The variance also displays two regimes, with a turning point at $\Delta =1$ s for conductivity (Figure 17) and temperature (Figure 18). Respectively, the exponents of $\Delta$ for conductivity are $\gamma =1.010\pm 0.001$ for the first regime, and $\gamma =0.213\pm 0.001$ for the second regime (

Table 5. Table of the $\gamma$ exponents for (a) conductivity and (b) temperature variance of Track 7, 8, 9, and the total of those.

(a) Conductivity
	Track 7	Track 8	Track 9	Total (7,8,9)
Regime 1	$0.954\pm 0.003$	$1.250\pm 0.015$	$1.301\pm 0.001$	$1.010\pm 0.001$
Regime 2	$0.211\pm 0.001$	$0.101\pm 0.001$	$0.292\pm 0.001$	$0.213\pm 0.001$
(b) Temperature
	Track 7	Track 8	Track 9	Total (7, 8, 9)
Regime 1	$1.172\pm 0.055$	$1.401\pm 0.082$	$1.264\pm 0.071$	$1.202\pm 0.071$
Regime 2	$0.171\pm 0.001$	$0.092\pm 0.001$	$0.274\pm 0.001$	$0.171\pm 0.001$

). Meanwhile, for temperature, the first regime offers an exponent $\gamma =1.202\pm 0.071$ and the second $\gamma =0.171\pm 0.001$ (Table 5). As per the moments, they appear to have two regimes as well, with turning points at $\Delta =1$ s for conductivity (Figure 19) and temperature (Figure 20). The structure function for the total of three hours appears to have a similar behavior with Track 7 in both conductivity and temperature, as seen from the parameters H and C in

Table 6. Values of parameters H and C obtained from Equation (4) for conductivity and temperature.

	Conductivity		Temperature
	$\mathbf{H}$	$\mathbf{C}$	$\mathbf{H}$	$\mathbf{C}$
Track 7	$0.162\pm 0.004$	$0.052\pm 0.002$	$0.137\pm 0.012$	$0.041\pm 0.005$
Track 8	$0.077\pm 0.008$	$0.022\pm 0.004$	$0.084\pm 0.014$	$0.027\pm 0.006$
Track 9	$0.189\pm 0.007$	$0.041\pm 0.003$	$0.189\pm 0.013$	$0.046\pm 0.006$
Tracks 7, 8, 9	$0.156\pm 0.004$	$0.049\pm 0.002$	$0.138\pm 0.013$	$0.042\pm 0.006$

for conductivity and temperature. These parameters were obtained from Equation (4).

4. Discussion

In the present work, vertical CTD profiles are presented for the first time with a full range of temperature and conductivity time series from the submarine caldera of Nisyros. The hydrothermal vent activity is relatively low and cannot be seen in the visual images captured by the cameras on board the ROV, as opposed to analogous findings from the Kolumbo submarine volcano that same year, where the activity was evident and the cameras captured ongoing plumes [15,17]. In 2010, Kolumbo underwent a period of unrest (Kolumbo2010), while in 2011, it was at rest (Kolumbo2011) [49]. These findings make it clear that hydrothermal venting occurs in different ways across periods of different activity. Understanding the data patterns as well as the underlying mechanisms stresses the importance of monitoring submarine volcanoes.

Comparing the data between the periods of rest and unrest within Kolumbo with our findings in Track 7 (Avyssos2010), we see that our results differ from both cases in Kolumbo. In

Table 7. 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			Temperature
	Kolumbo 2010	Kolumbo2011	Avyssos2010	Kolumbo2010	Kolumbo2011	Avyssos2010
H	$0.185\pm 0.002$	$0.371\pm 0.002$	$0.162\pm 0.004$	$0.051\pm 0.001$	0	$0.137\pm 0.012$
C	$0.030\pm 0.004$	$0.134\pm 0.00$	$0.052\pm 0.002$	$0.189\pm 0.007$	0	$0.041\pm 0.005$

, we compare the values of H and C from Avyssos2010 to the data published in Bakalis et al., 2017 (Kolumbo2010, Kolumbo2011). Note that all tracks presented in the table had a time-span of an hour and the missions in Kolumbo and Nisyros shared the exact same equipment and survey methodology; as such, they can be confidently used to distinguish different states of unrest in their hydrothermal vent fields. Kolumbo is considered to be one of the most dangerous submarine volcanoes in the Mediterranean [59]. Our results in Avyssos2010 differ from the results in Kolumbo in a way that makes it of lesser risk than Kolumbo2010, but also more prominently active than Kolumbo2011. During the resting period of Kolumbo, the values of temperature parameters are null, which is not the case in Avyssos.

On the other hand, the variance vs. lag time graph appears to have three regimes for conductivity in the data from both Kolumbo2010 and Kolumbo2011, while in Avyssos2010, we only have two regimes. Temperature in Kolumbo2010 has three regimes as well, while in Kolumbo2011, temperature variance is parallel to the time axis, which results in a $\gamma$ exponent equivalent to zero. In Avyssos2010, no parameters are null. The morphological structure of the two volcanoes are different, which might result in the difference in behavior after a time period, thus creating three regimes for Kolumbo and two regimes for Avyssos. The two submarine volcanoes differ in shape and size as Kolumbo has a circular cone with abrupt steep slopes, while Avyssos is a 3 km wide volcano [6].These geomorphological characteristics make Kolumbo a closed thermodynamic system, and Nisyros an open one. However, the thermodynamic state cannot account for the null values of temperature parameters during the resting period of Kolumbo, considering the homogeneity of hydrothermal fluid outflow.

The statistical moments in Kolumbo2010 appear to have two distinct regimes in the moments vs. lag time graph, with a turning point at $\Delta =5$ s. In Kolumbo2011 however, the turning point is at $\Delta =1$ s, the same as in Avyssos2010. This can be interpreted as an indicator for similar dynamics to Kolumbo2011 despite some structural differences in the variance function.

Last, the structure function has a convex shape in Kolumbo2010 and in Avyssos2010. As for Kolumbo 2011, the structure function of conductivity appears to have a similar behavior, whereas for temperature, $z\left(q\right)$ is 0 for all values of q.

Regarding the time series in Avyssos, the first results from the survey [47] originate from a simple data set gathered in 10 October 2010, during a total recording time interval of 1 h and 30 min. That time series preceded the data in the present work and were recorded over a specific location while the ROV was slowly hovering. According to those results, when the volcano is at rest, conductivity seems to obey a stationary process, and the crater of the volcano is in direct contact with the overlying water environment, thus defining an open system. In this research, the same methodology was used to extend the analysis of the entire data set gathered from the deepest parts of the caldera (10 and 11 October), covering approximately the full area under investigation.

We consider anomalies in the time series to be spikes that deviate from the baseline. The baseline in the case of Track 7 (Figure 13) would be a line almost parallel to the time axis at about 4.756 S/m for conductivity and at approximately 14.91 ${}^{\circ }$C for temperature. Both parameters would be otherwise stationary and would evolve in time in a manner that would prove irrelevant to a baseline, where there no activity (as is the case for Kolumbo2011). We use these values as a reference to evaluate the anomalies. For conductivity, the most prominent spike diverges at $0.9\%$ from the reference baseline, while for temperature, the respective percentage is $10.5\%$. The above observations are similar for all tracks presented in here (such as Track No 4, etc.). In Kolumbo2010, conductivity spikes scale at a maximum of $130\%$ (spike = 5.5 mS/m, baseline = 4.2 mS/m), while temperature anomalies spike at $480\%$ (spike = 21 ${}^{\circ }$C, baseline = 16.2 ${}^{\circ }$C), which shows the near-explosive state of the volcano. The case of Avyssos2010 clearly differs to the case of Kolumbo2010. Conductivity and temperature parameters in Avyssos2010 spike at a much smaller scale in terms of absolute values, as well as percentage, compared to Kolumbo2010. In Kolumbo2011, spikes appear to be random and constant with no particular baseline to reference any anomalies (see Figure 1c,d in [49]).

Activity is evident in both time series and vertical profiles, and seems to be consistent among conductivity and temperature, which are the main points of focus. The time scale peaks appear at the exact same point in time (Figure 13), whereas the vertical profiles appear at a similar depth for both the parameters (Figure 9 and Figure 10). In the case of the vertical profiles in Track 7, we notice spikes that skew at $1.6\%$ for conductivity and $11.7\%$ for temperature, respectively.

An additional piece of evidence is the comparison of the present vertical CTD profiles with reference CTD data. As such, a vertical profile between 20 and 250 m recorded with the exact same methodology and instruments during the same mission (but earlier time) in an open sea location, where no hydrothermal activity exists, has been considered and is shown in Figure 21a. The inset shows the full CTD profile, while the main panel shows a magnified region with the same x- and y-axis ranges as with Figure 21b which shows Track No 4 for comparison (see also the caption of the figure for further comments). From this comparison, it becomes evident that the vertical CTD profiles in Avyssos differ completely from respective (“reference”) CTD profiles in open sea, where effects such as mixing, thermal layers, turbulence, etc. can affect their shape. In this case, no abrupt changes are observed, rather step-like changes occur at certain depths. The situation can be easily generalized by comparing the open-sea vertical profile with any one shown in Figure 10.

Despite the activity of the volcanic field in Avyssos being weak and almost homogeneous in its entirety, as shown in the overall ROV surveys, the values of the exponents $\gamma$ are scattered for both conductivity and temperature in the first regime (Table 2 and Table 3). On the one hand, conductivity shows a normal to a slightly super-normal behavior, with the only exception in Track 1, where points are either uncorrelated or display persistent trends during reduced time intervals. On the other hand, the scaling of all tracks for the temperature indicates persistence. In the second regime, both fields present a strong anti-persistent motion, which mirrors on their slow changes.

These values reflect the existence of an active vent field, the behavior of which is not persistent. In addition, the value of the Hurst exponent indicates some delay between two similar value pairs in a time series. This situation reveals the existence of weak, homogeneous sources of hydrothermal fluid outflows, which are found throughout the entire examined area.

The exponents of $\Delta$ appear to have similar behaviors across various locations inside the caldera, as well as among conductivity and temperature. The similarity in behavior across parameters is also verified in their time series and their profiles. In Appendix A, the full list of maps, CTD profiles, and CTD time series for all Tracks mentioned in Table 1 are presented for Avyssos.

The analysis of the recorded time series of temperature and conductivity just above the Avyssos submarine volcano shows that both parameters seem to present an initial multiplicative behavior reflecting the mixing of at least two different random processes, which then turn to be quasi-stationary.

In the case of non-persistent hydrothermal activity, we expect to find dissimilar results when analyzing data from a longer period of time while the ROV remained at the deepest locations possible. As the analysis was performed while the ROV was hovering and slowly moving over a wider area at the same time, a direct comparison with the previous work in Nisyros [47] cannot be made. This seems to be distinctly reflected on the behavior of the conductivity time series, which have been found to be completely different from the single time series reported in [47]. The conductivity time series reported here present similarities with the temperature time series, which is not the case in the earlier work reporting results over a specific location in the vent field.

Comparing the results of the hourly intervals to the total of three hours (Tracks 7, 8, and 9), we find differences in the behavior of the moments ${\rho }_{(q,\Delta )}$. These differences are possibly attributed to a number of factors. One factor may be the particular geomorphological features of the caldera along these tracks, as the caldera bottom is not smooth. In addition to that, the ROV hovered over the seafloor in a slow motion, instead of being absolutely still. Finally, the weak variations in the underlying hydrothermal activity across the caldera can also potentially lead to these differences.

Although generally considered inactive, the island of Nisyros exhibits the basic properties of active volcanic regions in quiescent periods [40]. From time to time, a period that may range from months to years [40], Nisyros enters a period of unrest which leads to seismic activity, significant temperature increase in the hydrothermal system, gas detonations, steam blasts, mudflows and ground deformation [42,43,60,61,62,63,64,65]. The possibility of Nisyros entering the next periodic unrest in a few years cannot be fully excluded as suggested by the existence of weak hydrothermal sources found within the caldera.

Moreover, the present work confirms that the application of the GMM in geosciences seems to find a clear-cut way, at least in terms of analyzing CTD time series, to investigate the underlying mechanisms governing the activity of hydrothermal vent fields in submarine volcanic areas. As the GMM has been used in econometrics, commercial packages offering GMM analysis tools in geosciences are still not widely available, to the best of our knowledge. Based on our experience working with GMM to analyze CTD time series, we argue that such computer codes could significantly reduce the time needed to reach conclusions, something particularly useful for cases where the risk of natural hazards should be assessed and mitigated promptly. It is clear that monitoring natural hazards associated to submarine hydrothermal vent fields with state-of-the-art instrumentation is the spear to mitigate those risks. However, the continuous increase of computing power facilitates the introduction of advanced tools, such as Artificial Intelligence and Machine Learning, and can offer the chance to statistical methods, such as the GMM, to play an enhanced role as a complementary approach. An additional advantage seen in the present research is that the combination of a mathematical model (GMM) with fast, state-of-the-art sensors (ROV-based CTD) that collect data can provide a framework for understanding dynamic processes in the marine environment that could completely remain unknown by using traditional techniques of exploration (such as visual inspection or use of samplers).