Eight gas samples were collected from a fumarolic vent at Agia Kyriaki (AK) and 4 submarine gas emissions at Paleochori Bay (PA6, PA7, PA8, and PA9; Figures 1B, D) to analyze their chemical composition. These gases were collected in pre-evacuated glass bottles equipped with two Thorion stopcocks and partly filled with a solution of 4N NaOH. To sample the fumarole gas, we used a titanium tube inserted into the vent and connected to a Pyrex glass sampling line. The submarine gas emissions were instead sampled through an inverted funnel placed on the vent and connected to the glass bottles through silicon tubing. The gas samples were then analyzed at the laboratory of the INGV-Napoli, Osservatorio Vesuviano, using gas chromatography (He, H₂, Ar, O₂, N₂, and CH₄), ionic chromatography (H₂S as SO₄ ²⁻), and acidimetric titration (CO₂ as CO₃ ²⁻). These sampling techniques and analytical methods are discussed in detail in Bini et al. (2022) and Caliro et al. (2015).

To better understand the physicochemical processes that might drive hydrothermal explosions, we also used both geothermal well fluid and submarine thermal water data available from the literature. In particular, we used the chemical composition of the liquid sampled from the geothermal well Milos 2 [Table 1; data from Karabelas et al. (1989), Liakopoulos et al. (1991), Koutroupis (1992), and Dotsika et al. (2009)], and many pore waters from the Paleochori Bay recently published by Khimasia et al. (2021). This latter dataset consists of 152 analyses of thermal (T ≤ 122°C) and cold waters, and 13 analyses of seawater (Figure 1D). Nine water samples were excluded from our treatment since the electroneutrality parameter was >15% [according to Appelo and Postma (1993)]. Submarine thermal vent data from Valsami-Jones (2005), Wu et al. (2012), Price et al. (2013), and Roberts et al. (2021) were also used in this study.

In order to investigate the thermal-water origin of the Paleochori Bay and the eventual mixing processes between different sources, the concentration of Na and K in the porewater samples [from Khimasia et al. (2021)] are plotted against those of Mg in Figure 3. Samples are grouped according to the classification used in the original work of Khimasia et al. (2021) in “High Cl,” “Low Cl,” and “Background” samples (red, blue, and grey circles, respectively, in Figure 3 and the following Figures). Magnesium in a geothermal liquid equilibrates with the hosting-aquifer rocks and its concentration decreases with temperature according to water-rock equilibria (Giggenbach, 1988). At relatively high temperatures, such as >300°C, the equilibrium Mg content of a geothermal liquid approaches 0, if compared to that of seawater (∼1,330 mg·L⁻¹ at the Paleochori Bay). Figure 3 shows that Na and K are negatively correlated with Mg, suggesting a mixing between seawater and a deep geothermal component more saline than seawater. However, many samples are scattered (e.g., the “Low Cl” samples) likely due to the occurrence of other processes that overlap with those produced by a simple mixing. A process that explains this behavior is the boiling of ascending thermal waters, which experience steam separation and condensation of the separated vapor into the overlying liquids. Steam separation increases the concentration of the solutes in the liquid, while steam condensation would produce the opposite effect. However, hydrothermal steam separation and vapor condensation do not modify the ratios among the considered solutes, since “pure” water is either added or removed (see Eq. 5). In fact, the Na/K ratio plotted against the Mg/K precisely distributes along a mixing line between seawater and a Mg-poor geothermal component (PGL, Paleochori Geothermal Liquid in Figure 3C). Assuming the Mg content of PGL equals to 0, a Na/K ratio of 3.78 is extrapolated at the intercept of the best fit regression line at Mg/K = 0 (Figure 3C). It is worth noting that the inferred PGL Na/K ratio agrees with both that measured in the geothermal liquids encountered by the deep wells at Milos (3.81 ± 0.31; Table 1), and that extrapolated through linear regression from data of hydrothermal fluids discharged from the Paleochori Bay vents (3.67–4.03; Figure 3D).

