Analysis of Volcanic Thermohaline Fluctuations of Tagoro Submarine Volcano (El Hierro Island, Canary Islands, Spain)

: Temperature and conductivity ﬂuctuations caused by the hydrothermal emissions released during the degasiﬁcation stage of the Tagoro submarine volcano (Canary Islands, Spain) have been analysed as a robust proxy for characterising and forecasting the activity of the system. A total of 21 conductivity-temperature-depth time series were gathered on a regular high-resolution grid over the main crater of Tagoro volcano. Temperature and conductivity time series, as manifestations of stochastic events, were investigated in terms of variance and analysed by the Generalised Moments Method (GMM). GMM provides the statistical moments, the structure functions of a process whose shape is an indicator of the underlying stochastic mechanisms and the state of activity of the submarine volcano. Our ﬁndings conﬁrm an active hydrothermal process in the submarine volcano with a sub-normal behaviour resulting from anti-persistent ﬂuctuations in time. Its hydrothermal emissions are classiﬁed as multifractal processes whose structure functions present a crossover between two time scales. In the shorter time scale, ﬁndings point to the multiplicative action of two random processes, hydrothermal vents, which carries those ﬂuctuations driving the circulation over the crater, and the overlying aquatic environment. Given that both temperature and conductivity ﬂuctuations are nonstationary, Tagoro submarine volcano can be characterised as an open system exchanging energy to its surroundings.


