Oceanographic Research in the Thermaikos Gulf: A Review over Five Decades

Abstract

The Thermaikos Gulf (TG) is a semi-enclosed, river-influenced, marine system situated in the eastern Mediterranean Sea, sustaining both urban coastal regions and ecologically preserved natural areas. Facing a plethora of environmental and anthropogenic pressures, the TG serves as a critical nexus where human activities intersect with marine ecosystems. The quality and health of the TG’s marine environment are tightly linked to the socioeconomic activities of the coastal communities comprising approximately 1.5 million inhabitants. The main features of the TG’s environmental dynamics and ecological status have been scrutinized by dedicated research endeavors during the last 50 years. This review synthesizes the seminal findings of these investigations, offering an evaluation of their contribution to research, their present collective impact, and their trajectory toward the future. A severe deterioration of the TG’s environmental quality was detected in the 1970s and 1980s when the treatment of urban wastewater was completely absent. A steady trend of recovery was observed after the 1990s; however, so far, the goal of a “good environmental state” mandated by national legislation and European directives has not been achieved. A clear reduction in research was detected after 2010, associated with the recession of the Greek economy, following the “golden period” for research in the TG from the mid-1990s until the late 2000s. The most important research gaps and uncertainties are discussed, while specific targeted recommendations for the improvement of monitoring and understanding of the physical, biochemical, and ecological state of the gulf are provided: (i) increase in permanent observational stations (temporally and spatially); (ii) inclusion of all major environmental parameters; (iii) monitoring of the quantity and quality of all land-originated freshwater discharges; (iv) monitoring and management of important aspects of the marine environment that have received minimal attention in the past (e.g., coastal erosion, plastic pollution); (v) development of accurate prediction tools (e.g., numerical techniques) to support first-level responders and efficient management; (vi) establishment of a supervising public entity that would support the holistic overview and management of the entire TG. These suggestions are directed at overcoming the existing uncertainties in the knowledge of the TG, safeguarding its ecological integrity and its role as a crucial link to marine biodiversity and sustainability in the Mediterranean basin.

1. Introduction

The Thermaikos Gulf (TG) is a semi-enclosed microtidal coastal basin of the mid-latitude east-central Mediterranean Sea, influenced by multiple (transboundary and local) riverine inflows, and connected to the North Aegean Sea through its southern boundary. Thessaloniki city is located at the gulf’s northern tip and is the second largest metropolitan area of Greece with a population of almost 1.5 million residents. The TG is of particular importance, and regional, national, and international interest due to the large number of inhabitants, the various socioeconomic activities, and the existence and preservation of several ecologically protected areas located along its coastal zone [1]. It is a unique system of complex, natural processes that is severely influenced by human-related stressors, especially those involving dissolved pollutants or nutrients and particulate phase matter within the water column [2].

Six Special Areas of Conservation (SACs) of the Natura 2000 network have been assigned to the TG, including coastal and transitional waters, lagoons, and wetlands, highlighting its ecological importance. Highly valuable and unique marine ecosystems, such as seagrass meadows, stony and soft coral beds, saltpans, and saltmarshes, are deployed across the TG, hosting increased biodiversity, including strictly endemic species. A large number of algae, invertebrates, fishes, birds, and mammals have been reported in research studies considering TG marine life. This fact places the TG among the productive marine gulfs of Greece; owing to the extensive riverine system inflows, fisheries, aquaculture, and mussel farms are intensively developed. Posidonia oceanica seagrass meadows, which are among the most valuable coastal ecosystems, are broadly identified in the TG [3]. P. oceanica is one of the priority habitats of the European Union’s Habitats Directive [4], and it is protected by the Barcelona Convention [5]. However, cultural eutrophication is a major cause of seagrass decline in the Inner TG. Phytoplankton blooms and red tides affecting Posidonia are recurrent phenomena in the northern parts of the TG [6].

The first published research about the monitoring of the TG’s marine environment dates back to 1981 [7] and describes the occurrence of pollutant trace metals in sediments. In the following years, many numerical and observational studies were performed to investigate the physical, biological, geological, and chemical conditions, covering either the entire gulf or focusing on specific coastal sites of environmental and socioeconomic interest. Herein, we review the environmental state of the TG as derived from approximately five decades of field and numerical studies (1970s–2020s). The motivation for this retrospect study is to gather the main findings of previous research efforts and highlight research gaps over a socially and ecologically important Mediterranean marine area. Another important goal is to provide grounded recommendations for future monitoring and research toward reducing deficiencies, improving the ecological quality of the gulf, and designing optimal and holistic management plans.

A general description of the study area is presented in Section 2. Section 3 provides a statistical analysis of the observations dataset (monitoring activities, numerical simulations, satellite-derived observations, and laboratory experiments) comprising published research conducted from the 1970s to the present (approximately five decades). The most important findings associated with the oceanographic disciplines analyzed (physical, geological, chemical, biological) are described in Section 4, while the overall environmental quality state of the TG, the potential impacts of climate change, and the main research gaps and recommendations are discussed in Section 5. Additionally, an acronym list is provided in Glossary, the main characteristics used in the literature review are given in Appendix A, and Appendix B includes the typology of the TG’s benthic communities.

2. Description of the Study Area

The TG covers an aquatic area of approximately 3300 km² and extends from Thessaloniki city to Cape Possidi (Chalkidiki peninsula) in the south. It consists of three main sub-basins: the Inner TG, extending north from Cape Megalo Emvolo (22.82° E, 40.48° N); the Central TG, located between Cape Megalo Emvolo and Cape Epanomi (22.88° E, 40.37° N); and the Outer TG, southwards from Cape Epanomi (Figure 1a). Moreover, the Inner TG can be divided into two smaller basins, north (Thessaloniki Bay) and south (Central Gulf of Thessaloniki) of Cape Mikro Emvolo (22.93° E, 40.58° N; Figure 1b). The Inner TG is characterized by depths lower than 40 m, while the Central TG is also relatively shallow with deeper zones around 50 m, along its southern boundary. Finally, the Outer TG gradually reaches depths of 200 m towards its continental shelf break with the Aegean Sea and the deep basin of Sporades (Cape Possidi; Figure 1a).

The hydrography and the seawater quality of the TG are both affected by several factors, both environmental (e.g., meteorological conditions, freshwater discharges, and water mass exchanges with the open sea) and anthropogenic (e.g., encumbering human activities, marine and fluvial pollution, and coastal erosion due to waterfront infrastructure) [1,8,9,10,11,12,13,14,15]. The regional climate can be characterized as maritime temperate [16], typical of Mediterranean coastal inlets, where annual air temperatures vary between 9 °C and 17.5 °C and annual rainfall ranges from 400 mm to 1300 mm [8]. The prevailing Aeolian regime is characterized by northerly winds, which are quite frequent during each month of the year, but reveal their highest intensity during the winter (namely “Vardaris” wind; [17]) and summer months (named “Etesians”; [18]).

Figure 1. Bathymetry of (a) the Thermaikos Gulf (TG), divided into three basins (Inner TG, Central TG, and Outer TG), and (b) the Inner TG, divided into two basins (Thessaloniki Bay and Central Gulf of Thessaloniki). The main topographic characteristics (e.g., rivers, capes, Natura 2000 areas) and anthropogenic pressures (e.g., urban, industrial, agriculture) are marked [13]. The insert map in (a) shows the location of the TG in the Mediterranean Sea.

Figure 1. Bathymetry of (a) the Thermaikos Gulf (TG), divided into three basins (Inner TG, Central TG, and Outer TG), and (b) the Inner TG, divided into two basins (Thessaloniki Bay and Central Gulf of Thessaloniki). The main topographic characteristics (e.g., rivers, capes, Natura 2000 areas) and anthropogenic pressures (e.g., urban, industrial, agriculture) are marked [13]. The insert map in (a) shows the location of the TG in the Mediterranean Sea.

Five rivers, four (two large ones, Axios and Aliakmonas, and two smaller ones, Loudias and Gallikos) discharging along the western coast in the Inner and Central TG (Figure 1b) and one (Pinios) discharging in the Outer TG, are the main freshwater and nutrient input to the marine basin of Thermaikos. The natural fluvial inflow is further enhanced by smaller intermittent rivers and streams (e.g., Anthemountas, Figure 1b) and is connected to irrigation canals and a network of drainage trenches (in agricultural areas, Figure 1a). The influx of pollutants through the river system and the evolution of the brackish plume dynamics in the semi-enclosed parts of the Inner and Central TG significantly control the overall TG marine quality state and hydrodynamic processes [14,15]. The two largest rivers, Aliakmonas and Axios, accompanied by Loudias in between, form a complex integrated multi-channel deltaic system, with adjacent protected wetlands¹, i.e., several Natura 2000 sites, a Ramsar Site, a National Park², a coastal lagoon, and many wildlife refuges [19]. North of the integrated deltaic area, Gallikos discharges in Thessaloniki Bay. Its slow flow forms the Gallikos estuary, a vibrant ecosystem of three neighboring ponds, home to common and rare species of migratory birds³. Axios supplies more than 50% of the total freshwater [12] in the TG. Seasonality is not very strong in outflow rates because the daily releases are controlled by the operations of hydroelectric power dams along the two larger rivers’ routes [20]. The daily discharge rates range between 20 m³/s and 300 m³/s, respectively [15,21,22]. Freshwater outflow rates have shown a significant decrease during the past few decades [23], associated either with increased irrigation needs or with reduced precipitation rates [24,25].

The densely populated urban areas along the coastal zone accommodate several infrastructures and human activities, such as farming, treated and untreated sewage discharges, industrial work, touristic operations, aquaculture and intense fishing, seaport operations, and heavy maritime traffic that cause significant environmental pressures related to pollution and contamination of the marine environment (Figure 1a). The activities are mainly concentrated in the coastal zone of the Inner TG, especially in the more polluted Thessaloniki Bay which hosts the main population and industrial infrastructure (Figure 1b). Despite its environmental protection status (National Park; Natura 2000 and Ramsar site; Figure 1a), the western coastal zone is under pressure due to intense agricultural and aquacultural activities. The eastern part, along the coastline from the Central Gulf of Thessaloniki to the Outer TG, is a heavily exploited touristic area. Aquaculture installations, mainly mussel farming units, are located along the western deltas, while intensive fishing takes place over the entire region. Together with the riverine freshwater inputs that add nutrient pollution into the gulf, the anthropogenic stresses may increase the concentrations of organic and inorganic pollutants, heavy metals, hydrocarbons [26], and marine litter [27], downgrading the quality state of this ecologically sensitive marine environment.

3. Collection of and Statistical Analysis of the Studies about the Thermaikos Gulf

3.1. Literature Collection

We have gathered and examined all available literature regarding the marine state of the TG during the last 5 decades. The literature search was performed using all available search engines for scientific documents (e.g., Google Scholar, Scopus, Web of Science). Moreover, a large number of articles and reports about the marine environment of the TG during the last 40 years were produced with the participation of members of the current study’s authors’ team (Aristotle University of Thessaloniki: AUTh; Hellenic Centre of Marine Research: HCMR). Scientists from AUTh (e.g., [28,29,30,31,32,33]) and HCMR (e.g., [34,35,36,37,38]) participated in numerous studies conducted in the TG on different fields of research, using almost all available methodologies (e.g., field and lab work, numerical modeling, satellite imagery) and publishing a large number of articles and official technical reports for several basic and applied research projects since the 1990s. Additional articles about the TG, published in peer-reviewed journals, were also collected based on specific keywords related to all oceanographic disciplines (physical and chemical oceanography, hydrodynamic circulation, biogeochemical oceanography, ecological state, marine pollution, seawater quality, benthic state, sediment transport, marine geology, algae, invertebrates, plankton, etc.) and methodological approaches (field monitoring, in situ measurements, seawater and sediment sampling, chemical analysis, lab analyses, satellite observations, and numerical simulations) (see Appendix A for more details). The keywords represent a large number of specific variables and processes that are considered important for the marine environment and are associated with physical parameters (e.g., temperature, salinity, density) and processes (e.g., circulation patterns, eddies, vertical advection, stratification structures, water mass exchanges, river plume dynamics), geological and sediment characteristics (e.g., sedimentation rates, particulate matter concentrations, shoreline changes), chemical conditions and pollution (e.g., nutrients, heavy metals, organic pollutants, organic matter, gasses), and biological elements in the water column (e.g., chlorophyll-a (chl-a), plankton, eutrophication levels) and in benthic state (e.g., benthic communities and diversity, invertebrate species, ecological quality index).