We here estimate the PGL fraction and the amount of the vapor added or lost for each sample of the Paleochori Bay pore waters by applying the mass balance described in Section 2.2. Eqs 6, 7 can in fact be solved with respect to X and Y (X_PGL and Y_vap hereafter), given the concentrations of the 2 species i and j, namely, Na and K, in the endmembers. The seawater concentrations of Na and K were assumed as the mean of the 13 seawater samples from the Paleochori Bay (Khimasia et al., 2021). The Na concentration of PGL was assumed at 26,300 mg·L⁻¹, i.e., that of the liquid of Milos 2 well at reservoir conditions (Table 1). The K concentration was computed at 6,960 mg·L⁻¹ dividing the Na concentration by the estimated Na/K ratio of 3.78 (Figure 3C).

The validity of this model, i.e., the PGL equals to the Milos 2 liquid at reservoir conditions and the reliability of the computed X_PGL and Y_vap, was then tested with the following procedure. First, we corrected each porewater composition by eliminating the effect of the steam added or removed through the computed Y_vap [i.e., the corrected concentration C_i,corrected = C_i,measured × (1+ Y_vap)]. Then, we plotted the corrected concentrations versus X_PGL (Figure 4), and we estimated the original concentrations in the pure geothermal component for each solute at X_PGL = 1. The estimated concentrations of Cl, Ca, Mg, Sr, and Li for PGL are practically equal to those measured in the deep well Milos 2 (Figure 4). Both the high correlations and overlaps between estimated and measured concentrations of the geothermal component strongly support the reliability of the method, i.e., the used mass balance approach gives reliable X_PGL and Y_vap values. Furthermore, Figure 4 shows that the unmixed geothermal component at the Paleochori Bay is very similar to that of Milos 2 well and that Na, K, Cl, Ca, Mg, Sr, and Li behave as mobile species. It is worth nothing that, contrarily, SO₄ and SiO₂ are not conservative species. The estimated SO₄ concentration of 608 mg·L⁻¹ is in fact much larger than that of Milos 2 (10–18 mg·L⁻¹), while SiO₂ is much lower (Figure 4). This is of utmost importance for our aim of evaluating the conditions that favor hydrothermal explosions, as it highlights the extensive occurrence of two processes that contribute to form the impermeable cover of the hydrothermal system: (i) silica precipitation and (ii) addition of SO₄ from H₂S oxidation to sulfuric acid. The interaction of acidic solutions with the hosting rocks indeed causes the formation of the so-called advanced argillic alteration (kaolinite, alunite, gypsum, and so forth; Steiner, 1977). The development of these silica and clay caps enables fluid accumulation and pressure buildup in the subsurface.

The very similar composition of PGL with that of the liquid of Milos 2 is also suggested by the stable isotopes of water, which are plotted against the X_PGL values in Figure 5. To minimize the effects due to steam addition-separation in Figure 5, we only report the samples with −0.1 < Y_vap < 0.1. These two independent variables also highlight the mixing between seawater (local Mediterranean water, δD = 9.5 ± 1.2; δ¹⁸O = 1.65 ± 0.15) and PGL, for which we estimated, from the best-fit regressions at X_PGL = 1, a δD = −1.6‰ and a δ¹⁸O = 2.9‰. These values are very close to those reported for the fluids of the geothermal well Milos 2 (δD = −2.3‰ and δ¹⁸O = 3.2‰; Dotsika et al., 2009 and references therein).

The original temperature of PGL is investigated using the graphical approach based on the Na-K-Mg contents in Figure 6 (Giggenbach, 1988; Section 2.3). The Paleochori Bay pore waters plot along a mixing line between seawater and the PGL, whose Na and K concentrations are estimated at 26,300 and 6,960 mg·L⁻¹, respectively (see Section 3.1.2). Figure 6 also shows mixing lines constructed using a PGL with different salinities (i.e., different Na contents). The Paleochori pore waters roughly overlap the curve with Na content of 26,300 mg·L⁻¹, giving independent evidence of the validity of the Na content of Milos 2 as representative of that of the PGL. The Mg content, which in Figure 3C was assumed equal to 0, is here inferred through the binary mixing hyperbola method (Bini et al., 2022 and references therein) in the Na/Mg-Mg/K space (inset in Figure 6), using the High Cl waters (i.e., the most representative of the thermal component). Practically, we computed the Mg/K ratio of PGL as the vertical asymptote (0.0014 ± 0.0017) of the hyperbolic binary mixing between PGL and seawater (Figure 6). The Mg content is then estimated at 9.7 mg·L⁻¹ by multiplying 0.0014 by the K content.

