Oceanographic Research in Thermaikos Gulf: A Review over Five Decades

1. Introduction

Thermaikos Gulf (TG) is a semi-enclosed microtidal coastal basin of the mid-latitude east-central Mediterranean Sea, nurtured by multiple (transboundary and local) riverine inflows, and it is connected to the North Aegean Sea through its southern boundary (Figure 1a). 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. TG covers approximately 3,300 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^oE, 40.48^oN), the Central TG, located between Cape Megalo Emvolo and Cape Epanomi (22.88^oE, 40.37^oN) and the Outer TG southwards from Cape Epanomi (Figure 1a). Moreover, the Inner TG can be divided in two smaller basins, north (Thessaloniki Bay) and south (Central Gulf of Thessaloniki) of the Cape Mikro Emvolo (22.93^oE, 40.58^oN; Figure 1b). The Inner TG is characterized by depths lower than 40 m, while the Central TG is also relative 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). TG is of regional, national and international interest and importance 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 (Kaberi et al., 2023). It is a unique system of complex, natural processes severely, influenced by human-related stressors, especially those involving dissolved pollutants or nutrients and particulate phase matter within the water column (Price et al., 2005). The aim of this review is to compile the research findings conducted over the TG and integrate the knowledge on the oceanographic processes through the last five decades.

The hydrography and quality of the TG are both affected by several factors, both environmental (e.g., meteorological conditions, freshwater discharges, water mass exchanges with the open sea) and anthropogenic (e.g., encumbering human activities, marine and fluvial pollution, coastal erosion due to waterfront infrastructure; Kaberi et al., 2023). The regional climate can be characterized as maritime temperate (Flocas and Arseni-Papadimitriou, 1974), typical of Mediterranean coastal inlets, where annual air temperatures vary between 9^oC and 17.5^oC and annual rainfall range from 400 mm to 1300 mm (Poulos et al., 2000). 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 winter (namely “Vardaris” wind; Koletsis et al., 2016) and summer months (named “Etesians”; Anagnostopoulou et al., 2014). Five large rivers, four of which (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 flows’ input is further enhanced by smaller intermittent rivers and streams (e.g., Anthemountas, Figure 1b), and is conditionally 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 (Androulidakis et al., 2021; 2023a). The two largest rivers, Aliakmonas and Axios, accompanied by Loudias in between, form a complex integrated multi-channel deltaic system, with adjacent protected wetlands 1, i.e. several Natura 2000 sites, a Ramsar Site, a National Park 2 a coastal lagoon and many wildlife refuges (Vokou et al., 2018). Near the eastern boundary of the integrated deltaic area, Gallikos (a former intermittent flow river) acts as an overflow canal branch of Aliakmonas, during the recent decades. Its slow flow forms the Gallikos estuary, a vibrant ecosystem of three neighbouring ponds, home to common and rare species of migratory birds 3. Axios supplies more than 50% of the total freshwater (Karageorgis et al., 2003) in the TG. Seasonality is not very strong in outflow rates because the daily releases are controlled by the operations of hydroelectric power dams (Poulos and Chronis, 1997). The discharge rates range between minimum and maximum daily mean values are of the order of 20 m³/s and 300 m³/s, respectively (Karageorgis et al. 2001; 2005a; Androulidakis et al., 2023a). Freshwater outflow rates, have shown significant decrease during the past decades (Baltas and Mimikou, 2005), associated either with the increased irrigation needs or with reduced precipitation rates (Feidas et al., 2006; Pakalidou and Karacosta, 2018).

Figure 1. Bathymetry of (a) Thermaikos Gulf (TG), divided in three basins (Inner TG, Central TG and Outer TG), and (b) Inner TG, divided in 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 etc.) are marked (Krestenitis et al., 2012). The insert map in (a) shows the location of TG in Mediterranean Sea.

Figure 1. Bathymetry of (a) Thermaikos Gulf (TG), divided in three basins (Inner TG, Central TG and Outer TG), and (b) Inner TG, divided in 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 etc.) are marked (Krestenitis et al., 2012). The insert map in (a) shows the location of TG in Mediterranean Sea.

The densely populated urban areas along the coastal zone, accommodate several infrastructures and human activities, such as farming, treated and untreated sewage discharges, industrial enterprise, 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 on the coastal zone of the Inner TG, and especially in the more polluted Thessaloniki Bay that contains the main population and industrial infrastructure (Figure 1b). In spite of 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 (Zafirakou et al., 2015), and marine litter (Politikos et al., 2017), downgrading the quality state of this ecologically sensitive marine environment.

The first published research about the marine environment of TG dates back in 1981 (Chester and Voutsinou, 1981) and describes the occurrence of pollutant trace metals in sediments. 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 focused at specific coastal sites of environmental and socioeconomic interest. Herein, we review the environmental state of TG as derived by approximately 5 decades of field and numerical studies (1970s-2020s). The motivation for this retrospect study is to gather the main findings of all previous research efforts, and highlight the gaps of the conducted research over a socially and environmentally important Mediterranean marine area. The final goal is to provide robust recommendations for future monitoring and research that will reduce the deficiencies and contribute to the improvement of the gulf’s quality, assisting the adaptation of optimal and holistic management solutions. The statistical analysis of all monitoring activities, numerical simulations, satellite-derived observations, and laboratory experiments was based on published research conducted from the 1970s up to date (approximately five decades; Section 2). The most important findings, associated with all oceanographic disciplines (i.e., physical, geological, chemical, biological) are described in Section 3. The discussion about the overall TG’s environmental quality state, the impact of climate change, and the most profound gaps and recommendations are included in Section 4. Appendix A presents all the acronyms used in the study.

2. Statistical Analysis of the Studies About Thermaikos Gulf

We have gathered and examined all available literature regarding the marine state of TG during the last 5 decades. The total sample includes 208 papers, published either as articles in peer-reviewed journals (94%; Table 1) or as official technical reports of monitoring surveys (6%). 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). Only a small percentage of articles (2%) were literature reviews of previous studies.