The total sample includes 215 papers, published either as articles in peer-reviewed journals (94%; Appendix A) or as official technical reports of monitoring surveys (6%). Only a small percentage of articles (2%) were literature reviews of previous studies. The sample covers all types of methodological approaches (e.g., field and laboratory work, numerical simulations, satellite-derived observations) and oceanographic disciplines (e.g., physical, chemical, biological, geological). The potential limitations, regarding the research methodologies and accuracies, can be found in the respective published articles. The possible limitations of the present review are discussed in Section 5.

3.2. Evolution and Main Characteristics of the Research Activities through the Years

The evolution of the annual number of papers published about the TG during the last 5 decades is presented in Figure 2 (black line; derived from the publication dates). The corresponding number of studies conducted within each year (monitoring and simulation periods) is also shown (red line in Figure 2). The two time series show a similar distribution with lower numbers at the beginning and the end of the 50-year period. The highest Pearson correlation coefficient (R_P = 0.56 for Mann–Kendall [39,40]; 99% test of statistical significance: p_value < 0.0001) between the two series was derived after shifting the publications’ time series backward by 4 years, indicating an overall mean 4-year delay between the period of the conducted research and the respective publication date. For example, the highest peak number of conducted research activities was in 2001, while the respective peak of the published work occurred by 2005. Generally, the period with the most intense research efforts was between 1994 and 2008 (hereafter “golden period”; Figure 2) associated with continuous monitoring projects supported by both national and European funds. A strong drop was detected after 2009 and remained low during the 2010–2020 decade when the number of studies in the TG was reduced to almost half (18%; Appendix A) compared to the previous decade (39% in 2000–2009). The 2010–2020 decade showed even fewer research activities than the 1990s (20%), even though the research staff and investigation tools are expected to be more advanced and qualified as the years go by. It is noted that the Greek financial crisis, which started in 2009, had various socioeconomic consequences that affected all aspects of academic and research activity [41], although the environmental pressures were intensified due to the economic recession [42], increasing the need for environmental monitoring and holistic management even more (see Section 5.1). The small number of studies and the sharp reduction during the last five years (after 2018) could also be related to the delay that is usually required for research to be published (~4 years; Figure 2). Still, the lack of oceanographic research in the TG is profound after 2010 and especially between 2009 and 2015, which is probably related to the aforementioned socioeconomic factors and is not justifiable by the expected delays in the publication process.

Most of the research includes field observations (86%), while only 19% concerns numerical applications. Satellite-derived data (5%: Appendix A) were mainly used after 2009 (Figure 3a), likely related to improved data accessibility (e.g., via Copernicus Marine Service) and advances in filtering algorithms and remote sensing capabilities in the coastal ocean. Only 11% of the studies combined more than one methodological approach (multi-platform). During the “golden period” of intense research activity, from the mid-1990s until the end of the 2000s (Figure 2), observational monitoring was the main method for conducting research over the TG, peaking in 2001; there is a reduction after 2010 and a slight increase in the mid-2010s (Figure 3a). After 2010, fieldwork was partially replaced by satellite-derived data, which, however, cannot provide information for the entire water column and for all oceanographic parameters. In terms of scientific discipline, the highest percentages of research cover biological (55%) and chemical (50%) processes, associated with marine pollution and eutrophication events (Appendix A). Investigation of physical processes and hydrography was less than 40%, while geological and sediment processes were the least examined among all disciplines (25%). There were several years in which geological/sedimentary studies were completely absent (Figure 3b). It is noted that more than half of the studies (53%) covered several disciplines of the marine environment simultaneously, especially during the period from 1997 to 2003 (Figure 3b).

The majority of the studies focused on the Inner TG (Appendix A), with 77% conducting measurements or simulations over the area north of Cape Megalo Emvolo due to the increased monitoring needs of the more enclosed and heavily impacted Inner TG. The Central TG is investigated in 63% of the studies, while less attention was given to the southern Outer TG (42%). Only 25% of the research efforts approached the TG as a whole; the annual frequency of such studies is low throughout the 50-year period with relatively high numbers only during the late 1990s (Figure 3c). The interest in the Inner TG was enhanced compared to the Outer TG during the last decade even more (Figure 3c). The distribution of monitoring frequency is sparse in the Outer TG, especially in its central and eastern parts (>22.9° E; Figure 4). A large number of observational locations are concentrated in the most polluted area of the northern Inner TG. Thessaloniki Bay, and especially the urban coastal area of Thessaloniki (north and east coasts), contains the highest number of monitoring stations (20–30; Figure 4). Many field studies have also been conducted at the entrance of Thessaloniki Bay (Cape Mikro Emvolo) and on the east side of the Inner TG entrance (Cape Megalo Emvolo). Most of the studies focus on the Central Gulf of Thessaloniki (63%; Appendix A), where intense aquaculture and agriculture activities take place along its western coastal zone (Figure 1b); many stations have been placed around the Axios River delta (Figure 4). The strong research interest for these two sites, Thessaloniki Bay and Axios, is related to the environmental impact issues along the urban coastal zone of Thessaloniki and the ecological importance of the deltaic area, respectively. The simultaneous coverage of both the Bay and Central Gulf of Thessaloniki is included in half of the studies about the Inner TG, showing quite frequent field surveys from the 1970s until 2008 (Figure 3d). The low number of research efforts after 2009 (Figure 2) is also depicted in the lower spatial coverage of the Thermaikos sub-domains, and especially the implementation of studies that cover the entire study area (Figure 3c), or at least the entire Inner TG (Figure 3d). As far as the vertical coverage is concerned, most of the studies describe the water column state (73%), while only half of them refer to processes on the seabed (Appendix A). Benthic processes and properties (e.g., related to benthic fauna, ecological status, sediment characteristics) have been less investigated in comparison with other properties of the water column. However, as benthic processes are inherently more slowly evolving than the seawater conditions in the marine environment, lower monitoring frequency can be sufficient. The same applies to the spatial extent of seabed field observations (only a few sporadic samples can be representative of the benthic health condition over extended areas).

The temporal characteristics (duration and step) of the research studies are presented in Figure 3e,f. The short-term studies that cover periods of less than one year comprise 48% of the total (Appendix A). Study periods of one year or more (long-term) are less frequent (43%). The short-term studies mainly covered the years from 1995 to 2009, with a second but smaller increase after 2015 (Figure 3e). The temporal step of analysis in most studies exceeds one month, usually corresponding to 3-month (seasonal) monitoring (42% Appendix A). Monthly measurements are significantly less frequent (16%), while a similar percentage of studies refer to unique (conducted once) measurements and experiments (17%). However, most of the time, such one-off experiments are not included in the survey programs (Figure 3f). More frequent monitoring (<month) is also rare (16%), included mainly in studies between 1990 and 2010 (Figure 3f). Besides monitoring activities and realistic numerical simulations, a small number of studies (Figure 3e) refer to past-climatic periods (2%; e.g., decadal or centennial analyses) or process-oriented numerical experiments under idealized conditions (7%). The latter are usually related to numerical simulations that reproduce hypothetical hydrodynamic conditions. Studies based on future climatic scenarios are completely absent.

4. Main Research Findings

4.1. Hydrography and Circulation

The first publication simulating the hydrographic conditions of the TG was authored by Ganoulis and Krestenitis [43], presenting the development of a numerical model to simulate the transport, dilution, and degradation of the pollutant discharges from the sewerage system of Thessaloniki city. However, the earliest monitoring effort, with hydrographic and current measurements that covered the Inner and Central TG, took place in the mid-1970s [44]. This was the first campaign that described all the major physical processes of the gulf: the dominance of riverine waters in the western zone due to northerly winds and the Coriolis effect, the meteorological effect on the vertical structure, the main circulation patterns, and the water mass exchanges with the Aegean Sea. Since then, various oceanographic campaigns and numerical studies have been undertaken (analyzed in the following sections), leading to a deeper understanding of the dynamics and dominant forcing factors in the area.

4.1.1. Main Circulation Patterns

The distribution of physical properties and the general circulation of the TG are controlled by four major environmental factors: (i) the air–sea interactions and especially the wind forcing [11,14,15], (ii) the exchanges with the northern Aegean Sea through the southern boundary [10], (iii) the river plume dynamics [9], and (iv) the morphological setting of the enclosed basin [8,10]. The tidal amplitudes are small in the TG, reaching a few tens of centimeters [8,11,45] and leading to very weak corresponding tidal currents (<0.02 m/s) [10]. The river plume waters affect both the physical characteristics (density currents) and the biochemical properties of the TG coastal waters, due to their high concentrations of nutrients ([21]; see Section 4.3). The formation of the river plumes that control the expansion of the riverine discharge loads is also influenced by the prevailing wind-induced circulation and the proximity of the TG’s east and west coasts [46].

The morphological separation of the TG into three main sub-basins, the Inner, Central, and Outer TG (Figure 1a), and the further separation of the Inner TG into two smaller coastal areas (Figure 1b) control the hydrodynamic circulation conditions. Distinct circulation patterns may form in each sub-basin, controlling the connectivity pathways, the transport across the sub-basins’ entrances, and the renewal of the water masses. Wind forcing may induce different circulation patterns among the gulf’s sub-basins. Winds from the north-northwestern sector (Mistral type, “Vardaris” and “Hortiatis” local ravine winds) are the most dominant (combining strength and frequency) over the region, as derived by a 30-year investigation of the ERA5 Reanalysis database (Figure 5a), in agreement with older studies (e.g., [8]). This wind regime may form a cyclonic ocean circulation pattern in the Central Gulf of Thessaloniki ([15]; Figure 6a). These winds also enhance the southward transport of near-surface waters towards the Outer TG [47], moving the riverine brackish masses away from the confined northern parts of the gulf [10,13] (Figure 6a). This type of wind-induced circulation is commonly evidenced during the winter months, following the general cyclonic circulation along the northern Aegean coasts [47,48]. It also allows the northward cyclonic intrusion of denser Aegean Sea waters (ASWs) in the deeper layers [49] that facilitates the renewal of the Inner TG with clearer waters [13,14]. Such wind conditions can prevail during the entire annual cycle and enhance the renewal of the Inner TG, limiting the eutrophication events [14], a process that is crucial for the quality of the Thermaikos water masses. During these periods, a two-layer flow is formed mainly along the eastern coasts and across the entrances of the Inner TG (Cape Megalo Emvolo; Figure 1b) and Thessaloniki Bay (Cape Mikro Emvolo; Figure 1b). The two-layer flow consists of southward currents over the near-surface layer and northward streaming over the near-bottom layer [15]. A clear southward connectivity pathway, between the more polluted western Thessaloniki Bay ([50]; see Section 4.3) and the Central Gulf of Thessaloniki (Figure 6a), also affects the usually clearer waters along the eastern coasts under northerly winds [15]. The cyclonic general circulation along the Aegean coasts may further evolve along the western coasts and becomes more intense under northerly winds, thus connecting the northern parts of the TG with the southern coasts and the Aegean Sea [9,11]. The hyposaline signal of the Black Sea waters has also been detected in the TG, following the same cyclonic circulation of the northern Aegean Sea, with a peak inflow during the summer months [11,47,49].

Contrastingly, the southeast winds (Sirocco type; “Eurus” in Greek shipping notation), which are usually weaker but also very frequent over the TG (Figure 5a), may reduce the renewal capacity of the gulf [14], spread the freshwater plume over the Inner TG [13] and induce an anticyclonic circulation pattern over the Central Gulf of Thessaloniki [15] (Figure 6b). An anticyclonic pattern, formed in front of the Axios and Aliakmonas deltas over the Central TG (Figure 1a), is associated with the Coriolis force effect on the river plume (forming an “anticyclonic bulge”; [9,46]). The narrowness of the basin near the freshwater discharge sources (Cape Megalo Emvolo; Figure 1b) allows the offshore expansion of the plume towards the eastern coasts [9], especially under the influence of southerly winds [13,15]. Southerly winds may also spread river plume waters towards the Inner TG, even reaching the Thessaloniki Bay [13]. The ASW intrusion is weaker along the eastern coasts, evolving mainly over the Outer and Central TG. The two-layer flow is also weaker and can even become reversed in comparison to circulation under northerlies (e.g., Cape Mikro Emvolo; Figure 6b).

4.1.2. Thermohaline Variability and Stratification

