Thermohaline Temporal Variability of the SE Mediterranean Coastal Waters (Israel) – Long-Term Trends, Seasonality, and Connectivity

The high variability of coastal waters together with the growing need for assessing the state of the marine coastal ecosystem, require continuous monitoring at exceptional resolution and quality, especially during the Anthropocene changing seas. We perform a comprehensive analysis of a decadal (March 2011 to June 2021) thermohaline variability of the East Levantine Basin (LB) coastal waters (continuous measurements), its predominating temporal trends and their linkage with atmospheric forcing and advection. We identify statistically significant long-term warming and salinification trends with yearly rates of 0.048°C and 0.006, respectively. Through the use of the X11-ARIMA method temperature and salinity inter-annual trends are examined and associated with previously published open ocean dynamics as well as model reanalysis. We study the linkage between Northern and Southern coastal locations, and identify the along shore northward current as a primary cause of positive temperature anomalies arriving from the south. The coastal salinity long-term trend demonstrates a connection to local precipitation. A less coherent seasonal sequence is found with a bimodal behavior, where, salinity values drop in August on several summers. This drop is attributed to the intensification of the along shore current in the period of June-July, potentially advecting more Atlantic Water. The observations presented here emphasize the relatively strong coupling between coastal water and the open ocean, the influence of the general surface circulation of the LB on the coastal zone and the faster response time and higher sensitivity of the coastal environment to atmospheric forcing.


INTRODUCTION
The Mediterranean Sea (MS), being a marginal sea, relatively small in size, with a restricted connection to the Atlantic Ocean has a faster response to atmospheric forcing compared to the open ocean. It has therefore been repeatedly acknowledged as a prime location for climate change research (Giorgi, 2006;Giorgi and Lionello, 2008;Schroeder et al., 2012Schroeder et al., , 2016Pastor et al., 2018, UNEP/MAP andPlan Bleu, 2020). Atlantic Water (AW) entering the MS at the surface through the Gibraltar Strait undergo continuous transformation due to strong radiation and intensive evaporation exceeding runoff and precipitation (Robinson et al., 1992;Hamad et al., 2005;Malanotte-Rizzoli et al., 2014;Estournel et al., 2021). In summer, the Levantine Surface Water (LSW) occupies the upper part of the water column with salinity (S) as high as ∼39.7 and temperature (T) of 30 • C (Lascaratos et al., 1999;Gertman and Hecht, 2002). The coastal waters of Israel, at the eastern end of the LB, occupy the continental shelf and present some of the highest T and S values of the region. These uniquely high thermohaline values, coupled with extreme oligotrophy (Herut et al., 2000;Kress et al., 2014), make the marine ecosystem of the eastern LB highly sensitive to thermohaline changes, as identified for the global coastal zone (Gissi et al., 2021). Sub-mesoscale dynamics may also have a crucial influence on the physical and biogeochemical characteristics of a confined body of Robinson et al. (1992) coastal water. Such a formation of localized niches separated by transport barriers, was presented by Efrati et al. (2013) and was suggested to have a key role in the dynamics of planktonic communities. In addition to global and local climate induced processes, the coastal marine environment of the MS is affected by a variety of local anthropogenic stresses such as economic development, population increase, and changes in land use patterns (Micheli et al., 2013;Ramírez et al., 2018;UNEP/MAP and Plan Bleu, 2020).
Many authors have identified a SST warming trend of the Mediterranean basin or sub basins, either using in situ measurements (Rixen et al., 2005;Belkin, 2009;Ozer et al., 2017) or longer SST series originating from remote sensing data (Criado-Aldeanueva et al., 2008;Nykjaer, 2009;Skliris et al., 2012;Shaltout and Omstedt, 2014;Pastor et al., 2018Pastor et al., , 2020El-Geziry, 2021;Pastor, 2021). The reported warming rates vary based on the investigated area, period and water mass and are in the range of 0.024 to 0.12 • C year −1 ( Table 1). A salinification trend of the MS has also been repeatedly reported with yearly rates ranging from 0.0014 to 0.008 (Krahmann and Schott, 1998;Vargas-Yanez et al., 2010;Ozer et al., 2017;Schroeder et al., 2017;Skliris et al., 2018) (Table 1). Skliris et al. (2018) stated that the MS salinification trend is a result of changes in the regional water cycle including a long-term increase in evaporation coupled with a decrease in precipitation (Mariotti et al., 2002;Mariotti, 2010;Romanou et al., 2010;Skliris et al., 2012) and the reduction of river runoff both through climate change and the damming of major rivers (Rohling and Bryden, 1992;Skliris and Lascaratos, 2004;Skliris et al., 2007;Ludwig et al., 2009). We note that the highest warming and salinification trends (of 0.12 • C year −1 and 0.008, respectively) were identified in the southeastern MS LSW water (Ozer et al., 2017). Tanhua et al. (2019) pointed out to the gaps in ocean observing coverage and the ever-growing use of the marine environment for both open ocean and coastal applications and the need to supply these users with operational services. Sustainable, long-term in situ monitoring programs are required for the improvement of coastal forecasts (Brenner et al., 2006), the calibration of remote sensing data and the enhancement of knowledge of the coastal water dynamics and their effects on the coastal ecological system. Pinardi and Coppini (2010) articulated: "The high variability of ocean currents and the need for the assessment of the state of the marine environment require a continuous monitoring of the ocean environment at unprecedented resolution and quality. On the other hand, scientific discoveries require continuous observations of the essential state variables of the system in order to detect new basic mechanisms and processes." In this study, we perform a comprehensive analysis of the thermohaline variability of regional coastal waters, in order to better characterize the coastal marine environment, its predominating temporal trends and their linkage with atmospheric forcing and advection. Using the time series of the Israeli coastal monitoring stations, unprecedented for this region in terms of temporal resolution, we investigate coastal variability and dynamics in the eastern MS, and elucidate on the issues of: (1) long term and inter annual trends, including a comparison to open ocean findings, (2) the predominance of the seasonal cycle, and (3) the connectivity of coastal waters along the Israeli coast.

