Applicability of Landsat 8 Thermal Infrared Sensor to Identify Submarine Groundwater Discharge Springs in the Mediterranean Sea Basin

Submarine groundwater discharge (SGD) has received increasing attention over the past two decades as a source of nutrients, trace elements and pollutants to the ocean that may alter coastal biogeochemical cycles. Assessing submarine groundwater flows and their impacts on coastal marine 15 environments is a difficult task since it is not easy to identify and measure these water flows discharging into the sea. The aim of this study is to prove the great usefulness of the freely-available thermal infrared (TIR) imagery of the Landsat 8 thermal infrared sensor (TIRS) as an exploratory tool to identify SGD springs worldwide, from local to regional scales, for long-term analysis. The use of satellite thermal data as a technique to identify SGD springs in seawater is based on the identification of 20 thermally-anomalous plumes obtained from the thermal contrasts between groundwater and sea surface water. We propose a conceptual framework to apply this technique worldwide and also discuss the limitations of using this technique in SGD studies. The study was developed on a regional scale in karstic coastal aquifers in the Mediterranean Sea basin at different seasons and diverse meteorological conditions. Although this study demonstrates that the freely-available satellite TIR remote sensing is a 25 useful method to identify coastal springs in karst aquifers both locally and regionally, the limiting factors include technical limitations, geological/hydrogeological characteristics, environmental and marine conditions and coastal geomorphology.

Correspondence to: Sònia Jou-Claus (soniajouclaus@gemail.com) Abstract. Submarine groundwater discharge (SGD) has received increasing attention over the past two decades as a source of nutrients, trace elements and pollutants to the ocean that may alter coastal biogeochemical cycles. Assessing submarine groundwater flows and their impacts on coastal marine 15 environments is a difficult task since it is not easy to identify and measure these water flows discharging into the sea. The aim of this study is to prove the great usefulness of the freely-available thermal infrared (TIR) imagery of the Landsat 8 thermal infrared sensor (TIRS) as an exploratory tool to identify SGD springs worldwide, from local to regional scales, for long-term analysis. The use of satellite thermal data as a technique to identify SGD springs in seawater is based on the identification of 20 thermally-anomalous plumes obtained from the thermal contrasts between groundwater and sea surface water. We propose a conceptual framework to apply this technique worldwide and also discuss the limitations of using this technique in SGD studies. The study was developed on a regional scale in karstic coastal aquifers in the Mediterranean Sea basin at different seasons and diverse meteorological conditions. Although this study demonstrates that the freely-available satellite TIR remote sensing is a 25 useful method to identify coastal springs in karst aquifers both locally and regionally, the limiting factors include technical limitations, geological/hydrogeological characteristics, environmental and marine conditions and coastal geomorphology.