It is worth noting that this Na-K-Mg composition of the PGL plots on the full equilibrium curve (Figure 6), showing a temperature of 324°C. Notably, this high temperature estimated for PGL is supported by both direct temperature measurements (up to 323°C at ∼1,000 m depth in Milos 1; Figure 2; Mendrinos and Sullivan, 1987) and other geothermometers (Table 1). For example, the Na/Li geothermometer (Fouillac and Michard, 1981) returns a temperature of 321°C for PGL (Table 1).

Maps of temperature, fraction of geothermal component (X_PGL), and steam added/removed from pore waters (Y_vap; Figure 7), constructed using sequential Gaussian simulations (sGs; see Section 2.4), enable us to spatially constrain the ascent and decompressional boiling of geothermal fluids. All these three variables are spatially correlated (Figure 7), showing high values particularly in the western and central coastal sectors of the Paleochori Bay. The upflow of PGL predominantly occurs from the western part of the bay (Figure 7C), characterized by values of X_PGL > 0.2 (Figure 7C). Notably, steam condensation processes, that is high Y_vap values, focus on a few spots aligned along the NW-SE direction in the center of the bay, in correspondence of a known fault (Figure 7B).

The chemical composition of the eight gas samples collected from both the Agia Kyriaki fumarolic vent (AK) and the submarine gas emissions at high temperatures (>100°C) in the Paleochori Bay (PA6, PA7, PA8, and PA9; Figure 1D) is reported in Table 2, together with a geothermal vapor of Milos (MIL) available from the literature (Giggenbach, 1997). The gas composition is dominated by CO₂ with relatively high contents of H₂S (2.4%–4%) and H₂ (up to 2%). These high H₂ contents in the Paleochori Bay submerged vents were also reported by other authors [up to 3% in Daskalopoulou et al. (2018) and Dando et al. (1995)]. Notably, the H₂/CO₂ ratio of PA7 and PA8 (0.01–0.02) does not significantly differ from that measured at AK (0.022) and MIL (0.018), where the steam (99.3% and 99.6%, respectively) is the main component. These high H₂ contents, much higher than those normally measured in CO₂-dominated gas emissions from geothermal fields (such as in Italy and Greece; Chiodini and Cioni, 1989; Daskalopoulou et al., 2018), suggest a steam-dominated gas phase separated through boiling, whose vapor is quickly removed in the “cold” submarine environment. This is further evidence of a boiling process at the Paleochori Bay, as previously shown by the many hydrothermal vents close to or above the boiling point (inset in Figure 2) and the composition of the pore waters (see Section 3.1).

Mixing between different gas sources is supported by the chemical composition of these gases, particularly by their content of chemically inert species, such as N₂, He, and Ar. The relative composition of these species is reported in the N₂-He-Ar triangular diagram (Giggenbach, 1996) of Figure 8 and compared with those of the air and air saturated water (ASW), or seawater. The submerged gas emissions plot between the composition of geothermal vapors (AK and MIL) and that of ASW. This trend reflects a thermal water variably mixed with seawater, which supplies to the gas phase, separated through boiling, N₂, Ar, and He originally dissolved in the marine solutions (i.e., the ASW endmember). The gas samples collected in October 2021 thus represent the incondensable fraction of the vapor generated during the boiling of solutions formed by the mixing of seawater with a geothermal component. Notably, this binary mixing is also supported by the equilibrium temperature of the gases (T_H2-Ar; Table 2; Figure 8), estimated through the H₂-Ar geothermometer (Giggenbach, 1991) for a gas phase dissolved in liquid water. High T_H2-Ar of 328°C–339°C, similar to those measured in the geothermal wells and estimated for PGL through the Na-K-Mg content of pore waters (Section 3.1.3), characterize the AK fumarole and the geothermal vapor MIL. On the other hand, submerged gas emissions show decreasing T_H2-Ar down to 200°C, consistently with the addition of seawater (Figure 8).