The vertical structure of the water column reveals significant spatial and temporal variability, with periods of stratified or homogenized conditions. The characteristics of the surface layer vary considerably, although the lowest-salinity waters are usually concentrated in the western gulf [10]. Increased riverine outflows hinder vertical mixing, resulting in haline stratification [52], especially from December to April [10]. Strong cold northerly winds enhance vertical mixing and the homogenization of the water column, especially in the shallower Thessaloniki Bay [13], contributing to the formation of a southward-flowing, dense, near-bed current ([10]; analyzed further in Section 4.1.3). The seasonal pycnocline is stronger during summer (thermal stratification) and starts to recede in early autumn. Generally, the formation of the pycnocline starts in spring, and the highest density gradients occur in summer. It is shown that the winter thermocline in the Inner TG develops at greater depths (~15–20 m) than the rest of the domain, forming wider mixed layers, mainly due to the lower depths in the semi-enclosed northern part of the gulf [13]. The vertical homogenization during winter months results in a prevailing barotropic circulation, in contrast with the summer conditions, when the water circulation is mainly baroclinic [11,53]. Strong pycnoclines (combining thermoclines and haloclines) may also form in the near-bottom layer due to the intrusion of colder and more saline ASW; they either remain within the Central TG under southerly winds or enter the Inner TG under the effect of northerly winds [8,22]. It has been established that salinity presents a strong interannual variability, with intense short-term fluctuations. Field measurements after 2020 showed significantly higher salinities (often related to exceptionally high temperatures) at all depths of the Inner TG [15] than values measured during previous monitoring campaigns (e.g., 1994–2007; [13]). Pertaining to this, the sea temperature of the TG has also revealed significantly high peaks, especially during the last decade, and strengthened the formation of marine heatwaves (MHWs), designating the TG as an MHW “hot spot” ([54]; see Section 5.2).

4.1.3. Vertical Water Fluxes

The Inner TG also appears amongst the regions of dense water formation (DWF) in the Eastern Mediterranean Sea, a process taking place under intense cooling by northerly gales [13,53,55]. The phenomenon can be triggered by density gradients due to temperature, salinity, or the contribution of both [56]. Despite the existence of the deltaic system in the southwestern part of the Inner TG, which could inhibit such a process, the northern TG is a favorable area for DWF. Low air temperatures, strong northerlies, and weak river plume spreading in the Inner TG contribute to the formation of denser waters (Figure 6a). The shallow depths of Thessaloniki Bay (10–20 m) and the topography of the area allow rapid cooling of the water masses since the shallow bathymetry favors strong heat losses and fast mixing of the waters, while the morphology reduces the interactions with the ‘open sea’. Thus, under favorable forcing conditions, a semi-secluded mass that is gradually becoming cooler and denser can be formed in the gulf. After a DWF in the Inner TG, a southward density current can be formed towards the southern gulf. The denser water masses that cascade over the TG’s continental shelf are capable of renewing the North Aegean Trough’s deep water and especially the Sporades deep basin [57]. They can reach the southern boundary of the Outer TG approximately 2 weeks after their formation in the Inner TG [53].

Upwelling in the coastal region (close to land) takes place when the wind blows parallel to the coast, especially along the western and eastern coasts of the Central and Outer TG [58]. When winds are favorable to upwelling for the west Thermaikos coast (southerlies), large amounts of high-salinity water enter the surface (Figure 6b), as the buoyancy-driven circulation is reversed [9]. Downwelling may evolve periodically on a diurnal time scale, especially in summer in the northwestern Gulf (Figure 6b) during periods of onshore southeasterly winds [9,59]; note that such winds are very frequent in the TG (Figure 5a). The prevailing and usually stronger northwesterly winds (Figure 5a) may reduce the temperature at the surface after the coastal upwelling of colder waters over the eastern coasts of the Central and Outer TG ([47]; Figure 6a).

4.2. Sediment Transport and Modern Sedimentation

4.2.1. Historic Changes in River Deltas and Sediment Supplies

The deltaic complex of Gallikos, Axios, Loudias, and Aliakmonas is a prominent feature along the northern coast of the TG (Figure 1b). Sediment and freshwater inflows have changed significantly during the last century due to direct and indirect human interventions. The Axios River had two mouths that were periodically active: the presently active mouth at Cape Kavoura (current delta) and the now abandoned Palaiomana mouth, located around 15 km to the NE (22.86° E, 40.59° N; Figure 1b). The active mouth was at Cape Kavoura prior to 1900, after which the main flow shifted to Palaiomana, which remained as the active mouth until the artificial diversion of the flow back to the main tributary (outflowing to Cape Kavoura) in 1934, with discharges from the abandoned Palaiomana mouth thereafter only during floods [60]. This diversion was part of a major reclamation project in the deltaic plain (1925–1933), during which the routes of Axios, Aliakmonas, Gallikos, and Loudias were rearranged, and artificial banks were constructed along their straightened (after the works) distributary channels [61].

Historical data indicate significant accretion in the deltaic zone of Axios and Aliakmonas during the 20th century, with net gains of 175 km² and 139 km² between 1850 and 1987 [61]. Kapsimalis et al. [60] identified three main phases in the deltaic system evolution since 1850: (a) natural evolution between 1850 and 1916, with the system delivering a net sediment supply of 130 × 10⁶ m³ to the TG; (b) accretive phase from 1916 to 1956, during which the artificial realignment of the main river channels and drainage and reclamation projects in the coastal plain increased the sediment outflow to 900 × 10⁶ m³ (28 × 10⁶ t/yr) and led to coastal progradation at the active mouths of Axios and Aliakmonas; and (c) erosional phase after 1956 (1956–2000), due to human-induced (dam construction) reduction in sediment supplies (by 2.5 × 10⁶ m³/yr or 4 × 10⁶ t/yr), leading to coastline retreat in the deltas. Coastal erosion was tackled locally with hard coastal protection measures (1960–1974), like the construction of a 22 km long and 2.5 m high seawall along the abandoned Palaiomana mouth of Axios, equipped with five pumping stations to drain the agricultural plain [60]. The limitation of freshwater and sediment discharge has led to erosion of the protected deltas and shoreline retreat in more recent years. Petropoulos et al. [62] estimated losses of 1.26 km² in Axios and 0.63 km² in Aliakmonas between 1984 and 2009, while the observed subsidence rate in the deltaic plane since the 1960s is around 10 cm/yr [60,61]. Considering that these areas are low-lying, the risk of intensification of the coastal retreat and flooding under rising sea levels in the future is very high and could entail severe environmental and socioeconomic impacts [15]. Regarding the Pinios River, numerous paleochannels can be identified across the flood plain of the delta, with two major former river mouths north of the presently active (since 1955) mouth, and two more south from it [63]. Recent coastline changes shift from progradation around the presently active mouth to erosion over the coastal stretch of around 10 km (around 50% of the deltaic coastline) along the abandoned river mouths, with rates that range from +7.5 m/yr to −3.5 m/yr between 1955 and 2013 [63].

The present-day morphometric characteristics of the deltas indicate the dominance of fluvial processes in the north (Axios and Aliakmonas) and wave dominance in the south (Pinios) [64]. Axios, Aliakmonas, and Pinios provide the bulk of freshwater and fine sediment discharge in the gulf, while the Loudias River has a sporadic contribution [65,66]. While historic data estimate the average total riverine sediment influx to Thermaikos at 3–5 × 10⁶ t/yr [8,67], by the end of the 20th century, the total supply of suspended solids by the rivers had decreased to around 0.63 × 10⁶ t/yr (10.7% from Axios, 1.1% from Aliakmonas, 1.0% from Loudias, and 87.2% from Pinios; [65]).

4.2.2. Seasonal Patterns of Suspended Sediments

Riverine sediments entering the TG marine environment comprise very fine silt to clay, with mean diameters between 2 and 4 μm [46]. Due to the low settling rates, riverine sediments show high suspension times in the domain, often confined above the pycnocline and near the surface due to the stratification of the water column [65,66]. Elevated suspended particulate matter concentrations (PMCs) near the bed are attributed to associated current activity, either inhibiting settling or causing resuspension [8].

The distribution of suspended PMCs over the TG was investigated by various oceanographic campaigns. Among these, campaigns that covered the entire Central and Outer TG, extending southwards until the Sporades basin, were conducted in the framework of the European projects EURECOMARGE (one survey: June 1987; [68]), Metro-Med (four surveys: May 1997, July 1997, February 1998, and September 1998; [65]) and INTERPOL (three surveys: September 2001, October 2001, and February 2002; [53]). More recent monthly monitoring campaigns [35,36,69] were restricted to the Inner and Central TG. Satellite imagery has also been used to sense suspended PMCs in the surface waters, using Landsat [70,71] and MERIS [72,73]. The main conclusions of these studies regarding seasonal patterns of suspended particulate matter in the area are discussed in unison below.

Elevated PMCs are often concentrated near the surface and the seabed, forming surface and bottom nepheloid layers (SNLs and BNLs, respectively). The seasonal patterns of SNL formation in the gulf are mainly controlled by suspended sediment influx, stratification, and surface circulation, while the BNL structure is due to the resuspension of loosely consolidated sediments by near-bottom currents [65]. PMCs are typically lower in the rest of the column, with less pronounced, or even absent, intermediate nepheloid layers (INLs). Overall, under non-stratified conditions, INLs follow the distribution and magnitude of the SNL, while under stratified conditions, PMCs of INLs are significantly reduced. Winter conditions, with a high influx of solids from the rivers and strong haloclines, are typically characterized by well-defined SNLs (PMCs > 0.5 mg/L) that extend southwards from the Central TG and the deltaic zone in the north along the western coast of the Outer TG (Figure 7a), driven by the predominant cyclonic circulation (Figure 6a). INLs, when present, show lower PMC maxima and are confined closer to the north and west coasts. BNLs during this period (Figure 7b) can be extensive, with high suspended PMCs (>0.5 mg/L) covering most of the Outer TG, apart from its eastern part along the coast of the Chalkidiki peninsula. Higher PMCs (>1 mg/L) within the BNL can be found throughout the west coastal stretch, and locally elevated concentrations can also be present near the river deltas. Extraordinary dynamics and transport conditions near the bed, like density currents, can influence the BNL structure. Such conditions might be observed during winter when dense water masses can be formed in the Inner TG and slowly travel toward the Aegean Sea (Figure 6a; see Section 4.1.3), altering the BNL structure in their wake [52,55]. This has been observed in the field, with a peak in PMC of the winter BNL offshore the Pinios mouth (measured in February 2002) linked to a DWF event that had taken place 2 weeks prior [53]. Spring conditions can be associated with a very prominent SNL, with elevated maximum PMCs near all river mouths (3–7 mg/L), as well as throughout the Central TG and along the western coasts and the Inner TG (>1 mg/L; dashed lines in Figure 7a).

A reversal of the surface circulation near the southwestern part of the Outer TG can observed due to the occasional shift of the Sporades eddy from predominantly cyclonic in winter to anticyclonic in other seasons [11,47], which can lead to northward expansion of the Pinios SNL plume. High PMCs (>0.5 mg/L) in the spring BNL can cover the entire domain, with even higher concentrations (>1 mg/L) along the western half of the Outer TG (dashed lines in Figure 7b). Suspended PMCs are reduced during summer due to reduced influx from the rivers; under southerly winds, the sediment can be confined to the north and transported toward the Inner TG, with high PMCs (>1 mg/L) in the SNL northwest from Cape Megalo Emvolo (Figure 7c) and with very low values in the rest of the domain. The BNL can be more pronounced, with high PMCs (>1 mg/L) over the western part of the shelf and near the river mouths, as well as the Inner TG (Figure 7d). Autumn conditions (dashed lines in Figure 7c) are mostly characterized by higher outflows than summer and cyclonic circulation, leading to a more typical SNL formation in the domain, with higher PMCs (0.5–1 mg/L) along the river outflows in the north and with the sediment-laden plumes extending southwards along the west coast. Autumn surveys show that the BNL can often extend eastwards from the river mouths, with high PMCs (>1 mg/L) throughout the entire Central TG and extending along the east coast well south of Cape Epanomi (dashed lines in Figure 7d).

4.2.3. Near-Bed Processes and Sedimentation Patterns