1.Introduction
Submarine groundwater discharge (SGD) is an important component of the hydrological cycle and has 30 been commonly defined as any flow of water across the continental margin in the ocean-aquifer interface, regardless of fluid composition or driving force with spatial scale lengths of meters to kilometers (Burnett & Dulaiova, 2003;Moore, 2010;Taniguchi et al., 2019). This definition includes meteoric fresh groundwater resulting from inland recharge, but also seawater circulated through the sediments of coastal aquifers (Burnett & Dulaiova, 2003). Both water flows mix in coastal aquifers, 35 where biogeochemical reactions may occur when this groundwater interacts with the geological matrix (Moore, 1999). This dynamic mixing zone influences the transfer of chemical compounds such as nutrients, trace metals and other contaminants to coastal waters (Alorda-Kleinglass et al., 2019;Boehm et al., 2004;Rodellas et al., 2015;Trezzi et al., 2016). SGD-derived inputs of chemical compounds can highly impact coastal ecosystems, by influencing productivity, biomass, species composition and 40 sonification (Andrisoa et al., 2019;Garcés et al., 2011;Garcia-Orellana et al., 2016;Krest et al., 2000). Depending on its geological background SGD occurs in three different ways: coastal onshore springs (Garcia-Solsona et al., 2010;Mejías et al., 2012), diffuse submarine springs (Rodellas et al., 2012) and focused submarine springs (Bakalowicz, 2015;Fleury et al., 2007;Katz et al., 2009;Yobbi, 1992). Identifying and mapping groundwater discharge areas is challenging, despite the number of traditional 45 methods available to locate the main groundwater discharge locations and quantify their flow rates. These methods include easy procedures such as traditional local knowledge, visual observation, changes in vegetation, changes in water temperature and salinity, seepage meters or radioactive isotope tracers (Mejías et al., 2012;Schubert et al., 2014). Apart from these methods, several authors have suggested Thermal Infrared Remote Sensing (TIR-RS) as an alternative methodology to identify potential SGD 50 spring sites, since it enables the screening and study of inaccessible areas that have scarce hydrogeological information (Wilson & Rocha, 2012). Temperature has been used successfully to study SGD by using the relatively constant temperature of groundwater with that of surface seawaters, which fluctuates seasonally (Dale & Miller, 2007). In general, groundwater maintains a constant temperature between 5 m and 100 m depth, approximately 1 -2 ºC higher than the mean annual air temperature 55 (Anderson, 2005). The detection of SGD springs via TIR-RS is possible in any environment where there is thermal contrast between the discharging fluid and the receiving surface water body (Kelly et al., 2013). TIR images have the potential to identify the location of major SGD springs, as well as to study their spatial variability, by exploring the temperature difference between coastal seawater and brackish groundwater discharges. 60 There are two types of platforms to obtain thermal infrared information: airborne TIR-RS (airplane, helicopter and drone) and satellite TIR-RS (Modis, Aster, Landsat). Airborne TIR has been used for different applications; for example Shaban et al. (2005) and Akawwi et al., (2008), conducting aerial TIR surveys along the Mediterranean and Dead Sea coastlines to identify potential SGD sites. Kelly et 65 al. (2013) used TIRS images from localized point-source SGD to demonstrate that groundwater plume areas are linearly and highly correlated to in situ groundwater fluxes. Airborne TIR-RS has been also applied in combination with other methods, not only for qualitative SGD recognition, but also for quantifying groundwater flows from freshwater springs (Danielescu et al., 2009;Mejias et al., 2012). In the coastal carbonate aquifer of El Maestrazgo (Iberian Peninsula), a combination of complementary 70 techniques was used to locate submarine springs with airborne high-resolution thermal infrared, radon measurements and physical-chemical anomalies and to quantify the groundwater discharge by direct quantification with flowmeters and Ra isotopes (Mejias et al., 2012). Tamborsky et al. (2015) combined airborne TIR overflights with coastal radionuclide surveys to investigate the significance of SGD along the north shore of Long Island (NY, US) to provide quantitative evidence for TIR-RS as a tool to 75 remotely identify and measure SGD. Finally, Danilescu et al. (2009) assessed total freshwater discharge in two small nutrient-sensitive estuaries in Prince Edward Island (Canada) using a combination of TIR images, direct discharge measurements and numerical simulations.
Compared to airplane, helicopter and drone platforms, satellite TIR-RS has some characteristics that 80 limit its use for coastal water observation. The temporal resolution is fixed and varies depending on the satellite with a minimum daily revisit frequency and a maximum frequency of every 16 days for a specific area. Furthermore, the spatial resolution, which varies between 30 and 1000 meters, is lower than airborne resolution. This results in the fact, that small thermal anomalies, induced by small flows of SGD, are likely not to be detected. Additionally, satellite TIRS images are affected by atmospheric 85 conditions (i.e., clouds and shadows). However, satellite TIR-RS imagery has the great advantages of being free of charge (Landsat), easily accessible, globally available, multi-temporal, and covering a regional scale instantaneously. These advantages turn satellite TIR-based approaches into a viable and promising option to detect SGD worldwide. 90 There are several satellite missions capable of measuring sea surface temperatures (SST) with a moderate spatial resolution and acquisition: Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER), and the Landsat satellite among others, which provide appropriate spatial and temporal resolution for large-scale SGD monitoring. For those satellite suitable for SGD research, Landsat is the one with longest TIR data provision (1982 until today and already planned to beyond 95 2030 (Wulder et al., 2019)) and the one mostly applied for SGD research (Mejías et al., 2014;Wilson and Rocha, 2012, among others) The use of satellite imagery in SGD studies has evolved in parallel with the launch of new sensors that feature spatial resolution improvements over previous sensors. However, the application of satellite TIR 100 images is neither extensive nor widespread compared to the airborne TIR images. Several SGD studies used Landsat 7 to locate groundwater discharge areas (e.g. Wang et al., 2008 ), to detect known but previously unmapped SGD locations (Varma et al., 2010), to determine the spatial extent and scale of SGD-derived temperature anomalies (Wilson & Rocha, 2012), and to infer a SGD temporal variation using long-term thermal anomaly size variations . More recent studies used data 105 obtained by Landsat 8 TIRS to identify and characterize SGD sites using the sensors' technical improvements. For example, McCaul et al. (2016) proposed a multiapproach methodology to understand submarine and intertidal groundwater discharge patterns. Xing et al. (2016) evaluated the capability of satellite remote sensing methods (Landsat 7 and 8) to detect thermal anomalies related to SGD as a possible index of the presence of offshore low-salinity groundwater storage at local scale. To 110 the best of our knowledge, there is not a single study that thoroughly compares known SGD locations to satellite-data derived thermal anomaly indications over larger spatial scales, to analyze the appropriateness of satellite TIR-RS data for SGD research in a statistically robust manner.
The aim of this study is to prove the great usefulness of satellite TIR images in different sites, covering 115 a large scale and at different seasons to assess whether Landsat 8 TIR-RS can be used as an exploratory tool to identify SGD springs worldwide, from local to regional scales, for long-term analysis. The second aim is to define a conceptual application framework, identifying and discussing the principal limitations of satellite TIR-RS in the study of SGD at the local and regional level. 120 The study was carried out on the coastal karstic aquifers of the Mediterranean Sea basin, where there are many local studies that describe the discharge processes due to the great connectivity between this coastal aquifer type and the sea (Bakalowicz, 2005;Barberá and Andreo, 2015;Worthington, 1999). In this hydrogeological context, SGD takes place mainly through submarine or aerial springs (point source). Although groundwater discharge from submarine springs represents a negligible fraction of the 125 global SGD (Luijendijk et al., 2020), in some areas, such as the Mediterranean Sea, this fraction can be locally important, strongly influencing marine ecosystems and serving as a water resource for the population (Rodellas et al., 2015;Moosdorf & Oehler, 2017 year and under what hydrogeological, marine and environmental conditions these SGD springs are observable with a satellite.

Study area
The SGD contribution to the Mediterranean ranges from (3 -50)·10 11 m 3 ·yr -1 , where fresh groundwater 135 inputs represent 1 -25% of the total SGD inputs (Rodellas et al., 2015). SGD has been described and studied in several locations along the Mediterranean coast (e.g., Bakalowicz, 2015;Mejías et al., 2012;Tulipano et al., 2005;Bejannin et al., 2017). The Mediterranean basin is characterized by 46% of its coastline formed by karstic aquifers (Bakalowicz, 2015;Fleury et al., 2007;Trezzi et al., 2016). Its narrow continental shelves prevent large tidal amplification along the coast; tidal amplitude is usually 140 less than 0.2 m (Werner et al., 2013). The Mediterranean climate is seasonal, characterized by windy, mild, wet winters and by relatively calm, warm and dry summers. Strong local winds, such as the cold and dry Tramontana, Mistral and Bora, from the north, and the hot and dry Sirocco, from the south, are typical of the region. These strong regional and seasonal wind regimes provide a substantial amount of vertical mixing in the seawater column (GROUP, 1970). In general, the main rainfall season is during 145 fall and spring, with an average annual precipitation of 500 mm y -1 (Andreo & Carrasco, 1993); the sea temperature is approximately between 18 -27 and 14 -17 ºC in the summer and winter, respectively. In this study, we selected 54 springs located in the Mediterranean Sea basin (Supplementary material 1) based on published scientific literature (Bakalowicz, 2018;Basterretxea et al., 2010;Fleury et al., 2007;Mejías et al., 2012) where groundwater discharge is known to take place and the hydrogeological 150 context has been previously described. These springs are located in Spain, France, Italy, Croatia, Greece, Turkey, Syria, Lebanon and Libya ( Figure 1). The mean flow rates of the studied springs range between 0.009 and 50 m 3 ·s -1 . The depth and/or elevation of the discharge area with respect to sea level ranges from -30 m to 9 m. The distance of the discharge ranges between 0 and 3000 m from the shoreline. The type of coastal karst aquifer has been 155 defined using the same classification as in Tulipano et al. (2005) for Mediterranean coastal karst aquifers. The classification includes three types: (1) systems with poorly-developed, but highly fractured karstification; (2) systems with well-developed karstification below sea level and open to the sea and (3) systems with well-developed karstification below sea level, partially or totally closed to the sea. 160

Landsat 8 TIRS data acquisition
To determine the optimal time period for SGD detection using remote sensing, the SST of a series of images covering all seasons was compared. The TIRS instrument of Landsat 8 is a thermal imager two thermal infrared bands centred at 10.8 µm and 12.0 µm and a ground sampling distance (GSD) of 100 165 m. However, all thermal images are resampled using a cubic convolution to 30 m (Roy et al., 2014).
Only band 10 has been used to study temperature differences in water to detect SGD sites, because data collected in band 11 of the TIRS has some large calibration uncertainties (U.S. Geological Survey, 2014b). 170 A total of 25 path/row combinations were analyzed with the Landsat 8 TIR images of the Mediterranean coast between January 2017 and December 2018 to cover all 54 known SGD sites. To that end, a total of 1536 images were acquired from the U.S. Geological Survey (USGS) with cloud coverage between 0% and 90% in each image. A manual inspection of all images resulted in finer selection of 365 images with cloud-free conditions above the areas of interest. The finer selection of cloud-free images were 175 used for subsequent steps.

Deriving SST values from Landsat 8 TIRS data
Data processing included the conversion of digital numbers to SST, including an atmospheric correction of each image following the methodology presented by Chander et al. (2009). 180 Image processing began with radiometric correction, which was performed by converting the digital number (DN) to sensor spectral radiance through band-specific rescaling gain and bias factors according to Equation (1).
A web-based atmospheric correction tool developed by Barsi et al. (2003) based on MODTRAN was used to obtain values for atmospheric transmissivity and the upwelling and downwelling radiances of the atmosphere. The applied value for water emissivity in band 10 was 0.9904 (Wen-Yao et al., 1987). 205 Finally, to obtain the sea surface temperature (SST), the corrected radiances were introduced into Equation (3).
The resulting atmospherically-corrected SST data represent temperature with an error less than 1.3 K for the temperature range of 270 -330 K. This temperature represents the skin temperature of the water 215 (<1 mm of the upper most water layer), which differs from the bulk temperature below it by about 0.1 K due to sensible heat fluxes, evaporative heat loss, and long wave radiation (Donlon et al., 2002;Wloczyk et al., 2006).

Site inter-comparison between single image and multiple images 220
Temperature maps of coastal waters were created from temperature data to assess the significance of the SST anomalies. The identification of SGD spring sites was based on the assumption that temperatures of discharging groundwater may be different than seawater and less variable than seawater temperatures throughout the year. SGD spring sites were analyzed using two different procedures. First, single images were used to identify SGD springs by means of the water temperature anomalies. As a second 225 step of the single image approach, the change between images along the study period was also evaluated. This analysis allowed identifying variations in the morphology of the discharge plume. As a second approach, called multiple imaging, SGD spring sites were detected by evaluating the pixel-bypixel standard deviation (SD) between all image sets. Lower values of SD were used as indicators for groundwater discharge using sea surface temperature (SST) data series. This statistical parameter has 230 been previously applied in semi-arid areas to study groundwater and surface water interactions and identify spring discharge into lakes or enclosed seas Tcherepanov et al., 2005). It was assumed that groundwater tends to be less variable than surface water, which varies seasonally and daily. The applied multi-temporal thermal remote sensing approach was based on a variable number of Landsat 8 TIRS images. The images used to calculate the standard deviation varied between 5 and 17 235 images, depending on the number of images without clouds available for each studied site. The thermal mapping results were combined with orthophoto maps (only the land part) within a GIS to locate the SGD springs at the identified sites.

