Evaluating the Applicability of Thermal Infrared Remote Sensing in Estimating Water Potential of the Karst Aquifer: A Case Study in North Adriatic, Croatia

: One of the most prominent tourist destinations in the Adriatic coast, the city of Opatija, is facing a problem concerning seasonal drinking water shortages. The existing water resources are no longer sufﬁcient, and attention is being given to alternative resources such as the underlying karstic aquifer and several coastal springs in the city itself. However, the water potential of the area still cannot be estimated due to the insufﬁcient hydrological data. The goal of this research was to evaluate the use of thermal infrared (TIR) remote sensing as the source of valuable information that will improve our understanding of the groundwater discharge dynamics. Ten Landsat ETM+ (enhanced thematic mapper plus) and two Landsat TM (thematic mapper) images of the north Adriatic, recorded during 1999–2004 at the same time as the ﬁeld discharge measurements, were used to derive sea surface temperature (SST) and to analyze freshwater outﬂows seen as the thermal anomaly in the TIR images. The approach is based on ﬁnding the functional relationship between the size of the freshwater thermal signatures and the measured discharge data, and to estimate the water potential of the underlying aquifer. It also involved analyzing the possible connection between the adjusted size of the spring’s thermal signatures and groundwater level ﬂuctuations in the deeper karst hinterland. The proposed methodology resulted in realistic discharge estimates, as well as a good ﬁt between thermal anomalies with measured discharges and the groundwater level. It should be emphasized that the results are site speciﬁc and based on a limited data set. However, they conﬁrm that the proposed method can provide additional information on groundwater outﬂow dynamics and coastal springs’ freshwater quantiﬁcation.


Introduction
The increase in frequency and intensity of extreme hydrological events (e.g., droughts as well as floods) over the past few decades has emphasized the lack of freshwater availability throughout the world [1][2][3][4]. Pressure on water resources in terms of water quality and water demand is tremendous. The limitations of existing water resources encouraged researchers to turn their attention to new, unconventional water resources such as coastal and submarine springs. Investigations conducted over the last decades emphasized the importance of coastal and submarine springs worldwide (e.g., [5][6][7][8][9]), not only in terms of their availability to meet the freshwater demand but also for their importance for the coastal ecology. Coastal and submarine springs are common along the Mediterranean coastline [10], including the Adriatic Sea (e.g., [11]).
Discharge processes of the coastal and submarine springs are complex. Assessing the discharged quantities by means of direct measurements is often not possible, so finding an adequate approach to estimate water availability in dynamic contact conditions between fresh water and sea water is a challenging task. Various methods have been proposed to quantify coastal and submarine springs including techniques based on the natural or artificial tracers, water balance and hydrogeological modelling, direct measurements via the seepage meters, and eddy correlation approaches [7,12]. Development of thermal infrared (TIR) remote sensing (aerial high resolution [13][14][15][16][17] and satellite mid-resolution based [18,19]) enabled the detection of the springs, rather than estimating outflow quantities.
The research presented in this paper is a step forward in attempting to estimate outflows of the non-monitored springs using the size of the thermal anomalies created on the sea surface caused by the springs' water and the measured discharge data. Discharge and groundwater level data that was used was collected twenty years ago since the measurements were drastically reduced over the last two decades. Hydrological indicators for the springs' outflow regimes were monitored only during short term water exploration works twenty years ago, albeit without adequate continuity and processing results. Limited information on outflow characteristics imposed rough approximations, which increased the uncertainty of these estimations. However, the intention of the paper was not to reliably quantify the springs' outflows, but to explore the possibility of using remotely sensed TIR data in estimating a region's water potential.
With the goal of providing the hydrological data needed to estimate a region's water potential, continuous monitoring is planned to be restored by mid-2021. Apart from systematic hydrometric measurements conducted on the coastal springs, monitoring should also utilize the information derived from the TIR images.