The evolution of the annual number of papers published about 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 timeseries show a similar distribution with lower numbers in the beginning and the end of the 50-year period. The highest Pearson correlation coefficient (R_P=0.56; serves the Mann-Kendall (Mann, 1945; Kendall, 1975; 99% test of statistical significance: p_value<0.0001) between the two series was derived after shifting the publications’ timeseries backwards 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 the conducted research activities was in 2001, while the respective peak of the published work occurred by 2005. Generally, the period with 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 TG was reduced by almost half (18%; Table 1) compared to the previous decade (39% in 2000-2009; Table 1). The 2010-2020 decade showed even less research activities than the 1990s (20%), even though the research staff and investigation tools are expected to be more advanced and qualified, as years go by. It is noted that the Greek financial crisis, that started in 2009, had various socioeconomic consequences that affected all aspects of the academic and research activity (Koulouris et al., 2014), although the environmental pressures were intensified due to the economic recession (Troumbis and Zevgolis, 2020), increasing the need of environmental monitoring and holistic management even more (see Section 4.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 TG is profound after 2010, and especially between 2009 and 2015, probably related to the aforementioned socioeconomic factors and it is not justifiable by the expected delays in the publication process.

The majority of the research includes field observations (86%), while only 19% concerns numerical applications. Satellite-derived data (5%; Table 1) were mainly used to cover the period 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; Table 1). During the “golden period” of intense research activity, from the mid-1990s until the end of 2000s (Table 1; Figure 2), observational monitoring was the main method to conduct research over TG, peaking in 2001; there is a reduction after 2010 and a slight increase in the mid-2010s (Figure 3a). Post 2010, field-work 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 (Table 1). Investigation of physical processes and hydrography was less than 40%, while geological and sediment processes were the least examined among all disciplines (25%). There are several years that 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 examined research focused on the Inner TG (Table 1), with 77% of all studies 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 against the Outer TG during the last decade even more (Figure 3c). The distribution of monitoring frequency is scanty in the Outer TG, and especially in its central and eastern parts (>22.9^oE; Figure 4) where observational locations are very sparsely distributed. A large number of observational locations is 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 higher number of monitoring stations (20-30; Figure 4). A large number of field studies has also been conducted in the entrance of Thessaloniki Bay (Cape Mikro Emvolo) and in the east side of the Inner TG entrance (Cape Megalo Emvolo). Most of the studies focus on the Central Gulf of Thessaloniki (63%; Table 1), where intense aquaculture and agriculture activities take place along its western coastal zone (Figure 1b); many stations have been placed around 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 smaller 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). The research surveys of the Inner TG were equally divided across its sub-domains throughout the study period although the total number of them decreased after 2009 (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 (Table 1). Even fewer studies cover both the water column and the seabed (21%). Benthic processes and properties related to, e.g., benthic fauna, ecological status, sediment characteristics, etc. have been investigated less in comparison with other properties of the water column. However, note that benthic processes are inherently more slowly evolving through time than the seawater conditions in the marine environment, hence the generally lower frequency of respectively recorded field data. The same applies to the spatial extent of seabed field observations (only a few sporadic samples can be representative of the benthic health condition).

The temporal characteristics (duration and step) of the research studies are presented in Figures 3e and 3f. The short-term studies that cover periods of less than one year, comprise 48% of the total (Table 1). Study periods of one year or more (long-term) are less frequent (43%). The short-term studies mainly covered the years after 1995 until 2009, with a second, however smaller increase after 2015 (Figure 3e). Only long-term studies occurred from the late 1980s until mid-1990s (Figure 3e). The temporal step of analysis in most studies exceeds one month, usually corresponding to 3-month (seasonal) monitoring (42%; Table 1). Monthly measurements are significantly less frequent (16%; Table 1), 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%; Table 1), included mainly in studies between 1990 and 2010 (Figure 3f). Contrastingly, coarser temporal steps (e.g., bimonthly, seasonal, or semi-annual) were usually chosen to investigate the marine environment of TG in the 1970s and during the decade after 2011 (Figure 3f). Besides monitoring activities and realistic numerical simulations, a small number of studies (Table 1; Figure 3e) refer to past-climatic periods (2%; e.g., decadal or centurial 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.

3. Main Research Findings

3.1. Hydrography and Circulation

The first publication simulating the hydrographic conditions of TG was authored by Ganoulis and Krestenitis (1982), 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 (Voutsinou and Balopoulos, 1989). 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 were undertaken (analysed in the following Sections), leading to a deeper understanding of the dynamics and dominant forcing factors in the area.

3.1.1. Main Circulation Patterns

The distribution of physical properties and the general circulation of TG are controlled by four major environmental factors: (i) the air-sea interactions and especially the wind forcing (Kontoyiannis et al., 2003; Androulidakis et al., 2021; 2023a), (ii) the exchanges with the northern Aegean Sea through the southern boundary (Hyder et al., 2002), (iii) the river plume dynamics (Kourafalou, 2001), and (iv) the morphological setting of the enclosed basin (Kourafalou, 2001; 2004; Kontoyiannis et al., 2003). The tidal amplitudes are small in the TG, ranging up to a few tens of centimetres (Poulos et al., 2000; Tsimplis et al., 1995; Kontoyiannis et al., 2003), leading to very weak corresponding tidal currents (<0.02 m/sec; Hyder et al., 2002). The river plume waters affect both the physical characteristics (density currents) and the biochemical properties of TG coastal waters, due to their high concentrations of nutrients (Karageorgis et al., 2005a; see Section 3.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 (Kourafalou et al., 2004).

The morphological separation of TG into three main sub-basins of 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 at each sub-basin, controlling the connectivity pathways, the transport across the sub-basins’ entrances, and finally 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-years investigation of the ERA5 Reanalysis database (Figure 5a), in agreement with older studies (e.g., Poulos et al., 2000). This wind regime may form a cyclonic ocean circulation pattern in the Central Gulf of Thessaloniki (Androulidakis et al., 2023a; Figure 6a). They also enhance the southward transport of near-surface waters towards the Outer TG (Kourafalou and Tsiaras, 2007), removing the riverine brackish masses away from the confined northern parts of the gulf (Hyder et al. 2002; Krestenitis et al., 2012) (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 (Zervakis and Georgopoulos, 2002; Kourafalou et al., 2004). It also allows the northward cyclonic intrusion of denser Aegean Sea Waters (ASW) in the deeper layers (Karageorgis and Anagnostou, 2001) that facilitates the renewal of the Inner TG with clearer waters (Krestenitis et al., 2012; Androulidakis et al., 2021). Such wind conditions can prevail during the entire annual cycle and enhance the renewal of the Inner TG, limiting the eutrophication events (Androulidakis et al., 2021), 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 (Androulidakis et al., 2023a). A clear southward connectivity pathway, between the more polluted western Thessaloniki Bay (Christophoridis et al., 2009; see Section 3.3) and the Central Gulf of Thessaloniki (Figure 6a), also affects the usually clearer waters along the eastern coasts under northerly winds (Androulidakis et al., 2023a). 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 TG with the southern coasts and the Aegean Sea (Kourafalou et al., 2001; Kontoyiannis et al., 2003). 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 summer months (Kontoyiannis et al., 2003; Androulidakis and Kourafalou, 2011; Kourafalou and Tsiaras, 2007).

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 (Androulidakis et al., 2021), spread the freshwater plume over the Inner TG (Krestenitis et al., 2012) and induce an anticyclonic circulation pattern over the Central Gulf of Thessaloniki (Androulidakis et al., 2023a) (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”; Kourafalou, 2001; Kourafalou et al., 2004). 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 (Kourafalou, 2001), especially under the influence of southerly winds (Krestenitis et al., 2012; Androulidakis et al., 2023a). Southerly winds may also spread river plume waters towards the Inner TG, even reaching Thessaloniki Bay (Krestenitis et al., 2012). The ASW intrusion is weaker along the eastern coasts, evolving mainly over the Outer and Central TG. The two-layer flow is weaker and can even become reversed in comparison to the northerlies’ case (e.g., Cape Mikro Emvolo; Figure 6b).

Figure 6. Schematic 3-d representations of general circulation features (Cyclonic and Anticyclonic patterns) of TG under the prevailing (a) Northerly (cold period) and (b) Southerly (warm period) winds. River plume spreading, 2-layer flows, Aegean Sea Water (ASW) inflow, Deep Water Formation (DWF) under buoyancy loss conditions, and wind-induced (Northwesterlies: NW, Southerlies: S, Southeasterlies: SE) upwelling and downwelling processes are also marked.

Figure 6. Schematic 3-d representations of general circulation features (Cyclonic and Anticyclonic patterns) of TG under the prevailing (a) Northerly (cold period) and (b) Southerly (warm period) winds. River plume spreading, 2-layer flows, Aegean Sea Water (ASW) inflow, Deep Water Formation (DWF) under buoyancy loss conditions, and wind-induced (Northwesterlies: NW, Southerlies: S, Southeasterlies: SE) upwelling and downwelling processes are also marked.

3.1.2. Thermohaline Variability and Stratification

The vertical structure of the water column reveals significant spatial and temporal variability, with periods of stratified or homogenised conditions. The characteristics of the surface layer vary considerably, although the lowest-salinity waters are usually concentrated in the western gulf (Hyder et al., 2002). Increased riverine outflows hinder vertical mixing, resulting to haline stratification (Tragou et al., 2005), especially from December to April (Hyder et al., 2002). Strong cold northerly winds, enhance vertical mixing and the homogenisation of the water column, especially in the shallower Thessaloniki Bay (Krestenitis et al., 2012), contributing on the formation of a southward-flowing, dense, near-bed current (Hyder et al., 2002; see Section 3.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 sea depths in the semi-enclosed northern part of the gulf (Krestenitis et al., 2012). The vertical homogenization during winter months results to a prevailing barotropic circulation, in contrast with the summer conditions, when the water circulation is mainly baroclinic (Kontoyiannis et al., 2003; Zervakis et al., 2005). Strong pycnonclines (combining diverse thermoclines and haloclines) may also form in the near-bottom layer due to the intrusion of colder and more saline ASW, that either remain within the Central TG under southerly winds, or enter into the Inner TG under the effect of northerly winds (Poulos et al., 2000; Karageorgis and Anagnostou, 2001). 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 (Androulidakis et al., 2023a) than values measured during previous monitoring campaigns (e.g., 1994-2007; Krestenitis et al., 2012). Pertaining to this, sea temperature of the TG has also revealed significantly high peaks, especially during the last decade, and strengthened the formation of Marine Heatwaves (MHW), designating the TG as a MHW “hot spot” (Androulidakis and Krestenitis, 2022; see Section 4.2).

3.1.3. Vertical Water Fluxes

The Inner TG also appears amongst the regions of Dense Water Formation (DWF) in the Eastern Mediterranean Sea, process taking place under intense cooling by northerly gales (Estournel et al., 2005; Zervakis et al., 2005; Krestenitis et al., 2012). The phenomenon can be triggered by density gradients due to temperature, due to salinity or under the contribution of both (Shapiro et al., 2003). 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 favoured 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) that typically cascade along the continental over the Outer TG and towards the Sporades basin (Estournel et al., 2005). 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 favours strong heat losses and fast mixing of the waters, while the morphology reduces the interactions with the ‘open sea’. Thus, under favourable 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 to renew the North Aegean Trough’s deep water and especially the Sporades deep basin (Valioulis and Krestenitis, 1994). They can reach the southern boundary of the Outer TG, approximately 2 weeks after their formation in the Inner TG (Zervakis et al., 2005).

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 (Dodou et al., 2002). When winds are favourable 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 (Kourafalou, 2001). Downwelling may evolve periodically on a diurnal time scale, especially in summer at the northwestern Gulf (Figure 6b) during onshore southeasterly winds (Kourafalou, 2001; Savvidis et al., 2019); note that such winds are very frequent in TG (Figure 5a). The prevailing and usually stronger northwesterly winds (Figure 5a) may reduce temperature at the surface after the coastal upwelling of colder waters over the eastern coasts of the Central and Outer TG (Kourafalou and Tsiaras, 2007; Figure 6a).

3.2. Sediment Transport and Modern Sedimentation

3.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 TG (Figure 1b). Sediment and freshwater inflows have changed significantly during the last century due to direct and indirect human interventions. 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^oE, 40.59^oN; 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 (Kapsimalis et al., 2005). 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 (Poulos et al., 1994).

Historical data indicate significant accretion in the deltaic zone of Axios and Aliakmonas during the 20^th century, with net gains of 175 km² and 139 km² between 1850 and 1987 (Poulos et al., 1994). Kapsimalis et al. (2005) 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·106 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·106 m³ (28·106 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 of sediment supplies (by 2.5·106 m³/yr or 4·106 t/yr), leading to coastline retreat at 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 (Kapsimalis et al., 2005). 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. (2015) 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 (Kapsimalis et al., 2005; Poulos et al., 1994). Considering that these areas are low-lying, the risk of intensification of the coastal retreat and flooding under rising sea level in the future is very high and could entail severe environmental and socioeconomical impacts (Androulidakis et al., 2023b). Regarding 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 (Karymbalis et al., 2016). 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 (Karymbalis et al., 2016).

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) (Poulos et al., 1993). Axios, Aliakmonas and Pinios provide the bulk of freshwater and fine sediment discharge in the gulf, while Loudias River has a sporadic contribution (Karageorgis and Anagnostou, 2003; Kombiadou and Krestenitis, 2012). While historic data estimate the average total riverine sediment influx to Thermaikos at 3-5·106 t/y (Lykousis and Chronis, 1989; Poulos et al., 2000), by the end of the 20^th century the total supply of suspended solids by the rivers had decreased to around 0.63·106 t/y (10.7% from Axios, 1.1% from Aliakmonas, 1.0% from Loudias and 87.2% from Pinios; Karageorgis and Anagnostou, 2003).

3.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 (Kourafalou et al., 2004). 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 (Karageorgis and Anagnostou, 2003; Kombiadou and Krestenitis, 2012). Elevated suspended Particulate Matter Concentrations (PMCs) near the bed are attributed to associated current activity, either inhibiting settling or causing resuspension (Poulos et al., 2000).

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 projects EURECOMARGE (1 survey: June 1987; Durrieu de Madron et al., 1992), Metro-Med (4 surveys: May 1997, Jul. 1997, Feb. 1998 and Sep. 1998; Karageorgis and Anagnostou 2001; 2003) and INTERPOL (3 surveys: Sep. 2001, Oct. 2001 and Feb. 2002; Zervakis et al., 2005). More recent, monthly monitoring campaigns (Thermaikos 2004, 2006; Tsompanoglou et al., 2017), were restricted to the Inner and Central TG. Satellite imagery have also been used to sense suspended PMCs in the surface waters, using Landsat (Karageorgis et al., 2000; Balopoulos et al 1986) and MERIS (Monachou et al., 2014; Alexandridis et al, 2015). 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 (SNL and BNL, 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 resuspension of loosely consolidated sediments by near-bottom currents (Karageorgis and Anagnostou, 2003). 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 under non-stratified conditions, while under stratified conditions PMCs of INLs are significantly reduced. Winter conditions, with high influx of solids by the rivers and strong haloclines, are typically characterised 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 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 3.1.3), altering the BNL structure in its wake (Estournel et al., 2005; Tragou et al., 2005). 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 (Zervakis et al., 2005). 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 the occasional shift of the Sporades eddy from predominantly cyclonic in winter to anticyclonic in other seasons (Kontoyiannis et al., 2003; Kourafalou and Tsiaras, 2007), 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 by 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 from Cape Epanomi (dashed lines in Figure 7d).

3.2.3. Near-Bed Processes and Sedimentation Patterns

TG acts as a sediment trap, with the bulk of the riverine sediment (~90%) depositing within the bounds of the domain (Kombiadou and Krestenitis, 2012; Lykousis and Chronis, 1989). Sedimentation in the deltaic zone of 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 (Poulos et al., 1996). 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 (Karageorgis et al., 2005a; Zervakis et al., 2005). Numerical modeling and oceanographic observations indicate that nearshore currents are not able to resuspend sediment (Poulos et al., 1994), while wave-induced resuspension can be an important factor over low depths near the coast (Paphitis and Collins, 2005). The trawling activity taking place in the area has been shown to contribute to the BNL concentrations in the Outer TG (Tragou et al., 2005; Zervakis et al., 2005). In fact, trawling was estimated to increase BNL concentrations over the shelf by 2 to 3 times (Price et al., 2005). Similarly, using numerical modeling, Kombiadou and Krestenitis (2011) estimated that the sediment entering suspension due to trawling over the trawling season can supersede the annual influx by the rivers. Resuspended particles are removed rapidly from the system by sinking and/or horizontal advective processes (Muir et al., 2005). Overall, mechanically eroded sediment moves in the vicinity of the bed (3-6 m), contributing to the BNLs in the area, and resettle onto the bed shortly after (1-5 days; Kombiadou and Krestenitis, 2011).

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 (Lykousis et al., 1981; Lykousis and Chronis, 1989; Poulos et al., 2000). Sedimentation over the western shelf is dominated by the riverine water and sediment outflows, while 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 (Poulos et al., 2000). These coarser sediments are interpreted as relict sands that remain uncovered (or partly mixed) by modern sedimentation, due to the predominant cyclonic transport (Karageorgis and Anagnostou, 2001). Sedimentation rates near the river mouths were measured at 8.8 mm/yr for Axios and 1.8 mm/yr for (Karageorgis et al., 2005a). Corresponding modeled values are 18, 7 and 3 mm/yr near the mouths of Aliakmonas, Axios and Pinios, respectively (Kombiadou and Krestenitis, 2012). 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.

3.3. Marine Pollution and Biogeochemical Properties

TG is influenced by human activities (sewage and industrial discharges, agricultural runoff, shipping, airborne litter, hydrocarbon refueling, and harbour operations) while the five large rivers discharge substantial loads of nutrients, heavy metals and other compounds, after drainage resulting from human activities within their watersheds. Axios River has been recognized as the major supplier of suspended particulate matter even during the dry summer season (Tsompanoglou et al., 2017; see Section 3.2) and a source of heavy metals to the coastal environment of the gulf (Karageorgis et al., 2003). 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 TG. In Thessaloniki Bay and 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).

3.3.1. Heavy Metals

Metal enrichment is mainly detected in the northernmost part of the TG (Inner TG), and to a lesser extent close to the river mouths; then, the heavy-metal content gradually decreases towards the open sea (Karageorgis et al., 2005b). 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. Vasilikiotis et al., 1982; Fytianos and Vasilikiotis, 1983; Voutsinou-Taliadouri and Satsmadjis, 1983; Voutsinou-Taliadouri and Varnavas, 1995; Zabetoglou et al., 2002; Christophoridis et al., 2009; Violintzis et al., 2009; Kapsimalis et al., 2010; Christophoridis et al., 2019).

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, 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 3.2.3) controls the sediments’ geochemistry causing a substantial dilution of heavy metals concentration (Karageorgis et al., 2005b).