Overall identifications
As an example of the Landsat 8 TIR images analyzed in the Mediterranean Sea basin, the thermal images of the Almyros SGD spring in Agios Nikolaos (Greece) throughout the year 2017 are shown in Figure 2. In general, a thermal anomaly plume is observed in several of the 23 images (one of the July images is missing) occurring near and perpendicular to the coastline, reflecting the continuous discharge 245 of groundwater. Since a satellite image provides the sea surface temperature (SST) of the first millimeter of seawater (Donlon et al., 2002;Wloczyk et al., 2006), the thermal contrast is because the groundwater is located above the seawater because of its lower density due to its different temperature and salinity (Wilson & Rocha, 2012). The images series of Almyros Agios Nikolaos ( Figure 2) shows how thermal contrast caused by the 250 groundwater discharge cannot be observed throughout the whole year. The thermal contrast is more identifiable from the second half of April until the end of October, but the best thermal plume observations were from June to October ( Figure 2). Conversely, this spring cannot be identified from November to March. The groundwater discharge was not identified in February and December due to the absence of thermal contrast, and in the other months (January, March and November) because 255 clouds made identification difficult. Thus, the overall percentage of time for optimal SGD spring identification was 37.5% in 2017.
If the same approach is applied for the total of 54 SGD springs studied, only 23 springs were identified in individual images, representing a 44.4% success rate of the technique over the entire study period 260 (2017 and 2018). The success percentage for identifying the 23 SGD springs according to the period of the year is shown in supplementary material 2. The highest success percentage for SGD identification was during the summer, specifically from June to September (Figure 3), which corresponds to an average of 21% of the images analyzed. Conversely, the SGD springs were not identified in the remaining 79% of the images. The winter period, from December to February (Figure 3), had the lowest 265 percentage of SGD spring identifications.