Research Area and Environmental Settings
The city of Opatija, situated in the northern part of the Adriatic Sea, experiences a significant seasonal rise in water demand. As the number of inhabitants increases considerably during summer, water supply authorities are facing the difficult task of overcoming unfavorable hydrological conditions and insuring a continuous water supply.
Numerous coastal and submarine springs (vruljas) surround the city of Opatija ( Figure 1). The influence of the sea on the springs' water salinity is significant and the springs have not been seriously considered in the past as an alternative water source until recently. The considerable size of the contributing area (several hundreds of square kilometers), large potential storage capacity of the karst aquifer, and high amounts of rainfall (annual averages vary from 1500 to 2000 mm) triggered the attention of the researcher and water authorities to consider the springs as a possible alternative for meeting the increasing water demand.
Continuous hydrological monitoring on springs was never fully established. Sporadic measurements and analyses of the outflow regime in the springs' zone were conducted in the past as part of the occasional water exploration works. They focused primarily on mitigating adverse consequences of flooding rather than exploring the possibility of using the springs' water for water supply purposes and resulted only in rough estimations of the capacity of larger springs. Establishing sustainable management of freshwater resources and finding the right measure of using the water from coastal springs requires primarily a good understanding of the water quantities available for our use. In that respect, with the goal of providing the hydrological data needed to estimate region's water potential, continuous monitoring is planned to be restored by mid-2021. It will also utilize historic hydrological data and emphasize finding new approaches (e.g., methods based on integrating remote sensing and field measurements) to fill the data gaps.
Kristal and nearby coastal springs drain the catchment area of 365.80 km 2 ( Figure 1). More than 80% of the catchment (~300 km 2 ) is located on the Croatian territory while the rest is in Slovenia. With altitudes up to 1396 m a.s.l., it stretches in the northwest-southeast direction from the Croatian-Slovenian border across theĆićarija and Učka mountain ranges to the coast and has a dominant influence on the climate of the region. The average annual rainfall varies from 1500 mm on the coast to 2500 mm in higher altitudes [20], while the average annual air temperatures reach 15 • C on the coast and decreases as the altitude is increasing down to of 6 • C [21]. The minimum rainfall is recorded in summer, especially July, and the maximum rainfall is usually reached in autumn (i.e., November). Kristal and nearby coastal springs drain the catchment area of 365.80 km 2 ( Figure 1). More than 80% of the catchment (~300 km 2 ) is located on the Croatian territory while the rest is in Slovenia. With altitudes up to 1396 m a.s.l., it stretches in the northwest-southeast direction from the Croatian-Slovenian border across the Ćićarija and Učka mountain ranges to the coast and has a dominant influence on the climate of the region. The average annual rainfall varies from 1500 mm on the coast to 2500 mm in higher altitudes [20], while the average annual air temperatures reach 15 °C on the coast and decreases as the altitude is increasing down to of 6 °C [21]. The minimum rainfall is recorded in summer, especially July, and the maximum rainfall is usually reached in autumn (i.e., November).

of 15
Hydrographic network is limited to several short and steep torrential watercourses. Most of the precipitation infiltrates directly into the karstic aquifer and flows into the sea through several coastal and submarine springs. Due to immediate contact with sea, the salinity of the springs is rather high: the chloride concentration varies from 570 to 11,000 mg L -1 . Previously mentioned sporadic hydrological measurements in the Opatija's springs zone were conducted over the last two decades, from 1983 to 2004, mostly during dry periods. The springs' low flows varied between 0.4 and 1.0 m 3 s −1 [22].
Apart from the Opatija region, groundwater outflows are also evident in the City of Rijeka and in the nearby Rječina River catchment area. Zvir Spring located in the down- Hydrographic network is limited to several short and steep torrential watercourses. Most of the precipitation infiltrates directly into the karstic aquifer and flows into the sea through several coastal and submarine springs. Due to immediate contact with sea, the salinity of the springs is rather high: the chloride concentration varies from 570 to 11,000 mg L -1 . Previously mentioned sporadic hydrological measurements in the Opatija's springs zone were conducted over the last two decades, from 1983 to 2004, mostly during dry periods. The springs' low flows varied between 0.4 and 1.0 m 3 s −1 [22].
Apart from the Opatija region, groundwater outflows are also evident in the City of Rijeka and in the nearby Rječina River catchment area. Zvir Spring located in the downstream section of the Rječina River is the most significant one not influenced by the sea salinity and is therefore crucial to the city of Rijeka's water supply.