The pollution status of marine sediments in the coastal area of Thessaloniki Bay (Figure 1b) was evaluated by employing enrichment factors and comparing with sediment quality guidelines (Violintzis et al., 2009; Christophoridis et al., 2009; 2019). For Zn, Cu, Pb, As and Ag significant contribution from human-related sources was found, particularly in the inner part of the Bay (Violintzis et al., 2009). 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 (Christophoridis et al., 2019). 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 (Violintzis et al., 2009).

The dissolved fraction of heavy metals in 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 (Dassenakis et al., 2005). Zeri and Hatzianestis (2005) measured the dissolved Cu and Ni and reported that the maximum concentrations of both metals were encountered in front of the Axios and Aliakmonas rivers deltas. Stronger binding of both Cu and Ni with organic ligands of terrestrial origin was also observed than 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 at 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 (Tzempelikou et al., 2021).

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 PMC 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 (Tsompanoglou et al., 2017). 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 (Price et al., 2005). 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 Aliakmon 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, derived from partially treated domestic and industrial effluents (Tsompanoglou et al., 2017). Similarly, Pb and Cd were attributed to anthropogenic sources, derived from the rivers and the city of Thessaloniki (Price et al., 2005).

3.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 the 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 (Hatzianestis et al., 2000; 2001). These results confirm the findings of previous relevant studies (Fytianos et al., 1985; Larsen and Fytianos, 1989). 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 one, suggesting the absence of “fresh” DDT inputs in the gulf. The vertical distribution of ΣDDTs 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 Inner TG (Hatzianestis et al., 2000; 2001). Despite the ATR ban within the European Union in 2004, ATR and its TPs were still detected in seawater and sediments of TG collected during sampling campaigns in 2010 and 2011 (Nödler et al., 2013), although at relatively lower concentrations than in the period before the ban of ATR use (Albanis et al., 1992; Readman et al., 1993). 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 Gallikos river (Nödler et al., 2013).

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 (Dimou et al., 2006). Arditsoglou and Voutsa (2012) also conducted a study to determine the presence of phenolic and steroid Endocrine-Disrupting Compounds (EDCs), in all elements of the marine environment of Inner and Central TG (seawater, suspended particulate matter, sediments, 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 (Pothitou and Voutsa, 2008; Arditsoglou and Voutsa, 2010). 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 in 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 (Pantelaki and Voutsa, 2021).

Hatzianestis et al. (2000; 2001) 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 ΣPCBs 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 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 (Hatzianestis et al., 2000; 2001). 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 (Hatzianestis et al., 2000). Polycyclic Aromatic Hydrocarbons (PAH) 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 PAH (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 (Hatzianestis et al., 2001).

The concentrations of different organic contaminants (i.e. PAHs, PCBs, organochlorine pesticides and polybrominated diphenylethers) were measured at three sites of Inner and Central TG not only in surface seawater but also in the air in summer (July) 2012 (Lammel et al., 2015). 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 by a factor of 2–10 higher than in the Cretan sea (Southern Aegean). Furthermore, the spatial variation of seawater pollution across the three sites was remarkable with the highest PAHs levels recorded at a residential coastal site (Michaniona), the highest levels of organochlorine pesticides measured at Loudias River estuary implying influence from agricultural activities, whereas the highest levels of PCBs and polybrominated diphenylethers were recorded off-shore 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.

3.3.3. Nutrient Pressures

TG is one of the most disturbed Hellenic coastal areas regarding anthropogenic nutrient enrichment (Pagou, 2005). 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 TG that exhibit significant seasonal and annual variations (e.g., Samanidou et al. 1987; Balopoulos and Friligos, 1994; Pavlidou, 2012; Zachioti et al., 2022). More specifically the Inner and Central TG receives 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., Friligos and Koussouris, 1984; Friligos et al., 1997). From 1992 and on, Level I wastewater treatment was applied, while since 2001, Level II wastewater treatment is implemented (Krestenitis et al., 2012), with a 92% performance degree, in order to reduce the contribution of urban effluents to eutrophication events (e.g., Moncheva et al., 2001; see Section 3.4).

In general, the nutrients 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 the aquacultures, as well as in Thessaloniki Bay (port and sewage outfall area) (Pavlidou et al., 2005; Pagou, 2005). A timeseries (1995-2007) analysis of nutrients 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 (Zachioti et al., 2022). 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 of 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 (Androulidakis et al., 2021). 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 3.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 3.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 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 (Skoulikidis, 2009). In view of these differences, Aliakmon 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 (Skoulikidis, 2009). 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 since then gradually diminished (Skoulikidis, 2009; Tsiaras et al., 2014). 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 (Karageorgis et al., 2005c; Skoulikidis, 2009; Tsiaras et al., 2014). Accordingly, phosphate loads are unbalanced in relation to nitrate loads resulting to 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 nutrients 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, whereas 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 (Pagou et al., 2003; Pagou, 2005; Pavlidou, 2012). The above results were confirmed by Zachioti et al. (2022) 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. Nikolaidis et al., 2006; Pavlidou, 2012; Pavlidou et al., 2015; Simboura et al., 2016). 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 (Pavlidou, 2012). 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 (Genitsaris et al. 2019). Similarly, the N:P ratio was high (19.5) during the years 2020 -2022 (Kourkoutmani et al., 2023) indicating nitrogen pollution.