Conceptual framework to asses factors influencing the identification of SGD springs
The analysis of the images of the Mediterranean coast obtained with Landsat 8 TIRS during 2017 and 2018 was successful for several studied springs but did not identify SGD springs in all the images 270 analyzed. Thus, identifying SGD springs using Landsat 8 TIRS has some limitations, as the success rate was slightly less than 50%. The following potential limitations have been previously reported in the literature in some local studied areas: only information for the first millimeters of seawater is available (Donlon et al., 2002;Wloczyk et al., 2006), spatial resolution (Wilson & Rocha, 2012), period of the year (Bayari & Kurttas 2002;Wilson & Rocha 2012;Xing et al., 2016), the results are highly dependent 275 on atmospheric temperature, seawater currents, wind speed and direction, sea surface effects (Kelly et al., 2013) and cloud cover, and the need for specialist knowledge to convert the data into accessible (visualized) information (McCaul, 2016). Given the limitations, we propose a conceptual framework to assess the following types of factors that should be considered: (1) technical limitations, (2) geological/hydrogeological characteristics, (3)  280 environmental and marine conditions, (4) coastal geomorphology and (5) anthropogenic sources ( Figure  4).

Technical limitation factors
Some of the main limitations of the technique are related to the temporality of obtaining images, the spatial resolution, the availability of images in the desired period and the atmospheric conditions during 285 the image capture. Each satellite has an image capture spatial and temporal resolution, and therefore the results are subject to these pre-established conditions. For Landsat 8 TIRS, the temporal resolution is two weeks, and the spatial resolution is 100 m resampled to 30 m, so a smaller thermal anomaly plume than this produced by an SGD spring will not be identified. In the Mediterranean basin, the Alcalfar spring in Menorca (Spain) (Garcia-Solsona et al., 2010), which discharges in a semi-enclosed area with 290 an approximate width of 20 meters, was not identified. The availability of satellite images depends, in some cases, on technical problems that the satellite experiences during image collection, which implies that there are annual series of images with some missing images. In the example used in Figure 2 to identify the Almyros of Agios Nikolaos (Greece) SGD spring, only 23 of the 24 images that should have been produced during 2017 were obtained. The striping noise that it can be seen in the images 295 ( Figure 2) which affect TIR bands especially is another technical limitation. Sometimes the difference is so large that it would hind a proper SGD detection (Gerace and Montanaro, 2017).
Other important limitations related to the technique include the atmospheric conditions when the images are captured. Clouds and clouds shadows change radiometric information leaving the sea surface and 300 prevent a correct analysis of the images. For the Mediterranean Sea basin, of the 80% of the images in which SGD springs were not identified, 60% were due to the high presence of clouds. Thus, clouds are the main factor limiting the identification of SGD springs. The presence of clouds is higher in winter than in summer; therefore, warmer months are much better for identifying SGD springs in the Mediterranean Sea. However, in at least 20% of the 305 cloudless images, it was not possible to observe and locate SGD springs described in the literature. Therefore, a detailed analysis of the cloudless imagery is necessary to confirm the optimal conditions for locating SGD springs in the area of interest (Anderson, 2005).