Materials and Methods
Having temperature lower than the surrounding sea, the springs' outflows create a thermal anomaly on the sea surface that can be spotted on the TIR imagery, depending on the size of the springs' outflow and the spatial resolution of the TIR sensor. The goal of the research was to analyze the possible relationship between the extent of the thermal anomaly created by the springs' outflows and the discharge quantities, based on the existing hydrological data, and to use that relationship in estimating outflows at the nonmonitored locations. Apart from the occasional measurements on larger springs (Kristal Spring has constant outflow throughout the year while the springs Admiral and Slatina dry out during dry season), hydrological monitoring is established only on watercourses such as the Rječina River ( Figure 1). Hydrological conditions in the Rječina River catchment area are similar to those of the catchment drained through the springs in the Opatija's discharge zone. The dominant input into the Rječina's water balance are numerous karstic springs within the catchment. The minimum water temperature recorded in the river was 6 • C in the Rječina's spring located at the altitude of 325.24 m a.s.l. and the maximum did not exceed 10 • C [23,24]. This indicates that water remains in the karst aquifer for a long time and is not influenced by highly variable air temperature. Groundwater temperatures are almost the same at both analyzed locations; the Rječina River is fed by the numerous springs within the catchment and, even though the river upstream of spring Zvir often dries out during dry season, overflow from the spring continue to feed the river's discharge. The length of the river is 18.63 km [25] and the most downstream gauging station, Sušak tvornica, is located~450 m upstream from the river mouth. Most recent research shows that the Rječina River's catchment area is 381 km 2 [26]. The Kristal Spring and the Rječina River are geographically relatively close hence similar mixing conditions can be assumed, based on the existing knowledge of the surface conditions (waves and wind).
The approach to estimating the non-monitored coastal spring discharges assumes that the size of the thermal anomaly created by the spring's outflow on the sea surface is closely related to the outflow quantities. Detection and delineation of the anomalies is based on the sea surface temperature (SST) maps derived from the Landsat TM (Thematic Mapper) and ETM+ (Enhanced Thematic Mapper Plus) images. Conversion of a digital number to spectral radiance and then to at-satellite temperatures followed the standard procedure explained in detail in [27].
With surface temperatures lower than the surrounding sea, springs' outflows are easy to detect on the SST maps. The boundary of the area of freshwater influence is defined as the greatest change in temperature, relative to the distance as proposed by Johnson [16]. SST distribution on the maps varies. Instead of selecting single SST, for each outflow point a unique SST is defined to account for the differences in surrounding SST values. The greatest influence of the freshwater is at the outlet point where SST is equal to the freshwater temperature. As the distance from the outlet point increases, the freshwater mixes with the sea and the SST increases. The mixing process between the two fluids of different temperatures and densities is complex [28] as well as the heat transfer between them. Heat transfer occurs through processes of convection and radiation, whereas oceanographic elements (e.g., thermal stratification, currents) have a great influence. However, the data needed to thoroughly analyze the mixing process between colder groundwater and warmer sea water was not available. Instead, an adjustment was made by introducing the adjusted area of influence (i.e., the area that is somewhat smaller than the area delineated) on the SST maps: where A adjusted is the adjusted area of influence (m 2 ) and A i is the size of a pixel in SST map (m 2 ). Weighting coefficients w i in the Equation (1) are used to compensate for the decreasing influence of fresh water as the distance from the outflow increases. In such a formulation, estimation of the weighting coefficients w i depends on the remotely sensed SST. Fuzziness in this definition was dealt with using a linear fuzzy membership function: weighting coefficients are expressed on a continuous scale [29] from 1 (full membership) to 0 (no membership). A suitable fuzzy membership function, µ A was then defined for each fuzzy set (each location), depending on the boundary conditions, i.e., minimum, and maximum SST [30]: where SST is the sea surface temperature (K) of the pixel and SST MAX and SST MIN are the maximum and minimum SST at the observed outlet point. Adjusted area of influence was then compared to the observed discharges to analyse relationship between the two parameters.