The Israeli Coastal Monitoring Stations
In the framework of the National Monitoring Program of Israel's Mediterranean waters (INMoP), Israel Oceanographic and Limnological Research (IOLR) operates two continuous observation posts with near-real-time data transfer to IOLR server 1 . These stations are located at the westernmost edge of the coal terminals in Hadera (north) and Ashkelon (south), 2.2 kms offshore at water depths of 26 m (Figure 1). Hadera station (MedGloss # 80) was established toward the end of 1991, upgraded through the years and currently includes a bottom mounted, up looking, Teledyne RD Instruments 600 KHz Workhorse Monitor Waves ADCP located ∼100 m westward from the terminal end at water depth of 26 m. A Paroscientific Inc. Intelligent Digiquartz 8800 sea level sensor (PAROS) and a SBE16plusV2 Sea-Bird CTD (added to Hadera station in 2011) are installed on the westernmost terminal sporting pole at a depth of 11-12 m. The CTD is equipped with standard T/C duct and pressure sensors together with auxiliary sensors of dissolved oxygen (SBE63 Optode, Sea-Bird electronics) and fluorescence and turbidity sensors (FLNYU, WetLabs). A twin station at Ashkelon was established in mid-2012. Continuous CTD Raw data of T, S, Pressure, Oxygen, Fluorescence and Turbidity is collected at a high temporal sampling interval of 10 min. In order to ensure the systems durability against fouling it incorporates several antifouling measures including toxin emitting capsules, Bio-Wipers and an in-house developed underwater case minimizing the systems exposure to the surroundings. A routine maintenance procedure is preformed once every 2 months (for each of the stations) which includes the replacement of the CTD to a freshly authenticated instrument and the collection of reference measurements for data validation. Recovered CTDs go through a robust cleaning procedure (following Sea-Bird electronics guidelines) and are examined for biases both pre-cleaning (in comparison to the collected reference data) and post-cleaning in lab-performed comparisons. CTDs are sent to the manufacturer for calibration regularly every 2 years or earlier in the case of an identified drift of one of the sensors (i.e., over three times of the instruments initial accuracy of ±0.001 • C for T and ±0.004 for S). The Hadera station CTD dataset encompasses the period of March 2011 to June 2021. Due to various infrastructure difficulties, Ashkelon station experienced large data breaks since its establishment in 2012. For this reason, Ashkelon CTD measurements will only be used here as a basis of comparison to Hadera station data. Hadera station is located offshore (2.2 km) from the discharge of effluents at the shoreline of the nearby power plant and desalination facility. In order to assess the possible effect of these anthropogenic activities on our measurements, we examined the thermohaline marine monitoring report conducted in the area over the period of 2013 till 2020 (Cohen, 2021). T and S anomalies presented in the report, both for surface and near bottom values, are generally confined to a distance of ∼500 m from the shoreline and have a consistent southward spreading pathway, not effecting the station. Wood et al. (2020) used numerical simulations to study the distribution of desalination brines along the Israeli coast and the spreading of annual mean S anomalies in Hadera limited to under 1 km from the discharge location. To further confirm that Hadera monitoring station can be considered a good representative of northern Israeli coastal waters, a comparison of the daily averaged T time series to Achziv and SST datasets (see below) was performed. This comparison presented a statistically significant fit with a RMSE of 0.46 and Pearson R of 0.99 for the Achziv dataset and a RMSE = 0.66 and Pearson R of 0.98 for the SST dataset (Supplementary Figure 1).