Geology and hydrogeology limitation factors 310
The characteristics of coastal aquifers and the different driving forces that allow groundwater discharge to the sea are other factors that can limit the identification of SGD springs. According to Taniguchi et al. (2002), SGD is the combination of fresh groundwater discharge from the land and saline groundwater discharge from the seawater recirculation through the coastal and continental margin. Thus, the thermal signal of both types of groundwater will define the intensity and shape of the plume 315 that can be observed in the sea.
For coastal aquifer SGD springs, geological factors such as the lithology (e.g., type of rock, degree of karstification, presence of faults or fractures, etc.) determine the hydraulic properties that greatly influence the groundwater flow that discharges into the sea (Sander at al., 1996;Brunner et al., 2007;320 Edet et al., 1998). Another aspect to consider is the connectivity of the aquifer with the sea. Coastal aquifers may be totally or partially hydraulically connected to the sea, depending on several factors such as the geological structure of the aquifer and the seabed, as well as the hydraulic properties of both onshore and offshore geological formations. In addition, the amount of SGD into the sea depends on the aquifer water budget, which produces variations in groundwater discharge flow. The main factors that 325 influence this budget are the natural variations in recharge at different temporal scales (rainfall events, seasonality, inter-annual changes) and the abstraction of groundwater, which can reduce the discharge of the fresh groundwater component of the SGD, or even induce seawater intrusion into the coastal aquifer. Unfortunately, these factors are linked since in most cases when recharge is low, abstractions tend to be higher. 330 The present study was performed in a specific geological context -coastal limestone with different levels of karstification -with SGD occurring through springs (point source discharge). To analyze the influence of geological and hydrological component on the 54 studied karstic SGD springs, hydrogeological characteristics were analyzed only when detailed information was available, that is, for 33 springs representing 62% of the total (Supplementary material 1). According to the hydrogeological 335 characteristics, 15% of the springs correspond to a karstic system characterized by well-developed karstification connected to the sea (e.g., Moraig, Port Miou, Bestouan, Almyros of Iraklio, Almiros of Agios Nikolaos, Cephalonia, Ain Zayana and Chekka). Of these well-developed springs, most were identified; only the Ain Zayana and Chekka springs were not observed. 9% of the springs 9% of the springs for which karst type information was available, were systems with poorly-developed but highly 340 fractured karstification. These karst systems included 3 different fractured systems in which 1) faults dissect the aquifer such as in the Gokova (4 springs) (Bayari & Kurttas, 2002) and Ovacik spring, where the faults are located in the underlying beds that extend towards the sea (Elhatip, 2003); 2) groundwater flows along the zones of cracks, fractures and karst hollows, such as in the Donnalucata spring (Povinec et al., 2006); and 3) groundwater flows through stratification joints, such as in the Mortola spring 345 (Fleury et al., 2007). Of all these types of SGD spring located in fractured systems, thermal anomalies were only identified in the first type, in which the faults dissect the aquifer. The last type of defined karst system is one with well-developed karstification but low connectivity with the sea. This group is represented by only two springs: Kiveri Anavalos in Greece and Vise in France, of which only the first was identified. Therefore, the best geological and hydrogeological settings for identifying SGD springs 350 are karst systems with well-developed karstification and that are well-connected to the sea, as well as those systems with poorly-developed but highly fractured karstification. This is because these systems can drain large areas of recharge; therefore, SGD springs often have high flow rates, making it is easy to detect their thermal plumes. 355 SGD spring identification can also be limited by the groundwater flow rate, because the thermal anomaly generating the SGD thermal plume may by generated by the head transported by fresh groundwater flow (convection) and/or by the offshore movement of seawater intrusion (due to conduction and/or convection), which is also affected by groundwater flow rate. Therefore, low discharge water volumes it is unlikely be identified using this technique. In the studied case of the 360 Mediterranean Sea basin, the identified springs have an average flow rate between 0.75 m 3 ·s -1 and 50 m 3 ·s -1 . The identified springs with high flow rates were Almyros of Iraklio in Greece (with values between 2 and 30 m 3 ·s -1 ) and Jurjevska in Croatia (with values between 0.8 -50 m 3 ·s -1 ). The identified SGD springs with low flow rates were Ovacik in Turkey (with values of 0.75 m 3 ·s -1 ), Gokova in Turkey (with values between 0.02 -0.10 m 3 ·s -1 ), las Fuentes in Spain (between 0.03-0.17 m 3 ·s -1 ) and Moraig in 365 Spain (between 0.3 -9.0 m 3 ·s -1 ). However, there were springs with high flow rates, such as the Ses Ufanes spring in Spain (with values between 20 -30 m 3 ·s -1 ), and Font Dame-Font Estramar (with values between 1 -15 m 3 ·s -1 ), that were not identified, indicating that other factors affect the failure to identify these springs. 370 Interannual and annual changes, such as humid-dry periods/seasons, prolonged alteration to the groundwater systems (overexploitation) and extreme events such as strong rainfall storms, etc., can also modify SGD flows. Thus, the research period in which the images were collected is also important because seasonal variations (Lee et al., 2016) can lead to SGD thermal plume variations, not only due to the temperature difference between groundwater and seawater (Lee et al., 2016), but also annual SGD 375 flow rate variations (Michael et al., 2005) that allow their identification.
Regarding hydrogeological conditions in the Mediterranean Sea basin, it was expected that the best seasons for identifying SGD springs using satellite TIR-RS would be spring and autumn, when rainfall is higher and therefore higher discharge flow rate is expected. However, summer was the best season 380 for identifying SGD springs, indicating that other factors influence the identification of SGD springs using this technique.