3.3.4. Other Biogeochemical Characteristics

Particulate organic matter (POM) in 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 3.4.2). The distribution of particulate nitrogen and phosphorus was similar to that of POC (Tsompanoglou et al., 2017). 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 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 resulting 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 (Strogyloudi et al., 2023), low ratios of POC:chl-a were recorded in the mussel-farms area close to the rivers estuaries the high participation of chl-a to POC that 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% (Zervoudaki et al., 2014).

Inner TG characterized by high loads of organic and mineral particles, as well as inorganic nitrogen compounds, provided ideal conditions either for 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 (Marty et al., 2001 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 CH4), 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 Inner TG could act as a significant source of atmospheric N₂O and CH₄ (Marty et al., 2001).

In a study performed during May 1997 data of the CO₂-carbonate system of the Inner TG was obtained for the first time (Krasakopoulou et al., 2006). 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 Utilisation (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₂ (fCO2), 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^-2d^-1, depending on which formula for the gas transfer velocity was used, and in accordance to recent reports regarding other temperate continental shelves.

3.4. Biological Elements of the Water Column - Plankton

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

3.4.1. Bacterioplankton

Although the view on ecology of marine plankton has changed remarkably since the 1980’s (Azam et al., 1983, Fenchel, 1988) 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 (Mihalatou and Moustaka-Gouni, 2002) indicated high numbers of bacteria (from 0.9 to 6.8·10⁶ cells mL-1), reaching maximums 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^-3h^-1) indicating strong trophic coupling of bacterioplankton and phytoplankton. Another work, recently published by Genitsaris et al. (2023), aimed at the response of bacterioplankton communities from the Inner and the Outer TG to PAHs’ exposure under experimental conditions (mesocosms). No negative effects of the contaminants were apparent in the experiments and 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 TG. The Outer TG contained bacterioplankton communities with higher resistance sustaining bottom-up ecosystem stability.