The TG acts as a sediment trap, with the bulk of the riverine sediment (~90%) depositing within the bounds of the domain [66,67]. Sedimentation in the deltaic zone of the TG has been associated with seaward dispersion of the river outflow and settling of the suspended material, dominated by differential settling near the river mouth and by flocculation further offshore [74]. Sediment cores indicate that fine-grained sediments are frequently resuspended on the Thermaikos shelf, by natural (storm) or anthropogenic (intense trawling activities) forcing, while there is no evidence of erosion of the bed [21,53]. Numerical modeling and oceanographic observations indicate that nearshore currents are not able to resuspend sediment [1], while wave-induced resuspension can be an important factor over low depths near the coast [75]. The trawling activity taking place in the area has been shown to contribute to the BNL concentrations in the Outer TG [52,53]. In fact, trawling was estimated to increase BNL concentrations over the shelf by 2 to 3 times [2]. Similarly, using numerical modeling, Kombiadou and Krestenitis [76] estimated that the sediment entering suspension due to trawling over the trawling season can supersede the annual influx from the rivers. Resuspended particles are removed rapidly from the system by sinking and/or horizontal advective processes [77]. Overall, mechanically eroded sediment moves in the vicinity of the bed (3–6 m), contributing to the BNLs in the area, and resettles onto the bed shortly after (1–5 days; [76]).

Surface sediments over the Inner TG and the western half of the shelf (Central and Outer TG) comprise muds and sandy muds that transition to muddy sands in the eastern part of the shelf [8,67,78]. Sedimentation over the western shelf is dominated by the riverine water and sediment outflows, while the sand fraction is high (>50%) over the eastern shelf, with values that reach 60–85% along a relatively narrow band along the eastern shoreline of the Outer TG [8]. These coarser sediments are interpreted as relict sands that remain uncovered (or partly mixed) by modern sedimentation, due to the predominant cyclonic transport [22]. Sedimentation rates near the river mouths were measured at 8.8 mm/yr for Axios and 1.8 mm/yr for Aliakmonas [21]. Corresponding modeled values are 18, 7 and 3 mm/yr near the mouths of Aliakmonas, Axios and Pinios, respectively [66]. The same study indicated the following sedimentation patterns: Axios is the main sediment supplier for the Inner TG; Axios and Aliakmonas contribute equally to the Central TG, with the sediments from the former showing higher dispersal in the domain and sediment from the latter depositing mainly southwards along the western coastline; matter from Pinios settles along the west coast of the Outer TG, with a higher sediment accumulation north from the delta.

4.3. Marine Pollution and Biogeochemical Properties

The TG is influenced by human activities (sewage and industrial discharges, agricultural runoff, shipping, airborne litter, hydrocarbon refueling, and harbor operations) while the five large rivers discharge substantial loads of nutrients, heavy metals, and other organic compounds, after drainage resulting from human activities within their watersheds. The Axios River has been recognized as the major supplier of suspended particulate matter even during the dry summer season ([69]; see Section 4.2) and a source of heavy metals to the coastal environment of the gulf [12]. The influence of the rivers, the Thessaloniki metropolitan area, and the diverse sedimentological background contribute to the formation of a spatially varying and patchy distribution of major and minor elements in the surface sediments of the TG. In Thessaloniki Bay and the Central Gulf of Thessaloniki, the geochemical characteristics of the sediments are substantially influenced by the domestic and industrial effluents that were released partly treated for several decades resulting in high levels of organic carbon and contaminant elements (Cu, Zn, As, and Pb).

4.3.1. Heavy Metals

Metal enrichment is detected mainly in the northernmost part of the TG (Inner TG) and to a lesser extent close to the river mouths; from there, the heavy metal content gradually decreases towards the open sea [79]. Several other studies have also reported increased concentrations of heavy metals in sediments adjacent to the rivers, the industrial zone, the port of Thessaloniki, and the sewage outfall of the wastewater treatment plant (WTP) (Figure 1a; e.g., [49,80,81,82,83,84,85,86,87]).

The sediments in the western part of the TG’s continental shelf are dominated by major and minor elements of terrigenous origin that are transported by the rivers, reflecting the lithology of their catchment areas where abundant metal-rich formations, mainly ophiolites, are present. In addition, the Axios River carries suspended matter enriched in metals, due to mining and industrial activities in the upper catchment, while the Aliakmonas and the Pinios rivers have no significant industrial activity in their drainage basins. The geochemical signature of each river is identified by the supply of certain heavy metals: (a) Zn and Pb for the Axios River; (b) Cr, Co, Ni, Cu, and As for the Aliakmonas River; and (c) V, Co, Ni and Cu for the Pinios River. In the central and eastern parts of the shelf, the presence of relict sands (see Section 4.2.3) controls the sediments’ geochemistry, causing a substantial dilution of heavy metal concentrations [79].

The pollution status of marine sediments in the coastal area of Thessaloniki Bay (Figure 1b) was evaluated by employing enrichment factors and comparing observations with sediment quality guidelines [49,85,87]. For Zn, Cu, Pb, As, and Ag, a significant contribution from human-related sources was found, particularly in the inner part of the Bay [85]. The sediments from the area near the Thessaloniki port exhibited the highest impact of anthropogenic heavy metals, revealing the cumulative effect of anthropogenic load from the industrial and shipping activities operating for many decades in the area [87]. Most heavy metal concentrations were found at levels that could adversely affect marine biota, while sites at the inner part of the Bay could be considered of medium-low and medium-high priority in terms of toxicity potential [85].

The dissolved fraction of heavy metals in the TG’s water column corresponds to 33% up to 90% of the total metal load. The ratio of dissolved trace metal concentrations in the Inner TG, based on data collected from 1995 to 2003, showed three-fold enrichment in Cd, Cu, and Ni compared to the North Aegean; for Pb, the enrichment was seven times higher, whereas for Mn, the ratio was less than 1 [88]. Zeri and Hatzianestis [89] measured the dissolved Cu and Ni and reported that the maximum concentrations of both metals were encountered in front of the Axios and Aliakmonas river deltas. Stronger binding of both Cu and Ni with organic ligands of terrestrial origin was also observed in comparison with those of marine origin, underlining the important role of riverine dissolved organic matter in metal complexation at least during winter. Data on dissolved trace metals (Cd, Co, Cu, Hg, Ni, Pb, Zn), collected from 2012 to 2019 in the framework of the national monitoring for the implementation of the European Union Water Framework Directive⁴ (WFD), have shown no systematic metal enrichment. However, the adjacent estuaries and lagoons of the TG were classified as above-background-concentration water bodies for several metals, indicating that metal contamination is localized within these areas and does not extend extensively across the broader expanse of the gulf [90].

The mineralogical composition of particulate matter has been found to be fairly homogeneous over the Inner TG, consisting of clay minerals, quartz, plagioclase, K-feldspar, and calcite. The geochemical analyses of the PMCs showed the predominance of Al, Si, Fe, Ti, K, V, Mg, and Ba, indicating their terrigenous origin, i.e., aluminosilicate detrital material. The increased Al concentration at the surface suggested that river inputs constitute the main source of Al in the study area. Si and Ca are of terrigenous origin but also originate from autochthonous biogenic fraction [69]. Fe is predominantly held in the structure of aluminosilicate minerals, while Mn is mainly controlled by the resuspension of the surficial sediments, by redox recycling, as a result of organic matter degradation and biogeochemical processes [2]. The increased concentrations of particulate Co, Ni and Cr are attributed to riverine inputs produced from the weathering of the mafic and ultramafic rocks in the drainage basin of Axios and Aliakmonas watersheds. Enhanced Cu and Zn concentrations in Thessaloniki Bay and the southwestern part of the Central Gulf of Thessaloniki (Figure 1b) imply that a portion of the two trace metals is derived from partially treated domestic and industrial effluents [69]. Similarly, Pb and Cd were attributed to anthropogenic sources, derived from the rivers and the city of Thessaloniki [2].

4.3.2. Organic Pollutants

Pesticides, polychlorinated biphenyls (PCBs), and hydrocarbons are among the organic pollutants commonly found in the marine environment, entering the sea both indirectly and directly through a range of diffuse and point sources including sewage and industrial discharges, agricultural and urban runoff, ship transportation, and dry deposition from the atmosphere.

Agricultural activities in the river catchment areas seem to be the source of various pesticides like the herbicide triazine atrazine (ATR) and its transformation products (TPs) and organochlorine insecticides like DDT and its metabolites, lindane and dieldrin. The major organochlorine pesticides identified in surface sediments collected in 1995 at Inner TG were p,p’-DDT and its metabolites p,p’-DDE and p,p’-DDD, while the rest of the organochlorine compounds, namely lindane and dieldrin, were found at quite low concentration levels [91,92]. These results confirm the findings of previous relevant studies [93,94]. The maximum value of ΣDDTs (sum of concentrations of DDT, DDE, and DDD) was found inside the Thessaloniki port, while high levels of ΣDDTs were also measured close to the WTP sewage outfall and near the city of Thessaloniki. Moreover, the ratio DDE/DDT was consistently higher than 1, suggesting the absence of “fresh” DDT inputs in the gulf. The vertical distribution of ΣDDT concentrations in a sediment core showed the highest value 30 cm down-core, probably attributed to the banning of DDT use in Greece since 1972 and corroborating the absence of new DDT inputs in the Inner TG [91,92]. Despite the ATR ban within the European Union in 2004, ATR and its TPs were still detected in seawater and sediments of the TG collected during sampling campaigns in 2010 and 2011 [95], although at relatively lower concentrations than in the period before the ban of ATR use [96,97]. However, the detection frequency of the compound was 100%, probably due to leaching from agricultural soils where significant amounts were applied in the past and its slow degradation in the marine ecosystem. In parallel, terbuthylazine, which is an active substance replacing ATR, was detected at rather low concentrations except for one exceptionally high concentration that was observed in the western part of Thessaloniki Bay close to the estuary of the Gallikos River [95].

The occurrence of phenolic compounds (nitro- and chlorophenols) in the seawater was studied along the western coastline of the Outer TG for one year (October 2003–September 2004). Considering the dominant activities in the wider area, urban, industrial, and agricultural sources, as well as natural loading, are the main contributors to the presence of these compounds. The detected concentrations were higher than the levels set by the European Union for bathing waters, particularly for pentachlorophenol and 2,4-dichlorophenol, but do not appear to pose a threat to the aquatic environment of the area [98]. Arditsoglou and Voutsa [99] also conducted a study to determine the presence of phenolic and steroid endocrine-disrupting compounds (EDCs) in all elements of the marine environment of the Inner and Central TG (seawater, suspended particulate matter, sediments, and mussels), based on findings of previous studies that showed the presence of EDCs in inland waters (rivers and streams), sewage effluents, and industrial wastewaters discharged into the gulf [100,101]. Phenolic compounds were detected in all samples, with nonylphenol and its ethoxylates being the dominant pollutants. Nonylphenol occurred in sediments in relatively high concentrations, posing a significant risk to biota, while mussels exhibited low bioconcentration factors for nonylphenol and octylphenol.

Organophosphate esters (OPEs), a group of chemical compounds used as additives into plastic materials/products, were determined in samples originating from the coastal area of the Inner TG, as well as from rivers and streams outflowing into the gulf during 2019–2020. OPEs were ubiquitous in the aquatic environment, with higher concentration levels in streams. Alkyl OPEs were the dominant group of the dissolved fraction, contributing more than 50% in all samples from rivers and seawater, followed by aryl OPEs and chlorinated OPEs. Risk assessment did not show a potential threat due to OPEs to aquatic species, although certain streams exhibited relatively higher risk [102].

Hatzianestis et al. [91,92], using the same sediment samples analyzed for DDT and its metabolites, conducted a parallel study to determine seven PCB congeners. The sum of the concentrations of the seven PCB congeners (ΣPCBs) presented a spatial distribution similar to that found for DDTs with maximum ΣPCB values in the vicinity of the city of Thessaloniki. The congener distributions were generally dominated by hexachloro- and heptachloro-substituted compounds, as in the commercial products Arochlor and Clophen. Interestingly, the less volatile heptachloro-substituted compounds were present at very low levels or even absent in the more distant stations, indicating that they are more easily removed from the atmosphere and are partitioned close to the pollution sources.