Environmental and marine conditions
Other aspects that could affect the thermal contrast between the groundwater plume and the seawater 385 include environmental and marine conditions. Both factors can make the identification of SGD difficult, because they cause seawater mixing with the discharging groundwater, reducing the thermal contrast between them and modifying the sea surface temperature. The main environmental condition that can affect SGD identification is the action of wind, which can mix the first millimeters of the sea surface water, limiting the identification of SGD springs using remote sensing techniques. Similarly, marine 390 conditions can affect SGD identification if they reduce the thermal contrast between the groundwater plume and the seawater. These marine conditions are basically reduced to the coastal water hydrodynamic conditions that can be affected by processes such as the influence of tides or coastal currents, the formation of a pycnocline in surface seawater or the fetch due to wind action. Tides, coastal currents and fetch generate seawater movement and mix groundwater with seawater, causing a 395 thermal contrast attenuation. Similarly, the presence of a pycnocline can result in less vertical mixing of the water column. In subtropical areas with cold winters and hot summers, such as the Mediterranean Sea, coastal waters often develop a pycnocline during the summer months, as high temperatures increase the evaporation of seawater, generating an increase in salinity and therefore water density. This effect causes cold, fresh groundwater to flow over salty and dense seawater, generating an SGD layer 400 on the sea surface.
In the Mediterranean Sea, which has low tides and maximum thermal contrast between groundwater and seawater during summer and winter, the SGD spring visualizations in cloud-free images decrease significantly in winter months compared to warmer months ( Figure 4). Only 2 of 21 SGD springs were 405 identified during winter in 2017 (Anavalos Kiveri, Greece and Torre Badum, Spain) and 0 of 20 in 2018. The influence of marine factors in both cases was minimal; identification is therefore plausible during winter months. Since the number of SGD spring visualizations is much higher in summer than in winter, the environmental and marine conditions during winter months are unfavorable for the location of coastal springs. 410 Consequently, areas where the influence of various marine factors such as tides, coastal currents or fetch is small are ideal for identifying SGD springs, because such factors attenuate the thermal contrast needed to identify SGD springs. Similarly, in subtropical areas such as the Mediterranean Sea, it is easier to identify SGD springs because coastal waters often develop a pycnocline during the summer 415 months, resulting in less vertical mixing of the water column and therefore better identification of SGD springs. Therefore, in the Mediterranean basin environmental and marine conditions in summer months are much more favorable for the identification of coastal springs than in winter months.