Hydrological Data
Discharges of the coastal springs are significant, but the measurements are scarce. The available data are principally recorded at the most downstream gauging station Sušak tvornica ( Figure 1). Since 1999, the gauging station is recording only discharges and not water temperatures. However, water temperature measurements were conducted during the individual water quality samplings in the Rječina River and several coastal springs including Kristal (database of the Croatian water management legal entity Hrvatske vode).
Apart from the several discharge measurements in the past two decades (usually during the dry periods), the coastal springs in the Opatija region are not monitored. The most frequent measurements were conducted at the Kristal Spring: depending on the hydrological conditions in dry seasons, measured discharges vary up to 1.0 m 3 s −1 . In wet season discharges from the Kristal Spring as well as the other springs in the vicinity are far greater (few tens of cubic meters per second).
With insufficient hydrological data but with available measurements of the meteorological parameters [20] as well as the geological structure and the catchment boundaries that resulted from numerous hydrogeological tracings [22], Horvat and Rubinić [31] estimated water potential of the catchment (Table 1) analyzing the size of the catchment and the average annual surface runoff derived from the spatial distribution of the average annual rainfall and air temperatures. The average annual runoff was estimated to be 13.6 m 3 s −1 . When compared to the results of the previous research (e.g., [32]) estimated discharges fit into the hydrological characteristics of wider coastal area of the Dinaric karst as well as the Aegean karst [33]. The water potential is high, approximately 10.0 m 3 s −1 or higher, apart from several recorded droughts (such as the one in 1994). High salinity of the coastal and submarine spring water, however, makes this resource difficult to use for water supply. With few water intakes located in higher altitudes on Učka mountain (above 700 m a.s.l.), only part of this water potential is being used. Average intake is estimated to be 0.070 m 3 s −1 [34], which is less than 0.1% of the total potential. Several deep boreholes (deeper than the sea level) that were drilled on Učka mountain, in Opatija's hinterland for the purposes of the coal exploitation long time ago, have been restored and repurposed to monitor groundwater levels [35], including the piezometer Tu 43/83 (Figure 1)

Landsat Images
Coastal springs and the plume from the river Rječina outflow have temperatures lower than the sea and they are usually visible in the TIR images as they appear darker Remote Sens. 2021, 13, 3737 6 of 14 than the surrounding sea. To be detectable on images, the size of the plume, in general, must be equal to or larger than the spatial resolution of the TIR images analyzed, in this case 1200 m 2 for the images recorded with the ETM+ sensor and 14,400 m 2 for the images recorded with the TM sensor. Smaller plumes may also be detected if their reflectance dominates in the cell. Furthermore, if the data derived from the images is to be used for estimating outflow quantities, it is necessary that the images are recorded at the time of the field discharge measurements. In that respect, twelve images were chosen ( Table 2): two were recorded with TM sensor at 120 m spatial resolution of the TIR band and ten with ETM+ sensor with 60 m spatial resolution.

Results
SST maps were created for each Landsat image ( Figure 2). Twelve images showing the spatial distribution of SST indicate variation of the SST depending on the season. Lower SST values occur in March and September, while the maximum SST values are related to the summer months (i.e., June to August). Depending on the difference between the surrounding SST and the temperature at the outlet, freshwater plumes are visible on the images. Freshwater plumes smaller than the spatial resolution of the images could not be detected. In that respect, thermal anomaly of the Kristal Spring was too small to be detected by the TIR sensor in five images: 15  Even though the plumes are not visible in the image, the discharges may still be present, but the size of the anomaly is too small to be detected by the sensor's resolution. Such is the case with the image recorded on 27 September 2003: the plume is not visible but the discharge, although small, was measured. To avoid misinterpretation, only the images with visible anomalies were used to estimate the adjusted area of influence. For each of them, minimum SST was derived overlying the SST map with the outlet points (Table 3). To estimate maximum SST a curve was fitted through the scatter plot SST data vs. distance from the outlet ( Figure 3) and SST MAX was set at the curve's inflection point [16].    With boundary conditions known, fuzzy membership functions (weighted coefficients) were calculated separately for both locations using Equation (2). Area of each pixel influenced by the freshwater was then multiplied with the corresponding weighting coefficient, resulting in the spatial distribution of the weighted pixel areas (Figure 4). Adjusted area of influence was then calculated using Equation (1)    With boundary conditions known, fuzzy membership functions (weighted coefficients) were calculated separately for both locations using Equation (2). Area of each pixel influenced by the freshwater was then multiplied with the corresponding weighting coefficient, resulting in the spatial distribution of the weighted pixel areas (Figure 4). Adjusted area of influence was then calculated using Equation (1) ( Table 4).