The same sediment samples were also analyzed to determine hydrocarbons [91,92]. In surface sediments, the aliphatic hydrocarbon composition was dominated by petroleum-related constituents, and according to the observed concentration levels, the Inner TG can be characterized as moderately polluted with petroleum compounds. The evaluation of the gas chromatograms of the aliphatic hydrocarbons indicated the widespread presence of important petroleum-related residues in sediments with a distinctive decreasing trend offshore [91]. Polycyclic aromatic hydrocarbons (PAHs) were also determined in the upper 60 cm of one core. PAH concentrations in the sediment core were relatively low; however, the proportion of petrogenic PAHs (compounds found mainly in petroleum) over the other fractions was increasing in the deeper parts of the core, implying a long-lasting petroleum contamination history in the area [92].

The concentrations of different organic contaminants (i.e., PAHs, PCBs, organochlorine pesticides, and polybrominated diphenylethers) were measured at three sites of the Inner and Central TG not only in surface seawater but also in the air in summer (July) 2012 [103]. The results showed that the concentrations of PAHs, PCBs, and some organochlorine pesticides (hexachlorocyclohexane-HCHs and DDT and its degradation products) in surface seawater were higher than those in the Cretan Sea (Southern Aegean) by a factor of 2–10. Furthermore, the spatial variation of seawater pollution across the three sites was remarkable, with the highest PAH levels recorded at a residential coastal site (Michaniona); the highest levels of organochlorine pesticides measured at the Loudias River estuary, implying influence from agricultural activities; and the highest levels of PCBs and polybrominated diphenylethers recorded offshore, possibly indicating a shift of contaminant partitioning from particle-bound to dissolved, corresponding to a negative gradient of suspended particulate matter concentration, or influence from shipping activities.

4.3.3. Nutrient Pressures

The TG is one of the most disturbed Hellenic coastal areas regarding anthropogenic nutrient enrichment [104]. Environmental pressures of human origin, from urban, agricultural, industrial, commercial, marine, and aquaculture activities, have led to increased nutrient concentrations in the water column of the TG that exhibit significant seasonal and annual variations (e.g., [105,106,107,108]). More specifically, the Inner and Central TG receive domestic, stormwater, agricultural, and industrial effluents not only through the four rivers and their tributaries, but also from the treated wastewater discharges and the sewerage network overflows of the Thessaloniki metropolitan area. Large volumes of domestic and industrial wastewater have been directed for decades to the TG and especially to the northernmost part, Thessaloniki Bay, promoting the occurrence of eutrophication conditions (e.g., [109,110]). From 1992 on, Level I wastewater treatment was applied, and since 2001, Level II wastewater treatment has been implemented [13], with a 92% performance degree, in order to reduce the contribution of urban effluents to eutrophication events (e.g., [111]; see Section 4.4).

In general, the nutrient distribution shows significant spatial variability, with high concentrations of nutrients (e.g., nitrates, phosphates, nitrites, ammonia) near the west coastline of the Inner TG, in the proximity of the river estuaries and aquaculture installations, as well as in Thessaloniki Bay (port and sewage outfall area) [104,112]. A time series (1995–2007) analysis of nutrient data from Thessaloniki Bay and the area near the Axios River delta showed that nutrient concentrations were generally lower in Thessaloniki Bay than those near the river estuary, except for ammonia during the period 1995–2000 [108]. The seasonal variability of nutrient concentrations is controlled by several factors acting synergistically, such as the seasonality of the quantity and quality of both the treated and untreated wastewater effluents, the fluctuation in the flow rate of the main rivers discharging into the study area, and the mesoscale and sub-mesoscale ocean circulation in conjunction with the prevailing meteorological conditions [14]. This latter study has shown that localized-scale ocean dynamics may reduce the renewal of the gulf with clear waters of Aegean origin (see Section 4.2) or induce advection of the nutrient-rich waters from the rivers to the northern part of the gulf, enhancing the eutrophic characteristics of the gulf (see Section 4.4).

Among the rivers, Axios is the most polluted in terms of nutrient concentrations, characterized by a clear seasonal signal of high flow rate in spring and low in late summer, contrary to the Aliakmonas River, where the existence of dams along its waterway determines the temporal variability of the outflow with high to maximum discharge in summer due to the peak of hydropower production [113]. In view of these differences, Aliakmonas shows an increase in nitrate concentration with discharge, indicating the prevalence of agricultural land flashing, while Axios shows a nitrate decrease with discharge, indicating the prevalence of dilution processes [113]. Regarding long-term nutrient variation, Axios presents an increase in nitrate concentration during the 1980s that can be attributed to agricultural intensification. After a subsequent decrease, nitrate concentration increased again to reach the multi-year maximum in 1997–98 and has gradually diminished since then [113,114]. On the other hand, in Axios, phosphate concentration remained particularly high and slightly increased until 2000 due to increased inputs from industrial and urban waste point sources in the Former Yugoslavian Republic of Macedonia (North Macedonia now) which includes almost 80% of the river’s drainage basin [113,114,115]. Accordingly, phosphate loads are unbalanced in relation to nitrate loads, resulting in decreased values of N/P ratios, lower than the theoretical Redfield stoichiometry (N/P = 16:1), implying a possible nitrogen deficiency, in the coastal river-influenced environment.

Analysis of nutrient data time-series (1995–2002) in the Inner TG (including Thessaloniki Bay) has revealed that nitrate, nitrite, and ammonium concentrations showed a statistically non-significant decrease, nitrate showed a statistically non-significant increase, phosphate did not follow any trend, and ammonium showed a statistically significant (p_value < 0.05) decrease. This has resulted in a decrease in the N/P ratio which may cause indirect effects of eutrophication, like the blooming of certain harmful dinoflagellate species [104,107,116]. The above results were confirmed by Zachioti et al. [108], who found that the N/P ratio was constantly below 10, highlighting the strong N limitation in the inner part of the gulf, also demonstrated in other previous studies (e.g., [112,117,118,119]). However, based on weekly nutrient measurements in Thessaloniki Bay during March 2017–February 2018, the mean inorganic nutrient ratios were 25.1 for N/P, 0.70 for Si/N, and 18.4 for Si/P. The average N/P ratio in 2017–2018 was much higher (mean 25) compared to the N/P ratio (6.40) of the period 1995–2007 [107]. The high N/P ratio in combination with the extreme NH₄ concentrations measured during 2017–2018 may be linked to the relatively high contribution of dinoflagellates, such as Noctiluca scintillans and Spatulodinium pseudonoctiluca, to plankton community biomass and increased red tides [6]. Similarly, the N/P ratio was high (19.5) during the years 2020–2022 [120], indicating increased nitrogen loadings.

4.3.4. Other Biogeochemical Characteristics

Particulate organic matter (POM) in the TG is a mixture of terrigenous organic material, autochthonous biomass derived from planktonic and benthic organisms, and inputs from anthropogenic activities. In the period from June 2004 to June 2005, the highest particulate organic carbon (POC) concentrations were recorded at the surface of Thessaloniki Bay and Central Gulf of Thessaloniki during spring–early summer, linked to the development and decline of the spring phytoplankton bloom (see Section 4.4.2). The distribution of particulate nitrogen and phosphorus was similar to that of POC [69]. Furthermore, POC and the nitrogen content of the particles were significantly correlated, and the C/N ratio of the POM ranged between 5 and 12, with high values systematically observed over the river deltas and in Thessaloniki Bay and the Central Gulf of Thessaloniki. C/N ratios > 7 that have been calculated at different depth layers suggest that carbonaceous compounds were selectively accumulated in bulk POM and resulted in carbon enrichment. Nevertheless, the C/N ratio was close to the theoretical stoichiometric analogy of marine plankton (C/N = 106:16 = 6.6), implying that autochthonous biomass is the main component of particulate organic matter. In another study conducted from April 2000 to May 2001 [121], low ratios of POC/chl-a were recorded in the mussel farm area close to the river estuaries; the high participation of chl-a in POC was also supported by the strong positive correlation between total chl-a and POC in the Inner TG’s stations. In addition, a short-term deployment of sediment traps in the spring of 2003 showed that the proportion of phytoplankton carbon in the POC vertical flux was up to 45% [122].

The Inner TG, characterized by high loads of organic and mineral particles, as well as inorganic nitrogen compounds, provided ideal conditions for either the production of nitrous oxide (N₂O) by nitrification and denitrification or the production of methane (CH₄) by methanogenesis that consequently may be released to the atmosphere ([123] and references therein). Thus, the concentrations of the two greenhouse gases and the bacterial processes involved in their production (nitrification and denitrification for N₂O and methanogenesis for CH₄) were determined in surface waters in April 1998. High concentrations of both gases were recorded in the surface waters of the Inner TG. However, no direct relationship was found between the concentration and the production of the biogases, as they may also be produced in deep water or bottom sediment in shallow areas or derived from anthropogenic activity or ship contamination in polluted areas. Whatever their origin, the presence of extremely high concentrations of these two greenhouse gases in surface seawater suggests that the Inner TG could act as a significant source of atmospheric N₂O and CH₄ [123].

In a study performed in May 1997, data of the CO₂–carbonate system of the Inner TG were obtained for the first time [124]. Surface concentrations of total dissolved inorganic carbon (DIC) were lower than those recorded close to the bottom. The positive relatively good correlation between DIC and both the apparent oxygen utilization (AOU) and phosphate at the last sampling depth confirmed the regenerative origin of a large proportion of DIC. The decomposition of the organic material, imported directly from external sources (rivers, sewage) and produced in situ in the gulf, increased the concentration of dissolved inorganic carbon, particularly in the bottom layer. The correlations between surface fugacity of CO₂ (fCO₂), calculated from the total alkalinity (AT) and DIC measurements, and chl-a as well as AOU revealed that the carbon dioxide fixation through biological activity was the principal factor modulating the variability of fCO₂. A rough first estimate of the magnitude and direction of the air–sea CO₂ exchange showed that in May 1997, the Inner TG acted as a weak sink for atmospheric CO₂ at a rate between −0.60 and −1.43 mmolm⁻²d⁻¹, depending on which formula for the gas transfer velocity was used, and in accordance with recent reports regarding other temperate continental shelves.

4.4. Biological Elements of the Water Column—Plankton

Plankton—named from the Greek “πλαγκτόν”, meaning to “wander”—are a vital component of the marine environment. Bacterioplankton, comprising the highest abundance of pelagic biota, have key roles in global processes, influencing ecosystem structure and functioning as well as biogeochemical cycles. Marine phytoplankton contribute approximately 50% of global primary production and play a crucial role in removing CO₂ from the elevated atmospheric levels via carbon pumping to the deep sea. Protozooplankton are highly important for carbon and energy transfer in the microbial loop and plankton food web. Metazoan zooplankton are an important component of the pelagic food web, contributing to many ecosystem functions, such as the transfer of primary production to higher organisms [125].

4.4.1. Bacterioplankton

Although the view on the ecology of marine plankton has changed remarkably since the 1980s [126,127], recognizing the significant role of bacteria in plankton food webs (microbial loop), there is a lack of publications on bacterioplankton structure and dynamics in the TG. The first relevant work in the Inner TG [128] indicated high numbers of bacteria (from 0.9 to 6.8 × 10⁶ cells mL⁻¹), reaching maxima in May, within the water layer between 50% and 20% of the incident photosynthetic active radiation (PAR). This bacterial maximum coincided with high primary productivity (13.0 mg C m⁻³h⁻¹) indicating strong trophic coupling of bacterioplankton and phytoplankton. Another work, recently published by Genitsaris et al. [129], aimed at the response of bacterioplankton communities from the Inner and the Outer TG to PAH exposure under experimental conditions (mesocosms). No negative effects of the contaminants were apparent in the experiments, and the net growth rates of both bacterioplankton communities were always positive. The bacterioplankton community composition was also examined for the first time. High operational taxonomic unit (OTU) numbers were observed in both communities. In total, 841 and 693 OTUs were retrieved from the Inner and Outer TG, respectively. A high functional-response diversity of PAHs degrading OTUs reflected a high potential for acute mitigation responses of bacteria in the TG. The Outer TG contained bacterioplankton communities with higher resistance sustaining bottom-up ecosystem stability.

4.4.2. Phytoplankton