Achziv Islands
Over a period of 13 months, from May 2013 to May 2014 Several HOBO T loggers were maintained at two sites in Achziv Islands, located at the northern end of the Israeli Mediterranean coast (Figure 1). The averaged data of two loggers deployed at a water depth of 10 m (north and south Islands,

Haifa Section Cruises
Israel Oceanographic and Limnological Research conducts the biannual Haifa section cruise covering stations at water depths of 15-1900 m, extending 90 km off the coast of Israel, as a part of INMoP. A Sea-Bird SBE911plus CTD, interfaced to a SBE Carousel is used to collect continuous profiles of pressure, T, S, dissolved oxygen and fluorescence. Equipment, procedures, and calibrations are detailed in Ozer et al. (2017Ozer et al. ( , 2019. We compare our results to the updated analysis of long-term and inter-annual trends as presented in Ozer et al. (2017) adding the years of 2015 to 2020.

Meteorology
Air T, wind, and precipitation data from four chosen stations in the vicinity of Hadera (Zikhron Ya'akov, Hadera port, Ein HaHoresh, and Ein Carmel) was averaged to attain a daily time series for the studied period of 2011 to 2021. The use of multiple stations covering ∼30 km along the Israeli coastal plain allowed to reduce data breaks, while maintaining the climatic traits of the region. All data originates from the Israel Meteorological Service (IMS) archive available online 2 . 2 https://ims.data.gov.il/ims/7

Copernicus Monitoring Environment Marine Service Data
We use two products provided by the E.U. Copernicus Monitoring Environment Marine Service (CMEMS). (1) Daily gap-free remotely sensed L4 SST produced and distributed in near-real time by the Consiglio Nazionale delle Ricerche (Buongiorno Nardelli et al., 2013). We extract the T time series from the high spatial resolution product (HR 0.0625 • ) in the closest position to Hadera station (32.5 • N, 34.81 • E).
(2) The Mediterranean Sea Physical Reanalysis generated by a numerical system composed of a hydrodynamic model, supplied by the Nucleous for European Modeling of the Ocean and a variational data assimilation scheme (Escudier et al., 2020). Here we extract the daily time series of T and S for a depth of 10.5 m from two grid points of the reanalysis. The first representing coastal water closest to Hadera station (referred hereafter as CRA, 32.48 • N, 34.83 • E) and the second from an open water location (Z > 1200 m) northwest of Hadera station (referred hereafter as OWRA, 32.69 • N, 34.5 • E).

DeepLev Mooring Station
In November 2016, the first moored station at the Eastern LB was deployed at 1500m water depth, ∼50km offshore Haifa, Israel (Katz et al., 2020). For this study, we use the T time series from a downward looking Teledyne RD Instruments 300 KHz workhorse ADCP. The ADCP is installed on the upmost buoy of the mooring system, at a depth range of 20-30 m (depending on deployment) and sampling at a 15-min interval. The ADCP integrated T sensor has a reported precision of ±0.4 • C and a resolution of 0.01 • C. Inter-calibrations with CTD casts taken on deployment and recovery events yielded a mean variance of 0.015 • C with a 95% confidence interval of ± 0.056 • C (n = 11).

Data Analysis
At a preliminary stage the various time series data sets described above were aggregated into a unified database and synchronized based on the time of measurement. Next daily and monthly averages were calculated for all available parameters. T hourly averages grouped by month were attained to infer the range of the diurnal cycle ( Supplementary  Figures 2, 3). In addition, daily standard deviation averages of the complete CTD dataset (over 400,000 measurements in total) were calculated and yielded 0.19 • C and 0.032 for T and S, respectively. The seasonal Mann-Kendall test was used to examine the T and S CTD datasets for statistically significant long-term trends (Hirsch et al., 1982;Hussain and Mahmud, 2019).

Ekman Transport Calculation
To ensure that the IMS wind dataset adequately represents off shore conditions, the daily averaged wind was compared with a time series of the ECMWF ERA-5 reanalysis at position 32.5 • N, 34.75 • E and presented a strong agreement with R 2 of 0.85 and 0.87 for U and V components, respectively (not shown). Next, following Cropper et al. (2014) across and along shore wind stress (τ) components were derived from averaged daily wind speeds.
where, W ac and W al are across and along shore winds components achieved by a 25 • rotation of the original data, according to the shoreline inclination. ρ a is the air density of 1.22 kg/m 3 and C d the typical dimensionless drag coefficient of 1.3 × 10 −3 . Next Ekman transport (Q) components were derived from the wind stress field.
where, ρ w is the density of sea water of 1025 kg/m 3 and f is the Coriolis parameter defined as f = 2 sin(ϕ) (where, is the earths angular velocity of 7.292 × 10 −5 s −1 and ϕ is the latitude).

Filling Temperature Data Gaps
In order to complete missing CTD T data in the Hadera dataset a multi parameter regression (MPR) was performed. The MPR procedure was based on T records of the ADCP and PAROS, located in the immediate vicinity of the CTD, as well as the close by (∼1.3 km) SST dataset (Supplementary Figure 4). As a preliminary step the similarity of these datasets was assessed where they overlapped. Daily differences between the CTD T and the additional T datasets were computed and yielded 95% confidence intervals of 0.024, 0.016 and 0.028 • C for ADCP, PAROS and SST, respectively. Pearson R correlation ranged from 0.97 to 0.99 (the lowest found between the ADCP and SST), thus all datasets were accepted for the MPR procedure. Depending on the availability of data from these datasets a matrix of coefficients was calculated and used to achieve a complete daily data set of Hadera station over the period of March 01, 2011 to June 08, 2021 (Supplementary Tables 1, 2). This method was not implemented for Ashkelon station as the data gaps there were considered too large to overcome. In the analysis of S variations, the monthly averaged data set was used, as no other reference was available in order to complete missing daily data.

Time Series Decomposition
For the determination of seasonal, trend and irregular components we follow Pezzulli et al. (2005) and use the filter based X11-ARIMA method (X11), a variant of the Census Method II Seasonal Adjustment Program, assuming an additive model.
where, X t is the original time series, T t is the trend component, S t is the net seasonal component and I t is the irregular variation. The X11 is a filter based method for seasonal adjustment which follows the 'ratio to moving average' procedure first described by Macaulay (1931) and commonly referred to as "Classical Decomposition." In principle X11 is comprised of three steps which are implemented iteratively: (1) Estimation of the trend by a moving average, (2) Removal of the trend from the time series leaving the seasonal and irregular components, and (3) Estimation of the seasonal component using moving averages to smooth out the irregulars. After three iterations, the irregular term is finally calculated by removing the final trend and seasonal terms from the original time series. Key features of the X11 decomposition are that the trend is void of the annual cycle and any higher frequencies and that the seasonal term is defined locally in time, thus allowing for inter-annual variations in the seasonal cycle (Pezzulli et al., 2005). For the T decomposition, we use the daily averaged T and apply the moving averages windows of 30 × 365 and 30 × 30 for the trend and seasonal estimations, respectively. Due to the data gaps in the daily S dataset (Figure 2E), we implement the X11 decomposition on monthly averaged S and apply the moving averages windows of 2 × 12 and 2 × 2 for the trend and seasonal estimations, respectively. Similarly, the X11 method was also implemented on the daily CMEMS model reanalysis datasets (CRA and OWRA) in order to compare of the trend terms.

Synoptical Anomalies Delineation
The daily T irregular term (acquired by means of the X11 method) was further investigated by applying a threshold filter, parsing through this dataset and marking positive and negative amplitudes and durations. In order to set the thresholds for significant synoptical anomalies, we initially examined the amplitude of the diurnal cycle, which ranged from 0.15 • C in winter to 0.3 • C in summer for T (Supplementary Figure 2) and from 0.001 in winter to 0.003 in summer for S ( Supplementary  Figure 3). Since the S diurnal amplitude was at the level of sensor accuracy, we opted to use the values of ±2 times the average daily standard deviation (0.38 • C for T and 0.064 for S) for the anomalies delineation. In every incident that the irregular term crossed the threshold value, the duration of the anomaly was recorded together with the maximal amplitude. Cases of missing data whilst within a given anomaly event were omitted as the anomaly period could not be established. For the S time series, where the X11 method was performed on monthly averages, the following approach was used to calculate daily anomalies. The interpolated sum of trend and seasonal terms (of the monthly X11 analysis) were subtracted from daily averaged CTD data where it was available. The resulting residuals time series was then screened following the same procedure described above for T.

Long-Term Trends and Inter-Annual Fluctuations
A statistically significant warming trend was identified over the observed period (2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020)(2021) in Hadera station (11-12 m depth) using the seasonal Mann-Kendall test, with a yearly rate of 0.048 • C (p = 0, z = 10.92) ( Figure 3A). While the coastal dataset has the handicap of covering a relatively short period of time, the achieved trend is in agreement with previous studies (Table 1) Figure 5), the Hadera trend is similar to that achieved in the analysis of Levantine Intermediate Water (LIW) (0.03 ± 0.01 • C year −1 ), yet lower than the trend observed in surface water (0.13 ± 0.05 • C year −1 ). We attribute the disparity between Hadera and LSW trends to the fact that the analysis of LSW was limited to casts performed during the peak of the warming period (end of summer). In contrary, the analysis of LIW was based on bi-annual data (i.e., summer and winter casts). Furthermore, we consider the LIW to be less affected by mesoscale and seasonal processes, thus it is more integrative and reflective of long-term and inter-annual variability. A significant positive S yearly trend of 0.006 (p = 4.06 × 10 −7 , z = 5.06) was identified in the Hadera CTD dataset (Figure 4A). This trend is at the higher end of previously published salinification rates for the MS (Table 1), similar to the rates attained in the analysis of LSW and LIW (Ozer et al., 2017 and Supplementary Figure 5) and in the western MS (Schroeder et al., 2017). Pastor et al. (2018) reported upper layer salinification trends exceeding 0.004 year −1 in the easternmost MS, in agreement with our results. This trend is consistent with the increase in latent heat loss in the eastern MS since the early 1970s (Mariotti, 2010;Skliris et al., 2012) and the large reduction in Nile river runoff since the early1960s (Skliris and Lascaratos, 2004;Suari and Brenner, 2015). The short duration of the dataset analyzed in this study does not support the appraisal of our findings within a broader framework of surface T and S natural variability of longer decadal and multi-decadal scales. These will need to be addressed in the future, as data accumulates from the coastal stations.
Both T and S trends achieved through the X11 method show inter-annual fluctuation with a clear maximum in the beginning of 2018, which seemingly corresponds with the previously established inter-annual cycle of the Levantine basin in connection with the BiOS mechanism (Civitarese et al., 2010a,b;Gacic et al., 2010;Ozer et al., 2017). Smaller, less prominent extremes are evident at the start of 2013 for T and 1 year later for S (Figures 3B, 4B). A close similarity is found between the T and S X11 trend terms of Hadera and CMEMS reanalysis datasets (Figure 5). The T Trends of Hadera and CRA are consistently higher than that of the OWRA, displaying the stronger response and higher sensitivity of coastal waters to the atmospheric forcing and their potential impact on the marine environment ( Figure 5A). The S extrema of 2014 coincides with low precipitation in the winter of 2013-2014 ( Figure 4B). Indeed, a correspondence of the coastal S seasonal minimum to the total winter precipitation can be observed throughout the examined period. To demonstrate the influence of precipitation on coastal S, we calculate the regression between S and cumulative precipitation in two winters seasons (covering the months of November to March) where salinity measurements were rarely interrupted. A strong relation is found with R 2 of 0.9324 and 0.885 for the winters of 2014-2015 and 2019-2020, respectively. On several occasions, drops of S values below the regression line are supported by a passing winter storm with several days of continuous high precipitation levels (Figure 6). These results point to a strong linkage of the surface S with   the local rainfall. As S trends of CRA and OWRA follow the same temporal patterns (Figure 5B) we believe this is a close representation of the response to annual changes in the LB regional water budget.
The overall agreement of T and S time series of the OWRA with Hadera station in terms of variance and lags (Table 2), as well as the similarities of the T and S X11 trends discussed above, indicates a relatively strong coupling between coastal water and the open ocean. In contrast, a large difference is found between T time series of Hadera and DeepLev mooring with a time lag of 43 days ( Table 2). We attribute this inconsistency of DeepLev data to the fact that the T in this dataset is sampled at a deeper values were calculated by subtraction of the specified dataset from Hadera data. Maximal pearson R correlation value presented, yielded at specified offset from the Hadera station dataset (positive offset values signify that Hadera time series leads the compared dataset). level (∼30 m), below the surface mixed layer in summer, thus the changing level of the vertical mixing has a significant effect on the sampled data.

Climatology and the Seasonal Cycle
A clear seasonal cycle of surface T is evident in our data with average peak values of 29.86 • C in August and 17.62 • C in February, representing an average offset of +7 • C in summer and −5 • C in winter from the overlaying long-term trend (Figures 7A, 3B,C). Extreme water T seasonal anomalies are in concordance with large amount of days with air T below 10 • C for negative anomalies and above 25 • C for positive anomalies (not shown). Lag correlations performed between the average daily air T and CTD thermohaline values confirms the expected seasonal effect of atmospheric forcing on coastal water characteristics, with the air T time series leading by 15 and 50 days before water T and S, respectively (Supplementary Figure 6). The S dataset presents a less coherent seasonal sequence with a bimodal behavior in the summer seasons of 2014, 2016, and 2017, which is apparent also in the climatological S (Figures 4C, 7B). The drop in salinity in August can be attributed to the intensification of the along shore current in the period of June-July (Figure 7C), potentially advecting more AW which are lower in their S value. Wind derived along shore Ekman transport is kept close to zero in the summer months ( Figure 7D), suggesting that the strong current values are mainly of geostrophic origin. Rosentraub and Brenner (2007) reported maximal northward currents in the summer season FIGURE 8 | Kernel density estimation of amplitudes (A,D) and durations (B,E) for anomalies of temperature and salinity, respectively. Averaged temperature (C) and salinity (F) amplitudes of positive (red) and negative (blue) anomalies grouped by their duration. Standard error bars are shown where anomalies of a specific duration occurred more than once. over the Israeli continental shelf confined to the upper layer as well as an along-slope baroclinic jet during summer and early winter. The authors related their findings to the general surface circulation of the LB introducing AW into the basin through several established pathways. Mauri et al. (2019) has shown an increase in average dynamic topography during the months of June to August in latitudes of 33-34 N in further support of our hypothesis. Estournel et al. (2021) demonstrated that in summer, the shallower current flows closer to the coastline while in winter it flows above the 1000 m isobath. The climatological extrema S values of 39.47 and 39.01 are reached at the months of October and March, respectively. An extreme seasonal cycle is observed in 2019 with minimal S of 38.8, which we relate to the high precipitation of the preceding winter as discussed earlier. However, this minimum can also be related to an enhanced modified AW flux into the LB, following the reasoning of the BiOS mechanism. An unusual fast recovery brings S as high as 39.63 by October 2019 and is not supported by the meteorological datasets.

Synoptical Anomalies Analysis
T anomaly amplitudes range from +2 to −2.5 • C with modal values reached at −0.67 • C and 0.69 • C for negative and positive anomalies, respectively ( Figure 8A). The range of T anomalies duration is between 0 and 16 days with a maximal Kernel density at 1.5 days ( Figure 8B). Similar anomaly behavior is observed for S, with a slight preference to negative anomalies and modal Kernel densities at −0.141 and 0.136 for negative and positive anomalies, respectively (Figures 8D,E). Averaged T anomaly amplitudes grouped by durations show a relation between the duration of the anomaly and its corresponding amplitude, the longer the anomaly lasts, the higher the amplitude. This apparent dependency is kept for both positive and negative anomalies, but is more coherent in the case of the FIGURE 9 | Correspondence between along shore currents and positive temperature anomalies. (A) The along shore current time series is filtered for increasing positive temperature anomaly thresholds (shown in the X axis). The resulting filtered along shore current dataset is averaged (blue line, vertical bars show standard error) and the percentage of incidents where northward currents coincide with the positive temperature anomaly is calculated (gray bars, total number of data points above the temperature threshold). (B) Lag correlation between temperature and along shore current residuals time series. Peak synchrony marked by red dashed line.
latter. A similar, yet less clear relation can be observed for S. For both parameters, the majority of anomaly events last up to 10 days (Figures 8C,F). No significant seasonal T and S anomalies were identified (Supplementary Figure 7). We additionally perform an analysis of heat waves following the definition of the IPCC, 2019: "marine heatwaves are defined when the daily sea surface temperature exceeds the local 99th percentile over the period 1982 to 2016." The 99th T percentile for the period of measurements in Hadera was calculated, resulting in 30.56 • C. Running this filter brings up several heat waves (Supplementary Table 3) with somewhat higher maximum T toward 2018, yet our dataset is too short to identify any patterns.

Along-Shore Connectivity
Variance statistics and lag correlations calculated for T and S comparing the Hadera station data to available datasets (Figure 2) are presented in Table 2. No offsets were found between the 3 stations (Hadera, Ashkelon and Achziv) demonstrating that the 180 km strip of Israeli coastal waters undergo seasonal and inter-annual changes in synchrony. Ashkelon station was found to be 0.36 • C warmer than Hadera station, a fact we relate to its southerly location (100 km south of Hadera station) and to extended continental shelf at the south of Israel (over 20 km in width). A convincing relation was found between the positive T irregular term (from the X11 analysis) and the along shore currents, where the majority of T irregular values are explained by northward currents, the higher the T anomaly, the stronger the current (Figure 9A). A zero time offset yielded a maximal correlation between the irregular T and along shore currents ( Figure 9B). These results suggest that the majority of positive T anomalies are a result of northward transport of water along the Israeli continental shelf, advecting warmer water from the south. No T gradient was found between Hadera and Achziv (distanced ∼70 km, Table 2), perhaps signifying that the northern shallow waters of Israel share a stronger connection maintained by coastal dynamics. T comparisons of CRA and Hadera T data yielded similar results with a mean difference of 0.15 • C, high correlations and an offset of 2 days. The CRA S was found to be in average 0.08 below Hadera station. As both CRA and SST datasets are originated from westerly, deeper positions (compared with Hadera station), we infer that atmospheric forcing is enhanced at shallower, more easterly positions and is communicated toward the west over the continental shelf within synoptical timeframe. The large offset of S, where Hadera leads the CRA by 16 days remains unclear. No significant variance in S was found between Hadera and Ashkelon stations ( Table 2).

CONCLUSION
The continued study of coastal ocean state, variability, health and climate change through environmental monitoring programs plays an important role in broadening our knowledge and understanding of the changing seas, raises the public awareness, informs policy makers and supports the management of the coastal areas. The importance of sustained ocean and coastal observation, aimed at advancing our knowledge of the ocean and its link with climate, promoting the development of early warning systems, and enabling decision-making based on the best available science has become a global consensus (Tanhua et al., 2019). In this study we make use of continuous observations providing exclusive high temporal resolution data in the Eastern coast of the LB and address the issues of (1) long term and inter annual trends, (2) the coastal seasonal cycle, and (3) The connectivity of coastal waters along the Israeli coast. Significant long-term positive trends of T and S are validated by use of the seasonal Mann-Kendall test and compared to previously published findings in other areas of the MS. The T trend of 0.048 • C identified in this study supports the exceptionally high warming rates reported in the eastern MS (Nykjaer, 2009;Skliris et al., 2012;El-Geziry, 2021) (Table 1). Filter based time series decomposition allows the identification of inter-annual T and S fluctuations associated with previously published open ocean dynamics (Ozer et al., 2017) as well as with CMEMS physical model reanalysis. Our results of higher coastal T trend term (Figure 5A), the impact of local winter rainfall on S values (Figure 6), and the sensitivity of T anomalies to the prevailing northward alongshore current (Figure 9) demonstrate the faster response time and higher sensitivity of the coastal environment to atmospheric forcing. Nonetheless, we infer a relatively strong coupling between coastal water and the open ocean in terms of long-term, inter-annual and seasonal evolution of the surface thermohaline properties. Finally, we stress the importance of the persistence of these monitoring efforts and the introduction of additional sustainable coastal waters stations, which will improve the spatial coverage enabling a better interpretation and understanding of the predominating processes.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: Israel Marine Data Center -https://isramar.ocean.org.il/isramar2009/.

AUTHOR CONTRIBUTIONS
TO led manuscript drafting with the help of BH, IG, and HG. TO and IG contributed to the data acquiring. HG and BH contributed to the interpretation. All authors approved the submitted version.

FUNDING
This study was supported by the Israel Ministries of Energy and Environmental Protection through the National Monitoring Program of Israel's Mediterranean waters. This study is in fulfillment of a Ph.D. thesis of Tal Ozer at the University of Haifa.