Discussion
With determination coefficient R 2 = 0.9647, correlation between the size of the adjusted area of the thermal anomaly at the location of the Rječina River's outflow and the measured discharges are very high ( Figure 5).

Discussion
With determination coefficient R 2 = 0.9647, correlation between the size of the adjusted area of the thermal anomaly at the location of the Rječina River's outflow and the measured discharges are very high ( Figure 5). This relationship, assuming the same temperature of groundwater outflows and temperature conditions in the sea at the location of the Opatija's springs, is used to estimate the springs discharges (Table 5). Kristal Spring's outflow is separated into two branches. Measurements showed that the main branch provides in average 70% of the total spring's discharge. In summer 2004 region experienced very dry and stable recessive hydrological conditions, causing the surrounding springs Slatina and Admiral to dry out. This relationship, assuming the same temperature of groundwater outflows and temperature conditions in the sea at the location of the Opatija's springs, is used to estimate the springs discharges (Table 5). Kristal Spring's outflow is separated into two branches. Measurements showed that the main branch provides in average 70% of the total spring's discharge. In summer 2004 region experienced very dry and stable recessive hydrological conditions, causing the surrounding springs Slatina and Admiral to dry out. measured at the location where outflow from the Kristal Spring is the most concentrated while the anomaly registered on the TIR image also includes the outflows from the nearby springs Admiral and Slatina (Figure 1). Surface runoff at those two springs and in~350 m length of the Opatija's coastline was not visible during the measurements. Additionally, measurements are under the influence of the diurnal sea level dynamics, which has the great impact on the diurnal outflow regime and hence diurnal exchange of the coastal springs' discharges. Reliability of estimates was tested with relative absolute error (RAE) and relative standard error (RSE): where a i denotes the areas of influence derived from the TIR images at the Rječina's outlet and p i values are the predicted discharges. Average discharges calculated from the existing data are denoted withā, while n is the number of samples. As shown in the Table 6, low RAE and RSE values suggest good prediction, but reliability of the results should still be taken with caution due to low number of samples. To further examine the results, discharges were correlated with the measured groundwater levels (Figure 6), and high determination coefficient (R 2 = 0.9601) proved the strong connection between the two parameters. Extrapolation of the groundwater levels to dry periods when thermal anomalies could not be detected by the sensor's spatial resolution resulted in small but present outflow. That was the case with the image recorded on 27 September 2003. Coarse spatial resolution (120 m) could not detect presence of the spring's outflows, but extrapolated groundwater level indicated discharge~0.138 m 3 s −1 . It is possible that with higher resolution (e.g., aerial high resolution TIR sensor) thermal anomalies will be present in the image.
The relationship between the measured groundwater levels and area of influence at the Kristal Spring outlet was tested with RAE and RSE as the measure of reliability, using the Equations (3) and (4). Here, also, the RAE and RSE values are small (Table 6) suggesting the strong relationship between the two variables. However, this too should be taken with caution considering the small number of samples.

Limitations of the Proposed Method
Even though the RAE and RSE suggest reliable results of the presented research, the proposed method faces several limitations that should be taken into consideration when analyzing the results. The research relies greatly on the accuracy of the delineation of the thermal anomalies created by the freshwater outflow on the sea surface. To overcome the fuzziness in estimating the boundary and the size of the area of influence, an approach involving the fuzzy membership function definition was applied. Such an approach depends greatly on the researcher's knowledge and skills in choosing the appropriate function and defining the boundary conditions (i.e., minimum, and maximum SST). ens. 2021, 13, 3737 12 of 15