Several studies on phytoplankton in the TG have been carried out, although long time series are lacking. The existing studies provide important information on phytoplankton composition, abundance, and biomass; seasonal and spatial distributions; molecular diversity; and harmful algal blooms (HABs). The earlier studies on phytoplankton have focused on composition, abundance, and biomass. Diatoms and dinoflagellates were the richest groups in morphospecies of the diverse phytoplankton community of the Inner TG [130]. The phytoplankton seasonal succession was dominated by Leptocylindrus, Chaetoceros, Skeletonema, Pseudo-nitzschia, Thalassiosira, Cylindrotheca, and Cerataulina (diatoms) and Gymnodinium, Prorocentrum, Scrippsiella, Gonyaulax, and Heterocapsa (dinoflagellates) species forming phytoplankton biomass peaks up to 6.7 mgL⁻¹. Focusing on spatial phytoplankton distribution in the entire TG, Gotsis-Skretas and Friligos [131] showed a spatial variation of phytoplankton along the pollution gradient with a maximal total abundance of 7600 cells mL⁻¹. Most of the dominant phytoplankters of the Inner TG, mentioned above, also dominate the phytoplankton communities of the broader TG. New findings from works spanning the entire Aegean Sea (including the TG) indicate that circulation patterns leading to broad-scale compartmentalization might also have an important impact on the diversity and composition of phytoplankton assemblages in the coastal system of Thermaikos [132]. Androulidakis et al. [14] indicated that the phytoplankton species of the different habitat communities of the Inner, Central, and Outer TG were highly connected via circulation (see Section 4.1), while the largest species pool, including nutrient opportunists, was recorded in the phytoplanktonic community of a station located closest to the river deltas (see Section 4.3).

Genitsaris et al. [133], based on an eDNA analysis, showed that the temporal scale rather than the small spatial scale, was the main signal for shaping phytoplankton assemblages in the Inner TG. This metabarcoding analysis has detected phytoplankton groups frequently observed in the Inner TG, dominated by diatoms and dinoflagellates. The most abundant molecular species were Chaetoceros tenuissimus, Chaetoceros cf. wighamii, Skeletonema pseudocostatum, Thalassiosira sp., (diatoms), Scrippsiella trochoidea, Gonyaulax fragilis, Alexandrium margaelefii (dinoflagellates), the cryptophyte Teleaulax sp., the pico-chlorophyte Micromonas pusilla, and the haptophyte Haptolina sp. Dominant OTUs were closely related to species known to form harmful blooms and mucilage aggregates, together with rare taxa having potential negative impacts on human health not detectable with traditional microscopy. Among the potentially harmful species of the TG’s phytoplankton [6], Dinophysis species was found responsible for mussel intoxication in aquaculture areas [134,135]. Dinophysis bloom occurrence had been documented several decades ago [130,131,136,137]. Cysts of potentially toxic dinoflagellate species, such as Alexandrium cf. tamarense, A. cf. affine, A. cf. minutum, and Gymnodinium catenatum, were detected in the sediment of the Central and Outer TG [138].

Comprehensive long-term studies on red tides and algal blooms are scarce. The first study that examined the seasonal shift and described the responsible phytoplankton for the formation of these phenomena in the Inner, Central, and Outer TG was conducted by Petala et al. [139]. By applying traditional microscopy, they divided Thermaikos into two subregions based on phytoplankton abundance and, hence, the eutrophic status. The study estimated significantly higher abundance in the Inner and Central TG compared to the Outer TG. Algal bloom and red tide formation (Figure 8) by several mucilaginous diatoms (i.e., Chaetoceros sp., Cylindrotheca closterium, etc.) and dinoflagellates (i.e., Gonyaulax fragilis) including the heterotrophic dinoflagellate Noctiluca scintillans (Figure 9) and its relative Spatulodinium pseudonoctiluca were observed. Furthermore, the study found considerably lower plankton abundance in the Outer TG, demonstrating its oligotrophic state. A thorough study documenting the weekly phytoplankton species succession and identifying the responsible taxa for the formation of phytoplankton blooms and red tides in the urban marine environment of the Inner TG was conducted by Genitsaris et al. [6]. The study classified this subregion of Thermaikos as heavily eutrophic, indicating persistent phytoplankton blooms that were dominated by known mucilaginous species such as C. closterium, Chaetoceros spp., Leptocylindrus spp., and Skeletonema costatum (Figure 9). The accumulation of the organic mucilaginous material from the successive algal blooms and red tides resulted in the formation of an excessive mucilage aggregate phenomenon in June and July 2017 (Figure 8), extending from the Inner to the Outer TG [6], attributed to the same mucilage-producing diatoms, together with the mucilaginous haptophyte Phaeocystis sp. and the dinoflagellate G. fragilis.

The trophic conditions of the TG were assessed using the multiparametric Eutrophication Index (E.I.; [140]) by some studies. Based on nutrients and chl-a data of the period 1995–2007, an improvement in the trophic status of Thessaloniki Bay and of the area near the Axios River estuary was observed. The trophic status of both areas, which was classified as “Poor” until 2000, was improved and classified as “Moderate” in 2007 [108]. In addition, data collected in five samplings during 2012–2014 for the implementation of the WFD showed that the “Moderate” trophic status of the two areas is maintained [118].

Several studies in the prism of monitoring activities in the framework of the implementation of WFD 2000/60/EC and Marine Strategy Framework Directive⁵ (MSFD) 2008/56/EU evaluated the Inner TG, and particularly Thessaloniki Bay, with “Moderate” or “Poor” status [107,118,119]. Recently, a study on eutrophication assessment of the TG capitalized on current advances in remote sensing and employed satellite image analysis in order to evaluate the spatial distribution of chl-a levels [14]. In order to have a more accurate evaluation of chl-a levels, the authors combined microscopy with satellite ocean color image analysis and attempted to associate the chl-a with the “actual” plankton biomass during eutrophication events in the TG. It was found that chl-a concentrations evaluated using satellite imagery correlated well with phytoplankton biomass measured with time-consuming microscopy, revealing similar patterns. Thus, they concluded that by coupling the two methods, it is possible to determine the extent, as well as the responsible taxa, of the formation of these eutrophication phenomena. For instance, they found that during an annual cycle, chl-a peaked at 45 mgm⁻³ in December 2017 (Figure 8) under the prevalence of southerly winds (weak renewal; Section 4.1), when an extensive red tide caused by the photosynthetic ciliate Mesodinium rubrum expanded throughout the Inner and Central TG [14]. On the contrary, low chl-a levels were detected during periods with northerly winds that enhanced the renewal of the enclosed basins.

4.4.3. Protozooplankton

The protozooplankton community is underappreciated in plankton studies in the TG. The early study by Mihalatou and Moustaka-Gouni [128] focused on protozooplankton components (heterotrophic nanoflagellates, heterotrophic dinoflagellates, and ciliates) as grazers of bacteria. The ciliate M. rubrum as a member of protozooplankton was assigned to the photosynthetic functional category, supported by organelle robbery. Several heterotrophic dinoflagellates (e.g., species of the genera Protoperidinium, Gyrodinium) and even the red-tide-forming omnivorous feeder N. scintillans have been considered members of the phytoplankton community in earlier plankton studies [130,131]. Investigating the underexplored diversity of marine unicellular eukaryotes, Genitsaris et al. [133] identified phytoplankton heterotrophic unicellular eukaryotes frequently observed in the TG. In addition, they also revealed protozooplankton taxonomic groups previously undetected in the area (MALVs, MAST, and Cercozoa). The protozooplankton community showed temporal patterns rather than small-scale spatial separation responding to the variability of physical and chemical environmental factors. The specific functional characterization of protozooplankton, which had previously been taxonomically characterized by 18S rRNA gene amplicon sequencing, showed a diverse range of trophic preferences of picograzers, nanograzers, micrograzers, phagotrophs, decomposers, and parasites [133].

4.4.4. Metazooplankton

The metazooplankton community of the TG is also poorly known, resulting in a lack of data on trophic interactions occurring in the pelagic food web and the resource efficiency of zooplankton in relation to phytoplankton richness and harmful blooms. An early study on zooplankton [141] indicated that the spatial–temporal variability in the gulf is highly affected by riverine water inflow, pollution, and ocean circulation. Copepods and cladocerans were the most abundant groups, with high numbers of copepods in the Inner TG. The dominant species were Acartia clausi followed by Oithona nana and Podon polyphernoides. In the Central TG, the dominant species were Acartia clausi followed by Paracalanus parvus and Penilia avirostris. Appendicularians were well represented, in terms of species richness and number of individuals, in the Outer and Central TG. Another early study on metazooplankton focused on the annual cycle of neritic cladocera [142]. The distribution of mesozooplankton resting eggs was studied in the sediments of the TG and indicated a higher abundance of eggs at locations close to the river mouths, with a water column rich in zooplankters [143]. Recent research on metazooplankton reported the presence of the non-indigenous calanoid copepod Pseudodiaptomus marinus [144]. One of the very few studies dealing with marine rotifers across the Mediterranean Sea investigates the temporal distribution patterns of four coexisting Synchaeta species in the Inner TG. This work contributed to our understanding of the dynamics of this group of micrometazoans, previously underestimated among marine zooplankton, given that rotifers are lost when using larger mesh-size nets in sampling [144].

4.4.5. Planktonic Response to Climate Change and Pollutants

During the last decade, the Inner TG unicellular plankton community has been used in mesocosm experiments to explore the species response to multiple variables of climate change [145], and phytoplankton response to pollutants discharged in the seawater from exhaust gas cleaning systems on ships, known as scrubbers [1]. In regard to the climate change variables, the impacts of elevated temperature (see Section 5.2) in the form of an intense heat shock, together with changes in salinity, were examined on the natural phytoplankton community of the TG and demonstrated adverse effects on species diversity, resource use efficiency, and productivity [145]. The lower productivity and resource use function after the heat shock and the changes in salinity were supported by the lower inferred activity of the major phytoplankton taxonomic groups, as evidenced by the lower ratio of rRNA sequences to rDNA sequences (rRNA/rDNA) under these circumstances [146], which has been used as a proxy of metabolic activity in the literature. The major outcome of these studies was the synergistic negative effects on the phytoplankton community caused by the elevated temperature and increased salinity, which was also confirmed by measurements of physical properties by Androulidakis et al. [15]. Taking into consideration that climate change is associated with concurrent shifts in several variables besides temperature and salinity in the marine environment; it could also lead to more intensified and complex effects on primary producers, and this may modify the energy flow and consequently the biogeochemical cycles, therefore impacting ecosystem functioning. The potential impacts of scrubber effluents, containing several chemical compounds such as PAHs and metals, particularly vanadium and nickel, were examined in mesocosms using natural phytoplankton communities of the Inner and Outer TG. Adverse effects were observed only under high, 10% (v/v) dilutions of scrubber effluent exposure while no adverse effects were observed in communities exposed at 1% dilutions. Effluent discharge at 50 m along the ships’ wake is diluted 2000 times for vessels with established open-loop scrubbers sailing in the open sea [1].

4.5. Benthic Fauna: Diversity, Ecological Quality, and Environmental State

Despite the significance of zoobenthos in the assessment of the environmental status of marine ecosystems through biodiversity and food web descriptors, few studies have been conducted, so far, in the TG (see Section 3). Relevant attempts are almost exclusively limited to the research efforts of the Marine and Terrestrial Animal Diversity (LMTAD) and Ichthyology laboratories of the Biological Department of the AUTh, the Laboratory of Environmental Chemistry and Ecology of the Department of Environmental Engineering of the International Hellenic University (IHU), and the field studies of the HCMR. More specifically, the first research effort regarding the benthic fauna of Thermaikos Gulf was about the diversity of extant mollusks [147]. Quite a few years later, in the 1970s, the first scientific research program devoted to the benthic bionomy of the North Aegean Sea and funded by the United Nations Environment Programme/Mediterranean Action Plan (UNEP/MAP) was carried out by AUTh.

4.5.1. Benthic Communities

The seabed of the TG is dominated by various soft substrata types, which are locally covered by seagrass meadows of the species Posidonia oceanica and Cymodocea nodosa [148,149]. In the unvegetated sedimentary bottom, muddy sands constitute the predominate type [8]. The proportion of sand increases eastward and toward shallower depths, whereas the amounts of silt and clay are balanced in deeper waters and along the western coast, which is influenced by the transitional waters and terrigenous inputs of the Axios–Loudias–Aliakmonas deltaic system. Based on the European Nature Information System (EUNIS) typology (2022) of natural habitats⁶, 20 types of benthic communities occur in the TG, exhibiting a large variety of different facies, according to prevailing fauna (Appendix B). In many cases, the structure of the zoobenthos diverges from the typical one described in EUNIS, and this may be due to either the transitional stage of the studied communities, being under high pressure from various anthropogenic factors influencing the area, or the fact that the description was based on data derived from the western Mediterranean. In very shallow (around 3 m depth) waters, a narrow rocky platform exists, which together with the various artificial constructions (i.e., docks, marinas, aquaculture facilities, etc.) along the gulf represent the hard substrata type of the seabed. The hard substrata are covered by periwinkles, limpets, and barnacles in the supralittoral and mediolittoral zones [150], whereas in the sublittoral zone, four different facies of the photophilic algae community develop, depending on prevailing environmental factors [151,152,153,154].