Coastal morphology limitation factors 420
Another aspect that influences SGD thermal plume visualization is the coastal morphology; depending on its characteristics, the seawater residence time in the study zone may allow the formation of a thermal plume. In semi-enclosed areas, such as bays or coves, where the residence time of discharging water in the sea is in the range of one or several days (Tamborski et al., 2020), the formation of thermal plumes from SGD can be easily detected by satellite in comparison with open sea areas, where the 425 seawater residence time is shorter (less than a day), due to the effect of marine or weather factors such as waves, coastal currents or fetch. Groundwater discharge can occur below sea level, either at shallow (< 10 m) or at greater depths, right on the shoreline or a few meters from it (Zektser & Dzhamalov, 2006). The morphology of the coastal seabed combined with the geological characteristics of the aquifer (e.g., karstification degree and coastal hydraulic gradient) determine both the depth and offshore 430 distance of the groundwater discharge. These two characteristics are very important, since the SGD spring's location on the coast is critical to the correct formation of the plume and the thermal contrast between groundwater and seawater that must be recorded by satellite. While groundwater discharges produced at coastal level or several meters inland usually generate thermal contrasts easily detectable by satellite (Mejías et al., 2012), submarine springs located several meters deep often represent a challenge 435 for satellite detection. This is because when groundwater discharges below sea level, it rises to the sea surface generating a buoyant thermal layer of several millimeters, due to its lower density that arises from its temperature and salinity characteristics. Therefore, a greater discharge depth requires more time to reach the sea surface, adapting the SGD's thermal signal at the sea surface to the surrounding water temperature. Thus, marine factors may have a greater influence in balancing the seawater and 440 groundwater temperatures, preventing the recognition of thermal anomalies on the sea surface, as stated by Mallast et al. (2013).
However, the time required for groundwater to reach the sea surface depends not only on the depth of the discharge below sea level, but also on hydrogeological factors, such as the discharge flow. Thus, 445 SGD springs characterized by shallow depths and large flow rates will favor the detection of SGD induced thermal plumes on the sea surface, while deep areas with small flows will be undetectable.
Intermediate situations, such as shallow areas with small flows or deep areas with large flows, may be detected depending of the relative importance of the environmental and/or marine limiting factors. 450 For the karstic springs studied in the Mediterranean Sea basin, 28 of the 54 springs studied were located in semi-enclosed coastal areas. Of these, 16 were identified by satellite. These springs were Port Miou and Bestouan in France; Cephalonia, Anavalos Kiveri and Almyros of Agios Nikolaos in Greece; Martinscica, Perilo, Novijanska, Jurjevska and Patan in Croatia; and Gokova (four springs), Ovacik and Antalya (four springs) in Turkey. These results show that when discharge occurs in semi-enclosed 455 areas, where the seawater residence time is large, the springs can be identified until late autumn, even though the flow of groundwater is relatively small (0.75 m 3 ·s -1 was reported for the Ovacik spring in Turkey). Therefore, the success in identifying SGD springs in semi-enclosed areas also is higher throughout the year. 460 Our study showed that, of the 54 studied springs for which information is available, the seashore distances ranged from the shoreline to 1 km, and the discharge depths vary between 7 to -150 m with respect to sea level. The first group of springs discharge near the seashore and reach the sea through small streams; these karstic springs are located between 300 and 500 m inland and at elevations of 2 m, 3 m and 15 m above sea level for Patan in Croatia, Almyros of Iraklio in Greece and Maro in Spain, 465 respectively. Of them, only Patan and Almyros of Iraklio were identified. The second group of springs discharges in coastal lagoons at a distance of 100 m from shore and a depth of -4 m (Font Dame and Font Estramar in Salses-Laucate lagoon in France) and an unknown shore distance and -30 m (Vise in Thau lagoon). None of these springs were identified. The third group of springs is located between 0 and 10 m from the shoreline and in shallow waters between 0 and -7 m (Torre Badum, Las Fuentes, 470 Font de Dins in Spain, Ain Zayana in Libya, Agios Nikolaos, Cephalonia and Anavalos Kiveri in Greece and Ovacik and Gokova in Turkey). All of these springs were identified except the Ain Zayana spring in Libya. The fourth group of springs is located close to the shoreline, but at a water column depth of -12 m. The two springs of this group were identified (Moraig in Spain and Port Miou and Bestouan in France). 475 Neither of the two springs representing the fifth and last group, in which discharges occurs offshore between 100 m and 1 km with a water column depth between -35 and -150 m (Mortola in Italy and Chekka in Lebanon) were identified.
From the analysis of the images obtained by Landsat 8 TIRS and considering the hydrogeological data 480 of the various coastal springs, we conclude that groundwater discharges that occur at significant depths (> 12 m) it is unlikely be identified by this technique, probably because the thermal anomaly generated between groundwater and seawater does not reach the sea surface. This is the probable reason for the failure to identify the springs of Vise in France, Mortola in Italy and Chekka in Lebanon. For the Vise spring, which discharges at -30 m depth, the hydrogeological information indicates that it is a mixture 485 of karstic water, thermal water and seawater (Aquilina et al., 2002); most likely the discharging groundwater does not produce enough thermal contrast to be detected with a satellite sensor. On the other hand, SGD springs located in very shallow areas, such as the Salses-Laucate lagoon, which has an average maximum depth of -2 m (Bejannin et al., 2017), and where there are two well-known springs (Font Dame and Font Estramar) were not identified. A possible reason is that these types of shallow 490 coastal lagoons are highly influenced by atmospheric temperature, since the small water column has little capacity to buffer temperature fluctuations (Lee et al., 2016). However, Anavalos Kiveri and Torre Badum, identified during winter months, are shallow (0 and to -10 m) SGD springs. Therefore, although initially these types of springs should be easily detectable by satellite, atmospheric factors can greatly influence the formation of the thermal contrast needed to identify them. 495 The shore distance at which groundwater discharge occurs is another factor that affects the identification of groundwater discharge in coastal areas. SGD springs discharging further than 500 m from the seashore were not recognized, such as the Mortola spring with a distance of 800 m. This may be related to the fact that as the distance to the coast increases, the depth of the water column increases, 500 and the travel time of the fresher groundwater also increases, allowing a complete temperature equilibrium. Furthermore, these zones are more affected by ocean currents that increase water mixing, attenuating the thermal signal.
In some cases, there is not a single factor affecting SGD spring identification, but rather a combination 505 of several factors (flow rate, season, coastal morphology, etc.) which makes difficult to explain why an SGD spring was not identified.

Anthropogenic sources
There are some coastal outflows originated from anthropogenic sources may mask or modify the SGD 510 thermal signal or even generate a similar thermal signal such as wastewater treatment outflows, fish farms, ports, power plants, among others. For example, wastewater treatment facilities discharge directly to the coastal waters in relatively shallow waters. Even when discharged at depth during periods with high stratification, it is possible to detect its thermal effect on the coastal waters due to buoyant plumes to surface (DiGiacomo et al., 2004) with SST differences of at least 0.5 ºC identifiable with 515 TIR-RS (Gierach et al., 2017). In the case of inland aquaculture, pumped seawater temperature may be modified to reach optimal growth conditions for some species. For instance, water at 15ºC is used in European sea bass hatcheries (Navarro-Martín et al., 2009).

Application of the multi-temporal SST series method 520
In some studied areas, there may be low thermal contrast that can prevent the identification of SGD springs. In these cases, it is possible to deepen the analysis of the images by using the multi-temporal SST series proposed by Mallast et al. (2014) in the Dead Sea that minimizes the effect of the previously-indicated limiting factors. Unlike single images, the multi-temporal SST series method allows several images and thus several points in time to be integrated, accentuating the small thermal 525 anomalies for easy identification. An example of the usefulness of the multi-temporal SST series method is in the Serra d'Irta (Eastern Iberian Peninsula) (Mejias et al., 2012;Garcia-Solsona et al., 2010;Trezzi et al., 2016). Although the multi-temporal SST series method is better for identifying SGD springs because thermal contrast is enhanced, this method did not enable SGD spring identification in other places where the single image method did not identify them. Furthermore, with multi-temporal SST series, the temporal morphological information of the discharge plume is lost (Mallast et al., 2013). Conversely, single images allow the identification of morphological variations of the discharge plume and the temperature 540 variations along the plume. Therefore, if the study objective is to identify SGD springs, the multitemporal SST series method is the best option, but if the objective is to study the SGD plume variation, single images are better.