3.4.2. Phytoplankton

Several studies on phytoplankton in the TG have been carried out, although long timeseries 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 (Nikolaides and Moustaka-Gouni, 1990). The phytoplankton seasonal succession was dominated by Leptocylindrus, Chaetoceros, Skeletonema, Pseudo-nitzschia, Thalassiosira, Cylindrotheca, Cerataulina (diatoms) and Gymnodinium, Prorocentrum, Scrippsiella, Gonyaulax, Heterocapsa (dinoflagellates) species forming phytoplankton biomass peaks up to 6.7 mgL^-1. Focusing on spatial phytoplankton distribution in the entire TG, Gotsis-Skretas and Friligos (1990) showed a spatial variation of phytoplankton along the pollution gradient with maximal total abundance of 7,600 cells mL^-1. 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 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 (Spatharis et al., 2019). Androulidakis et al. (2021) indicated that the phytoplankton species of the different habitat communities of the Inner, Central and Outer TG were highly connected via circulation (see Section 3.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 3.3).

Genitsaris et al. (2020), 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 (Genitsaris et al., 2019), Dinophysis species was found responsible for mussel intoxication in aquaculture areas (Koukaras and Nikolaidis, 2004; Varkitzi et al., 2013). Dinophysis bloom occurrence had been documented several decades ago (Athanassopoulos, 1931; Gotsis-Skretas and Friligos, 1990; Nikolaides and Moustaka-Gouni, 1990; Nikolaidis and Evagelopoulos, 1997). Cysts of potentially toxic dinoflagellate species, such as Alexandrium cf. tamarense, A. cf. affine, A. cf. minutum, as well as Gymnodinium catenatum, were detected in the sediment of the Central and Outer TG (Giannakourou et al., 2005).

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. (2018). 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 tides 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 at the urban marine environment of the Inner TG was conducted by Genitsaris et al. (2019). 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 (Genitsaris et al., 2019), attributed to the same mucilage-producing diatoms, together with the mucilaginous haptophyte Phaeocystis sp. and the dinoflagellate G. fragilis.