Seven out of the twenty biocenoses found in natural habitats of the TG are considered threatened at the European level, according to the International Union for Conservation of Nature (IUCN) Red List assessment (EUNIS, 2022). More specifically, the photophilic algae community, the Cymodocea nodosa meadows on muddy sands in sheltered waters, and the littoral muddy sands and muds in lagoons and estuaries are assessed as Endangered (EN) in European waters. The Posidonia oceanica meadows, the sublittoral muddy detritic bottoms, the euryhaline/eurythermal muddy estuaries, and the mediolittoral sands are classified as vulnerable (VU), whereas the lower mediolittoral rock and the infralittoral mixed sediments are classified as near-threatened (NT). Apart from the above habitats that demand specific management plans to improve their conservation status, two facies of the stony coral Cladocora caespitosa on rocky and mixed sediments and one facies of the sea pen Crassophyllum thessalonicae should also be considered as threatened since both species are assessed as endangered according to the regional assessment of Mediterranean anthozoans [155].

4.5.2. Benthic Diversity and Species of Community Interest

Overall, 1368 benthic invertebrate species have been reported from the faunistic studies that have been carried out in the TG since 1957, with the vast majority belonging to the phylum Mollusca (57%). Annelida (18%) and Arthropoda (16%) are the next most speciose phyla (Figure 10). At the class or order level, Gastropoda prevail, followed by Polychaeta, Bivalvia, Peracarida, Eucarida, Demospongiae, and Anthozoa (Figure 11). This increased diversity of the benthic invertebrate fauna emerges from the relevant diversity in habitat types (see Section 4.5.1) and from the extensive taxonomic studies focusing in particular on the phylum Mollusca [156,157,158,159]. Among the reported invertebrate species, some are commercially exploited, such as the clams Venus verrucosa and Callista chione, the scallops Pecten jacobaeus and Flexopecten spp., the prawns Penaus kerathurus and Parapenaeus longirostris, the mud and mantis shrimps Upogebia tipica and Squilla mantis, the crabs Callinectes sapidus and Maja squinado, and the sea cucumbers Holothuria tubulosa and H. poli [160,161,162,163], raising concerns regarding the future viability of natural stocks and habitats.

Other species, however, face direct threats, such as the endangered corals Crassophyllum thessalonicae (Figure 12a) and Cladocora caespitosa (Figure 12b) and the critically endangered fan mussel Pinna nobilis (Figure 12c). The soft coral C. thessalonicae, originally described by Vafidis and Koukouras [164], lives exclusively in the TG, where it forms sparse populations (0.5 to 4 colonies/20 km) on the silty sediments of the Central TG. The habitat of this species is heavily disturbed by bottom trawling, urbanization, and pollution [155]. There are no active protection measures, so far, for the sea pen, whose population remains outside of the established SACs of the TG. The pillow coral Cladocora caespitosa forms moderately dense populations on coarse sedimentary and rocky habitats, especially along the eastern coast of the gulf [154]. Climate change and fisheries mainly threaten the species [155], whereas in areas influenced by sewage treatment plants, siltation constitutes an additional threat, as the polyps of the colonies are buried within the silty sediment. Cladocora caespitosa is an ecosystem-engineering species due to the provision of a biogenic habitat that constitutes an important biodiversity reservoir [165]. Although there are no specific conservation measures for the species, a part of its population in Thermaikos is protected within the SAC GR1220005. The Thermaikos Gulf hosted the densest population of the fan mussel P. nobilis in the Greek Seas [166] with densities reaching 1.04 individuals/m² [167]. Unfortunately, the species has been affected by Haplosporidia parasites and mycobacteria that caused mass mortality events (MMEs) that totally devastated the population [168]. Currently, only dead eroded shells heavily colonized by epiphytes and sessile invertebrates (Figure 12c) may be found in the TG, including the designated SACs in the region.

4.5.3. Ecological Quality Based on Benthic State

According to European Union directives MSFD and WFD, the monitoring of coastal and transitional marine waters is mainly assessed through biological quality elements under an integrative approach that uses physicochemical data as supporting variables. Among the four proposed biological quality elements (i.e., phytoplankton, phytobenthos, zoobenthos, and fish), the zoobenthos, which includes all animal organisms living in the benthic zone, is the most promising tool for developing appropriate indices to assess the ecological quality status and the seabed integrity. This is due to its widespread and abundant presence in both species’ richness and number of individuals within marine habitats, the rather limited natural variability in the synthesis of its biotic communities, and lastly its fundamental life traits. The zoobenthos species are relatively long-lived sessile or discretely motile organisms that do not perform large-scale migrations, and accordingly, their synthesis reflects the cumulative environmental conditions over time in a specific site. Moreover, the different species have different tolerance limits and can respond to various stress factors, such as eutrophication, organic pollution, mechanical disturbance of substrata, and heavy metal load. They also play a main role in food webs, nutrient cycling, and energy flow in benthic–pelagic coupling. Several biotic indices have been developed to assess the ecological quality of marine waters based on benthic invertebrates, in full compliance with the WFD. Among the various proposed indices, the BENTIX index originally developed by Simboura and Zenetos [169] has been successfully inter-calibrated between European countries to be used over Mediterranean coastal waters, and the M-AMBI index [170,171] has been developed for transitional waters [172].

The National Monitoring Water Network (NMWN⁷) estimates the ecological quality of coastal and transitional waters in the TG by analyzing the synthesis of benthic macroinvertebrates and applying the BENTIX and M-AMBI indices, respectively. The operational program has been, so far, implemented in three periods: an initial reference period (2000–2008), a first implementation period (2012–2014), and a second implementation period (2018–2022). According to relevant reports, the ecological quality in the Inner and Central TG has been assessed as moderate in the reference and the first implementation period and improved to good in the second one [173,174,175]. In the Outer TG, the ecological quality state, originally assigned as good [173], has degraded to moderate [174,175]; however, the latest assessment has not been based on biological samplings, but on an extrapolation of the previous evaluation. Considering the transitional waters, the state has been steadily assessed as moderate [174,175].

Based on the few studies that have applied the BENTIX index in the TG [28,29,30,31,32,33,38,152,172,176,177,178,179], the ecological quality state alternates from moderate to good according to the specific geographic area and the sampling period, revealing increased spatiotemporal variability in ecological quality status. A spatial trend of a decreased state in the Inner TG together with an improved state in the Outer TG emerges, whereas in the Central TG, the eastern shores appear to be in a better state that the western ones, which are influenced by the estuarine systems of Axios–Loudias–Aliakmonas. Seasonal assessment in three sampling stations in the Inner and the Central TG during the period 2017–2020 [28,29,30,31] revealed a gradual improvement in the quality of the western shores but a deterioration in the eastern ones, where the ecological quality state shifted from moderate to good and vice versa, respectively. In the transitional waters of Thermaikos, all relevant studies report moderate quality in Loudias and Aliakmonas estuarine systems [32,33] and good quality along the western coasts of the Outer TG [33].

5. Discussion and Concluding Remarks

Thermaikos Gulf (TG) is a productive, natural, coastal ecosystem with extensive estuarine systems that is highly impacted by anthropogenic pressures. Severe deterioration of its quality has been documented between the 1970s and 1990s; thereafter, a steady trend of recovery has been identified, especially after the operation of the two main WTPs in the late 1990s, though without achieving, so far, the goal of good environmental state, as outlined by the European Union’s Marine Strategy Framework Directive (MSFD) and Water Framework Directive (WFD).

5.1. Research Trends and Gaps

Most of the monitoring surveys were conducted from the mid-1990s until the late 2000s, during the “golden period” for research in the TG. The most investigated region is the northern part of the TG (Inner TG). Although Thessaloniki Bay is more polluted and its impact is more direct on the coastal area near the urban waterfront of Thessaloniki (red tides, dirty sea phenomena), most of the research has consistently focused on the southern part of the Inner TG (Central Gulf of Thessaloniki), and especially around the river deltas. Very few studies integrally include the entire TG (<25%), and these studies were mainly conducted during the 1990s. A clear reduction in research was detected after 2010, associated with the recession of the Greek economy. The global financial crisis that started in the late 2000s had several impacts on the socioeconomic environment of Greece, with strong implications for academia and, by extension, for environmental research and monitoring activities. The fiscal austerity years after 2009 are characterized by the brain drain of young scientists [180,181], budget cuts in applied research and especially environmental monitoring [182], decreasing human resources in research institutions and universities [41,183], and reduced funding for the management of environmentally protected areas (e.g., Natura 2000; [1]). These factors, combined with the complex bureaucracy imposed by national and European Union-funded research projects [184,185], the continuously changing Greek legislation regarding the management of funds [186], the over-taxation of those funds, and the endless delays, imposed a tremendous burden on research development [41]. The lack of research efforts is even more challenging, given the fact that fiscal austerity measures have further intensified the anthropogenic pressures and environmental impacts, especially in urban coastal areas such as the TG [187]. The misuse of biotic resources in protected areas reached unprecedented levels in Greece during the period of the economic recession [42]. Scientific research activity focusing on the TG marine ecosystem is still relatively limited, although the country’s economy shows signs of stabilization (or recovery in some financial sectors, like tourism, that may, on the other hand, negatively impact the coastal environment) after more than a decade of recession. The COVID-19 pandemic after 2020 and the successive lockdowns have impacted scientific research even more [188], especially in environmental monitoring that requires group cooperation and higher funding; safety restrictions and different funding priorities have created additional gaps in the TG’s monitoring.

The continuous operational monitoring of the TG is deficient and is based only on four stations operated by the NMWN, where only one in Thessaloniki Bay has collected measurements after 2020. The entire TG has been originally documented in a moderate state, whereas according to the recent NMWN results, the ecological quality has been improved in the Inner and Central TG’s coastal waters (good state) but remained moderate in the Outer TG. This outcome seems odd, as the Outer TG, influenced by waters of the open Aegean Sea that are classified as excellent (see Section 4.1), appears to be in lower quality than the Inner and Central parts of the gulf, which are highly impacted by eutrophication and organic pollution. This specious observation is probably due to the fact that the sole station from the Outer TG has not been monitored again by NMWN after its initial assessment. This highlights the severe limitations of the existing operational NMWN network, which is characterized by reduced temporal and spatial (only two or three stations) coverage that is not only unable to cover a wide and heterogeneous marine area impacted by various anthropogenic pressures, but may also lead to uncertain or erroneous conclusions about the quality state of the TG. The increase in observational stations, both temporally and spatially, to cover the entire domain of the gulf is crucial. The monitoring approach should include all the major environmental parameters and provide a continuous and holistic view of the TG’s environmental indicators. Continuous monitoring (e.g., with permanent monitoring stations/buoys) over the entire water column at the two entrances (Cape Mikro Emvolo and Cape Megalo Emvolo) would also provide vital real-time information about the renewal ability of the enclosed basins.

Necessarily and additionally to the monitoring of the marine environment, systematic measurements of supply and freshwater quality must also be included in the rivers and channel networks that discharge into the TG, something done circumstantially in the study area. Open access data are unavailable. The lack of measurements of the freshwater supply and quality leads to uncertainties and an inability to estimate potential environmental impacts due to the influx of nutrients and other pollutants to the coastal system.

While several public agencies oversee certain aspects of the TG environmental monitoring, the absence of a supervisory public entity, capable of organizing and overseeing all coastal research and protection endeavors, is a major hindrance. Such an organization would facilitate and harmonize monitoring initiatives across the entirety of the gulf’s expanse. This deficit may lead to managerial inconsistencies and “dispersion” of responsibilities, especially during unexpected hazardous events (e.g., pollution accidents; [15]) that require immediate actions by first-level responders and a focused operational plan. Another benefit of such an entity (organization or association or policy actor) would be the holistic overview and management of the entire coastal environment (atmosphere, land, and ocean) based on coupled monitoring of air–land–sea interactions that are currently not considered.