Introduction
Hydrothermal vents, found in submarine volcanic islands, release important emissions of hot fluids that rise in a buoyant turbulent plume over the source [1]. Those hydrothermal emissions range from intense high-temperature plumes to diffuse low temperature diluted discharges emanating directly from the seafloor [2]. The fluids mix with the surrounding seawater and rise to a level where their density matches the water outside the plume being able to alter the temperature and salinity fields specific to the vertical seawater structure [3]. These vertical motions generate a horizontal flow, resulting in an anticyclonic vortex of sufficient magnitude to trap water, minerals, tracers, and organisms in a local recirculation [4]. As a result, temperature and conductivity fields may fluctuate, influencing the local circulation. At the same time, fluctuations are the outcome of random perturbations acting on a field, noise sources caused by the short-term circulation over the volcano [5]. Decoding fluctuations with stochastic analysis tools can shed light and benefit understanding the complexity of hydrothermal emissions [6][7][8].
In the Mediterranean, a newly developed mathematical technique was applied to the time series of temperature and conductivity fluctuations, recorded by remotely operated vehicles (ROV), to understand the nature of the stochastic processes that drive the activity of volcanoes [6,7]. These studies suggest that temperature fluctuations can be used as an indicator to determine whether the volcano operates as an open or closed system, i.e., if it exchanges energy with the surroundings or not. Open stands for nonstationary temperature time series and closed for stationary ones. On the other hand, conductivity, the ability of a water sample to conduct an electrical current, fluctuations can provide information about its current state of activity. During unrest (active) degassing periods, conductivity and temperature follow a universal multifractal behaviour, i.e., a multiplicative action of random processes, with continuous energy dissipation. Otherwise, conductivity behaves as a universal multifractal but with a stationary pattern [6,7].
El Hierro is the youngest and southwesternmost island of the Canary Archipelago characterised by a three-armed rift system and originated by hotspot volcanism [9,10] ( Figure 1). On 10 October 2011, a submarine volcano eruption took place 1.8 km south of the coast of El Hierro Island, forming the shallow submarine volcano Tagoro, located at 27 • 37.1160 N, 017 • 59.4660 W [11]. When the eruption finished in March 2012, hydrothermal manifestations were characterised by the release of hot fluids, gases, inorganic nutrients and metals as a result of the new degasification phase of the post-eruptive stage of the Tagoro submarine volcano [12]. This degasification process is exclusively associated with hydrothermalism linked with shallow submarine volcanoes, resulting in water properties changes [13,14].
In the Mediterranean, a newly developed mathematical technique was applied to the time series of temperature and conductivity fluctuations, recorded by remotely operated vehicles (ROV), to understand the nature of the stochastic processes that drive the activity of volcanoes [6,7]. These studies suggest that temperature fluctuations can be used as an indicator to determine whether the volcano operates as an open or closed system, i.e., if it exchanges energy with the surroundings or not. Open stands for nonstationary temperature time series and closed for stationary ones. On the other hand, conductivity, the ability of a water sample to conduct an electrical current, fluctuations can provide information about its current state of activity. During unrest (active) degassing periods, conductivity and temperature follow a universal multifractal behaviour, i.e., a multiplicative action of random processes, with continuous energy dissipation. Otherwise, conductivity behaves as a universal multifractal but with a stationary pattern [6,7].
El Hierro is the youngest and southwesternmost island of the Canary Archipelago characterised by a three-armed rift system and originated by hotspot volcanism [9,10] ( Figure 1). On 10 October 2011, a submarine volcano eruption took place 1.8 km south of the coast of El Hierro Island, forming the shallow submarine volcano Tagoro, located at 27° 37.1160′ N, 017° 59.4660′ W [11]. When the eruption finished in March 2012, hydrothermal manifestations were characterised by the release of hot fluids, gases, inorganic nutrients and metals as a result of the new degasification phase of the post-eruptive stage of the Tagoro submarine volcano [12]. This degasification process is exclusively associated with hydrothermalism linked with shallow submarine volcanoes, resulting in water properties changes [13,14].  Since 2011, the Tagoro submarine volcano has been deeply monitored by more than 28 different multidisciplinary cruises in order to understand submarine volcanic processes better and how the volcano evolves in terms of different physical-chemical seawater parameters and the relationship of these with the surrounding marine ecosystem [8,11,12,[15][16][17][18][19][20][21][22][23][24][25][26][27].
However, there is a big gap in knowledge in terms of understanding how fluctuations in temperature and conductivity can provide information regarding the hydrothermal system and its state of activity.
The main goal of this manuscript is, using stochastic calculus tools, to assess and classify the state of activity (rest/unrest) and characterise the nature of stochastic mechanisms pertaining to around and over the main crater of the Tagoro submarine volcano in its degassing stage. This objective will be evaluated by applying the mathematical methodology proposed in [6] and using the temperature and conductivity measurements collected by a CTD sensor (Conductivity-Temperature-Depth) mounted on a rosette over the Tagoro submarine volcano. In addition, this will be the first time that a cost-effective oceanographical instrument (rosette) is implemented to register statically hydrothermal emissions in the Atlantic Ocean instead of the ROV, widely used in the Mediterranean Sea [6,26,28]. Although ROVs are a valuable tool for collecting in situ samples and recording HD images, compared with the rosette sampler, with a significantly higher measurement frequency, they are much more limited.

Materials and Methods
In October 2016, the multidisciplinary cruise VULCANO-II-1016 on board of R/V Ángeles Alvariño was carried out. One of the survey objectives consisted of deploying up to 21 oceanographic stations with a CTD sensor mounted on the rosette sampler in a regular 10.5-m high-resolution grid above the main crater, main crest and surroundings ( Figure 1d).
Each station records a unique time series of conductivity and temperature fluctuations, referred to here with the E c and θ parameters, between a time length from 9 to 13 min, with the exceptions of CTDs 41, 43 and 47, as a result of discarding some periods of fluctuations because the rosette was not static one meter above the seafloor (Table 1). Thus, each CTD time series was sampled with a frequency of 24 Hz or 0.0417 s distance between two consecutive measurements maintaining the rosette sampler, through the altimeter sensor, at approximately 1 m above the seabed. At the same time, the ship remained stationary on the surface through its Dynamic Positioning System. Data were acquired using an SBE 911-plus CTD equipped with dual temperature and conductivity sensors, calibrated at the SeaBird laboratory before and after the cruise, with accuracies of 0.001 • C and 0.0003 S/m, respectively. Table 1. The location, duration and date of the 21 CTD time series recorded during the multidisciplinary cruise VULCANO-II-1016. The initial and final times correspond to the entire time in which the oceanographic rosette sampler remained in the water. The duration at the seabed (in minutes and seconds) corresponds to the length of the CTD time series fluctuations without the descending and ascending track. All the minutes until the oceanographic rosette reach the volcano edifice and remain steady one meter above the seafloor were discarded, resulting in significatively shorter time series for CTDs 41, 43 and 47. The time series of temperature (θ) and conductivity (E c ) recorded at each location inside the grid are subject to fluctuations and have been analysed as the result of stochastic processes. In this sense, an initial evaluation of the dynamics of the hydrothermal emissions fluctuations driven by the activity of the volcano, x(t), can be made by examining the time dependence of its variance, W(t) = x 2 (t) − x (t) 2 . If the process, x(t), corresponds to uncorrelated random events (normal distribution), then the variance of x(t) grows linearly in time (Gaussian behaviour). Any departure from linearity, i.e., anomalous behaviour, might indicate the existence of anti-correlated/correlated events and/or systems where the environment fluctuates at similar time scales as the random variable x(t) [29]. In the anomalous regime, the behaviour is sub-normal/super-normal for growth slower/faster than linear, and the variance reads [30]:

Nº CTD
where K γ is a generalised coefficient expressed in proper units, e.g., for temperature the units are • C 2 s −γ , and Γ () is the Gamma function. According to the exponent γ, the stochastic processes can be classified as sub-normal (anti-persistent process) for values from 0 < γ < 1, Brownian or normal for γ = 1, super-normal (persistent process) for 1 < γ < 2, ballistic for γ = 2, and stationary for γ = 0. It is worth mentioning that Equation (1) does not apply to instantaneous changes of enormous amplitude (Lévy flights) where the first moment diverges. For discrete time series, Equation (1) is expressed as time average and reads: where T = N τ is the total time, N is the total number of measurements, and τ = 0.0417 s is the reciprocal of the maximum sampling rate. M represents the number of time series monitored under the same conditions, the first summation being omitted when M = 1. The parameter ∆ is the lag time and plays the role of the actual time in the analysis; it takes values in the range τ ≤ ∆ ≤ T/10. Temperature or conductivity time series are represented by x i (n) with n = 1,2, 3, . . . , N and i = θ or E c . Equation (2) provides an adequate description of when the process is normal, γ = 1. Any departure from normality requires the estimation of more moments, including fractional ones. We apply the Generalised Moments Methods (GMM) following these required steps: • Construction of N/10 new time series each of length T − ∆. The new time series, y n (∆), contains the absolute change between two values of the original series that are set apart by ∆, so that y n (∆) = |x(nτ + ∆) − x(nτ)| for n = 1, 2, . . . , ((T − ∆)/τ) and for ∆ = mτ with m = 1, 2, ..., N/10 where τ and Nτ/10 define the minimum and maximum time lags. • Estimation of the statistical moments, order of q, for each one of the new time series, y n (∆), is carried out according to: where only positive moments, q, are taken into account. Moments between 0 < q ≤ 2 are responsible for the core of the probability density function (pdf), while moments q > 2 contribute to the tails of the pdf (Bakalis et al., 2017 and the references therein).

•
The moments will scale according to: where z(q) is the structure function and its shape gives information about the stochastic mechanisms that govern the process. When the structure function fits a linear form, that is, z(q) = Hq, then there is a direct relation between the scaling exponent γ of Equation (1) and the Hurst exponent, H, γ = 2H, being the process characterised as monofractal. The parameter H is the mean fluctuating scaling exponent. For 0 < H < 0.5 the process is anti-persistent (sub-normal), normal for H = 0.5, and super-normal or persistent for 0.5 < H < 1. If z(q) presents any convex shape, then the process is multifractal. Among multifractals, universal multifractals are ubiquitous, and their structure function reads [31]: where the Lévy index, a, 0 ≤ a ≤ 2, indicates the class to which the probability distribution belongs and shows how fast the inhomogeneity rises. Inhomogeneity refers to the degree of mixing of various stochastic mechanisms shaping the final process. For a = 0, z(q) = Hq, the process is monofractal, and the field homogeneous. For a = 1, z(q) = Hq − Cqlog(q), the process draws steps from a Cauchy-Lorentz distribution. For a = 2, z(q) = Hq − C q 2 − q , the process is the multiplication of two random processes and is described by a log-normal or Kolmogorov distribution. For 1 < a < 2, the multifractal character is the multiplicative result of more than two random processes. Finally, C is the intermittency parameter that measures the mean inhomogeneity of the field and always takes positive values.

Results
Implementing a high-resolution grid of temporal series around and over the main crater reveals temperature and conductivity fluctuations in the whole domain, suggesting that the fields are disturbed by emissions from the volcano and may be far from a thermodynamic equilibrium. Using non-affected vertical profiles of θ ( • C) and E c (mS/cm) measured simultaneously on the submarine volcano surroundings (see Figure 1c), we obtained the average temperature and conductivity reference profiles and their standard deviation down to 127 m (shown in Figure 2). The high variations of both properties during the whole time series over the seabed, concerning the reference profile's standard deviation, confirm the significance of the fluctuations due to the continuous release of hydrothermal emissions. These fluctuations as a function of time and their correlations for the randomly selected CTD30, one of the 21 time series collected with a time interval of 11 min and 21 s, are illustrated in detail in Figure 3. The variances of temperature and conductivity have been obtained as a function of lag time for all CTD time series. Figure 4 shows the variance of both temperature ( Figure  4a) and conductivity (Figure 4b) of the CTD30 sampling point. Both properties present a quite similar scaling with a crossover point at Δ ≈ 3.5 s, defined as the time where the scaling exponent changes. The crossover point divides the time series into two time scales or regimes: the shorter time scale corresponds to the first regime and small lag times, and the longer time scale considers all the lag times and therefore, we will refer to it as the whole regime. We fitted each curve with Equation (1), W(Δ) = bΔ γ , obtaining the early (small lag times) and the long-time (all lag times) behaviour for θ and Ec. For small lag times, 0.042 < Δ ≤ 3.5 s, temperature and conductivity pose similar weak sub-normal behaviour, with scaling exponents 0.84 and 0.85, respectively. For the long-time behaviour, the whole lag times range has been considered. Again, both θ and Ec scale similarly, with values of exponents γ = 0.75 and γ = 0.77, respectively. For both small and large times, temperature and conductivity follow sub-normal behaviour, anti-persistent variations, which become stronger for longer times and lesser when scaling exponents. The absence of stationary behaviour (γ = 0) in temperature confirms the presence of a mechanism that causes the alterations (see also below), suggesting the existence of a remarkable hydrothermal vent field, consistent with what has been described so far in the area [8,21]. The variances of temperature and conductivity have been obtained as a function of lag time for all CTD time series. Figure 4 shows the variance of both temperature ( Figure 4a) and conductivity (Figure 4b) of the CTD30 sampling point. Both properties present a quite similar scaling with a crossover point at ∆ ≈ 3.5 s, defined as the time where the scaling exponent changes. The crossover point divides the time series into two time scales or regimes: the shorter time scale corresponds to the first regime and small lag times, and the longer time scale considers all the lag times and therefore, we will refer to it as the whole regime. We fitted each curve with Equation (1), W(∆) = b∆ γ , obtaining the early (small lag times) and the long-time (all lag times) behaviour for θ and E c . For small lag times, 0.042 < ∆ ≤ 3.5 s, temperature and conductivity pose similar weak sub-normal behaviour, with scaling exponents 0.84 and 0.85, respectively. For the long-time behaviour, the whole lag times range has been considered. Again, both θ and E c scale similarly, with values of exponents γ = 0.75 and γ = 0.77, respectively. For both small and large times, temperature and conductivity follow sub-normal behaviour, anti-persistent variations, which become stronger for longer times and lesser when scaling exponents. The absence of stationary behaviour (γ = 0) in temperature confirms the presence of a mechanism that causes the alterations (see also below), suggesting the existence of a remarkable hydrothermal vent field, consistent with what has been described so far in the area [8,21]. The similarity between both variable fluctuation patterns as well as the occasional distinct peaks is remarkable, and (c) their normalised cross-correlation coefficient rθ,Εc showing a strong correlation for short times (t < 3.5 s), which remains large even for long times. (d) The differentiation with respect to time of the temperature time series, Δθ, (e) the differentiation with respect to time of the conductivity time series, ΔEc, and (f) the normalised cross-correlation coefficient rΔθ,ΔEc of the differentiation with respect to time of temperature and conductivity time series going rapidly to zero, at the crossing point 0.6 s. In grey is indicated the 95% confidence bounds, where the time series are uncorrelated.
The rest of the CTD time-series exhibit crossover points whose positional values range from 1.8 to 5.04 s. Their variances for the whole regime are consistent with the CTD30 sub-normal behaviour, always with scaling exponent in the range 0 < γ < 1 ( Table  2). In the short time scale (first regime), although most of the sampling points show a similar pattern in terms of variance and can be defined as sub-normal, the temperature variances for CTDs 37, 38, and 50 and conductivity variances for CTDs 32, 38, 40, 48, and 50 suggest in those sampled locations processes that can be characterised as Brownian (γ ≈ 1), since 0.95 ≤ γ. The latter highlights the necessity of a closer investigation of the analysis of the recorded data which extends beyond the standard variance. The rest of the CTD time-series exhibit crossover points whose positional values range from 1.8 to 5.04 s. Their variances for the whole regime are consistent with the CTD30 sub-normal behaviour, always with scaling exponent in the range 0 < γ < 1 ( Table 2). In the short time scale (first regime), although most of the sampling points show a similar pattern in terms of variance and can be defined as sub-normal, the temperature variances for CTDs 37, 38, and 50 and conductivity variances for CTDs 32, 38, 40, 48, and 50 suggest in those sampled locations processes that can be characterised as Brownian (γ ≈ 1), since 0.95 ≤ γ. The latter highlights the necessity of a closer investigation of the analysis of the recorded data which extends beyond the standard variance.
Geosciences 2021, 11, x FOR PEER REVIEW 8 of 16 Figure 4. Logarithmic plots of the variance for (a) temperature, θ, (orange dots) and (b) conductivity, Ec, (dark green dots) with respect to the lag time for the CTD30 time series. The variance for both temperature and conductivity has been fitted by Equation (1). In each figure, the short black solid line represents the best fit for short times (first regime), and the long one the best fit for long times (whole regime). Notice that Δ (s) refers to the unit of time in seconds. Table 2. γ, α, H, and C parameters of the temperature and conductivity for both first and whole regimes (short and long times) and for each one of the recorded CTD time series. Notice that for short times the alpha index is 2. Generalised moments for both temperature and conductivity measurements up to the fourth order were determined by using Equation (3) and the moments for the CTD30 time series are illustrated in Figure 5. A vertical solid line corresponds to the crossover  (1). In each figure, the short black solid line represents the best fit for short times (first regime), and the long one the best fit for long times (whole regime). Notice that ∆ (s) refers to the unit of time in seconds. Table 2. γ, α, H, and C parameters of the temperature and conductivity for both first and whole regimes (short and long times) and for each one of the recorded CTD time series. Notice that for short times the alpha index is 2. Generalised moments for both temperature and conductivity measurements up to the fourth order were determined by using Equation (3) and the moments for the CTD30 time series are illustrated in Figure 5. A vertical solid line corresponds to the crossover point dividing the two time scales, and its choice is based on the turnover point of the variance. For both θ and E c , each moment shows dependence on the time lags.

Temperature θ Conductivity Ec
Each moment has been fitted by Equation (4), ρ(q, ∆) ≈ ∆ z(q) , and the obtained exponent is equal to the value of the structure function for this specific moment. In this way, the values of z(q) for the moments defined in the range 0.25 ≤ q ≤ 4 will be fitted later by Equation (5) for the estimate of the structure function. For the CTD30, the structure function for each property is represented in Figure 6. The value z(q) has a convex shape as a function of the order of the moment, q, for both short and long times. This convexity indicates the multifractality of temperature and conductivity time series. For the first regime (short times), the Lévy stability index, a, of the structure function is equal to 2 and points to the existence of a log-normal or Kolmogorov distribution, which describes the multiplicative effect of two random processes. In the simplest scenario, the first process is likely the volcano hydrothermal vents that would drive the second one, a short-term circulation over the crater that contributes to the rise of the temperature and conductivity variables. For longer times, the process becomes even more complicated. The stability index takes values a = 1.61 and a = 1.58 for temperature and conductivity, respectively, highlighting an increase in complexity as the observation time increases. And yet, this increase in complexity may reflect the existence of more than the previously mentioned two random processes that shape the overall multifractal character. The same is true for all CTD time series analysed in the vicinity of the crater because the structure function for short times always presents a log-normal distribution with Lévy stability index equal to 2 and for long times corresponds to a universal multifractal, Equation (5), whose Lévy stability index takes value in the range 1.04 ≤ a < 1.90, see Table 2.
Geosciences 2021, 11, x FOR PEER REVIEW 9 of 16 point dividing the two time scales, and its choice is based on the turnover point of the variance. For both θ and Ec, each moment shows dependence on the time lags. Each moment has been fitted by Equation (4), ( , ) ≈ ( ) , and the obtained exponent is equal to the value of the structure function for this specific moment. In this way, the values of ( ) for the moments defined in the range 0.25 ≤ ≤ 4 will be fitted later by Equation (5) for the estimate of the structure function. For the CTD30, the structure function for each property is represented in Figure 6. The value ( ) has a convex shape as a function of the order of the moment, q, for both short and long times. This convexity indicates the multifractality of temperature and conductivity time series. For the first regime (short times), the Lévy stability index, , of the structure function is equal to 2 and points to the existence of a log-normal or Kolmogorov distribution, which describes the multiplicative effect of two random processes. In the simplest scenario, the first process is likely the volcano hydrothermal vents that would drive the second one, a short-term circulation over the crater that contributes to the rise of the temperature and conductivity variables. For longer times, the process becomes even more complicated. The stability index takes values = 1.61 and = 1.58 for temperature and conductivity, respectively, highlighting an increase in complexity as the observation time increases. And yet, this increase in complexity may reflect the existence of more than the previously mentioned two random processes that shape the overall multifractal character. The same is true for all CTD time series analysed in the vicinity of the crater because the structure function for short times always presents a log-normal distribution with Lévy stability index equal to 2 and for long times corresponds to a universal multifractal, Equation (5), whose Lévy stability index takes value in the range 1.04 ≤ < 1.90, see Table 2.  For short times (Figure 6), the H and C parameters of the structure function of CTD30 take values of 0.662 and 0.055 for θ and 0.669 and 0.054 for Ec, respectively. Considering the structure functions obtained for all the CTD time series for the first regime ( = 2), listed in Table 2, the H exponent is varied between 0.540 and 0.695, implying super-normal processes if the intermittency parameter, parameter C, were zero. However, the latter is not true, and the non-zero intermittent effects establish conductivity and temperature mean fields as non-conservative. In the long-time limit, the structure-function of CTD30 provides values of H and C of 0.424 and 0.059 for θ and 0.434 and 0.060, for Ec, respectively. For long times, the Hurst exponents are lower than 0.5, implying stronger sub-normal processes and an increasing complexity reflected in the value of the alpha-stable index, which is different now for temperature and conductivity. For short times, the structure functions of both temperature and conductivity have a high similarity indicating a strong correlation of one another. Indeed, the normalised cross-correlation coefficient of the two, rθ,Ec, is higher than 0.95 for short times and decreases as the time increases (Figure 3c). Strong correlation does not necessarily imply the existence of causality between temperature and conductivity since this is the result of a broader mechanism affecting both properties equally and simultaneously, as the hydrothermal field does. Actually, by differentiating, with respect to time, the time series of temperature, Δθ, and conductivity, ΔEc, (Figure 3d,e), respectively, the normalised cross-correlation coefficient of the two, rΔθ, ΔEc, goes rapidly to zero (it is not exactly delta correlated) and the crossing of the abscissa is at 0.6 s. After that point and for a short window, before fluctuations start around zero (see Figure 3f), Δθ and ΔEc are anticorrelated to one another. These findings support the argument that the hydrothermal field is the reason for the strong correlation between temperature and conductivity since it operates as a vehicle carrying material to a sampling point.
The findings show that for both short and long time scales as defined above, the volcano area operates as an open system releasing a continuous flux of energy into the seawater, establishing thus temperature and conductivity to be non-conservative mean fields as a result of local swirls and/or vortexes. Alternatively, a non-conservative conductivity mean field might also be due to the presence of a magnetic field generated by the ions discharged from the crust of the volcanic hydrothermal chimneys to the water [6]. For short times (Figure 6), the H and C parameters of the structure function of CTD30 take values of 0.662 and 0.055 for θ and 0.669 and 0.054 for E c , respectively. Considering the structure functions obtained for all the CTD time series for the first regime (a = 2), listed in Table 2, the H exponent is varied between 0.540 and 0.695, implying super-normal processes if the intermittency parameter, parameter C, were zero. However, the latter is not true, and the non-zero intermittent effects establish conductivity and temperature mean fields as non-conservative. In the long-time limit, the structure-function of CTD30 provides values of H and C of 0.424 and 0.059 for θ and 0.434 and 0.060, for E c , respectively. For long times, the Hurst exponents are lower than 0.5, implying stronger sub-normal processes and an increasing complexity reflected in the value of the alpha-stable index, which is different now for temperature and conductivity. For short times, the structure functions of both temperature and conductivity have a high similarity indicating a strong correlation of one another. Indeed, the normalised cross-correlation coefficient of the two, r θ,Ec , is higher than 0.95 for short times and decreases as the time increases (Figure 3c). Strong correlation does not necessarily imply the existence of causality between temperature and conductivity since this is the result of a broader mechanism affecting both properties equally and simultaneously, as the hydrothermal field does. Actually, by differentiating, with respect to time, the time series of temperature, ∆θ, and conductivity, ∆E c , (Figure 3d,e), respectively, the normalised cross-correlation coefficient of the two, r ∆θ, ∆Ec , goes rapidly to zero (it is not exactly delta correlated) and the crossing of the abscissa is at 0.6 s. After that point and for a short window, before fluctuations start around zero (see Figure 3f), ∆θ and ∆E c are anticorrelated to one another. These findings support the argument that the hydrothermal field is the reason for the strong correlation between temperature and conductivity since it operates as a vehicle carrying material to a sampling point.
The findings show that for both short and long time scales as defined above, the volcano area operates as an open system releasing a continuous flux of energy into the seawater, establishing thus temperature and conductivity to be non-conservative mean fields as a result of local swirls and/or vortexes. Alternatively, a non-conservative conductivity mean field might also be due to the presence of a magnetic field generated by the ions discharged from the crust of the volcanic hydrothermal chimneys to the water [6].
To understand how the environment and seawater above Tagoro main crater are affected by those mechanisms governing the stochastic process and characterising the system, the spatial distribution of CTD time series according to its stochastic behaviour in the first regime (short times) is shown in Figure 7. There is not a significant trend in stochastic processes behaviour concerning the main crater distance, which is located over and between the CTDs 30, 34, 35 and 36. Analysing together the scaling of the mean field of temperature and its level of complexity, we found that the entire area, and not only the area above the crater, presents a higher activity and non-conservative field acting as an open system far from equilibrium. Near the crater, a less non-conservative region is found, characterised by an H smaller than 0.57. Additionally, in the middle of the high-resolution grid, the conductivity field can be considered an almost conservative one, C~0 (Figure 7). These findings likely reflect local structures which do not allow strong mixing with the overlying surrounding environment, or it might mirror a stronger discharge of the fluid near the crater that weakens the effect of the marine short-term circulation at the sampling point in the window of observation, or it can be a combination of the two. The overall result reduces the intermittency parameter. Located mainly in the grid edges, the CTDs 40, 48, and 50, whose short time behaviour in terms of variance points to Brownian motion, present a multifractal character as well and match with a lower intermittency parameter, C, shown in Table 2.
To understand how the environment and seawater above Tagoro main crater are affected by those mechanisms governing the stochastic process and characterising the system, the spatial distribution of CTD time series according to its stochastic behaviour in the first regime (short times) is shown in Figure 7. There is not a significant trend in stochastic processes behaviour concerning the main crater distance, which is located over and between the CTDs 30,34,35 and 36. Analysing together the scaling of the mean field of temperature and its level of complexity, we found that the entire area, and not only the area above the crater, presents a higher activity and non-conservative field acting as an open system far from equilibrium. Near the crater, a less non-conservative region is found, characterised by an H smaller than 0.57. Additionally, in the middle of the high-resolution grid, the conductivity field can be considered an almost conservative one, C ~ 0 (Figure 7). These findings likely reflect local structures which do not allow strong mixing with the overlying surrounding environment, or it might mirror a stronger discharge of the fluid near the crater that weakens the effect of the marine short-term circulation at the sampling point in the window of observation, or it can be a combination of the two. The overall result reduces the intermittency parameter. Located mainly in the grid edges, the CTDs 40, 48, and 50, whose short time behaviour in terms of variance points to Brownian motion, present a multifractal character as well and match with a lower intermittency parameter, C, shown in Table 2. In the extremes of the grid, the intermittency parameter of the conductivity field increases, highlighting higher complexity and displays the opposite behaviour for the temperature field. Possible contributions to the conductivity field at each sampling point can come directly from the volcano, and/or from the crust, and/or from diffusive ions from nearby points, and/or from the drifting of physical-chemical plumes, and/or from entrainment processes. If volcanic activity were the only source that affects the conductivity field, In the extremes of the grid, the intermittency parameter of the conductivity field increases, highlighting higher complexity and displays the opposite behaviour for the temperature field. Possible contributions to the conductivity field at each sampling point can come directly from the volcano, and/or from the crust, and/or from diffusive ions from nearby points, and/or from the drifting of physical-chemical plumes, and/or from entrainment processes. If volcanic activity were the only source that affects the conductivity field, higher values of intermittency, C, would describe a non-continuous outgassing process reminiscent of short volcano breaths followed by longer ones, whereas low or even zero values would describe a continuous flow of material whose H value provides the risk of eruption (high risk for H~1, ballistic behaviour, and low risk for H~0.5, random behaviour).
The hypothesis of diffusivity as a source that affects the conductivity field is also discarded because it is not able to wash out local inhomogeneities (e.g., a typical diffusion coefficient of 10 −9 m 2 /s for ions moving in water during our observations results only in a displacement of~1.3 mm) and to create a homogeneous one in such a reduced time of observation. Dynamic drifting is always present in open oceans, and its contributions might differ at different sampling points, thus mirroring the structure of the submarine volcano. The existence of drift is compatible with the variations of the temperature field. At sampling points where intermittency is almost zero for the conductivity field, its corresponding value for the temperature field is high indicating a strong mixing process. Certainly, there are points where intermittency takes high values for both fields mirroring strong mixing/competition of various random processes at these locations.

Discussion
The non-conservative mean fields are rather a consequence of the action of multiplicative noise, which can affect a physical system in many ways. To mention some of them, multiplicative noise can (i) stabilise a transient state [32], (ii) affect growth rates [33], and (iii) turn stationary distributions to non-Gaussian ones [34]. Furthermore, the strength of a multiplicative noise can modify the response of a system and unravel hidden periodicities, in the case that they exist, due to stochastic resonance [35]. Moreover, it is demonstrated that the effect of the buoyant seawater plume on the surrounding environment, together with the entrainment effect from the background waters and the general convective circulation over the hydrothermal sources, could potentially induce rapid changes in the temperature and conductivity mean fields [4]. This demonstrates that those physicalchemical anomalies or organisms remain near the source around the Tagoro system. Other physical oceanic processes, such as waves, tides, internal tides, etc., influence the variability of the temperature and conductivity mean fields at the Tagoro system on a larger time scale; however, a detailed study to understand such potential processes will be considered for futures studies.
The use of a high-resolution grid ( Figure 7) allows us (i) to unveil the sub-normal behaviour of both temperature and conductivity fields in the entire area, and not only inside the main crater whose alterations are detected even several meters above the seafloor, (ii) to designate potential mechanisms for the functionality of the volcano, such as intermittency describing the outgassing ions and drift terms redistributing this material in the entire area, and (iii) to characterise the volcano as an active one, since multifractal properties of the conductivity field show that the volcano continuously provides material to the marine environment, and that it yet operates as an open system because of the non-conservative temperature field. The results confirm that the Tagoro submarine volcano is an active vent system significantly altering the properties of the seawater surrounding the main crater. Such physical-chemical changes have often been documented as changes in the temperature, salinity, pH, and ORP [8] in the pH and carbonate system [12,36], in the emission of reduced species and O 2 [11], in nutrient enrichment [27], and also in biological changes in terms of species composition [11,16,26], mainly related with bacterial and planktonic groups [20,22]. In agreement with the high fluctuations observed in the CTD profiles and the emerging gasses and particles documented by ROV videos, our results show that the system is far from both equilibrium and homogeneity.
The observed multifractality needs a closer look at the causes that produce it. A system driven by noise sources, multiplicative and/or additive, is usually described by a general class of stochastic differential equations (see Equation (1) of [37]). Alternatively, complex systems can be described by (i) fractional calculus and generalised fractional Langevin equations [38,39], (ii) linear Langevin equation subject to multiplicative and/or to additive noise [34], and (iii) stochastic population models [40][41][42]. A stochastic approach and description of the present system should take into account the findings of fluctuation analysis. We leave this task for forthcoming work. It is worth mentioning that the method (GMM) analytically described in the present work is not the only one that can be used for fluctuations analysis. Methods such as detrended fluctuation analysis (DFA) [43], multifractal detrended fluctuation analysis (MFDFA) [44], diffusion entropy analysis (DEA) [45], or even the classical rescaled analysis (RA) [46] can be used if the corresponding time series is stationary.
For the first time, our study applied the methodology described by [6] to a shallow submarine volcano in the Atlantic Ocean. Previous studies had used time series data collected by a remotely operated vehicle (ROV) in Avyssos (Nisyros Island) [7], and in Kolumbo (Santorini Island) [6], a shallow submarine volcano, whose last eruption happened in 1650 [47]. The results of the study [6] showed that Kolumbo is a far-from-equilibrium system following sub-normal behaviour, similar to the results obtained in the present study for the Tagoro's hydrothermal system. For Kolumbo the crossover point, where the scaling exponent changes, is 5 s, similar to 3.5 s for the present study; however, a difference in the dynamics between the two (Kolumbo and Tagoro) might not be excluded due to closed/open sea conditions.
Overall, our results show that the GMM methodology can be applied to different systems and that the use of other oceanographical instruments for data collection is consistent with the effectiveness of the methodology. This can be useful since it allows choice of the instrument to be employed, although the ROV represent higher cost-effectiveness related to a rosette sampler since the utilisation of ROVs in oceanographical campaigns presents a notably higher cost.
The cost-effective methodology applied in this study may be exploited for forecasting future eruptions. By using this sampling methodology in time and using GMM for method analysis, we can characterise the level of activity (values of H and C for both temperature and conductivity) and monitor the evolution of such a physical system, predicting the natural hazard according to its potentiality. A combination of a high value of H and a value of C close to zero (almost homogeneous process, which means that the hydrothermal material dominates over the other processes) underlines an intense activity of the volcano, possibly according to a high hazard. On the other hand, for values of H close to 0.5 and non-vanishing values of intermittency (C parameter), the hazard is low since random environmental processes can affect the direct hydrothermal vent turning. The validity of the behaviour of H and C as prediction indexes should be tested in future experiments, which will monitor the dynamics of the El-Hierro submarine complex. These values of C and H are likely of importance for the existence of aquatic life in volcano surroundings, but extensive simulations in terms of fractional Langevin equations should be done to see how these parameters, especially C, affect the system. Such an analysis will deliver the memory kernel of the temperature and conductivity, examine possible repeated patterns, and help to study the future behaviour of Tagoro submarine volcano, see [6].

Conclusions
The Tagoro shallow submarine volcano, located south of El Hierro Island at the Canary Archipelago, has been deeply monitored in the last ten years since its origin. Its emissions, recorded with high sensitivity temperature and conductivity sensors mounted on an oceanographic rosette, exhibit important fluctuations in the surrounding of the main crater. The 21 θ-E c time series have been analysed by the GMM method in a short window of observation, delivering a wealth of findings. GMM returns structure functions with convex shapes indicating that both temperature and conductivity fluctuations correspond to multiplicative processes. This behaviour indicates the direct competition between the overlying marine short-term circulation and the degassing material of the volcano leading to long-lived non-equilibrium states. The nonstationary nature of the temperature field points to an open system that releases energy continuously to the environment. The estimated H and C parameters showed that the risk of eruption is low at the corresponding time window. However, continuous monitoring of the volcanic activity is necessary to further compare and classify the future volcano activity. The latter can be done by using the