The trophic conditions of TG were assessed using the multiparametric Eutrophication index (E.I.; Primpas et al., 2010) by some studies. Based on nutrients and chlorophyll-a (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” till 2000, was improved and classified as “Moderate” in 2007 (Zachioti et al., 2022). 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 (Pavlidou et al., 2015).

Figure 8. Eutrophication events in Thermaikos Gulf (TG) in July 2017 (Photo courtesy: Eurokinissi.gr), February 2017 (Photo courtesy: Thes.gr), March 2017, and December 2017. The plankton species associate with the red tide events in March and December 2017 are also shown.

Figure 8. Eutrophication events in Thermaikos Gulf (TG) in July 2017 (Photo courtesy: Eurokinissi.gr), February 2017 (Photo courtesy: Thes.gr), March 2017, and December 2017. The plankton species associate with the red tide events in March and December 2017 are also shown.

Several studies in the prism of monitoring activities in the framework of the implementation of WFD 2000/60/EC (WFD) and Marine Strategy Framework Directive 2008/56/EU (MSFD), evaluated the Inner TG, and particularly the Thessaloniki Bay, with “Moderate” or “Poor” status (Pavlidou et al., 2012; 2015; Simboura et al., 2016). Recently, a study on eutrophication assessment of TG capitalized on current advances in remote sensing, and employed satellite image analysis in order to evaluate the spatial distribution of chl-a levels (Androulidakis et al. 2021). In order to have 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 the 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, for the formation of these eutrophication phenomena. For instance, they found that during an annual cycle, chl-a peaked at 45 mgm^-3 in December 2017 (Figure 8) under the prevalence of southerly winds (weak renewal; Section 3.1), when an extensive red tide caused by the photosynthetic ciliate Mesodinium rubrum expanded throughout the Inner and Central TG (Androulidakis et al., 2021). On the contrary, low chl-a levels were detected during periods with northerly winds that enhanced the renewal of the enclosed basins.

Figure 9. Micrographs of common planktonic organisms of TG as seen by microscopy. Epifluorescence micrographs were taken by UV excitation for (A) DAPI-fluorescence cells or by green excitation (B) and blue excitation for (C) Chl-a red auto-fluorescence. Light - phase contrast (D, E, F, G) micrographs. (A) Heterotrophic bacteria. (B) Pico- cyanobacteria. (C) Pico-chlorophytes. (D) The diatoms Chaetoceros sp. and Leptocyllindrus sp. along with the dinoflagellate Ceratium cf. horridum. (E) The red-tide forming ciliate Mesodinium rubrum. (F) The bloom forming dinoflagellate Karenia with other plankton species. (G) Cells of the red-tide forming heterotrophic dinoflagellate Noctiluca scintillans along with the pelagic tunicate Oikopleura dioica. Scale bars in μm.

Figure 9. Micrographs of common planktonic organisms of TG as seen by microscopy. Epifluorescence micrographs were taken by UV excitation for (A) DAPI-fluorescence cells or by green excitation (B) and blue excitation for (C) Chl-a red auto-fluorescence. Light - phase contrast (D, E, F, G) micrographs. (A) Heterotrophic bacteria. (B) Pico- cyanobacteria. (C) Pico-chlorophytes. (D) The diatoms Chaetoceros sp. and Leptocyllindrus sp. along with the dinoflagellate Ceratium cf. horridum. (E) The red-tide forming ciliate Mesodinium rubrum. (F) The bloom forming dinoflagellate Karenia with other plankton species. (G) Cells of the red-tide forming heterotrophic dinoflagellate Noctiluca scintillans along with the pelagic tunicate Oikopleura dioica. Scale bars in μm.

3.4.3. Protozooplankton

Protozooplankton community is underappreciated in plankton studies in the TG. The early study by Mihalatou and Moustaka-Gouni (2002) focused on protozooplankton components (heterotrophic nanoflagellates, heterotrophic dinoflagellates and ciliates) as grazers of bacteria. The ciliate M. rubrum as 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) even the red-tide forming omnivorous feeder N. scintillans, have been considered members of phytoplankton community in earlier plankton studies (Nikolaides and Moustaka-Gouni, 1990; Gotsis-Skretas and Friligos, 1990). Investigating the underexplored diversity of marine unicellular eukaryotes, Genitsaris et al. (2020) identified phytoplankton heterotrophic unicellular eukaryotes frequently observed in 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, that had previously been taxonomically characterized by 18S rRNA gene amplicon sequencing, showed a diverse range of trophic preference of picograzers, nanograzers, micrograzers, phagotrophs, decomposers and parasites (Genitsaris et al., 2020).

3.4.4. Metazooplankton

The metazooplankton community of TG is also poorly known resulting to 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 (Siokou-Frangou and Papathanasiou, 1991) 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, 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 (Alvanou, 1999). The distribution of mesozooplankton resting eggs was studied in the sediments of TG and indicated higher abundance of eggs at locations close to the river mouths, with water column rich in zooplankters (Siokou-Frangou et al., 2005). Recent research on metazooplankton reported the presence of the non-indigenous calanoid copepod Pseudodiaptomus marinus (Kourkoutmani and Michaloudi, 2022). 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 (Kourkoutmani et al., 2023).

3.4.5. Planktonic Response to Climate Change and Pollutants

During the last decade, Inner TG unicellular plankton community has been used in mesocosm experiments to explore the species response to multiple variables of climate change (Stefanidou et al., 2018a), and phytoplankton response to pollutants discharged in the seawater from exhaust gas cleaning systems on ships, known as scrubbers (Genitsaris et al., 2023). In regard to the climate change variables, the impacts of elevated temperature (see Section 4.2) in the form of an intense heat shock, together with changes in salinity, were examined on the natural phytoplankton community of TG and demonstrated adverse effects on species diversity, resource use efficiency and productivity (Stefanidou et al., 2018a). 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 (Stefanidou et al., 2018b), 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 phytoplankton community caused by the elevated temperature and increased salinity, which was also confirmed by measurements of physical properties by Androulidakis et al. (2023). 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; 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, and particularly vanadium and nickel, was 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 open sea (Genitsaris et al., 2023).

3.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 TG (see Section 2). 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 Aristotle University of Thessaloniki (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 Hellenic Centre of Marine Research (HCMR). More specifically, the first research effort about benthic fauna of Thermaikos Gulf was about the diversity of extant mollusks (Sakellariou, 1957). 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) has been carried out by AUTh.

3.5.1. Benthic Communities