Finally, some aspects of the marine and coastal environment of the TG, although having been covered by mass media and having captured the public’s interest, have received minimal attention from authorities, with only occasional and inefficient initiatives for protective action and integrated coastal management during the last decades⁸. Although coastal erosion and shoreline retreat is a major problem, especially along the eastern coasts of the Central Gulf of Thessaloniki (e.g., Perea, Agia Triada) and the western coasts of the Central (e.g., river deltas) and Outer TG (e.g., Katerini beaches), assessment and monitoring of its evolution and, more importantly, testing and implementation of mitigation and protection measures are completely absent. An example of the Axios delta shoreline retreat between 2019 and 2024, based on the Normalized Difference Water Index (NDWI) derived from two Sentinel-2 L2A images (10 m resolution; Copernicus Dataspace Browser⁹), is presented in Figure 13. The estimated reduction in the delta area is approximately 11 km². This reduction is alarming, considering the ecological importance of the delta, and indicates the urgent need for coastal monitoring and protection in the area, especially in the environmentally fragile wetlands of the gulf.

Another human-induced environmental stressor that severely impacts both biodiversity and human health is the high concentration of marine plastic litter in the TG [27,189]. This stressor is also under-monitored and unfortunately not included in any current protection plans and management strategy for the region, although it is included in Descriptor 10 of MSFD. Microplastic pollution in the TG is mainly controlled by river discharges [190,191] which carry significant amounts of microplastics (40–50 times higher than the maximum ocean surface concentrations; [192]); monitoring and mitigation of this stressor in the rivers are also very limited.

Therefore, several limitations and challenges of TG-related environmental research can be identified as a result of our exhaustive review:

(a)

Temporal issues: Our review synthesizes findings from the last 50 years of research; while this provides a substantial historical perspective, it also implies that some recent developments or changes may not be fully captured. Research conducted in the early years might have different methodologies, monitoring instruments, calibration techniques, and finally data quality compared to contemporary studies.

(b)

Spatial issues: There is an ad hoc spatial diversification of monitoring density in several areas; the open sea waters are more difficult to cover compared to coastal areas. Delineating the spatial coverage of monitored areas can be complemented only by extended numerical modeling, which is referred to in our review, but only for physicochemical processing yet not biological and ecological modeling.

(c)

Data scarcity and quality: Potential limitations in funding (e.g., fiscal austerity) and available resources (e.g., pandemic lockdown) could have impacted the availability and quality of data collected during these periods. Moreover, recently or currently (<4 recent years) conducted research may have not been published yet, rendering the study region a scarce-data area at present.

(d)

Integrating disciplines: A further limitation is the poor status of past interdisciplinary research in the relevant literature covering socioeconomic issues, besides marine biology, environmental science, and hydrology/hydrodynamics, which are already extensively covered by the review.

(e)

Overarching suggestions: The suggested increasing, establishing, and maintaining of new monitoring infrastructure can be resource-intensive and logistically challenging. Moreover, recommendations for the establishment of a supervising public entity indicate a need for better governance and management of the TG. However, implementing such recommendations requires political will, stakeholder engagement, and institutional capacity, which may face barriers such as bureaucratic hurdles, conflicting interests, and limited resources.

5.2. Thermaikos Gulf under Climate Change

Evidence of a changing climate can be recognized in the recent evolution of physical properties in the seawater of the TG. Figure 14 shows the prevailing ocean warming conditions in the gulf as derived from daily satellite sea surface temperature (SST; Copernicus Reprocessed Mediterranean L4 dataset; 0.05° resolution grid) between 1982 and 2023. The annual mean SST, averaged over the entire TG, reveals a strongly increasing rate (Sen’s slope of 0.52 °C/decade; Figure 14a), higher than the general trend of the broader Aegean, Ionian, and Cretan seas (0.49 °C/decade; [54]). The TG is characterized as a “hot spot” area for MHWs, with events that not only occur near the sea surface, but also extend through the entire water column, especially in the shallower Inner TG [15]. A statistically significant (p_value < 0.0001) SST shift of around 1.2 °C was computed after 2006, based on the Pettitt homogeneity test [51], while the mean SST increased from around 17.5 °C in 1982 to 19.5 °C in 2023. The SST trend follows the respective air temperature interannual increase (0.43 °C/decade) that also showed a clear shift after 2006 (Figure 5b). The highest SST peaks of the maximum values (99th percentiles) were observed in 2021 and 2023 (>29 °C; Figure 14a). The summers of these two years were the warmest of the entire period, affecting several environmental and socioeconomic factors of the area (e.g., aquaculture; [54]). The related trend slope of the maximum levels is even steeper than the annual mean (0.67 °C/decade), increasing by approximately 2.5 °C during the 1982–2023 period. The respective MHWs were computed based on the satellite-derived SSTs and the methodology introduced by Hobday et al. [193]; the annual MHW cumulative intensity (°C × days), averaged over the entire TG, is presented in Figure 14b. The events were defined using the 90th percentile climatology (1982–2023) as the threshold baseline and a 5-day minimum duration for each event. Very high MHW intensities occurred successively after 2018, confirming the intensification of warming in the TG. Intensity peaks were computed for 2012 (related mainly to winter MWHs as derived from the high mean SST and relatively low 99th percentiles; Figure 14a) and 2021. Especially in 2021, the TG’s MWHs had significant biochemical implications, affecting the marine environment and causing mass mortality events of prominent biota with considerable impacts (e.g., on aquaculture, fishing, tourism, human health; [54,165,194]). The importance of MHWs raises the need for accurate prediction tools for both short-term (e.g., a few days based on operational hydrodynamic modeling) and interim-future periods (e.g., forecasts of weeks to months based on artificial intelligence techniques; [195]).

An increase has also been observed in salinity, which reached the highest measured levels after 2020, associated with elevated temperatures and possible changes in precipitation rates and river discharge outflows [15]. Temperature increases together with changes in salinity may also affect species diversity, resource use efficiency, and productivity [145]. Experiments in the TG [145] showed how changes in the species inventory driven by an MHW would influence the ability of the resultant phytoplankton communities to cope with salinity changes, especially in the context of climate change, where short-term temperature changes during heatwaves are expected to increase in frequency and magnitude during the 21st century [196].

Temperature also greatly affects the equilibrium constants of the CO₂–carbonate system and particularly the solubility coefficient of CO₂; higher temperature and salinity decrease the solubility of CO₂ (e.g., for a temperature increase of 1 °C, pCO2 rises by ~4%, [197]), which may induce CO₂ oversaturation and change of the character of the area, potentially acting as a CO₂ source to the atmosphere. High-resolution simulations of the physical and biogeochemical state of the Mediterranean were used to study the climate-change-related impacts on the marine ecosystems of the basin under different future CO₂ emission scenarios [198,199]. The projections showed a DIC concentration increase as well as an acidification signal in the upper water column, as a result of increased respiration, but also because of the increase in atmospheric CO₂ and of the changes in CO₂ solubilization driven by changes in temperature and salinity, with stronger changes in the eastern basin than in the western one. Even though future projections are unfavorable and the TG is a coastal ecosystem characterized by complex biogeochemical processes that are directly linked to the CO₂–carbonate system, there is a lack of regular monitoring of related parameters.

The limitation of freshwater and sediment discharge has led to erosion of the protected deltas and coastline retreat in more recent years [8]. In this regard, the risk of intensification of coastal retreat and flooding under rising sea levels in the future may have severe environmental and socioeconomic impacts. The interannual variability of the sea level anomaly over the TG is also characterized by a strong increasing trend of approximately 4 cm/decade during the last 30 years, resulting in increasing flooding frequencies over the low-lying coastal areas [200]. Climate change is also associated with concurrent shifts in several variables of the marine environment and not only in physical parameters (temperature, salinity, sea level), leading to more complex effects on primary producers and impacts on the biochemical cycles and ecosystem functioning. Despite these alarming signs, it is noted that climatic studies analyzing long-term future projections of physicochemical and biological features, based on different climatic scenarios combining atmospheric inputs from ensemble modeling, are absent for the TG, while the research over past climate-scale periods forming the present-day situation corresponds to only 2% of the entire research activity in the gulf.

5.3. Overall Environmental and Ecological State

The major supplier of pollutants (e.g., nutrients, metals) and sediments in the TG is the Axios River in the western coastal zone, as it traverses extended agricultural plains with intensely fertilized farmlands. Metal concentrations were mainly detected in the estuaries and in the northernmost part of the gulf (Thessaloniki Bay), where renewal conditions are weaker and numerous human activities have been concentrated for decades (Thessaloniki metropolitan area). Water renewal in this enclosed basin is mainly controlled by the prevalence of northerly winds and two-layer transport across the entrance of each sub-basin. Contrastingly, equally prevalent southerlies may weaken the renewal and constitute a favorable pre-condition for the formation of eutrophication events. Although the nutrient pollution of the TG was reduced during the last two decades, after the strong reduction in untreated urban wastewater discharges in the marine environment, eutrophication events are still very common, especially in the shallower parts of the Thessaloniki Bay’s coastal waters, where the ability of hydrodynamic circulation for renewal is limited. Organic chemicals (e.g., organochlorine pesticides, PCBs) are also major pollutants across the entire TG and have been detected in higher concentrations in Thessaloniki Bay and in the western Outer TG that may even exceed the environmental quality standards set by the European Union for priority hazardous substances in surface waters [201]. The Inner TG is also impacted by hydrocarbons and can be characterized as moderately polluted with petroleum compounds. Hydrocarbon leakages by tankers occurring around the large-ship anchorage and the supply pipeline to the local oil refinery industrial infrastructure, located near the port of Thessaloniki, enhance coastal water pollution locally. In addition, scrubber effluents, which are usually warmer and more acidic than the surrounding water, may increase acidity and make toxic heavy metals more bioavailable to wildlife.

The multi-parametric Eutrophication Index, derived from measured nutrients and chl-a concentrations, confirmed an improvement in the trophic state of Thessaloniki Bay and Axios delta after 2000 (from “Poor” to “Moderate”) without, however, reaching a “Good” quality level, despite the operation of two main WTPs. The distribution of phytoplankton species is strongly controlled by the prevailing circulation and the concentration of nutrients, with large species pools of nutrient opportunists near the river deltas. Significantly higher phytoplankton abundance was observed in the Inner and Central TG compared to the Outer TG, demonstrating the oligotrophic state of the latter. The Inner TG is still highly eutrophic, with persistent blooms leading to episodes of accumulation of organic mucilaginous material.

The TG is a marine area of high importance for the benthic invertebrate fauna of the Greek Seas, due to its increased biodiversity at both the taxonomic and ecological levels. Several endemic, protected, threatened, and fishery-exploited species have dense and widespread populations in the gulf, as do invasive ones, which seem to spread via navigation. The area hosts a unique species of sea pen, C. Thessalonicae, which is threatened and not yet protected. Since the 1970s, the benthic state has diachronically degraded in the Inner TG (Thessaloniki Bay) and on the western shores of the Central TG, whereas the relevant state has improved on the eastern shore and in the Outer TG. However, this pattern seems to be currently reversed, as the ecological quality is becoming worse on the eastern and better on the western shores. The ecological status and especially the seabed integrity can be assessed based on the evaluation of the zoobenthos. In the Outer TG, the ecological quality, based on benthic-derived indexes, has degraded from “Good” to “Moderate” during the last decade. Transitional waters have also been characterized as “Moderate”. Despite the limited time series of available data, all relevant studies reveal organic pollution as the main threat to the benthic biota. Infectious epidemics and MHWs should be added to the main threats, as they have caused mass mortality of the fan mussel P. nobilis since 2019 and of other shellfish, including the cultivated Mediterranean mussel Mytilus galloprovincialis that collapsed in 2021 due to intense MHWs. Other threats to benthic communities are related to excessive fishing and especially bottom trawling.

In conclusion, the TG emerges as a resilient yet fragile coastal ecosystem, facing intense natural stresses and human-induced pressures. Despite a discernible trajectory of recovery following a period of degradation, achieving the environmental objectives outlined by regulatory frameworks remains elusive. Key recommendations involve the establishment of a broad monitoring network with both permanent and frequently revisited stations, identification and observation of anthropogenic stressors, bolstered protection measures for vulnerable benthic communities, and integration of advanced forecasting technologies. The establishment of a supervisory body is advocated in order to assist governance and take action toward sustainable solutions for mitigating anthropogenic impacts and fostering eco-friendly practices which are paramount for the resilience of the TG and the maintenance of its ecosystem services in the future.