Valentine Puzenat
email: valentine.puzenat@gmail.com

Javier Escartin

Jean-Emmanuel Martelat

Thibaut Barreyre

Sven Le Moine Bauer

Paraskevi Nomikou

Nuno Gracias

Pascal Allemand

Varvara Antoniou

Omer Coskun

Rafael Garcia

Philippe Grandjean

Steffen Leth Jørgensen

Lluís Magí

Manolis Mandalakis

William Orsi

Paraskevi Polymenakou

Anders Schouw

Guillem Vallicrosa

Othonas Vlasopoulos

Shallow-water hydrothermalism at Milos (Greece): Nature, distribution, heat fluxes and impact on ecosystems

Keywords: Milos, shallow-water, hydrothermalism, AUV, drone, photomosaics

Submarine hydrothermal activity is responsible for heat and chemical exchanges through the seafloor. Shallow-water hydrothermal systems (SWHS), while identified around the globe, are often studied in a way that is less comprehensive than their deep-ocean counterparts (e.g., along ridges), where systematic optical and acoustic mapping is more prevalent and coupled to in situ observations and sampling. Using aerial drones, an AUV, and temperature measurements at 10-40 cm subseafloor, we investigated in 2019 one of the most extensive SWHS known to date, in Paleochori (south of Milos, Greece). Hydrothermal venting, found from the shore to water depths of almost 500 m, shows emissions of gases and high-temperature fluids, often associated with bacterial mats and/or hydrothermal mineral precipitates. This study provides extensive drone mapping coupled with local AUV surveys for seafloor characterization and ground-truthing from the shore to ~20 m water depth. Seafloor photomosaics also provide a detailed context to samples, measurements and observations carried in situ. We interpret the photomosaics to define distinct seafloor types, linked to this hydrothermal activity. White hydrothermal patches (WHPs) often show a clear polygonal organization, together with outflow areas that are both more dispersed and distributed. Polygonal patterns likely result from fluid convection in a sandy porous medium heated from below. These WHPs display elevated subseafloor temperatures, typically >50 °C, with maximum values of ~75 °C. Photomosaics also display textures of biological origin, including seagrass and bioturbation patterns. Widespread bioturbation by burrowing shrimps is often associated with WHPs, bounding them, but also occurs on sandy seafloor away from hydrothermal patterns. Subseafloor temperatures at these bioturbated areas are of ~30-40 °C, and are thus transitional between hot WHPs and sedimented seafloor unaffected by hydrothermal activity (~24 °C). In addition to linking subseafloor temperature data and interpreted seafloor 4 photomosaics, our results provide a comprehensive general overview of this SWHS, of the organization of its hydrothermal outflow through the seafloor, and of the underlying subseafloor fluid circulation. This paper also gives the first perspectives on the heat fluxes of the system, and constitutes a background for other studies on the nature and distribution of microbial communities, controlled by this hydrothermal activity.

Introduction

The circulation of fluids in the Earth's lithosphere is often driven by heat sources at depth, and much of this activity occurs underwater. Hydrothermal circulation drives heat and chemical exchanges through the seafloor, is associated with water-rock reactions at depth and sustains numerous ecosystems. To constrain the setting and overall structure of these hydrothermal systems, detailed optical surveys of the seafloor are needed for a comprehensive geological mapping. These data are also crucial to understand the distribution, geometry and nature of hydrothermal outflow through the seafloor, and to extrapolate local observations (e.g., temperature and chemical measurements) in order to evaluate the associated mass, thermal and chemical fluxes (e.g., [START_REF] Price | Processes influencing extreme As enrichment in shallow-sea hydrothermal fluids of Milos island, Greece[END_REF][START_REF] Callac | Modes of carbon fixation in an arsenic and CO2-rich shallow hydrothermal ecosystem[END_REF]. Systematic and comprehensive mapping also provides a baseline for temporal studies, providing a context for in situ measurements, and establishing the links between hydrothermal activity, ecological and biological processes. To date, this is done partially in deep water settings (e.g., [START_REF] Maki | Large-area visual mapping of an underwater vent field using the AUV "Tri-Dog 1[END_REF][START_REF] Barreyre | Structure, temporal evolution, and heat flux estimates from the Lucky strike deep-sea hydrothermal field derived from seafloor image mosaics[END_REF][START_REF] Marcon | Megafaunal distribution and assessment of total methane and sulfide consumption by mussel beds at Menez Gwen hydrothermal vent, based on geo-referenced photomosaics[END_REF], but such mapping studies are less developed in shallow-water hydrothermal systems. While abundant and particularly ubiquitous at volcanically active areas (e.g., [START_REF] Tarasov | Deep-sea and shallowwater hydrothermal vent communities: Two different phenomena[END_REF][START_REF] Price | A review of the geochemistry and microbiology of marine shallow-water hydrothermal vents[END_REF][START_REF] Martelat | Terrestrial shallow water hydrothermal outflow characterized from out of space[END_REF], their geometry, temporal evolution, and associated fluxes are yet not well constrained. Systematic and comprehensive mapping of shallow-water systems to date is therefore partial and limited to a few sites (e.g., [START_REF] Khimasia | Hydrothermal areas, microbial mats and sea grass in Paleochori Bay, Milos, Greece[END_REF][START_REF] Khimasia | Hydrothermal flux and porewater geochemistry in Paleochori Bay, Milos, Greece[END_REF].