Limitations of the Proposed Method
Even though the RAE and RSE suggest reliable results of the presented research, the proposed method faces several limitations that should be taken into consideration when analyzing the results. The research relies greatly on the accuracy of the delineation of the thermal anomalies created by the freshwater outflow on the sea surface. To overcome the fuzziness in estimating the boundary and the size of the area of influence, an approach involving the fuzzy membership function definition was applied. Such an approach depends greatly on the researcher's knowledge and skills in choosing the appropriate function and defining the boundary conditions (i.e., minimum, and maximum SST).
The accuracy of the proposed method depends greatly on the differences between the freshwater temperature and the temperature of the surrounding sea. Smaller difference affects the precision of defining the boundary of the thermal anomalies and estimating the area of influence. In the frame of the proposed method, it is a fortunate circumstance that low flows occur at the end of the dry period when the temperatures differences are the most significant, which is in favor of detecting the size of temperature anomalies.
An analysis of the relationship between the discharges with the size of the thermal anomaly, derived from the TIR images, must be based on temporal alignment of the discharge measurement and the satellite overpass. It cannot always be the case, since the visibility of the surface scanned by the sensor depends greatly on the weather conditions (i.e., passive sensors cannot penetrate clouds). Without sufficient number of discharge measurements that coincide with the satellite images, it is difficult to derive satisfactory results. The accuracy of the proposed method depends greatly on the differences between the freshwater temperature and the temperature of the surrounding sea. Smaller difference affects the precision of defining the boundary of the thermal anomalies and estimating the area of influence. In the frame of the proposed method, it is a fortunate circumstance that low flows occur at the end of the dry period when the temperatures differences are the most significant, which is in favor of detecting the size of temperature anomalies.
An analysis of the relationship between the discharges with the size of the thermal anomaly, derived from the TIR images, must be based on temporal alignment of the discharge measurement and the satellite overpass. It cannot always be the case, since the visibility of the surface scanned by the sensor depends greatly on the weather conditions (i.e., passive sensors cannot penetrate clouds). Without sufficient number of discharge measurements that coincide with the satellite images, it is difficult to derive satisfactory results.

Conclusions
Coastal and submarine springs have been analyzed as valuable water resources for decades using different methods and approaches, but they still represent an unknown as a potential resource to meet the continuously growing demand for fresh water. Based on the research that has been conducted so far on Opatija's springs, we can assume that underlying karst aquifers can provide enough water to meet those demands. However, further analyses based on more detailed data will give better insight into the water potential of the area.
Differences in water temperature between freshwater outflows and the surrounding sea enabled the use of TIR images to estimate outflow quantities of the coastal water resources. Coastal and submarine springs are phenomenon easy to detect in TIR image but the accuracy of the estimated size of their thermal signature depends greatly on the ratio of the discharge and the spatial resolution of the sensor. Thematic mapper's spatial resolution of 120 m can comply with the accuracy requirements in case of larger outflows, but when the discharges are low, as was the case with the image recorded on 27 September 2003, the anomaly cannot be detected by the sensor.
The first step was to analyze relationship between the size of the River Rječina's thermal anomalies and the discharges measured at the time of the satellite overpass. This relationship was then used to estimate outflows of the springs in Opatija. Unfortunately, apart for the image recorded on 29 September 2004 when the measured discharge was 0.504 m 3 s −1 , for the rest of the images the time of the discharge measurements did not coincide with the satellite overpass. As the additional control of the outflow characteristics and the size of the thermal anomalies, thermal anomalies were correlated with the groundwater levels. Undoubtfully, the size of the thermal anomalies is influenced by the connection between the groundwater and the sea, their mixing process as well as the temperature changes.
Described approach could not be based on a thorough deterministic analysis of mixing process of fresh and sea water (i.e., two fluids with different temperatures and their thermal signatures) since it would require detailed input data which, in this case, was not available. However, despite simplifications, the proposed method provided realistic estimation of the discharges. Further development should put an emphasis on planned field measurements aligned with satellite overpasses. That will certainly improve understanding of the coastal springs' outflow regime, the physics of the mixing process, and the heat transfer, as well as the method's accuracy. The estimation of the discharges at the non-monitored springs gives a rough approximation of the water potential in the area, but it is a good starting point for further investigations of the non-monitored coastal springs.  Data Availability Statement: Publicly available Landsat datasets were analyzed in this study. This data can be found here: https://earthexplorer.usgs.gov/; The hydrological data presented in this study is available on request from The Croatian Meteorological and Hydrological Service.