Identification of new SGD springs in the Mediterranean basin 545
One of the main objectives of this study is to demonstrate the great usefulness of satellite TIR imagery at local and regional scales, identifying SGD springs not previously described in the scientific literature. An example of this usefulness is the identification of 21 SGD springs undescribed in the bibliography along the Croatian coast using single images ( Figure 6). This demonstrates that the analysis of single images obtained from Landsat 8 TIRS allows the identification of new SGD springs. Therefore, this 550 economical technique is very useful in inaccessible areas (Gumma & Pavelic, 2013;Edet et al., 1998;Sander et al., 1996).

Challenge and recommendations for future SGD studies and application in other areas
This study has shown the usefulness of satellite images for identifying coastal springs in limestone 555 karstified aquifers, both locally and regionally. Therefore, in the design of SGD studies at both local and regional levels, it is recommended to first screen the coastal water temperatures for thermal anomalies using the presented satellite based TIR approach that will help narrow the sampling surveys. Although the study of the thermal images during a single year should be sufficient to identify potential coastal springs, it is highly advisable to analyze thermal images over several years, since there may be one or 560 several factors (technical, marine, environmental, hydrogeological, etc.) that may alter the SGD thermal signal.
Although this study focused on the discharge of groundwater from springs located in karst aquifers, the study of satellite TIR images could be extrapolated to other types of aquifers where the discharge is 565 more diffused and where the proportion of meteoric water is lower. To identify thermal plumes, water discharging into the sea, whether of meteoric origin or marine circulation, must have acquired sufficient thermal contrast in the coastal aquifer before discharge. This thermal contrast can also occur in places where coastal aquifers are volcanic rocks or fractured granites, and even in sedimentary formations. However, in each case, the specific geological context (aquifer matrix, hydraulic parameter distribution, 570 etc.) should be considered in the analysis.
Several studies have shown the possibility of quantifying SGD through the study of thermal plumes obtained by aircrafts (Tamborski et al., 2015;Danielescu et al., 2009). By determining the thermal plume area, which is often directly related to the discharge volume, SGD may be quantitatively 575 identified. However, there are several factors that can alter the thermal plume shape that could result in error. Furthermore, because the TIR image only allows a observation of the less than 1 mm of thermal contrast on the sea surface prevents, a determination of the volume beneath the thermal plume is missing. Moreover, quantification by means of images represents a major challenge because satellite images (e.g., Landsat 8 TIRS) can only be obtained once every 16 days, and many factors can alter the 580 shape of the plume independent of the actual discharge. Therefore, unless the discharge occurs in locations where the environmental and marine conditions are well-controlled, quantification by satellite images represents a major challenge that should be further investigated.
Although this study focused regionally on the Mediterranean Sea basin, it can be extrapolated to other 585 parts of the world, in places where the SGD has sufficient thermal contrast when discharging into the sea. This technique has been successfully used in regional areas such as Ireland (Wilson & Rocha, 2012;McCaul et al., 2016) or the Laizhou Bay in China (Xing et al., 2016). However, in areas where the thermal contrast between the sea and the aquifer is low and/or temperature does not fluctuate during the year, such as in tropical zones, the use of satellite images represents a challenge that must be 590 explored in detail. Furthermore, it is important to perform an initial study to attempt to identify the previously-described main factors that can limit SGD spring identification in the study area, to try to control/avoid these effects.

Conclusions
This study demonstrates that satellite remote sensing is a useful tool for the identification of coastal 595 SGD springs in karst aquifers, both locally and regionally, by performing an initial screening of the coastal water temperature images to identify possible thermal anomalies that can help narrow SGD study sampling surveys. However, this study highlights limiting factors that should be considered: (1) technical limitations, (2) geological/hydrogeological characteristics, (3) environmental and marine conditions and (4) coastal geomorphology. 600 In the karstic coastal aquifer of the Mediterranean Sea basin, the highest success percentage of SGD visualizations was during June, July, August and September. The best geological and hydrogeological settings for identifying SGD springs are karst systems with well-developed karstification that are wellconnected to the sea, as well as those systems with poorly-developed, but highly fractured karstification. 605 Environmental and marine conditions in the summer months are much more favorable for the identification of coastal springs than in the winter months. Semi-enclosed areas, where the seawater residence time is large, favor the identification of SGD springs throughout the year. Groundwater discharges that occurred deeper than 12 m or in very shallow waters were not identified. Also, SGD springs that discharged further than 500 m from the seashore were not recognized. 610 Multi-temporal SST series are better for identifying SGD springs in coastal zones, but the information of the discharge plume changes is lost, while single images are more suited to study the morphological variations in the discharge plume and temperature variations along the plume. 615 Although this study focused regionally on the Mediterranean Sea basin, it can be extrapolated to other parts of the world and other aquifer types, in places where the SGD has sufficient thermal contrast when discharging into the sea and where the specific geological contexts (e.g., aquifer matrix and hydraulic parameters distribution) are considered. Long time series are better for identifying SGD spring areas, since there may be one or several factors (technical, marine, environmental, hydrogeological, etc.) that 620 can alter the SGD's thermal signal. It is recommended to first screen the coastal water temperature images obtained by satellite to identify possible thermal anomalies that will help narrow the sampling surveys. This technique allows the identification and quantification of SGD springs in zones without hydrogeological information. 625

Data availability
All data are available from the corresponding author upon request.