Here we present a comprehensive study on the organization of hydrothermal outflow at the shallow water hydrothermal system off the coast of Milos (Greece). Numerous sites around Milos extend to depths of >200 m, and several studies have investigated those on the south-eastern coast of the island over the last few decades (e.g., Aliani et al., 1998a;[START_REF] Aliani | Periodicities in sediment temperature time-series at a marine shallow water hydrothermal vent in Milos Island (Aegean Volcanic arc, Eastern Mediterranean)[END_REF] Journal Pre-proof [START_REF] Valsami-Jones | The geochemistry of fluids from an active shallow submarine hydrothermal system: Milos island, Hellenic Volcanic Arc[END_REF][START_REF] Price | Processes influencing extreme As enrichment in shallow-sea hydrothermal fluids of Milos island, Greece[END_REF]. These sites are most visible in sedimented, sandy seafloor, where hydrothermal outflow is associated with bacterial mats and hydrothermal precipitates. Therefore, prior studies have focused primarily on the associated geochemical and microbiological processes (e.g., [START_REF] Callac | Modes of carbon fixation in an arsenic and CO2-rich shallow hydrothermal ecosystem[END_REF][START_REF] Price | A review of the geochemistry and microbiology of marine shallow-water hydrothermal vents[END_REF][START_REF] Houghton | Spatially and temporally variable sulfur cycling in shallow-sea hydrothermal vents, Milos, Greece[END_REF][START_REF] Megalovasilis | Geochemistry of Hydrothermal Particles in Shallow Submarine Hydrothermal Vents on Milos Island, Aegean Sea East Mediterranean[END_REF]. While earlier studies have provided maps with distribution of active venting areas extending from the shore to deep regions (e.g., [START_REF] Dando | Gas venting rates from submarine hydrothermal areas around the island of Milos[END_REF], the actual geometry and distribution of sites and venting areas are poorly constrained. Only recently we have obtained detailed and accurate mapping of the Milos hydrothermal system to depths of ~20 m using satellite imagery [START_REF] Martelat | Terrestrial shallow water hydrothermal outflow characterized from out of space[END_REF], and drone aerial imagery coupled with both temperature measurements and fluid chemistry [START_REF] Khimasia | Hydrothermal areas, microbial mats and sea grass in Paleochori Bay, Milos, Greece[END_REF][START_REF] Khimasia | Hydrothermal flux and porewater geochemistry in Paleochori Bay, Milos, Greece[END_REF]. Here we use underwater surveys to adequately ground truth interpretation of drone imagery, and to provide a detailed context where temperature measurements can be precisely linked to the observed seafloor patterns. This approach, that differs from that of prior studies, allow us to correlate directly spatial patterns of seafloor hydrothermal activity and subseafloor temperature.

For this study we have thus conducted extensive aerial drone surveys at the Milos nearshore sites of Agia Kyriaki and Spathi Bay, in addition to Paleochori Bay, the main study site. The resulting aerial photomosaics provide a regional overview of the distribution of hydrothermal activity at these sites, complementing the prior study of [START_REF] Khimasia | Hydrothermal areas, microbial mats and sea grass in Paleochori Bay, Milos, Greece[END_REF] at Paleochori bay, and extending it to adjacent areas. Drone surveys are coupled with Autonomous Underwater Vehicle (AUV) surveys to generate high-resolution seafloor photomosaics that are used for seafloor characterization and ground-truthing of the broader aerial drone seafloor imagery. Drone and AUV data reveal the detailed organization of hydrothermal activity at these sites and their link to associated ecosystems. We completed our study with subseafloor temperature measurements, Journal Pre-proof sampling (sediment cores), and in situ observations carried out by divers at Paleochori Bay.

These observations and measurements are thus carefully located relative to the structure and organization of hydrothermal outflow at specific sites imaged by drone or AUV.

With the analysis and interpretation of drone and AUV imagery, temperature measurements, and in situ observations, we discuss (a) seafloor hydrothermal patterns and their links with subseafloor hydrothermal circulation, (b) the correlation of temperature records with these patterns, that can subsequently be used in hydrothermal heat flux quantification, and (c) the study of seafloor patterns and temperature data to provide information on the spatial associations of certain ecosystems with hydrothermal activity. While the in situ observations, measurements, and samples were acquired at Paleochori Bay, they allow us to propose a general, regional overview of the shallow water hydrothermal system of Milos, and to discuss the organization of hydrothermal outflow, its nature and spatial distribution. Our results may also constrain the organization of hydrothermal outflow at other systems in sedimented areas, and provide a context to adequately constrain heat fluxes.

Journal Pre-proof

Geological background

Milos (Greece) is located at the central part of the Hellenic Volcanic Arc, where the African plate subducts beneath the Aegean micro-plate (Figure 1). In this particular volcano-tectonic setting, numerous volcanic areas often show intense hydrothermal activity with underwater manifestations (e.g., Methana, Santorini, Milos, Kos; [START_REF] Dando | Hydrothermal studies in the Aegean Sea[END_REF]. Among these, Milos hosts one of the most important and extensive shallow-water hydrothermal system known to date [START_REF] Fytikas | Volcanology and petrology of volcanic products from the island of Milos and neighbouring islets[END_REF][START_REF] Papanikolaou | Correlation on neotectonic structures with the geodynamic activity on Milos during the earthquakes of March 1992[END_REF]Dando et al., 1995a;[START_REF] Aliani | Periodicities in sediment temperature time-series at a marine shallow water hydrothermal vent in Milos Island (Aegean Volcanic arc, Eastern Mediterranean)[END_REF][START_REF] Nomikou | Submarine volcanoes along the aegean volcanic arc[END_REF]. At Paleochori Bay (SE of Milos island, Figure 1), known hydrothermal venting is distributed over a wide region from the shore to deeper areas, as identified both during ROV dives [START_REF] Nomikou | Submarine volcanoes along the aegean volcanic arc[END_REF], and as gas flares at water depths of at least 300 m [START_REF] Dando | Gas venting rates from submarine hydrothermal areas around the island of Milos[END_REF][START_REF] Dando | Hydrothermal studies in the Aegean Sea[END_REF]. Additional venting sites have been identified around the island [START_REF] Martelat | Terrestrial shallow water hydrothermal outflow characterized from out of space[END_REF] but with no detailed image surveys (e.g. using aerial drone in the shallow water, or underwater imagery).

Hydrothermal venting at Paleochori displays emissions of high temperature fluids of low pH and varying chlorinity, involving free gases (CO 2 , H 2 , H 2 S,CH 4 ) that are variably mixed with seawater [START_REF] Dando | Gas venting rates from submarine hydrothermal areas around the island of Milos[END_REF][START_REF] Valsami-Jones | The geochemistry of fluids from an active shallow submarine hydrothermal system: Milos island, Hellenic Volcanic Arc[END_REF]. Fluids are saline and sulphidic, depleted in Mg 2+ and SO 4 2-and particularly enriched in As and CO 2 , in addition to Na + , Ca 2+ , K + , Cl -, SiO 2 and other dissolved gases [START_REF] Fitzsimons | Submarine hydrothermal brine seeps off Milos, Greece: observations and geochemistry[END_REF][START_REF] Valsami-Jones | The geochemistry of fluids from an active shallow submarine hydrothermal system: Milos island, Hellenic Volcanic Arc[END_REF]Price et al., 2013). Hydrothermal venting is often associated with orange and white precipitates of bacterial and hydrothermal mineral origin, which are respectively composed of As-sulfide and a mix of silica and sulfur [START_REF] Price | Processes influencing extreme As enrichment in shallow-sea hydrothermal fluids of Milos island, Greece[END_REF], Yücel et al., 2013;[START_REF] Godelitsas | Amorphous As-sulfide precipitates from the shallow-water hydrothermal vents off Milos Island (Greece)[END_REF][START_REF] Callac | Modes of carbon fixation in an arsenic and CO2-rich shallow hydrothermal ecosystem[END_REF]Durán-Toro et al., 2019;[START_REF] Houghton | Spatially and temporally variable sulfur cycling in shallow-sea hydrothermal vents, Milos, Greece[END_REF][START_REF] Khimasia | Hydrothermal flux and porewater geochemistry in Paleochori Bay, Milos, Greece[END_REF]. Since these white seafloor areas are highly reflective, hydrothermal seafloor textures can be mapped subaerially Journal Pre-proof using satellites [START_REF] Martelat | Terrestrial shallow water hydrothermal outflow characterized from out of space[END_REF], or drones [START_REF] Khimasia | Hydrothermal areas, microbial mats and sea grass in Paleochori Bay, Milos, Greece[END_REF]. A recent satellite study has thus documented hydrothermal activity at and around Paleochori and on the North-East coast of Milos island [START_REF] Martelat | Terrestrial shallow water hydrothermal outflow characterized from out of space[END_REF], but this study is limited to water depths of <20-30 m.

Venting at other sites is documented from water column plumes identified in echo sounding profiles [START_REF] Dando | Hydrothermal studies in the Aegean Sea[END_REF]. Satellite imagery also reveals an apparent stability of the geometry and extent of these sites over time scales of a few years to a decade, both at Paleochori and at other areas around Milos [START_REF] Martelat | Terrestrial shallow water hydrothermal outflow characterized from out of space[END_REF].

Heat, mass, and chemical fluxes are important to understand the dynamics of hydrothermal systems. Seafloor and subseafloor temperatures have been measured locally at Paleochori (e.g. Aliani et al., 1998b;[START_REF] Aliani | Periodicities in sediment temperature time-series at a marine shallow water hydrothermal vent in Milos Island (Aegean Volcanic arc, Eastern Mediterranean)[END_REF][START_REF] Valsami-Jones | The geochemistry of fluids from an active shallow submarine hydrothermal system: Milos island, Hellenic Volcanic Arc[END_REF], and recent papers report systematic temperature measurements coupled to geochemical observations (e.g., [START_REF] Khimasia | Hydrothermal areas, microbial mats and sea grass in Paleochori Bay, Milos, Greece[END_REF][START_REF] Khimasia | Hydrothermal flux and porewater geochemistry in Paleochori Bay, Milos, Greece[END_REF][START_REF] Megalovasilis | Geochemistry of Hydrothermal Particles in Shallow Submarine Hydrothermal Vents on Milos Island, Aegean Sea East Mediterranean[END_REF], with an associated mass flux for Paleochori estimated at ~11300 ± 1100 m3/day [START_REF] Khimasia | Hydrothermal flux and porewater geochemistry in Paleochori Bay, Milos, Greece[END_REF]. These measurement sites cannot be spatially collocated with seafloor structures (imaged by either drone or AUV), as required by our approach to establish a reliable flux estimate that takes into account the geometry of outflow at the seafloor.

Gas flow rates have been measured at individual gas vents [START_REF] Dando | Gas venting rates from submarine hydrothermal areas around the island of Milos[END_REF][START_REF] Botz | Hydrothermal gases offshore Milos Island, Greece[END_REF][START_REF] Megalovasilis | Geochemistry of Hydrothermal Particles in Shallow Submarine Hydrothermal Vents on Milos Island, Aegean Sea East Mediterranean[END_REF]. As available measurements are discrete and punctual, and lacking constraints on the density and distribution of degassing, we cannot extrapolate these values regionally.

Journal Pre-proof 3 Methods

The data used in this paper were acquired during two field work operations: Drone and Autonomous Underwater Vehicle (AUV) images were collected in June 2019, and completed in September 2019 with additional drone imagery, in situ temperature measurements, seafloor sampling as well as photos and videos from divers and GPS positions at the surface during all operations. We conducted studies over a wide range of scales, including satellite [START_REF] Martelat | Terrestrial shallow water hydrothermal outflow characterized from out of space[END_REF], drone, AUV, and in situ images (Figure 2), together with in situ measurements and samples. Study areas along the SW shore of Milos included Agia Kyriaki, Spathi and Paleochori Bays. At this latter site we conducted a focus study on a polygonal hydrothermal pattern, as detailed below. The sampling of sediments and associated microbial ecosystems will be the scope of another study focusing on the links between microbiology and the thermal and chemical environment controlled by this hydrothermal activity.

Seafloor mapping

The distribution, nature and geometry of the seafloor hydrothermal structures and features reveal the organization of subseafloor hydrothermal circulation. A systematic mapping of our drone and AUV imagery at different scales, groundtruthed by in situ images and observations (Figure 2), is thus required to characterize precisely the different seafloor textures and their distribution, linking them to hydrothermal processes.

Drone imagery

Drone surveys at Agia Kyriaki, Paleochori and Spathi Bays (Figure 1) were conducted using two drones: a Phantom Pro (weight 1.4 kg, Camera FC6310 8.8 mm 20 Mpix) and a Mavic Pro (weight 0.907 kg, Camera L1D-20c 10.26 mm 20 Mpix). Surveys were pre-programmed before

Journal Pre-proof take-off and flown on automatic mode. We processed and photomosaiced drone imagery using the Agisoft software to obtain orthophotos with a resolution of ~10 cm per pixel. Away from the shore some images failed to be aligned as we lacked clear seafloor markers for matching. In these cases, images were visually inspected and manually aligned. Best images were usually acquired at sunrise around 05:30 am, with no or limited sun reflection on the sea surface, and with low winds. Images taken later during the day displayed more contrast visually but sun reflections were present in ~30 % of the individual images. As for satellite imagery [START_REF] Martelat | Terrestrial shallow water hydrothermal outflow characterized from out of space[END_REF], drone image quality also depends on water turbidity. The areas covered in the drone surveys are shown in Figures 1 and S1 and a summary of information from each area is shown in Figure 3.

We visually identified and defined several seafloor types, based solely on the visual inspection of drone imagery. First, the most prominent pattern is seafloor with white patches of hydrothermal origin, named White Hydrothermal Patches (WHPs) from here on, as shown in Figures 2 and3.

Throughout the study area the intensity and homogeneity of these patches vary significantly, and hence their visual appearance also displays significant variability. We thus define three subcategories based on the reflectivity and whiteness of the seafloor (High, Medium, and Lowdensity WHPs, Figure 3). Second, we observe patches of darker material that seem to be closely associated to the WHPs, with similar patterns and geometrical arrangements, named Dark Hydrothermal Patches (DHPs) from here on (Figures 3e,3f). Associated to DHPs we repeatedly identify clear linear features, often bounding WHPs instead of forming broad and diffuse patches, that are here identified as Dark Lines (Figures 3a,3b). Lastly, we map dense meadows of seagrass that correspond to Posidonia oceanica (Figure 3), based both on prior studies (Aliani et al., 1998b[START_REF] Martelat | Terrestrial shallow water hydrothermal outflow characterized from out of space[END_REF] and in situ field observations. Both the shapefiles that we have

Journal Pre-proof generated and the drone mosaics used in this study are provided in open access (see Data Availability Statement below for details).

AUV imagery

We used the SPARUS II Autonomous Underwater Vehicle (AUV, [START_REF] Carreras | Sparus II AUV-A Hovering Vehicle for Seabed Inspection[END_REF] from the University of Girona (Spain) to image the seafloor at a resolution higher than that of drone imagery, and to gather direct seafloor images for ground-truthing. We conducted a total of 25 AUV surveys at both Paleochori and Spathi Bays (See Figures 1 andS2). The AUV was programmed to follow lawn-mower patterns, flying at an altitude of ~5 m above the seafloor. The combined total AUV track flown during the surveys was ~41 km.

Preliminary optical photomosaics were created in the field for data quality checks, using only the AUV navigation information. Final AUV photomosaics were created offline using image registration and pose-graph optimization [START_REF] Elibol | Augmented state-extended Kalman filter combined framework for topology estimation in large-area underwater mapping[END_REF]. This processing also uses the AUV navigation information [START_REF] Campos | Underwater Multi-Vehicle Trajectory Alignment and Mapping Using Acoustic and Optical Constraints[END_REF], and thus provides a general georeferencing of the mosaic. As the surveyed seafloor was generally flat with limited local relief, a 2D image reconstruction approach was preferred over of a full 3D reconstruction. This optimized the handling of low-overlap images that hinder 3D reconstructions [START_REF] Gracias | Application Challenges of Underwater Vision[END_REF]. Final mosaics were blended [START_REF] Prados | A novel blending technique for underwater gigamosaicing[END_REF] and rendered with a ground resolution ranging from 3 to 10 mm per pixel, compared to the 6-10 cm per pixel resolution of drone imagery.

As the AUV navigation is primarily inertial, there is an inherent vehicle drift, and as a result AUV photomosaics, while generally well-placed, display both spatial shifts and mild distortions relative to drone photomosaics. To work on a geographically consistent dataset, we shifted and Journal Pre-proof wrapped the AUV photomosaics relative to the drone ones using tie points (Figure S3), generating a new set of AUV seafloor mosaics co-registered with drone imagery.

We divide AUV photomosaics into three zones: Paleochori West, Paleochori East and Spathi.

Seafloor textures, digitized at a scale of 1:40, are identified solely based on the AUV imagery, and are briefly described here. First, as on drone mosaics, patches of white seafloor are highly reflective. We also refer to them as White Hydrothermal Patches (WHPs, Figure 4), defining three subcategories based on their density and reflectivity (Figure 4d). WHPs -High density are patches with highly reflective white material, with no sediments visible within. WHPs -Medium density are patches of less reflective white material, locally discontinuous, showing within sediments and associated textures (e.g. ripple marks). WHPs -Low density correspond to sediments showing clear ripples through discontinuous white materials. The boundaries of these textures, and particularly those of low density, are gradual on both AUV and scuba diver imagery.

To define and digitize the limits of these textures, we evaluated them visually throughout the different AUV photomosaics and at different scales. We also observe patches of darker material in Spathi Bay (Figures 4e,4f), here referred to as Dark Hydrothermal Patches (DHPs), as well as Dark Lines (Figures 4e,4f) that are similar to the patterns defined above from the drone photomosaics.

The high resolution of the AUV relative to drone imagery allows us to identify substrates of biological origin that are not observable in drone photomosaics. First, we observe continuous zones of sediment mounds of up to a few tens of cm in diameter. These mounds often surround the White Hydrothermal Patches but are also found on sandy seafloor lacking hydrothermal patterns. Here we refer to this seafloor texture as Bioturbation (Figure 4). Second, we observe meadows of Posidonia oceanica (Figures 4e,4f), also identifiable on drone photomosaics, and

Journal Pre-proof reported in prior studies of the area (Aliani et al., 1998b;[START_REF] Martelat | Terrestrial shallow water hydrothermal outflow characterized from out of space[END_REF]. Lastly, we can identify another type of seagrass that forms sparse meadows with visible sandy seafloor corresponding to Cymodocea nodosa (Figure 4c, 4d, Aliani et al., 1998b). These Cymodocea meadows show significant variations in density and often disappear gradually. Here we digitize the densest patches only, determining boundaries based on their optical appearance at multiple scales. The digitized shapefiles and the uncorrected and corrected AUV photomosaics (coregistered with the drone photomosaics) are provided as open-access data (see Data Availability

Statement below for details).

Surface quantification from drone and AUV imagery

The seafloor features described above were digitized as polygons or lines in QGIS. This digitization allows us to estimate the surfaces covered by the different seafloor types, and to investigate their spatial relationships. Resulting shapefiles were then used to create binary grids (resolution of 10 cm) for each substrate type on both drone and AUV data (0 being the absence of substrate, and 1 the occurrence of this substrate). These binary grids facilitate calculations and analyses of seafloor surfaces (performed in MATLAB ® ), which are reported in Tables 2 and3.

We also report in Table 4 surfaces from drone imagery calculated solely over AUV survey areas.

As discussed below, quantitative comparisons of AUV and drone interpretation are critical to validate the extrapolation of local AUV observations to regional scales, using drone imagery.

Temperature measurements and processing

We used different autonomous temperature sensors to measure subseafloor temperature. Sensors have a rod with a thermocouple at their tips, and a housing with the data logger and batteries (Figures 5b, 5c, Figure 6a). In this experiment we deployed two WHOI-MISO high-temperature

Journal Pre-proof sensors, one WHOI-MISO low-temperature sensor, four NKE S2T6000DH-Ti high-temperature sensors and three NKE S2T6000D low-temperature sensors (see details in Table 1).

We also combined instruments to measure subseafloor temperature at different depths and obtain vertical temperature gradients (Tables S1 andS2). These included a set of 3 NKE instruments, named "Poseidon" from here on (Figure 6a), and two sets of 2 NKE sensors named "Set 1" and "Set 2". We also deployed a thermal blanket to measure the heat flux at the center of a polygonal hydrothermal cell (Figure 11). The thermal blanket includes two high-precision Antares thermistors (-4 -50 °C, accuracy of r 0.001 °C, sampling rate of 2 s, see [START_REF] Salmi | Quantitative estimate of heat flow from a mid-ocean ridge axial valley, Raven field, Juan de Fuca Ridge: Observations and inferences[END_REF], at the bottom and top of the blanket, separated by the blanket material. As the thermal conductivity of the blanket material is known, heat flux F can be estimated using Fourier's Law (F=-k dT(z)/z,

where k is the thermal conductivity, T the temperature, and z the depth).

Sensors were deployed in September 2019, solely at Paleochori Bay (Figure 1), to obtain records ranging from a few seconds to days. Instruments were synchronized before and after their deployment, and the time was subsequently corrected assuming a linear time clock drift between the two synchronizations (see details in [START_REF] Barreyre | Temporal variability and tidalmodulation of hydrothermal exit-fluidtemperatures at the Lucky Strike deepsea ventfield, Mid-Atlantic Ridge[END_REF]. Deployments details are given in Tables S1 andS2, and examples of short-term deployments are illustrated in Figure 6b. The temperature data processing strategy described below depends on the record length, and is reported in Table S1 and Table S2.

The thermal blanket and the sensor WHOI-HiTT-018 were deployed over two days and their records are thus considered as "long-term deployments". The record from sensor WHOI-HiTT-sensors showed mild temporal variations at the end of the record (~0.8 °C in amplitude), and this record section was excluded to calculate the average reported in Table S1.

We visually inspected the temperature records identified as "mid and short deployments" (several seconds to several hours) to determine if they reach a stable temperature (i.e., if thermal equilibrium with the subseafloor is achieved), or if temperature is still evolving at the end of the record (i.e., thermal equilibrium not reached). When temperature reaches equilibrium, we report the average of the last minute (Figure 5a, Figure 6b). For deployments still showing a mild temperature evolution but approaching thermal equilibrium and a stable temperature, we report the average over a shorter period of time at the end of the record, that is defined in each case based on temperature evolution. The averaging time windows vary between 5 to 30 seconds, and the resulting temperatures should be close to but slightly lower than the actual subseafloor temperatures (see details and values in Table S1 and Table S2). For deployments with significant temperature evolution at the end of the record, we report the last temperature value instead, which therefore corresponds to a minimum temperature estimate of the equilibrium temperature.

In addition to the subseafloor temperatures and the processing used to obtain them (see above), we include in Table S1 and Table S2 other observables that are key to correlate these temperatures to the patterns of hydrothermal outflow.

First, we provide positions of temperature measurements obtained from surface GPS positions.

These have errors associated with offsets between surface and seafloor positions, in addition to inherent GPS errors. GPS readings were obtained with hand-held GPS units (Garmin GPSMAP62), operated by a diver at the surface over the instrument deployment site. Tables S1 and S2 also report temperature measurement positions relative to the drone and AUV

Journal Pre-proof

photomosaics presented in this study. This position is determined using QGIS after locating each deployment on the drone and AUV photomosaics, and based on field observations and scuba diver imagery. Divers acquired both video and still imagery at deployments, allowing us to precisely identify and correlate structures observed in situ with those on the drone and AUV photomosaics. With the positions extracted from the photomosaics, we obtain a consistent geographic co-registration between temperature measurements, geological interpretations from the photomosaics, and in situ observations and measurements performed by divers.

Second, we also assigned to each deployment a seafloor type as classified from AUV and/or in situ observations (e.g., Figure 5b,5c). Reported seafloor types correspond to those identified in AUV photomosaics, as described above, in addition to transitional ones for measurements at boundaries between seafloor patterns, and to other areas associated with venting that cannot be identified in the AUV or drone mapping. First, we acquired temperature records within the White Hydrothermal Patches (WHPs) and Bioturbation (BT) zones. The limit between WHPs and Bioturbation zones, when adjacent, is often gradual and not sharply defined. We thus establish a transitional category (WHPs/BT). We also installed temperature sensors in sedimented seafloor unaffected by hydrothermal activity, that is classified as Non-hydrothermalized Sand (NS) from here on. We assign measurements to Hydrothermal Sand (HS) when the seafloor shows significant coloration (ochre or grey primarily) from hydrothermal precipitation that is not imaged in the drone or AUV photomosaics. Lastly, we report an additional category, Fluids (FL), for all measurements conducted directly at bubble outflows or in brine seeps, and therefore in areas of active and focused hydrothermal outflow.

Third, the subseafloor depth of the temperature measurements is determined from markers on temperature sensor rods, placed at known distances from their tips. Using field observations and

Journal Pre-proof underwater imagery acquired during deployments, we determine the depth of temperature measurements during deployments.

Fourth, Table S1 and Table S2 report the vertical temperature gradients acquired using the thermal blanket and the sets of instruments Poseidon, Set 1, Set 2, that in fine are required to calculate local heat fluxes. Thermal gradients are calculated using a least-square fit to the reported temperatures and their depths. Table S1 and Table S2 also include the goodness-of-fit to identify thermal profiles that are non-linear or show no clear thermal gradients with depth (Figure 6c). Table 5 reports the ranges and averages of temperatures and thermal gradients for each seafloor type. To calculate heat fluxes from these gradients, and to understand the relationship between conductive and advective processes, among other parameters we require constraints on the thermal conductivity of subseafloor sediment, which was not measured during this experiment.

Temperature at Paleochori Bay was measured at different subseafloor depths, depending on probe length and its insertion by the divers. As temperature gradients in the study area can be high (in many cases >10 °C/m, and as high as 400 °C/m), we discuss and interpret subseafloor temperature patterns based on temperature data normalized to a subseafloor depth of 35 cm. For sites occupied by the combined set of instruments ('Poseidon', 'Set 1' and 'Set 2'), we report the interpolated temperature at a depth of 35 cm. For individual temperature probe deployments lacking temperature gradients we used instead the maximum thermal gradients measured from each seafloor type. We then provide a minimum and maximum temperature at 35 cm subseafloor, calculated using these gradients and the measured temperature (Table 5, Table S1, Table S2).

Journal Pre-proof 4 Results

Seafloor textures

Drone photomosaics

The clearest seafloor structures visible on drone imagery are those associated with hydrothermal activity, in addition to biological patterns such as Posidonia oceanica meadows (Figure 3, Figure 7). Systematic mapping of different seafloor features on the drone photomosaics at the three study (Paleochori Bay, Agia Kyriaki, and Spathi Bay, see Figure 7) constrains surfaces and proportions of the different seafloor patterns (Table 2).

The WHPs cover ~5 % of the drone-surveyed seafloor at Paleochori Bay (Figure 7a), ~8 % at Agia Kyriaki (Figure 7b), and <2 % at Spathi Bay (Figure 7c). The DHPs cover a much lower proportion of the surveyed seafloor, and account for <1 % at Spathi Bay. These estimates exclude the Dark Lines visible at the three bays, often crossing and/or bounding WHPs (see Figure 3).

The DHPs reported in 2020) study. This study also identifies 'Sea Grass', that correlate well with some of the main Posidonia meadows we identify, but in other areas [START_REF] Khimasia | Hydrothermal areas, microbial mats and sea grass in Paleochori Bay, Milos, Greece[END_REF] reports features that are absent in our photomosaic. While some differences in the interpretation between studies may be attributed to

Journal Pre-proof criteria used in the classification, other differences may reflect instead temporal changes, that are out of the scope of this paper. For comparison, the 'Bacterial mats' identified by [START_REF] Khimasia | Hydrothermal areas, microbial mats and sea grass in Paleochori Bay, Milos, Greece[END_REF] correspond to slightly less than 2% of the mapped seafloor, about half of the surface that we identify here as WHPs. Similarly, the seafloor this prior study identifies as 'Sea Grass' is ~2%, somewhat lower than the 3% of Posidonia oceanica that we interpret in our drone data.

AUV photomosaics

The high resolution of AUV photomosaics, acquired at a few meters above the seafloor, provide imagery of higher resolution than the one from the drone, without the water-column filter and sea surface illumination artifacts. These AUV data adequately image the seafloor textures described above, along with additional ones of biological origin not identified on drone photomosaics.

While AUV-surveyed surfaces are much smaller than those from drone (see Figure 1), the overlapping imagery allows ground-truthing and facilitates the extrapolation of local AUV observations to broader areas covered by drone photomosaics. AUV surveys were conducted East and West of Paleochori Bay, and at Spathi Bay, with no surveys at Agia Kyriaki (Figure 8).

In Paleochori West (Figure 8b), WHPs of High-, Medium-and Low-density cover ~34 % of the surveyed area (Table 3). Here we also mapped broad areas of bioturbated seafloor (~32 %, Table 3). This seafloor pattern lacks clear bacterial mats or hydrothermal precipitates but its geometry and distribution show that it is often associated with hydrothermal activity: Bioturbated areas

Journal Pre-proof

In Paleochori East (Figure 8a), WHPs of High-, Medium-and Low-density cover a significant proportion of the surveyed seafloor (~43 %, Table 3). Bioturbation covers ~35 % of the surveyed seafloor (Table 3). As at Paleochori West, these areas often develop over Low-density WHPs and bound the High-and Medium-density WHPs, particularly in the region showing clear polygonal patterns (northeast of the mosaic, Figure 8a). Bioturbation is also found away from any recognizable hydrothermal seafloor pattern. Finally Cymodocea nodosa patches are distributed throughout the area, often bounding WHPs, but covering much less surface than at Paleochori West (<2.5 %).

At Spathi Bay (Figure 8c), High-and Medium-density WHPs are limited (<4 %, Table 3). We clearly identify seafloor areas that correspond to Dark Hydrothermal Patches (~1.7 %; Table 3), in addition to Dark Lines that frame the WHPs. Bioturbated areas (~7.2 %), both bound WHPs, and develop away from any hydrothermal patterns, as observed in Paleochori. AUV surveys at Spathi Bay also mapped a broad meadow of the seagrass Posidonia, covering here ~26 % of the surveyed seafloor.

Drone vs AUV mapping

Quantitative comparison of imagery and interpretations from drone and AUV data is necessary to determine if the patterns identified in these different datasets are equivalent despite the differences in imaging conditions, or cannot be directly correlated. This is also required to extrapolate local AUV data to regional scales. Table 4 provides AUV imagery also provides detailed maps where sampling and temperature measurements can be precisely located relative to seafloor structures, therefore enabling their interpretation.

Subseafloor temperature patterns

Based on mapping at Paleochori Bay, we assign a geological setting for each temperature measurement using 6 seafloor categories, as reported in Tables 5, S1 and S2, and in Figure 10 (see their description in Section 3.2). Temperature measurements in Paleochori were made at different subseafloor depths, depending on probe rod lengths. As there are significant temperature gradients throughout the study area (see Tables 5,S1, and S2, Figure 10), we provide temperatures at 35 cm subseafloor, to minimize any depth bias in spatial temperature patterns, facilitating their interpretation (see methods section, and full data, including temperature gradients, in Tables S1 andS2). We first present results from a study site displaying a polygonal

Journal Pre-proof hydrothermal pattern (Hexagon study site, Figure 11), and then results from other measurements throughout Paleochori Bay (Figure 12).

Hexagon study site

To establish the association between the different seafloor types and subseafloor thermal regimes,

we conducted a focus study on an area at Paleochori that displays a clear polygonal organization of hydrothermal outflow (Figure 11). In particular, we selected a hexagonal structure with white hydrothermal seafloor at its edges, referred to as the "Hexagon" from here on. We measured temperatures along six profiles, crossing the edges of this hexagon, gathering measurements at The temperature sensor set 'Poseidon' was positioned at the center of each profile, thus providing six vertical temperature gradients at the Hexagon (Figure 10b, Table S1). The lowest gradient of the area is measured in bioturbated (BT) seafloor (50.85 °C/m), while the WHPs display the Journal Pre-proof maximum gradient (94 °C/m). These thermal gradients are completed by values acquired over a broader study area (Paleochori Bay), as discussed below.

Paleochori Bay

While polygonal patterns are common at Milos (e.g. Paleochori or Spathi Bay), we studied other areas of interest that show structures with varying sizes, shapes, textures, compositions and that lack polygonal organization (Figure 12). Field work at several sites throughout Paleochori allows us to test if the subseafloor thermal regimes at the Hexagon and elsewhere show similar patterns or if the differences in outflow geometry, in addition to the seafloor texture, also control subseafloor temperatures, and result in different thermal regimes for similar textures. We thus acquired subseafloor temperature values at White Hydrothermal Patches (WHPs), Bioturbations (BT), Normal Sand (NS), Hydrothermal Sand (HS) and Fluids (FL). We note that these last two categories correspond to measurements at or near actively venting sites (gas, fluids, or both),

where small-scale variations in seafloor nature and texture can only be identified by in situ observations, and are hence not documented on the AUV or drone photomosaics.

As in the Hexagon, we document the same lateral thermal variation from WHPs to NS (Figure 10a, red symbols). Indeed, the mean temperature at 35 cm subseafloor and under WHPs is ~65 °C (~67 °C in average at the Hexagon). Under NS terrain these temperatures tend to decrease, with maximum temperatures of ~32 °C, and in average near-seawater temperatures (~28 °C vs. 21 °C at the Hexagon). The bioturbated terrain and the transitions to both WHPs and to NS terrain display intermediate values.

Vertical temperature gradients were measured at every site, using one or several sets of instruments ('Poseidon', 'Set 1' and 'Set 2', Figure 10b). Gradients for WHPs, BT and NS are all Journal Pre-proof <100 °C/m, as seen on the hexagon, and while the data is sparse and shows significant variability, it suggests a decrease in observed thermal gradient from WHPs towards NS terrain, mimicking the lateral variations in subseafloor temperature.

In addition to seafloor types observed in the photomosaics, we measured temperatures at hydrothermalized sediments, and around gas and fluid vent sites (HS and FL, Table 5 and Figure 10). These sites show extreme variability in both the temperature ranges (up to 100 °C and >150 °C for HS and FL respectively, Figure 10a) and the associated thermal gradients (up to 150 °C/m and >400 °C/m, respectively, Figure 10b). While most of the measurements in HS seafloor throughout Paleochori bay show gradients <100 °C/m, we identify a single site with >100 °C/m (Figure 10b). While we cannot establish a direct comparison with the results of [START_REF] Khimasia | Hydrothermal flux and porewater geochemistry in Paleochori Bay, Milos, Greece[END_REF], we have calculated both the gradient and the subseafloor temperature based on their data (Supplementary Table 2 in [START_REF] Khimasia | Hydrothermal flux and porewater geochemistry in Paleochori Bay, Milos, Greece[END_REF], and selecting the sites that show a temperature increase with depth. We observe that a) the subseafloor temperatures obtained in this study, and reported at 35 cm subseafloor, range from ~25 to ~100 °C, with a single value exceeding 100 °C, and b) that the gradients are typically <100 °C/m, with only two sites reporting gradients >100 °C/m. These values are thus in good agreement with our results, and suggest that the associations that we observe between seafloor patterns and subseafloor temperature can be generalized to the study area.

Discussion

Distribution and quantification of hydrothermal activity

The interpreted seafloor patterns, derived from drone imagery and ground-truthed by AUV photomosaics and in situ observations (Figure 2) characterize in detail the organization of Journal Pre-proof hydrothermal outflow through the seafloor. In particular, we document actively venting areas, their regional distribution through the study area (Figure 1), and characterize in detail their nature and geometry. WHPs are also found as broad and extensive patches, often with a gradual transition to sedimented seafloor, and without showing a clear geometric organization (e.g. limits, internal structure). We interpret these patches as evidence of widespread and diffuse outflow through the seafloor, lacking the focusing inferred for areas displaying polygonal patterns. Hence, the variations in the density of these patches (qualitatively classified as High, Medium, and Low density WHPs, see Figure 4), may correspond to spatial changes in either the hydrothermal fluid mass flux, temperature variations, or both, that are gradual and distributed. Our temperature measurements within these WHPs are not dense enough to verify this, and additional studies with in situ measurements are required instead.

In addition to the white reflective areas, drone imagery allows us to identify and map darker areas (DHPs) and lineations, with geometries and patterns that are similar to those of the WHPs. AUV imagery shows that these darker seafloor textures are of several possible origins. For example, the dark patches (DHPs) and lineations identified at Spathi Bay are sandy and seem to be closely associated with the WHPs spatially, in some instances appearing on their continuity. Similar sharp or more transitional change from white areas toward darker material (brown, green brown or gray) have been described (Aliani et al., 1998b, Sievert et al., 1999[START_REF] Bayraktarov | The pH and pCO2 dependence of sulfate reduction in shallow-sea hydrothermal CO2-venting sediments (Milos Island, Greece)[END_REF] Journal Pre-proof

showing sometimes complex patterns [START_REF] Wenzhöfer | In situ microsensor studies of a shallow water hydrothermal vent at Milos, Greece[END_REF]. These transitional areas before normal Non-hydrothermalized Sand are interpreted as change in bacterial population and type of precipitate containing manganese and iron oxide [START_REF] Sievert | Spatial heterogeneity of bacterial populations along an environmental gradient at a shallow submarine hydrothermal vent near Milos Island (Greece)[END_REF][START_REF] Wenzhöfer | In situ microsensor studies of a shallow water hydrothermal vent at Milos, Greece[END_REF][START_REF] Callac | Modes of carbon fixation in an arsenic and CO2-rich shallow hydrothermal ecosystem[END_REF]. In other areas these dark patches and lineations appear in zones that do not display any WHPs. While we did not investigate these areas systematically with temperature sensor deployments, some of the measurements indicate elevated temperatures (up to 95 °C at 35 cm, see Table S1 and Table S2) suggesting that at least some of these zones are also active. This is also confirmed by recent in situ studies over a hydrothermal patch at Spathi Bay [START_REF] Callac | Modes of carbon fixation in an arsenic and CO2-rich shallow hydrothermal ecosystem[END_REF], suggesting that these are due to goethite precipitates. In situ observations also display seafloor areas with ochre hydrothermal precipitates (see Section 2), that may appear as dark seafloor in both AUV and drone imagery. Hence, such dark patterns may correspond to fossil, inactive areas of hydrothermal outflow, or to active areas with hydrothermal precipitates yielding a darker seafloor pattern, often adjacent to bacterial mats or white hydrothermal precipitates (e.g., [START_REF] Callac | Modes of carbon fixation in an arsenic and CO2-rich shallow hydrothermal ecosystem[END_REF]. These types of seafloor patterns, with both visually white seafloor areas and darker seafloor structures, have been identified in deep-sea hydrothermal systems, and interpreted as active and inactive hydrothermal outflow areas respectively (e.g., Lucky Strike, [START_REF] Barreyre | Structure, temporal evolution, and heat flux estimates from the Lucky strike deep-sea hydrothermal field derived from seafloor image mosaics[END_REF].

AUV imagery also reveals dark patches and lineations corresponding to the seagrass Cymodocea nodosa, at both Paleochori West and East (Figure 8). While we have limited AUV coverage, our observations suggest that these are closely associated with the WHPs, often surrounding them (e.g., Figure 8a andb). These observations are consistent with the spatial associations discussed by Aliani et al. (1998b) between the distribution of seagrass and their link to hydrothermal activity.

To discriminate between a biological origin (seagrass, or debris of seagrass, Figure S5) or a hydrothermal one (active or inactive hydrothermal outflow zones and associated recent or altered precipitates), and in the absence of direct observations or AUV surveys, drone photomosaics should be acquired at different times, as debris are mobile, and thus likely associated with transient patches. In contrast, we speculate that dark areas of hydrothermal origin may be more persistent over time, as suggested by the preserved geometry of the hydrothermal outflow in the area over several years, documented by the imaging of WHPs through satellite imagery [START_REF] Martelat | Terrestrial shallow water hydrothermal outflow characterized from out of space[END_REF], and by their association to white patches [START_REF] Callac | Modes of carbon fixation in an arsenic and CO2-rich shallow hydrothermal ecosystem[END_REF].

Quantitative comparisons of seafloor types visible on both drone and AUV photomosaics show that surface estimates are equivalent using both image datasets (Figure 9, Table 4). Minor differences observed in our surface estimates may have multiple causes. First, we may document temporal, local changes between the acquisition of the drone and AUV imagery. Second, imaging conditions differ significantly between the AUV (near-seafloor acquisition, with images acquired a few meters above seafloor), and the drone (imaging through the full water column). Drone imagery is strongly influenced by several external factors, such as sun illumination angle, sea surface reflections or water turbidity (e.g., [START_REF] Martelat | Terrestrial shallow water hydrothermal outflow characterized from out of space[END_REF]. AUV imagery may also be affected by flickering (light refraction by surface waves), among others. Finally, drone imagery is acquired through a water column that varies in thickness from near-zero (shore) to up to 15-20 m in the deepest parts. These factors may in part explain the differences observed between a prior survey [START_REF] Khimasia | Hydrothermal areas, microbial mats and sea grass in Paleochori Bay, Milos, Greece[END_REF], owing to the combination of differences in imaging conditions and temporal changes between the drone surveys. Hence, with drone data we expect a loss of resolution and definition that increases away from the shore. Despite these factors, our

Journal Pre-proof comparison indicates consistency across datasets, and that we can confidently use broader drone imagery and its interpretation to extrapolate local AUV observations.

Links between seafloor patterns and subseafloor temperatures

Subseafloor temperature measurements display spatial variations that we can directly link to seafloor hydrothermal patterns, and thus provide information on subseafloor organization of hydrothermal flow. Subseafloor temperatures and gradients reported in other areas of Paleochori [START_REF] Khimasia | Hydrothermal flux and porewater geochemistry in Paleochori Bay, Milos, Greece[END_REF], while they cannot be directly correlated to specific seafloor textures in our study, are consistent with our results, and suggest that these lateral temperature gradients, at scales of a few to ~10m, are well established and consistent throughout. Gradients are controlled by the temperature decrease from the WHPs (average of ~67 °C and reaching ~75 °C at some sites) to the sediments lacking hydrothermal alteration, NS (temperatures between those of seawater and ~30 °C), with transitional values in bioturbated areas, BT, that bound the WHPs (Figure 10). These temperature gradients away from individual WHPs, documented here at the "hexagon" and at other sites around Paleochori Bay (Figure 11, Figure 12), are also described at Paleochori Bay by prior studies, at subseafloor depths of up to 10 cm, and with maximum measured temperatures of ~100 °C instead (e.g. Figure 3 in [START_REF] Wenzhöfer | In situ microsensor studies of a shallow water hydrothermal vent at Milos, Greece[END_REF]. Such significant lateral gradients over short spatial scales require active fluid circulation, with focused hydrothermal flow and discharge at the WHPs seafloor, and a possible recharge at nonhydrothermalized sediments, either at the center of polygonal structures, or around the isolated and disconnected patches elsewhere in Paleochori. We anticipate that this fluid flow organization may be common at sedimented shallow-water systems showing polygonal patterns, such as Spathi and Voudia Bays in Milos, or Vulcano in Italy [START_REF] Martelat | Terrestrial shallow water hydrothermal outflow characterized from out of space[END_REF], even though no temperature data are available from these sites.

Journal Pre-proof

The hydrothermal patterns visible in AUV and drone imagery and described above are associated with areas showing gas venting (instrumented sites indicated as FL in Table S2), or sands

showing hydrothermal deposits of colors other than white and only identified with in situ scuba diver observations (areas categorized as HS). Clearly associated with hydrothermal outflow, these hydrothermalized deposits can only be identified by in situ divers and are only partially imaged by the drone and AUV. Subseafloor temperatures at these sites display significant variations (ranging from ~30 °C to ~95 °C) and cannot be correlated with the horizontal thermal gradients described above. The temperature measurements assigned to the 'Fluids' category regroup measurements in both degassing outflow and brine seeps. Focused venting likely results in temperature variations at scales similar to the venting sites (a few cm to <1 m), and may explain the extreme variability observed in our data. To better understand these structures, subseafloor temperatures should be acquired in a systematic manner with dense sampling to characterize and constrain the spatial variations at scales of <1 m.

Overall, the thermal gradients in the shallow subseafloor are extremely elevated, often above 100 °C/m below the WHPs. These values point to a thin conductive layer (up to 1-2 m) that corresponds to the cap of sediment saturated with hydrothermal fluids and/or gases at ~100 °C below WHPs, either in zones that are narrow and focused (e.g., small patches or edges of polygonal patterns), or over broader areas underlying the large patches instead. Furthermore these gradients suggest efficient conductive cooling near-seafloor of the upwelling high-temperature fluids, a mixture with cool seawater, or a combination of both. This is consistent with a porous flow of hydrothermal fluids at very slow speeds, with no visual expression at the seafloor, where there is conduction in the shallow subseafloor overimposed on an advective system. Figure 13 shows a schematic representation of subseafloor hydrothermal circulation, depicting the lateral Journal Pre-proof thermal gradients, with local, focused gas discharge to the seafloor, that is decoupled from the surrounding hydrothermal flow at shallow, near-seafloor levels.

Hydrothermal patterns impact local ecosystems

Submarine hydrothermal outflow sustains chemosynthetic ecosystems, such as microbial mats (e.g., [START_REF] Jannasch | Morphological survey of microbial mats near deep-sea thermal vents[END_REF][START_REF] Barreyre | Structure, temporal evolution, and heat flux estimates from the Lucky strike deep-sea hydrothermal field derived from seafloor image mosaics[END_REF], and associated macrofauna. At Milos, AUV imagery reveals a strong hydrothermal control on macrofauna distribution. As seen on the AUV photomosaics, bioturbation is abundant, clearly recognizable in the imagery, and attributed to the burrowing activity of the mud shrimp Callianassa truncata [START_REF] Ziebis | Cymodocea nodosa 11.50 % 2.33 % Posidonia oceanica 25.94 % Drone AUV Spathi -DHPs 456.11 m2 463.12 m2 Spathi -Posidonia[END_REF][START_REF] Robinson | Situ Investigations Of Shallow Water Journal Pre-proof Hydrothermal Vent Systems, Palaeochori Bay, Milos, Aegean Sea[END_REF]. This burrowing continuously remobilizes sediment, obliterating the ripple marks that develop at the seafloor. We identify two different settings for this bioturbation. The most common one corresponds to bands up to a few meters wide that frame WHPs, and is seen in almost all of the polygons and patches imaged by AUV (Figures 8 and13). The second one corresponds to broad areas extending over several tens of square meters that are not adjacent to

WHPs. Associations between the distribution of hydrothermal activity and other type of macrofauna have also been proposed, but could not be verified with our AUV mapping. This is the case for sessile epibenthos at Milos around vents in rocky areas [START_REF] Bianchi | Increased diversity of sessile epibenthos at subtidal hydrothermal vents: seven hypotheses based on observations at Milos Island[END_REF], and for fish and crabs (Aliani et al., 1998b;Gamenick et al., 1998).

Our subseafloor temperature measurements suggest that bioturbated zones may map a thermally optimal habitat for C. truncata, often at the transition between the WHPs (~70 °C) and the sedimented seafloor with no hydrothermal activity (NS, 25 °C). In addition, the nearby bacterial mats may provide a source of nutrient and carbon for the shrimps, as proposed for other organisms (e.g., [START_REF] Bianchi | Increased diversity of sessile epibenthos at subtidal hydrothermal vents: seven hypotheses based on observations at Milos Island[END_REF]. The bioturbation zones found in areas away from WHPs Journal Pre-proof may indicate diffuse outflow insufficient to support bacterial mat growth. There is to our knowledge no description of the upper thermal tolerance of C. truncata. However, it seems unlikely that it can survive the temperatures higher than 50 °C found in WHPs, suggesting a temperature threshold between the bioturbated areas and the WHPs. The temporal stability of the hydrothermal system is well documented by satellite imagery [START_REF] Martelat | Terrestrial shallow water hydrothermal outflow characterized from out of space[END_REF], and is also supported by a qualitative comparison of our drone data with that from an earlier survey [START_REF] Khimasia | Hydrothermal areas, microbial mats and sea grass in Paleochori Bay, Milos, Greece[END_REF]) that shows only minor and local differences in the distribution of hydrothermal patches and structures. However, new AUV surveys will be critical to verify if the zones burrowed by Callianassa truncata are persistent in the long term, or if these ecosystems show significant spatial variations over time that may be linked to changes in subseafloor temperature.

We also observe that Posidonia oceanica only grows south and outside of Paleochori Bay, while the seagrass Cymodocea nodosa is found in the vicinity of the hydrothermal vents and around patches. This is consistent with the proposed interpretation by Aliani et al. (1998b) suggesting that C. nodosa has a higher thermal tolerance than P. oceanica, and that both are absent in areas with elevated subseafloor temperatures (Aliani et al., 1998).

Bacterial mats at the seafloor are the most visible ecosystem, and the one used here to infer the geometry and distribution of active hydrothermal activity. The microbiology of Paleochori Bay has been the focus of several studies in the past decades (e.g. [START_REF] Brinkhoff | Distribution and diversity of sulfur-oxidizing Thiomicrospira spp. at a shallow-water hydrothermal vent in the Aegean Sea (Milos, Greece)[END_REF], Giovannelli et al., 2013, Price et al., 2013). These studies have shown the presence of multiple microbes involved in the sulfur cycle. However, most of these studies used now outdated methods, and no investigation using next-generation sequencing technology has been published yet. As a follow up of this study and the fieldwork conducted in 2019, there is a study in progress providing a

Journal Pre-proof detailed description of the microbial communities in relation to the different subseafloor thermal regimes and seafloor patterns described here.

Photomosaics from the Milos hydrothermal site document the distribution of ecosystems associated with hydrothermal circulation, and may be used for quantitative studies. For example, the areas bioturbated by C. truncata, mapped in AUV surveys, cover a surface that is commeasurable to that of WHPs at three different locations (Table 4). We thus expect that this terrain and the associated ecosystem is widespread throughout our study area, covering significant portions of the seafloor and surrounding WHPs areas identified by our drone mapping.

Hence, detailed studies of hydrothermal ecosystems, including their distribution, nature, and quantification, require systematic seafloor mapping with underwater imagery (i.e. AUVs) in addition to in situ observations (e.g., scuba divers) that cannot be used for quantitative studies.

Aerial surveys do provide a regional view of these systems, including predictive subseafloor temperature patterns and associated overall hydrothermal outflow, but fail to capture the details of these structures (e.g., certain macrofauna ecosystems, zones of focused fluid and gas venting, etc), in addition to being limited to near-shore environments (<20 m water depth).

Perspectives on heat fluxes

The shallow water hydrothermal system of Milos is complex, extending over a significant area to depths well below those studied here, and displaying different modes of hydrothermal outflow.

Systematic mapping and quantification of these systems is needed to estimate in fine associated fluxes (mass, heat, chemical output). While we focused our study on systematic seafloor imaging coupled with subseafloor temperature measurements, the results presented here provide the background for future studies and fieldwork to constrain the heat flux of this system. The total Journal Pre-proof heat flux of a hydrothermal system such as Milos is partitioned between conductive and advective fluxes (Figure 13). Pure conductive heat fluxes can be enhanced by the proximity of a heat source, and at Milos we anticipate a subseafloor saturated in high-temperature fluids at a shallow level, with a thin conductive layer, as required by the significant thermal gradients observed. Advective fluxes are partitioned into diffuse and focused flow (Figure 13), with both fluid and gas venting.

The WHPs mapped likely correspond to diffuse fluid percolation at small spatial scales (<1 m), while the focused flow is observed at gas venting sites and brine pools, where higher temperatures with extreme variability are observed (Figure 10). Fluxes associated with focused flow are difficult to evaluate, as the number of venting sites, their distribution throughout the study area, and the associated characteristics (temperature, mass flux, outflow velocity) are difficult to estimate systematically, as discussed below.

Heat flux from diffuse hydrothermal flow has been estimated in deep-ocean systems (e.g. [START_REF] Barreyre | Structure, temporal evolution, and heat flux estimates from the Lucky strike deep-sea hydrothermal field derived from seafloor image mosaics[END_REF], and is calculated from:

𝑄 = ∆𝑇 * 𝜌 * 𝜈 * 𝐶𝑝 * 𝑆
where ∆𝑇 is the temperature difference between outflow temperature (38 °C) and ambient seawater (24 °C), 𝜌 is the density of seawater, 𝜈 is the flow velocity (m/s), 𝐶𝑝 is the specific heat of the diffuse fluid and S is the considered area of diffuse outflow (m 2 ). Here we can estimate ∆𝑇 from the temperature dataset and S from the photomosaics (Table 2), while assuming values for 𝜌 (1025 kg.m 3 ) and 𝐶𝑝 (4.2 * 10 3 J.kg -1 .°C -1 ) (see [START_REF] Barreyre | Structure, temporal evolution, and heat flux estimates from the Lucky strike deep-sea hydrothermal field derived from seafloor image mosaics[END_REF]. However, despite local studies of gas flow rates (e.g. [START_REF] Botz | Hydrothermal gases offshore Milos Island, Greece[END_REF][START_REF] Megalovasilis | Geochemistry of Hydrothermal Particles in Shallow Submarine Hydrothermal Vents on Milos Island, Aegean Sea East Mediterranean[END_REF], we lack a precise and consistent estimate of the fluids velocity, their composition and nature (and hence density) and

Journal Pre-proof their variability, that could be extrapolated regionally. To obtain a flux estimate, we assume a diffuse velocity flow through the WHPs between a minimum and a maximum of 1 and 5 mm/s respectively, that is based on diffuse flow studies through bacterial mats elsewhere [START_REF] Sarrazin | A dual sensor device to estimate fluid flow velocity at diffuse hydrothermal vents[END_REF][START_REF] Barreyre | Structure, temporal evolution, and heat flux estimates from the Lucky strike deep-sea hydrothermal field derived from seafloor image mosaics[END_REF]. Table 6 provides values for each study area and WHPs type.

Focusing on Paleochori Bay, the best studied area, our preliminary diffuse heat flux estimates vary between ~2350 and 11760 MW (Table 6). We note that [START_REF] Khimasia | Hydrothermal flux and porewater geochemistry in Paleochori Bay, Milos, Greece[END_REF] Journal Pre-proof sites, nor to determine any temporal variability (e.g., tidal). Finally, our preliminary observations do not show a clear spatial association between the position of venting sites and seafloor structures mapped in photomosaics, as gas venting is observed in areas of sedimented seafloor (NS), zones of bioturbation, and in WHPs. While temperature measurements at a few gas venting sites are elevated (close to ~100°C), and with highly variable gradients (Table 5 and Figure 10), we cannot relate gas outflows to the systematic patterns of subseafloor temperature associated with bacterial mats.

Our observations and results suggest that in the shallow subseafloor, and at small spatial scales (a few 10s of m or less) gas flow is decoupled from diffuse fluid flow, and hence seafloor patterns cannot be used to predict or quantify gas outflow location and density. For example, [START_REF] Botz | Hydrothermal gases offshore Milos Island, Greece[END_REF] report surface gas flow rates of up to ~200 ml/m 2 /min while recently [START_REF] Megalovasilis | Geochemistry of Hydrothermal Particles in Shallow Submarine Hydrothermal Vents on Milos Island, Aegean Sea East Mediterranean[END_REF] reports gas and fluid flow rates at individual sites ranging from 4.5 to 63 l/h (gas) and 25 to 92 l/h (water), with temperatures ranging between 30 °C and 110 °C. We conducted an experiment at Paleochori at two sites, and obtained two gas flow rates values of 18 l/h and 360 l/h, the latter significantly higher than the values reported by [START_REF] Megalovasilis | Geochemistry of Hydrothermal Particles in Shallow Submarine Hydrothermal Vents on Milos Island, Aegean Sea East Mediterranean[END_REF]. 6). Heat flux associated with the gas outflow is negligible as the heat capacity of gases is two orders of magnitude lower than that of seawater.

Journal Pre-proof

As thermal gradients are significant in the shallow seafloor, part of the heat loss is also conductive, associated to an enhanced and thin conductive layer. With a mean thermal gradient of ~64 °C m -1 for the WHPs (Table 5), and a typical thermal conductivity for sediments (1.7 W.m - 1 °C-1 ), the conductive heat flux through all the WHPs at Paleochori would correspond to ~4 MW (109 W/m 2 ), which is significantly lower than our diffuse heat flow estimates.

While poorly constrained, these heat fluxes provide a first-order view into heat flux partitioning at Milos. Prior studies at deep-sea hydrothermal sites indicate that diffuse flow may account for at least ~50 % and up to more than 90 % (e.g., [START_REF] Veirs | Magnitude and variance of near-bottom horizontal heat flux at the Main Endeavour hydrothermal vent field[END_REF][START_REF] Ramondenc | The first measurements of hydrothermal heat output at 9°50′N, East Pacific Rise[END_REF][START_REF] Barreyre | Structure, temporal evolution, and heat flux estimates from the Lucky strike deep-sea hydrothermal field derived from seafloor image mosaics[END_REF][START_REF] Escartin | Hydrothermal activity along the slow-spreading Lucky Strike ridge segment (Mid-Atlantic Ridge): Distribution, heatflux, and geological controls[END_REF]. The preliminary heat flux estimates presented here lack additional constraints, such as temperature measurements, knowledge of physical properties of sediments (e.g., conductivity of sediments, nature and density of fluids, outflow velocities), or chemical indicators, such as such as patterns of Rd concentration in hydrothermal fluids and in the water column (e.g., [START_REF] Rudnicki | Theory applied to the Mid-Atlantic Ridge hydrothermal plumes: the finite-difference approach[END_REF]. As both the heat fluxes associated with focused flow (limited to the vents investigated by [START_REF] Megalovasilis | Geochemistry of Hydrothermal Particles in Shallow Submarine Hydrothermal Vents on Milos Island, Aegean Sea East Mediterranean[END_REF] and the conductive heat flow are only a fraction of the lowest estimate of diffuse heat flux, we postulate that heat loss at the Milos shallow water hydrothermal site is dominated by diffuse flow. This partition may also prevail at deeper water depths that have not been systematically surveyed yet [START_REF] Nomikou | Submarine volcanoes along the aegean volcanic arc[END_REF]. As we have systematically mapped only to water depths of ~20 m, and with hydrothermal activity to depths of 500 m or more, we anticipate that this is one of the most extensive and significant underwater hydrothermal systems identified to date. For comparison, heat fluxes at hydrothermal systems along the Mid-Atlantic Ridge such as Lucky Strike (~1000 MW, [START_REF] Barreyre | Structure, temporal evolution, and heat flux estimates from the Lucky strike deep-sea hydrothermal field derived from seafloor image mosaics[END_REF][START_REF] Escartin | Hydrothermal activity along the slow-spreading Lucky Strike ridge segment (Mid-Atlantic Ridge): Distribution, heatflux, and geological controls[END_REF] or Rainbow (1000-5000 MW, [START_REF] Thurnherr | Hydrography and high-temperature heat flux of the Rainbow hydrothermal site (36 14′ N, Mid-Atlantic Ridge)[END_REF] Journal Pre-proof are within the range of diffuse flow estimates at Paleochori Bay at depths <20 m, which is a small fraction of the whole Milos system (additional bays and hydrothermal activity at depths >20 m).

To document the full extent of this hydrothermal system and fully evaluate the fluxes, we require systematic mapping and surveying coupled to sampling, and conducted with underwater robots.

This mapping will also allow us to establish the links to ecosystems, and the controls on outflow geometry. Future studies will seek to constrain the thickness of the permeable layer hosting this fluid flow organization with subseafloor geophysical surveys, acoustic mapping of the water column to systematically map gas plumes, in addition to physical properties of sediments (thermal conductivity, porosity, permeability).

Conclusions

Drone imagery, ground-truthed by AUV photomosaics, provides a detailed geological context for in situ samples, measurements and observations. With these data and observations we establish a regional overview of the shallow-water hydrothermal system along the SE coast of Milos and to depths of ~20 m. We constrain not only the organization of its hydrothermal outflow, but also its subseafloor flow, the effects of this hydrothermal activity on both macrofauna and microbial ecosystems, and discuss the magnitude and partition of associated heat fluxes.

Photomosaics are interpreted to constrain the distribution and quantification of hydrothermal activity. At Milos, bacterial mats and hydrothermal mineral precipitates show up as White Hydrothermal Patches with elevated subseafloor temperatures (>50 °C and up to 75 °C) that are organized in polygonal patterns, isolated patches, or in more broad and dispersed areas. These patches are systematically bound by bioturbated areas, with subseafloor temperatures of (~30-40°C), that transition to sandy seafloor lacking any hydrothermal signature, and with subseafloor temperatures close to those of seawater (typically <32 °C). These gradients both in seafloor texture and subseafloor temperature document the position of persistent hydrothermal outflow, that is diffuse and associated with low outflow velocities through the seafloor.

This hydrothermal activity, and its organization, impacts both the distribution of bacterial mats at the seafloor, and of macrofauna. Indeed, the lateral thermal gradients are exploited by a burrowing shrimp, Callianassa truncata, as this narrow band of seafloor is an optimal environment thermally close to white hydrothermal areas, a likely source of nutrients. The burrowing activity of this shrimp continuously reshapes the seafloor, obliterating its otherwise rippled morphology, forming bioturbation, a seafloor texture that can be identified in nearseafloor AUV imagery.

Hence, our results show that drone photomosaics can be used to link seafloor hydrothermal patterns and subseafloor temperature data in order to predict subseafloor hydrothermal circulation regionally, to water depths of ~20-30 m. AUV imagery provides visual ground-truthing that is critical to identify structures and patterns that cannot be observed in aerial images, hindered by water column turbidity and attenuation and surface illumination artifacts such as sea surface reflections. This AUV imagery also allow us to establish the associations between ecosystems and hydrothermal flow, and to predict their distribution regionally based on drone imagery. This systematic mapping of the shallow-water hydrothermal system of Milos is also a first step to evaluate the thermal, mass, and chemical fluxes in the area and to evaluate the controls of hydrothermalism on microbial communities subseafloor, both spatially and in depth. Our preliminary analysis to establish heat fluxes also suggest that diffuse flow accounts for a significant portion of the total heat output in the surveyed area, but heat and mass flux quantification requires future dedicated experiments, with adequate constraints on the physical properties of the sediments hosting the hydrothermal flow. We postulate that Milos is one of the most active underwater hydrothermal systems on Earth, based on the fluxes that we estimate are associated with its shallow, near-shore portion, and its extension to deep waters (>500m) and throughout the island, in areas yet unsurveyed. S2). Red portions of the curve at each deployment correspond to the data used to calculate average temperatures reported in Table S2. The type of seafloor at each site is also 5, Table S1, and Table S2). WHPs: White Hydrothermal Patches. WHPs/BT: Transitional area. BT:

Journal Pre-proof

Bioturbation. NS: Non-hydrothermalized Sand. HS: Hydrothermalized Sand. FL: Fluids. Journal Pre-proof Table 1. Details and specifications on the temperature sensors used to acquire the data presented in this paper.

Table 2. Surfaces for all seafloor types identified on drone photomosaics.

Table 3. Surfaces for all seafloor types identified on AUV photomosaics.

  bound and outline High-and Medium-density WHPs, and often develop in association with Lowdensity WHPs. Additionally, light patches of the seagrass Cymodocea nodosa, covering ~11.5 % of the Paleochori West photomosaic, often bound WHPs.

  the comparison of surfaces for equivalent AUV and drone-derived seafloor patterns at Paleochori West, Paleochori East and Spathi. These reported surfaces correspond to those from the AUV surveys on the one hand, and to those interpreted on the aerial drone photomosaics but solely over the surfaces surveyed by the Journal Pre-proof AUV, on the other hand. In Paleochori West and East we compared mapped surfaces for WHPs with all density combined. In Spathi bay we compared the surfaces of DHPs and Posidonia oceanica, in addition to the WHPs. As shown in Figure 9, the surfaces of interpreted seafloor patterns from both AUV and drone show extremely consistent values, with a slight overestimation of Posidonia oceanica surface on the drone imagery relative to the AUV, and a slight underestimation of the WHPs at Paleochori East, also on the drone imagery. This excellent agreement between measured surfaces of seafloor textures identified both in drone and AUV imagery demonstrates that these seafloor patterns (WHPs, Posidonia, DHPs) are equivalent, and that drone imagery can be used for accurate mapping regionally. At the same time, AUV groundtruthing reveals other seafloor textures not imaged by the drone (Bioturbation, Cymodocea patches).

  WHPs, transition from WHPs to BT (WHPs/BT), Bioturbation (BT), and Non-hydrothermalized Sand (NS)(Figures 10 and 11

  ).As shown by the blue symbols in Figure10a, the highest temperatures are recorded below the WHPs, at the edges of the Hexagon (~67 °C in average, maximum temperatures of ~73 °C), with a progressive and gradual radial decrease towards the Non-hydrothermalized Sand, in adjacent polygons and at the center of the hexagon. These NS areas display seawater-like temperatures (21.82 °C in average), reaching up to ~33 °C in close proximity to the transition towards the bioturbated terrain. In addition to these lateral thermal gradients, we document a limited data spread of <~20 °C (standards deviations of ~5.25-4.10 °C). Subseafloor temperatures thus show consistent lateral variations around the edges of this polygon, that are well-correlated with seafloor types.

  Areas of active hydrothermal outflow are highlighted by White Hydrothermal Patches(WHPs), which are associated with elevated subseafloor temperatures (reaching up to ~75 °C). Mapped both on drone and AUV photomosaics, WHPs are necessarily associated with active hydrothermal flow of high-temperature fluids through the seafloor, required to sustain the white mineral precipitates and/or bacterial mats. Our mapping documents two clearly different organizations of the outflow, that we attribute to different modes in the subseafloor fluid circulation. First, some areas display clear polygonal patterns (e.g., Figure8a, Figure11), with spatial scales that vary from ~5 m to ~50 m. This type of organization develops in sandy, mobile seafloor, that blankets most of the seafloor at Paleochori, Spathi and Agia Kyriaki bays. The center of these polygons corresponds systematically to sediment showing no hydrothermal signature, with clear and pervasive ripple marks (Figure4, FigureS4). This type of outflow organization is clearly visible in Paleochori East (Figure8a) and in Spathi Bay, although with smaller polygons (<20 m, see Figure8c). In addition to the well-defined and continuous polygonal patterns, numerous isolated patches throughout Paleochori show alignments defining discontinuous polygon edges. This polygonal organization of hydrothermal outflow, together with the systematic temperature patterns observed across them (Figure10), clearly point to a well-structured subseafloor flow organization. These seafloor geometries are reminiscent of polygonal patterns recognized during fluid convection in sandy porous media heated from below (e.g.,[START_REF] Nield | Convection in Porous Media[END_REF]. While polygonal convection in the laboratory can result both in upwelling and downwelling at polygon Journal Pre-proof edges depending on experimental conditions, our temperature patterns suggest upwelling below WHPs, while cell centers are downwelling zones and hence at seawater temperatures. The geometry and size of the cells depend on both the permeability and thickness of the porous layer, and on hydrothermal fluids properties. However, to date, we lack constraints on the physical properties of the sediments in the area (porosity and permeability, in addition to thermal conductivity), and on the thickness of the sediment layer at Paleochori and other bays, likely underlain by rocky basement.

  report a mass flux of ~11300 ± 1100 m3/day, based on the extrapolation of vertical velocity of fluids based on a model of advection and cooling from subseafloor temperature measurements that we cannot colocate systematically with hydrothermal patterns at the seafloor, which yield values of 0.4-0.95 m/day (Fig. 11 in Khimasia et al., 2021), or .004-.01 mm/s. With these values, the advective heat flux would be negligible (a few tens of MW or less). This demonstrates that additional studies are required to better constrain fluid flow values in a systematic manner, and taking into account seafloor patterns of hydrothermal outflow. Heat flux from focused flow contributes further to the advective fluxes of the area, and depends on the fluid outflow temperature, the mass flux, and the abundance of venting sites (spatial distribution, size and mass fluxes). Constraints on this advective flux thus require a systematic mapping of the venting sites (or a representative subset), coupled to flux and temperature measurements and fluid sampling. For example, Figure 8a shows a detail of Paleochori displaying polygonal structures imaged by the AUV. Careful examination of the photomosaics allows us to identify gas venting sites as bubble trails emanating from the seafloor. When gas flow rate is low, trails are not captured, as individual bubbles imaged in the water column cannot be traced to a seafloor source. In situ observations by divers are limited to a few sites, and are insufficient to establish a quantitative map on distribution and abundance of gas and fluid venting

  Lacking a clear understanding of the density, distribution, flow rate, and nature (fluid and gas composition) of focused gas and fluid outflow (vents), we cannot extrapolate of these values regionally to infer well constrained fluxes. To evaluate qualitatively the relative importance of diffuse vs. focused heat flux, we use values provided by Megalovasilis (2020) for ~50 vents, with average outflow temperature of 93 °C (∆T=69 °C) and average flow rate of ~58 l/h. These surveyed sites yield a combined heat flux of ~0.2 MW, and thus a very small fraction of the diffuse heat flow estimate (~2350 -11760 MW, Table

Figure 1 .

 1 Figure 1. (A) Location map showing drone and AUV coverage. Orange square locates Figures 4e,

Figure 2 .

 2 Figure 2. Example of a hydrothermally active area at Paleochori Bay, with the same

Figure 3 .

 3 Figure 3. Details of the uninterpreted (left) and interpreted (right) drone photomosaics for Spathi

Figure 4 .Figure 5 .

 45 Figure 4. Details of the uninterpreted (top) and interpreted (bottom) AUV photomosaics for

  indicated attached to the deployment name (WHPs: White Hydrothermal Patches. BT: Bioturbation. HS: Hydrothermalized Sand). The records show clearly the probe insertion at the seafloor, followed by a sudden increase of temperature, and the extraction by sudden drops. (B) and (C) show the field image of the Poseidon instrument installed at the seafloor, in bioturbated seafloor (B) and in a white hydrothermal patch (C), showing that these correspond to a thin layer (bacterial mat, hydrothermal precipitates) covering the otherwise sandy seafloor.

Figure 6 .

 6 Figure 6. (A) Set of NKE temperature probes ("Poseidon"), deployed at the seafloor (scuba diver

Figure 7 .

 7 Figure 7. Drone photomosaics (background) with interpreted seafloor textures and features from

Figure 8 .Figure 9 .

 89 Figure 8. AUV photomosaics over the drone imagery (background image), and associated

Figure 10 .

 10 Figure 10. (A) Subseafloor temperature at 35 cm depth observed for the different seafloor types

Figure 11 .

 11 Figure 11. Location of temperature measurements on the hexagon area, with corresponding

Figure 12 .

 12 Figure 12. (A) Location of temperature measurements over the Paleochori Bay drone

Figure 13 .

 13 Figure 13. 3D interpretative sketch of the organization of the subseafloor thermal regime, with

Table 4 .

 4 Seafloor surfaces estimated from AUV surveys and from aerial drone surveys calculated over the areas surveyed by the AUV.

Table 5 .

 5 Summary of subseafloor temperatures (at 35 cm) and vertical thermal gradients at Paleochori Bay, including the Hexagon study site. Number of observations in parenthesis. WHPs: White Hydrothermal Patches. WHPs/BT: Transitional area. BT: Bioturbation. NS: Nonhydrothermalized Sand. HS: Hydrothermalized Sand. FL: Fluids.*Negative measured values are likely associated with advection and complex thermal regimes, not conductive regimes.

Table 6 .

 6 Estimates of heat flux in MW for each bay. Estimates are provided for both diffuse and conductive heat flow, based on the surfaces in Table 2 corresponding to WHPs surfaces. For diffuse flow the minimum and maximum values correspond to outflow velocities of 1 and 5 mm/s. See text for the parameters used in the calculations including thermal conductivity, heat capacity, and density. indicate measurements with no depth data available, and the nominal rod length is used instead. Light green cells indicate records where the reported temperature does not correspond to the average of the last minute, as reported and explained in the text. WHPs: White Hydrothermal Patches. WHPs/BT: Transition between WHPs and bioturbated areas (BT). NS: Nonhydrothermalized Sand.

	Journal Pre-proof

cells

Table S2 .

 S2 Master table of temperature data acquired throughout Paleochori Bay. Values are sorted and reported by dive, probe, and deployment. Absolute coordinates are based on GPS data. Relative coordinates are based on GIS maps after placing measurements on the photomosaics. Time in hours is in UTC, with time=0 corresponding to 25/09/2019 at 00:00. Start and end of deployments were handpicked to readjust uncertainties from GPS time that might influence the final temperature values. reported correspond to the mean of the last minute of deployment, when the record is most stable Red cells indicate measurements with no depth data available, and the nominal rod length is used instead. Light green cells indicate records where the reported temperature does not correspond to the average of the last minute, as reported and explained in the text (usually 30 seconds). Light red cells indicate records where the highest temperature is reported instead, and corresponding to short deployments where the temperature record was far from steady state, as indicated in the text. WHPs: White Hydrothermal Patches. WHPs/BT: Transition between WHPs and bioturbated areas (BT). NS: Non-hydrothermalized Sand. HS: Hydrothermalized Sand. FL: Fluids.

			Journal Pre-proof
	Temperature
	(°C)
		WHPs (n =	WHPs/BT (n	BT (n =	NS (n =	HS (n =	FL (n = 9)
		12)	=7)	21)	11)	19)
	Min	52.83	35.39	21.74	21.82	30.98	35.01
	Max	73.01	60.75	41.42	32.54	95.47	156.74
	Mean	66.81 ±	47.93 ± 8.26	28.75 ±	25.83 ±	60.26 ±	96.10 ±
		5.25		5.65	4.10	16.39	33.60
	Gradient
	(°C/m)
		WHPs (n =	WHPs/BT (n	BT (n =	NS (n =	HS (n =	FL (n = 9)
		5)	= 2)	9)	3)	19)
	Min	26.02	56.94	-1.30 *	7.90	1.06	13.00
	Max	94.00	59.00	50.85	24.20	133.75	415.00
	Mean	63.58 ±	57.97 ± 1.46	16.50 ±	13.90 ±	51.07 ±	156.67 ±
		28.30		18.70	8.96	43.75	127.85

Journal Pre-proof

Acknowledgments

We acknowledge the assistance of Cedric Hamelin during field operation in September 2019.

This project was partially funded by INSU-CNRS Tellus and Syster Projects to JE (2016) and JEM (2018), and by the RAMONES project (EU H2020 FET 101017808) to JE. Additional support was provided by internal funds from ENS, IPGP, U. of Lyon, and U. of Bergen. We also benefitted from the comments of two anonymous reviewers, that helped us clarify and improve this manuscript.

Data Availability Statement

The drone photomosaics used in this article can be found at https://doi.org/10.17882/77670, an open-source online data repository hosted at SEANOE [START_REF] Martelat | Milos Shallow Water Hydrothermal System: Drone seafloor photomosaics[END_REF].

The AUV photomosaics used in this article can be found at https://doi.org/10.17882/77669, an open-source online data repository hosted at SEANOE (Puzenat et al., 2019a).

The shapefiles for AUV and Drone photomosaic interpretations and instrument locations used in this article can be found at https://doi.org/10.17882/77668, an open-source online data repository hosted at SEANOE [START_REF] Puzenat | Milos Shallow Water Hydrothermal System: Shapefiles for AUV and Drone photomosaic interpretations and instrument locations[END_REF].

The temperature records used in this article can be found at https://doi.org/10.17882/77656, an open-source online data repository hosted at SEANOE [START_REF] Puzenat | Milos Shallow Water Hydrothermal System: Autonomous temperature records (Summer2019 fieldwork)[END_REF].

Table S1. Master table of temperature data acquired at the Hexagon study site. Values are sorted by profile and then by probe. Absolute coordinates are based on GPS data. Relative coordinates are based on the GIS mapping. Distance along profile is calculated using the Poseidon instrument as the origin (distance=0). Negative values are for the inside of the hexagon while positive values are for the outside. Time in hours is in UTC, with time=0 corresponding to 25/09/2019 at 00:00.

Start and end of deployments were handpicked to readjust uncertainties from GPS time that might influence the final temperature values. Unless otherwise specified, temperatures reported correspond to the mean of the last minute of deployment, when the record is most stable. Red

Journal Pre-proof

Declaration of interests

☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

☐The authors declare the following financial interests/personal rela�onships which may be considered as poten�al compe�ng interests:

HIGHLIGHTS

-Comprehensive study on the shallow-water hydrothermal system of Milos -Seafloor mapping using drone and AUV imagery