The seabed of TG is dominated by various soft substrata types, which are locally covered by seagrass meadows of the species Posidonia oceanica and Cymodocea nodosa (Haritonidis et al., 1990; Panayotidis et al., 2022). In the unvegetated sedimentary bottom, muddy sands constitute the predominate type (Poulos et al., 2000). The proportion of sand increases eastward and to shallower depth, whereas the amount of silt and clay balance 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 4, 20 types of benthic communities occur in 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 is covered by periwinkles, limpets and barnacles in the supralittoral and mediolittoral zone, respectively (Kitsos, 2003), whereas in the sublittoral zone four different facies of the photophilic algae community develop, depending on prevailing environmental factors (Chintiroglou et al., 2004; Antoniadou and Chintiroglou, 2009; Antoniadou et al., 2011; Ganias et al., 2023).

Seven out of the 20 biocenoses found in natural habitats of TG are considered as 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 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 (Otero et al., 2017).

3.5.2. Benthic Diversity and Species of Community Interest

Overall, 1,368 benthic invertebrate species have been reported from the faunistic studies that have been carried out in TG since the 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 3.5.1) and to the extensive taxonomic studies focusing in particular to the phylum Mollusca (Manousis et al., 2010; 2012; Manousis and Galinou-Mitsoudi, 2014; 2022). Among the reported invertebrate species, some are commercially exploited, such as the clams Venus verrucosa, Callista chione, the scallops Pecten jacobaeus and Flexopecten spp., the prawns Penaus kerathurus, 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 (Damianidis et al., 2016; Kevrekidis and Thesslaou-Legaki, 2011; Kevrekidis and Antoniadou, 2018; Vafidis and Antoniadou, 2023), raising concerns on the future viability of natural stocks and habitats.

Other species, however, face direct threats 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 in by Vafidis and Koukouras (1991), lives exclusively in TG, where it forms sparse populations (0.5 to 4 colonies/20 km) on the silty sediments of the Central TG. The species habitat is heavily disturbed by bottom trawling, urbanization, and pollution (Otero et al., 2017). There are no active protection measures, so far, for the sea pen, whose population remains out of the established Special Areas of Conservation (SACs) of TG. The pillow coral C. caespitosa forms moderately dense populations on coarse sedimentary and rocky habitats, especially along the eastern coast of the gulf (Ganias et al., 2023). Climate change and fisheries mainly threaten the species (Otero et al., 2017), whereas in areas influenced by sewage treatment plants, siltation constitute 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 biogenic habitat that constitutes an important biodiversity reservoir (Antoniadou et al., 2023). Although there are no specific conservation measures for the species, a part of its population in Thermaikos is protected within the SAC GR1220005. Thermaikos Gulf hosted the denser population of the fan-mussel P. nobilis in the Greek Seas (Vafidis et al., 2014) with densities reaching 1.04 individuals/m2 (Galinou-Mitsoudi et al., 2006). Unfortunately, the species have been affected by haplosporidia parasites and mycobecteria that caused Mass Mortality Events (MMEs) that totally devastated the population (Zotou et al., 2020). Until nowadays, only dead eroded shells heavily colonized by epiphytes and sessile invertebrates (Figure 12c) may be found in TG, including the designated SACs in the region.

3.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 to develop appropriate indices to assess the ecological quality status and the seabed integrity. This is due to their widespread and abound presence in both species’ richness and number of individuals within marine habitats, also to the rather limited natural variability in the synthesis of their biotic communities, and lastly to their fundamental life traits. They 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 are able to respond to various stress factors, such as eutrophication, organic pollution, mechanic disturbance of substrata, and heavy metal load. They also play a main role to 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 (2002) has been successfully inter-calibrated between European countries to be used over Mediterranean coastal waters and the M-AMBI index (Borja et al., 2004; Muxica et al., 2007) for transitional waters (Simboura and Reizopoulou 2008).

The National Monitoring Water Network (NMWN 5) estimates the ecological quality of coastal and transitional waters in 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 along three periods: an initial reference period (2000-2008), a first (2012-2014) and a second (2018-2022) implementation period. 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 (YPEKA, 2012; 2017; 2023). In the Outer TG, the ecological quality state -originally assigned as good (YPEKA, 2012)- has degraded to moderate (YPEKA 2017; 2023); however, the latest assessment has not been based on biological samplings, but to an extrapolation of the previous evaluation. Considering the transitional waters, the state has been steadily assessed as moderate (YPEKA, 2017; 2023).

Based on the few studies that have applied the BENTIX index in TG (Antoniadou et al., 2004; Antoniadou and Chintiroglou, 2005; Chintiroglou et al., 2006; 2007; Simboura and Reizopoulou, 2007; Antoniadou and Chintiroglou 2009; Simboura et al., 2014; HCMR, 2015; AUTh 2018a; 2018b; 2019; 2020a; 2020b; 2023), 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 decreased state in the Inner TG together with an improved state in the Outer TG emerges, whereas in the Central TG, the eastern shores appears in 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 (AUTh, 2018a; 2018b; 2019; 2020a) revealed a gradual improvement in the quality of the western shores but a deterioration in the eastern ones, where 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 (AUTh, 2020b; 2023) and good quality along the western coasts of the Outer TG (AUTh, 2023).

4. Discussion and Concluding Remarks

Thermaikos Gulf (TG) is a productive, natural, coastal ecosystem with extensive estuarine systems that is highly impacted by anthropogenic pressures. A severe deterioration of the quality has been documented between 1970s-1990s; thereafter a steady trend of recovery has been identified, especially after the operation of the two main WTPs in 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).

4.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 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 to 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%), mainly conducted during the 1990s. A clear reduction in research is 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 (fiscal austerity), with strong implications for academia and, in extension, for environmental research and monitoring activities. The fiscal austerity years after 2009 are characterized by brain drain of young scientists (Pellicia, 2013; Azaria et al., 2019), budget cuts in applied research and especially environmental monitoring (Centre for the Development of Educational Policy, 2015), decreasing human resources in research institutions and universities (Koulouris et al., 2014; Chalari, 2016), and reducing funding for the management of environmentally protected areas (e.g., Natura 2000; Paliogiannis et al., 2018). These factors, combined with the complex bureaucracy imposed by national and European Union funded research projects (Skoulikidis et al., 2021; Ladi et al., 2022), the continuously changing Greek legislation regarding the management of funds (Paliogiannis et al., 2019), the over-taxation of those funds and the endless delays, caused a tremendous burden on research development (Koulouris et al., 2014). The lack of research efforts is even more challenging, given the fact that the fiscal austerity measures has further intensified the anthropogenic pressures and environmental impacts, especially in urban coastal areas such as the TG (Petrakakis et al., 2013). The misuse of biotic resources in protected areas reached unprecedented levels in Greece during the period of the economic recession (Troumbis and Zevgolis, 2020). 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, yet usually with adverse effects on the coastal environment) after more than a decade of recession. The COVID-19 pandemic after 2020 and the successive lockdowns have impacted the research activities even more (Jackman et al., 2022), especially environmental monitoring that requires group cooperation and higher funding; safety restrictions and different funding priorities have created additional gaps in TG’s monitoring.

Currently, 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. A tendency for improvement in the ecological quality of the Inner and Central TG’s coastal waters was documented in the NMWN results. This outcome seems odd, as the Outer TG, influenced by waters of the open Aegean Sea that are classified as excellent (see Section 3.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 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 of observational stations, both temporarily 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. Centralized management, offering free (open-) access to the relevant data is minimal. The lack of measurements of the freshwater supply and quality leads to uncertainties and 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 facets of the TG environmental monitoring, the absence of a supervisory public entity, capable of organizing and supervising all coastal research and protection endeavours, 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; Androulidakis et al., 2023a) 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 by authorities, with only occasional and inefficient initiatives for protective action and integrated coastal management during the last decades 6. Although coastal erosion and shoreline retreat is a major problem, especially along the eastern coasts of the Central Gulf of Thessaloniki (e.g., Peraia, 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 retreat, based on the Normalized Difference Water Index (NDWI) derived from two Sentinel-2 L2A images (10 m resolution; Copernicus Dataspace Browser 7) in 2019 and 2024 is presented in Figure 13. The estimated reduction of the delta area is approximately 11 km².

Another human-induced environmental stressor, that severely impacts both biodiversity and human health, is the high concentration of marine plastic litter in the TG (Politikos et al., 2017; Kermenidou et al., 2023). 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 the Descriptor 10 of MSFD. Microplastic pollution in the TG is mainly controlled by river discharges (Zeri et al., 2020; Kaandorp et al., 2020) which carry significant amounts of microplastics (40-50 times higher than the maximum ocean surface’s concentrations; Meijer et al., 2020); monitoring and mitigation of this stressor in the rivers are also very limited.

4.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 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 Sen’s Slope (0.52°C/decade; Figure 14a), higher than the general trend of the broader Aegean, Ionian and Cretan Seas (0.49°C/decade; Androulidakis and Krestenitis, 2022). The TG is characterized as a “hot spot” area for MHWs, with events that do not only occur near the sea surface, but extend through the entire water column, especially in the shallower Inner TG (Androulidakis et al., 2023a). A statistically significant (p_value<0.0001) SST shift of around 1.2°C was computed after 2006, based on the Pettitt Homogeneity test (Pettitt, 1979), 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 (99^th 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; Androulidakis and Krestenitis, 2022). The related trend slope of the maximum levels is even steeper than the annual mean (0.67^oC/decade), increasing by approximately 2.5^oC during the 1982-2023 period. The respective MHWs were computed based on the satellite-derived SSTs and the methodology, introduced by Hobday et al. (2016); the annual MHW cumulative intensity (^oCxdays), averaged over the entire TG, is presented in Figure 14b. The events were defined using the 90^th percentile climatology (1982-2023) as threshold baseline, and a 5-days minimum duration for each event. Very high intensities successively occurred after 2018 confirming the intensification of the warming impact in the TG. Intensity peaks were computed for 2012 (related mainly to winter MWHs as derived from the high mean SST and relatively low 99^th percentiles; Figure 14a) and in 2021. Especially in 2021, the TG’s MWHs had significant biochemical implications, affecting the marine environment and imposing mass mortality events of prominent biota with considerable impacts (e.g., on aquaculture, fishing, tourism, human health, etc.; Androulidakis and Krestenitis, 2022; Lattos et al., 2022; Antoniadou et al. 2023). The importance of MHWs raises the need of accurate prediction tools for both short-term (e.g. a few days based on operational hydrodynamic modeling) and interim-future periods (e.g., weeks to months forecasts based on artificial intelligence techniques; Krestenitis et al., 2024).

An increase has also been observed in salinity that reached the highest measured levels after 2020, associated with elevated temperatures and possible changes in precipitation rates and river discharge outflows (Androulidakis et al., 2023a). Temperature increases together with changes in salinity may also affect species diversity, resource use efficiency and productivity (Stefanidou et al., 2018a). Experiments in the TG (Stefanidou et al., 2018a) showed how changes in the species inventory driven by a 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 21^th century (Darmaraki et al, 2019).

Figure 14. Annual (a) mean and 99^th Percentile Sea Surface Temperature (SST; ^oC) and (b) cumulative intensity of Marine Heatwaves (MHWs), derived from the satellite-derived daily fields and averaged over the TG during 1982-2023. The linear trends, the Sen’s Slopes, and the μ1 and μ2 levels of the Pettit homogeneity test for the mean values are also shown.

Figure 14. Annual (a) mean and 99^th Percentile Sea Surface Temperature (SST; ^oC) and (b) cumulative intensity of Marine Heatwaves (MHWs), derived from the satellite-derived daily fields and averaged over the TG during 1982-2023. The linear trends, the Sen’s Slopes, and the μ1 and μ2 levels of the Pettit homogeneity test for the mean values are also shown.

Temperature also greatly affects the equilibrium constants of the CO₂-carbonate system and particularly the solubility coefficient of CO₂; warming and increasing salinity decrease the solubility of CO₂ (e.g., for a temperature increase of 1°C, pCO2 rises by ~4%, Takahashi et al., 1993), which may induce CO₂ oversaturation in many cases and change of the character of the area, acting finally 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₂ emissions scenarios (Solidoro et al., 2022; Reale et al., 2022). 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 unfavourable and the TG is a coastal ecosystem that is already subject to climate change and 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 (Poulos et al., 2000). To 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, especially over the western coastal area of the TG. 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 to increasing flooding frequencies over the low-lying coastal areas (Androulidakis et al., 2023b). 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 analysing long-term future projections of physicochemical and biological features, based on diverse climatic scenarios combining atmospheric input from ensemble modeling, are absent for the TG, while the research over past climate-scale periods forming the present day situation corresponds to only the 2% of the entire research activity in the gulf.

4.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 favourable pre-condition for the formation of eutrophication events. Although the nutrient pollution of TG was reduced during the last two decades, after the strong reduction of 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., 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 levels set by the European Union for bathing waters. 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, can trigger seawater acidification, and enhance coastal water pollution locally. These effluent discharges are usually warmer and more acidic than the surrounding water, and therefore increasing acidity can 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 concertation 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 in both taxonomic and ecological level. Several endemic, protected, threatened and fishery exploited species have dense and widespread populations in the Gulf, as well as invasive ones, which seem to spread via navigation. The area hosts a unique species of a 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 the western shores of the Central TG, whereas the relevant state has improved in the eastern shore and the Outer TG. However, this pattern seems to be currently reversed, as the ecological quality is getting worse in the eastern and better in 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 timeseries of available data, all relevant studies reveal organic pollution as the main threat for 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 for benthic communities are related to excessive fishing and especially bottom trawling.

In conclusion, 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.