A High-Resolution Numerical Model of the North Aegean Sea Aimed at Climatological Studies

: A new, high-resolution model for the northern part of the Aegean Sea, aimed primarily at climatological research (relaxation and data assimilation-free climate simulations), is hereby presented, along with the results of a 28-year-long simulation covering the period from 1986 to 2013. The model applied is the Regional Ocean Modelling System (ROMS). A signiﬁcant improvement over previous models of the Aegean introduced in this work is the replacement of parameterizations of the Dardanelles exchange by a fully three-dimensional simulation of the ﬂow in the Strait. The incorporation of part of the Marmara Sea in the model domain enables the interaction with other regional climate simulations, thus allowing climatic variability of the exchange of the Mediterranean and Black Seas. An extensive validation is carried out comparing the model output with all the available observations from several different platforms, i.e., satellite sea surface temperature and height, T/S proﬁles from R/V ships, and HF radar surface currents velocity. We focus on the model’s ability to reproduce, to some extent, the distinct thermohaline features and circulation patterns that characterize this speciﬁc area of the Mediterranean Sea. Our ﬁndings, after comparing simulation results with all the available observations, revealed the model’s sufﬁciency to simulate very adequately the complex hydrology of the North Aegean Sea, and the model’s ability to reproduce incidents of deep-water formation that took place in the region in previous decades during the Eastern Mediterranean Transient (EMT).


Introduction
The North Aegean Sea is one of the most intriguing seas of the Mediterranean, being the first receptor of waters from the Black Sea through the so-called Turkish Strait System (TSS), comprising the Dardanelles and Bosphorus Straits and the Marmara Sea). While it is one of the northernmost extremities of the basin, and thus, an alternative source of Eastern Mediterranean Deep Waters (EMDW), along with the Adriatic Sea [1,2], the inflow of very light, low-salinity waters from the Black Sea through the TSS plays a critical role in controlling deep-water formation over the basin for the following reasons: (a) the lateral buoyancy inflow through the Dardanelles reduces surface density and increases stratification, hindering the formation of very dense waters, and (b) the Black Sea Waters (BSW) form a thin surface layer, colder than its underlying waters of Aegean/Levantine origin, a condition minimizing heat losses to the atmosphere in the winter [3]. This critical role of the Dardanelles outflow in controlling overturning processes in the North Aegean was proposed after it was recognised that the massive Eastern Mediterranean Transient (EMT) was triggered in 1987 in this basin [2,4], although the involved processes were not fully understood at the time. The above event shook the until-then established perception of the Eastern Mediterranean, overturning circulation and demonstrating that the Aegean Sea has the capacity to become a major dense-water formation site and an alternative producer of Eastern Mediterranean Deep Water (EMDW) [5].
The inflow of Black Sea waters in the North Aegean plays a crucial role in controlling not only overturning processes, but also circulation [6], as well as ecosystem functioning and productivity [7][8][9] of the basin, which also explains the fact that most Greek coastal fishery takes place in the North Aegean [10]. Thus, the development of a numerical model suitable for climatological studies of the North Aegean would require a simulation of the interannual variability of the Dardanelles inflow of BSW. Until now, most simulations of the circulation of the Aegean Sea have been aimed mainly for operational oceanography applications [11][12][13][14], for which parameterizations of the seasonal cycle of the Dardanelles exchange [15] are considered adequate.
More recent works have treated the Dardanelles exchange as an open boundary exhibiting only seasonal variability of the exchanged volume rates [16], while the temperature characteristics of the inflowing BSW have been allowed to vary interannually based on observational data [17]. In an effort to address interannual variability of the Dardanelles outflow, Androulidakis et al. [18] have exploited the water budget of the Black Sea throughout the period 2002-2003, while in a more recent work, Mavropoulou et al. [19] showed that the interannual variability of the Dardanelles outflow plays a little role in the general circulation of the North Aegean, although it can have an important impact in the water column structure.
The present work extends the domain of the numerical simulation of the North Aegean Sea into the Marmara Sea to enable the application of open boundary conditions there, obtained from any regional-scale simulations available for the period of interest. In the present application of the model, simulating the period 1986-2013, we have used results from a TSS analytical model covering that period [20].

Model Description
We have employed the Regional Ocean Modelling System (ROMS), a free-surface, hydrostatic, primitive equation model, widely used by the scientific community for a diverse range of applications [21,22]. A rectangular grid of 1.0 km resolution with 31 vertical sigma levels covering the North Aegean Sea area-black square area ( Figure 1)-was developed using General Bathymetric Charts of the Ocean (GEBCO) 2014 bathymetry. To ensure the hydrostatic consistency of the grid (minimizing perturbation), we have used a linear programming method [23] to smooth bathymetry and set the ROMS' bathymetric smoothing parameters rx0 [24] and rx1 [25] equal to 0.125 and 6.91, respectively.
For momentum advection, we have used the 3rd-order upstream split scheme [26][27][28] as implemented in the ROMS model. This choice came after several sensitivity test runs mainly focused on the Turkish Straight System (TSS). The latter scheme gave the smallest Bias and Root-Mean-Square-Error (RMSE) and mainly a more realistic representation of the steep T/S distribution in the vicinity of the Dardanelles exit, compared to the default 3rd-order upstream Bias scheme. For tracer advection, the Multidimensional Positive Definite Advection Transport Algorithm (MPDATA) scheme [29], with a weak harmonic operator along constant geopotential surfaces, was used. The vertical mixing scheme was the classical Mellor, Yamada 2.5 [30], using the default background values for eddy viscosity and diffusivity. For the bottom stresses, a quadratic bottom drag scheme was applied using the model default value (of 3 × 10 −3 ). The selection for bottom-stress type and its value came after several sensitivity-test runs and the previous combination gave the slightly smallest Root Mean Square Error (RMSE) for sea surface height when compared with altimetric data for sea level anomaly. Initial conditions for temperature and salinity came from a Mediterranean Sea hydrographic atlas [31]. The T/S data combined with ERA5 [32] winds were used to estimate the initial values for the baroclinic geostrophic and Ekman velocity components, setting a level of no motion at 600 m depth to determine the total velocity. Boundary conditions for the south side of the computational grid were provided by a basin scale model of the Mediterranean Sea distributed from CMEMS [33] on a daily resolution. Data for the eastern boundary, located at the western extremity of Marmara Sea, were derived from the work of V. Maderich [19], providing temperature, salinity, sea surface height and volume flux time series for both inflow and outflow layers on a daily basis. At all open boundaries, mixed radiation-nudging conditions [34]-were used for 3D fields (T, S, u, v), Chapman [35] for free surface (ζ) and Flather [36] for barotropic velocities (ubar, vbar). Riverine inputs from the 8 larger rivers outflowing in the Aegean were provided in the form of a climatological time series from the SMHI e-hype hydrological model [37], covering the period from 1980 to 2010. Atmospheric forcing of 1 4 -degree spatial resolution and 1 h temporal resolution came from the European Centre for Medium-range Weather Forecast (ECMWF) ERA5 reanalysis dataset and was incorporated in the model using the COARE bulk formulae [38].
The atmospheric forcing used consisted of wind velocity components at 10 m height, 2 m height air temperature, dewpoint temperature (to calculate relative humidity) at 2 m height, mean sea level pressure, total rainfall, net shortwave radiation, and downwards longwave radiation.
The simulation covers the period from 1 January 1985 to 1 January 2014. The model was initialized with a cold start on 1 January 1985. The first year (from 1 January 1985 until 31 December 1985) is considered the spin-up period and the results were excluded from our analysis. Finally, it should be noted that our experiment excluded any kind of relaxation (climatological or atmospheric).

Data for Model Validation
Various sources of oceanographic information have been used to validate the model. Sea-surface temperature information was provided by the Copernicus Marine Environment Monitoring Service (CMEMS-SST GLO SST L4 REP OBSERVATIONS 010 024, [39]) for the period 1986 to 2013, and an ultra-high resolution SST product [40] for the year 2010. Sea Level Anomaly (SLA) and Mean Dynamic Topography (MDT) were also provided by CMEMS (product names SEALEVEL MED PHY L4 REP OBSERVATIONS 008 051, SSALTO DUACS, SEALEVEL MED PHY MDT L4 STATIC 008 066). Historical data of temperature and salinity covering the 28-year period came from SeaDataNet database (Mediterranean Sea-Temperature and salinity Historical Data Collection SeaDataCloud V2 [41]). The validation of the model is extended to surface current information from the northeast Aegean, the immediate receptor of Black Sea Waters in flowing into the Aegean Sea through the Dardanelles Strait. This region is covered by the WERA-type High-Frequency radar "Dardanos", with information provided for the year 2010, January-June 2011 and July-August 2021. Of course, the ocean currents from the HF radar cover a much smaller geographical area than the numerical domain of the simulation, however, our validation with this source of information is important for several reasons, such as reproducing exchange fluxes between the two seas, simulating regional dynamics, as well as developing tools for coastal pollution management and search-and-rescue (SAR) operations.

Methods for Model Validation
A crucial post-process following the completion of the numerical simulations is the estimation of the model's statistical error. To that aim, we extensively validated our results with all available information; satellite data (Sea Surface Temperature and Sea Level Anomaly), in situ (T/S profiles from R/Vs) and HF radar surface velocity fields. The model's error for each run was estimated using standard statistical metrics such as Root Mean Square Error (RMSE), Bias skill, Pearson correlation coefficient and Murphy skill score [42] as a model skill index when compared with altimetry data. The significance of each statistical metric used is dependent on the nature of the investigated parameter: for example, in comparing modelled and observed T/S profiles, RMSE and Bias skill offer a more thorough understanding than correlation coefficient or Murphy Skill score.
One of the criteria used to evaluate the model's performance is the ability to reproduce the surface circulation. As the broadest spatiotemporal coverage for surface circulation is provided from satellite altimetry and the consequent calculation of geostrophic currents, this method was chosen for the comparison with simulated currents in the North Aegean.
The kinetic energy of the surface geostrophic currents is obtained from the model's sea surface height output (ζ) using a finite differences approach for the geostrophic approximation where u g and v g are the zonal and meridional components of the geostrophic velocity, g is the gravitational acceleration and f the Coriolis parameter. The same procedure is followed for the estimate of kinetic energy using sea surface heights derived from satellite observations: Absolute Dynamic Topography (ADT) was obtained as the sum of SLA plus the multiyear climatological MDT, covering the period from 1992 to 2012 for the Mediterranean Sea.

Sea-Surface Temperature
Initially we compared the sea surface temperature (SST) of our model to high resolution Level 4 reprocessed satellite SST data. The results from our analysis are shown in Figure 2. The region exhibiting the poorest performance of the model is the vicinity of the Dardanelles exit to the Aegean Sea, as revealed by both Bias and RMSE indices (Figure 2a,b). The complexity of the Turkish Straits System and, especially, the hydraulic control that occurs at the Nara Passage [17,43] combined with the hydrostatic nature of the model-although the number of vertical levels is high-pushes the vertical mixing scheme locally to its limits. This results in a mismatch between the modelled and the observed values, but not for all months, as can be seen from Figure 3, where the seasonal RMSE from the mean climatological year as calculated from the model's results are presented. The error seems to be systematic and tends to decrease during the year with higher values during winter months and lower during autumn. Figure 4 shows the interannual variability and linear trends of SSTs from satellite observations and model results. Trend values from model results are lower (0.038 • C/year) compared to the satellite data trend (0.054 • C/year). Although the numbers are of the same order of magnitude, the lower trend of the simulations seems to be due to the initial overestimation of the SST. Nonetheless, most of the variability remains within the confidence interval, as defined by one standard deviation of the satellite observations, even during the years of 1992 and 1993 when the EMT peaked.   From the above figures, we conclude that our model tends to overestimate the sea surface temperature with a positive bias, both interannually and seasonally, but for a 28year long simulation the magnitude of error RMSE and Bias is considered acceptable for a free run. Overall, the model satisfactorily reproduces the interannual and intra-annual SST variability, with a Murphy skill score approximately equal to 0.97 on the SST's evolution full signal and 0.92 after removing the seasonal cycle.
The southwest-northeast direction of the increased SST differences in the Dardanelles area (Figures 2 and 3) can be attributed to the misrepresentation of the southward extent and the path, followed by the BSW plume by the model, mainly affected by the uncertainty of volume inflow at the Dardanelles Strait and the influence of wind forcing.
In fact, we should stress that the area east of Lemnos (the Lemnos Plateau) hosts the strongest frontal zone in the Aegean Sea [44]. The inflow of BSW creates a thin, buoyant layer of brackish water in the Dardanelles area [17] which, owing to its small depth, responds rapidly to changes of wind stress [45]. Wind forcing also affects the position and extent of the Samothraki anticyclone [46], the south-eastern flank of which may contribute significantly to the circulation of the northern part of the Lemnos Plateau [47]. Under the influence of north winds, the path of the BSW bifurcates and its southern branch circulates around Lemnos Island. To the south, the BSW is bound by Levantine-originated waters, carried northward by the general cyclonic circulation of the Aegean Sea [48]. Except for during the summer season, the BSW can be identified in SST images by its low temperature signal [49]. During summer, the dominant north winds in the area cause upwelling in the eastern coasts of the Aegean Sea [50,51]. Two of the upwelling centres, easily spotted in the historical satellite SST images during summer, are in the vicinity of the Dardanelles area; one along the Turkish coast south of the Dardanelles exit and another at the west coasts of Lesvos Island [52]. Moreover, the intensity of wind curl and the response to Ekman transport are linked to the intrusion of the Levantine-originated waters from the south-eastern part of the Aegean Sea further to the north [53,54].
To investigate the increased differences between the model and satellite SSTs, we have used the available observations of surface currents from the HF radar during August Note that the wind conditions during August 2010 are markedly different from those during August 2021 (Figure 5c,g). A stronger Ekman transport during 2010 intensified the Samothraki anticyclone and moved its east flank to the west, between Lemnos and Imvros islands (Figure 5d,h). The cold upwelled waters of the eastern coasts of the central Aegean were transported to the southern part of the area measured by the radar (Figure 5b), where they formed a frontal region of strong convergence (Figure 5a). During 2021, the influence of the anticyclone of Samothraki was found between Imvros and the Dardanelles exit ( Figure 5h). Cold waters moving to the west from the upwelling centre to the southeast of Tenedos Island can be observed in Figure 5f. The convergence region has moved northwards, and its direction has changed with respect to that of August 2010 ( Figure 5e). As the region of the largest model SST error is aligned with the thermal front and drops moving away from it for both periods (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013), the model error increased for the period forced by the climatological field (Figure 5i-l), and taking into account that the interannual variability of the BSW is larger than the seasonal variability [20], we conclude that the largest part of the model SST error in this region is caused by the misrepresentation of the position of the plume of BSW, attributed to the combination of the uncertainty of BSW inflow, the low spatial resolution of wind forcing and the artificial numerical mixing at the Dardanelles Strait.

Sea Level Anomaly
The direct comparison of satellite sea level anomalies required the estimation of steric height as a post process, due to the model's hydrostatic approximation. Although non-Boussinesq steric effects have a negligible imprint on monthly or longer regional sea level patterns, they do impact global mean sea levels [55,56]. The mean sea level of the regional model [57], from which we extracted the boundary conditions, does not include the global mean steric effect. Thus, the steric height component was estimated using a maximum integration depth of 600 m and the Gibbs Seawater (GSW) Oceanographic toolbox [58]. The estimated daily steric height was added to the model's dynamic height (ζ). The results of the sea-level validation are presented in Figure 6 in the form of error maps (Bias and RMSE) and as basin-integrated time series for 21 years. We also validate our simulation of mean dynamic topography (MDT), which reflects the average strength of geostrophic currents. We used the CMEMS MDT 2020 product which provides an estimate of the mean sea surface height above the geoid on a 4.6 km resolution in both directions horizontally. It is based on altimetry data, in situ measurements of temperature and salinity and covers the period 1992-2012. Figure 7 shows a comparison between simulations and satellite observations, both for the general circulation (as expressed by Mean Dynamic Topography, panels a and b) and its variability (as expressed by the standard deviation of the Sea-Level Anomaly, panels c and d). The Mean Dynamic Topographies are expressed as anomalies with respect to the area-weighted regional mean, as estimated from observations and simulation.
The relatively crude grid-size of satellite altimetry products (1/8 of a degree) results in lower effective resolution compared to our model's. For that reason, small-scale features in the model results are difficult to evaluate. Moreover, land contamination leads to a lower accuracy of satellite altimetry in coastal regions [59], while corrections for the inverse barometer effect, tides, and high frequency winds [60], in satellite altimetry products lead to reduced accuracy in regions such as the North Aegean Sea, where these corrections are difficult to model.
Comparison of simulated and observed Mean Dynamic Topography reveals that the model reproduces very satisfactorily the anticyclonic flow around Samothraki Island, and the cyclonic circulation south of Lemnos Island (Figure 7a,b). The Sea-Level Anomaly standard deviation maps show highest variability in the area of the Samothraki anticyclone, and least in Skyros basin, with the model exhibiting higher overall variability than the observed values, a fact possibly related to the lower resolution of the satellite data, as discussed above.  Based on our calculations, the model underestimates the sea level by approximately 2 cm (Figure 6a,c) and has an overall error of 4.5 cm. Murphy's skill score for the full signal is approximately 0.75 and for the de-seasoned 0.76.

Hydrographic Properties-θ/S Profiles
The most important part of the validation is, in our opinion, the direct comparison of the model's results with θ/S profiles to evaluate its capability of reproducing the vertical structure and characteristics of the area of interest. The overall model skill for potential temperature, salinity, and potential density from the available casts (mainly R/V CDT data, complemented by XBTs and Bottles) for the period from January 1986 to December 2013, are presented in Figure 8. Positions are shown in Figure 1 and the number of acceptable data-QC flag "good"-is 5119. Unfortunately, only a small proportion of the available profiles-639 casts or 12.5% of the total number of casts-were recorded from the North Aegean's deep basins, and as seen from the histogram in Figure 8, a significant issue is the temporal inhomogeneity of observations during the last few decades. The vast majority (87.5%) covers depths from the surface to 450 m. A quick and obvious finding from the model performance positive Bias is that it tends to slightly overestimate both salinity and temperature. However, in terms of potential density, the results show the that the model is capable of reproducing the observed water types, having a positive Bias of 0.04 kg/m 3 and a good-higher than 0.5-performance index, according to Murphy's skill score for all fields (T,S,ρ). A more focused analysis of the hydrographic results for the surface and the intermediate water column is provided in Figure 9. The presented results are the outcome of the comparison between the model's output and the available θ/S profiles for the period from 1986 to 2013 in terms of hydrographic properties and error estimation. First, the mean profile from all observational data was estimated by averaging in time and space for the full extent of the North Aegean. Model data were interpolated in time (to accurately match the cast time of the observations because we kept only daily means for the prognostic fields of temperature and salinity) and averaged spatially. From the mean observed and modelled profiles, the Bias and the RMSE of the first 400 m were estimated, using a 25 m vertical step. Panels a and b (Figure 9) show the mean salinity and potential temperature profile as calculated from model's data and observations, as well as an estimate of the observed variation (±one standard deviation). Panels d, e, f, and g give the error estimation expressed using Bias and RMSE for salinity and temperature. Finally, the overall θ/S diagram for all profiles for the upper 400 m for the 28-year-long (panel c) exhibits a very good reproduction of the observed water-masses by the hindcast. The temperature overestimation extends from the sea-surface (as seen already in Section 3.1.1) to the upper 160-170 m (Figure 9g) having about the same magnitude (0.2 • C) as the SSTs. The Bias reverses sign to negative in the layer below and gradually decreases. This is in agreement with the RMSE vertical distribution (Figure 9f), where the best score-of almost 0.03 • C-is approximately at 180 m, decreasing below 300 m. A possible cause for the temperature positive Bias in the surface layer may be the Jerlov water type selected to replicate the optical properties of the seawater in the North Aegean. After several sensitivity runs, it was set as type 3 in ROMS nomenclature, i.e., IB Western Mediterranean [61], constant for the whole numerical grid.
For the case of salinity, a high RMSE appears at the first 100 m ( Figure 9d) and again, as in temperature, a positive Bias for the upper 170 m. Below that depth, both RMSE and Bias gradually decrease until 400 m, where they obtained the best scores. A part of the error comes from the southern open boundary, where data from a Mediterranean regional-scale model were used [33]. According to available quality information documents for the version providing the boundary conditions, the regional-scale model exhibited an overall positive Bias for salinity, notably for the Aegean Sea region. The contribution of the eastern open boundary to the overall salinity error, at first glance, is mainly local, as it mostly affects the Dardanelles exit several kilometres away from it but contributes to the overall salt budget. An analysis for the eastern open boundary behaviour and its effect are shown in Figure 10.  [20] for the eastern open boundary of our computational grid offered not only locally, as presented at Figure 10a,b, but overall (for both surface and the upper layers of intermediate waters, Figure 9) reasonable results concerning the characteristics/properties of the Black Sea waters that enter the North Aegean through Dardanelles straits.
Another contribution to the salinity positive Bias comes from the riverine contribution in the North Aegean Sea. The riverine importance in terms of nutrient load, as well as the thermohaline functioning of the area, is extensively discussed by Tsiaras et al. [62]. In our case, the lack of salinity data in the river outflow time series from E-HYPE forced us to use a constant value for the simulation. After several test runs, we set the salinity value equal to 14 in order to obtain stable solutions and mainly to remove unrealistic freshwater surface fronts with sharp vertical gradients at the first 10 m in the coastal areas close to the rivers exit. Numerical stability was eventually achieved, but with the aforementioned cost potentially added to the total salinity error.
The second part of our analysis for the θ/S profiles is focused on the deep basins of the North Aegean Sea and their variability. We have extensively compared the available hydrographic data for each of the basins with the model's output and the final results are presented at the following figures (Figures 11-13). Lemnos basin (Figure 1 cyan rectangular) is an oblong basin with a mean depth of approximately 500 m at the northeast part of the North Aegean. The deepest part is 1582 m and is located north of the island of Lemnos. The number of θ and S profiles is 184 but with significant temporal inconsistencies. Figure 11 shows the mean salinity and temperature profile from the surface to 1000 m depth and its errors estimation plus a θ/S diagram. The mean profiles were calculated with the same approach as described previously for surface and the intermediate depths.   The θ/S diagram (panel c) of the above figure exhibits a realistic representation of the Lemnos basin water mass structure, especially for the deeper-than-400 m layers. The mean simulated salinity profile is within one standard deviation of the observed (panel a), while the corresponding temperature profile is very close to observations (panel b). The model slightly underestimates salinity, with an almost constant with depth error (panels c and e). A similar picture is seen for temperature, with a decreasing error as we go deeper in the water column (panels f and g). The maximum error for both tracers is found in the upper surface layers. For temperature, although the overall Bias has a positive sign (0.16 • C), from the surface to a depth of approximately 60 m it becomes negative with the maximum absolute value for Bias to be equal to 0.32 • C. From that depth and below, Bias and RMSE gradually decrease until the depth of 250 m. The maximum positive Bias is found at the 400 m and is equal to 0.1 • C. Figure 12 shows the hydrographic characteristics of the Athos basin (Figure 1 blue rectangle) for the upper 1000 m (maximum depth 1561 m), where the number of the available observations is almost twice that of the Lemnos basin. The mean profiles (panels a and b) and the θ/S diagram (panel c) reveal the model's capability for this basin also its ability to replicate, in a realistic way, the structure of the water column from the surface to the deepest point. The errors for temperature and salinity (panels d-g) are smaller compared to north Lemnos basin and exhibit the same behaviour. An interesting finding is salinity's Bias value below 200 m, which is constant and even smaller than the north Lemnos basin, with a value equal to −0.019. Summarizing for the Athos basin, it can be concluded that the model has a very good performance in this part of the North Aegean.
The results for the north Skyros basin are shown at Figure 13. The number of available θ/S profiles was smaller than for the other basins (133) and had the biggest temporal gaps. The maximum depth of the basin is 1091 m. Both parameters exhibit a negative Bias which is almost constant for salinity and variable for the upper 400 m for temperature. The highest undershot is located in the 50 m where the RMSE is 0.32 • C and Bias −0.29 • C. As seen in the previous figures, the error is reduced below 400 m. Again, the agreement of the water masses in the θ/S diagram is acceptable as for the other two basins, despite the relative scarcity of observations.
As the validation results point to a very good performance of the model in the deep layers, which are considered the "clearest" recorders of climatic variability, it is possible to analyse the evolution of the deep layers, expressed through the interannual variability of potential density (σ θ0 ), integrated from 600 m to the maximum depth of its basin with a daily step (Figure 14). The major deep-water formation events that happened during the late 1980s and early 1990s-1987, 1992 and 1993-known as the EMT, are seen in the time series presented along with other (smaller in magnitude) events that took place during the simulated 28-year period. One more argument about the validity of the simulated results is provided in Figure 15, where the interannual variability of the north Aegean deep basins' potential density from the available observations is compared with the model output, integrated vertically below 600 m. The Bias sign in the three basins reveals that the model slightly underestimates the density of the deep waters. The total error (RMSE) is exceptionally small for the time span of the simulation.
The overall scores of the simulation show that the model's error is in an acceptable range for a free run and even comparable/adequate with an assimilative one [63]. Notably, the model underestimates the density of the deep-layer waters produced in the EMT period (1987)(1988)(1989)(1990)(1991)(1992)(1993), however it reproduces well the rate at which the density of the deep layers progressively decays in the three layers.

Surface Circulation
Regarding the surface circulation, the model reproduces to a high degree the general pattern, as described from previous works [18], and is in agreement with the satellite observations [64]. In the northern part of the domain (Figure 16b), the most prominent feature is the distinct pathway of the Black Sea waters that enter the North Aegean Sea through the Dardanelles and split east of Lemnos Island creating two branches, one going to the northern part and contributing to the Samothrace eddy. The second branch initially follows a south-westward route to the central Aegean through the Evia island jet and later westwards forming the Chios cyclone. Regarding the geostrophic kinetic energy distribution, the highest energy levels for the model (Figure 16d) are located mainly south of Lemnos Island and at the north-eastern part of the domain, in the region between the islands of Thasos, Samothrace and Lemnos, where an anticyclonic system dominates the local circulation.  Consequently, the observed highest kinetic energy is found at the same areas but with a smaller magnitude. An interesting finding is the low signal of the Sporades eddies in the satellite data results (Figure 16c). In both cases, the differences in magnitude have to do with the higher spatial resolution of the model (1 km) compared to the original satellite information (1/8 of a degree). Due to the low resolution of the observed sea-level fields, the high-wavenumber part of the spatial two-dimensional spectrum does not contribute to the total variance, which explains the mismatch.
Another validation of the model's ability to reproduce the surface circulation can be provided by the presence of the HF radar system "DARDANOS" covering the wider area of the Dardanelles outflow in the northeast Aegean [46]. The validation is held via comparison of the seasonal surface circulation at the Dardanelles strait exit for 2010, as derived from HF radar data and the model. The radar's reconstructed data have 30 min temporal resolution and 0.5 km spatial resolution. To compare them with the model's results, daily mean fields are calculated and interpolated linearly to match the model's resolution. Overall, the model reproduces well the local circulation ( Figure 17). The model has an error pattern similar to that of the seasonal sea surface temperature (Figure 3). The maximum deviation takes place during the winter months (January, February and March) and then gradually reduces as time evolves. A possible explanation beyond the parameterization applied for the Dardanelles outflow after 2009 (a climatological mean year from the period 1985 to 2009), can be the low spatial resolution of the wind field input forcing the simulation. This, combined with outflow mismatches due to the parameterization, may be the source of the observed deviations.

Discussion
A new model has been developed for the Aegean Sea, based on the ROMS platform, suitably configured to enable climatological studies of the basin. The replacement of seasonal-cycle parameterizations of the exchanges with the Black Sea by time-dependent fluxes, as determined by Maderich et al. [20], as well as the removal of all climatological relaxation terms from the numerical schemes, enables the simulation of changing climatic conditions. An evaluation of the model performance through the comparison of the 30-year-long hindcast to all available oceanographic information revealed that the model, despite the removal of relaxation terms and the lack of data assimilation, exhibits very small Bias and RMSE, and, especially important for long-term simulations, shows no sign of significant numerical drift (as far as can be deduced from the available observations). Furthermore, the model satisfactorily reproduces the surface circulation of the North Aegean Sea, both on a large scale (i.e., the whole basin) as well as over the critical region of the vicinity of the Dardanelles exit. Thus, the model provided a well-validated 30-year-long hindcast  of the hydrodynamics of the Aegean Sea, suitable for a range of further studies. Of special interest is that the model exhibits its best performance in the deep basins of the North Aegean Sea, repositories of climatic variability information, shedding light on the evolution of the deep layers throughout the whole hindcasting period. Moreover, the model hereby introduced exhibits very good performance, with small biases in the temperature and salinity profiles. Potential sources of the salinity Bias, such as the lateral boundary conditions and the riverine inputs, have already been discussed. An open question is the net E-P contribution, as estimated from the COARE [38] bulk formulae using ERA5 input (for the total accumulated precipitation and radiative fluxes).
The question is raised from the known issues that previous ERA datasets (ERA 40 and Interim) had in the overall rainfall rates globally and fluxes locally in the Mediterranean Sea [65,66]. Although the ERA5 model is significantly improved in total and over the global ocean, we do not have any knowledge from the available bibliography for its skill over the Mediterranean Sea and how it may impact the overall E-P-R budget.
Regarding the reproduction of the sea-level, the skill score index of our model results is lower compared to the SSTs score. However, the model adequately reproduces the seasonal, inter and intra-annual variability, and some of the deviation between the simulated and observed sea level should be attributed, in our opinion (and according to the available literature), to the temporal variability of the mass component of sea level (contribution from meltwater and salt, [67]). A noteworthy fact is that the model reproduces the large positive sea-level anomalies in 2010 and 2011, which is in agreement with findings by Landerer and Volkov [68] during the same period in the Mediterranean Sea.
The hereby produced 30-year hindcast, exhibiting the least Bias in the deep layers, fills the observational gaps in the evolution of the deep basins (Figures 14 and 15). Based on the slopes of the density time-series, three significantly different periods can be identified: The short and episodic dense-water formation events (namely in the winters of 1987, 1992, 1993 and smaller episodes in 2008 and 2012, and a gradual increase of density after 2008), the periods of stagnation (characterized by decreasing densities), and the steady-state periods (where the density remains more of less constant). The periods of stagnation correspond to a dynamical balance between the rate of buoyancy gain in the deep basin and the downward vertical buoyancy flux through the interface between intermediate and deep waters [69,70]. The periods of steady state usually take place after prolonged stagnation periods have led to gradual homogeneity between the intermediate and deep layers, and, thus, the turbulent fluxes across the boundary are minimal. Based on this work's result, which has filled gaps in the observational record due to lack of frequent of CTD casts, the most distinctive stagnation period in the North Aegean deep basins followed the EMT and lasted until about 2004, while the period between 2004 and 2008 can be characterized as a steady-state period.
Overall, north Skyros basin appears more sensitive to buoyancy forcing, exhibiting higher overall density variability and direct deep-water formation in the winters of 2006, 2008 and 2012, while the deep layers of Athos exhibit a gradual density increase, suggestive of a slow replenishment mechanism through lateral exchanges with neighbouring basins, possibly the Skyros basin.
Thus, the hereby introduced model is suitable for high-resolution climatological studies of the Aegean Sea. The high resolution of the numerical domain enables the full simulation of mesoscale geostrophic circulation, as well as the simulation of tidal propagation and the dynamics of the low-wavenumber part of the internal-wave spectrum. Despite the fact that the hereby-presented hindcast exhibits acceptable deviation from the observed Aegean Sea hydrographic properties and variability, the model set-up permits the planning and execution of more advanced numerical experiments, incorporating astronomical forcing, and thus simulating high-frequency processes, such as the generation, propagation and dissipation of internal baroclinic tides, the energy exchange between tidal frequencies and the internal-wave continuum, and the energy flow towards turbulence and mixing. Furthermore, the success of the simulation of the above small-scale processes can be assessed (as compared to the "basic" hindcast hereby presented) based on the same methodology and observations used in this work. Funding: This work has been supported by the Policy-oriented marine Environmental Research for the Southern European Seas (PERSEUS) funded by the EU under FP7 Theme "Oceans of Tomorrow" OCEAN.2011-3, Grant Agreement No. 287600. This work was also partially supported by the project "Hellenic Integrated Marine Inland water Observing, Forecasting and Offshore Technology System (HIMIOFoTS)", implemented in the framework of Operational Programme Competitiveness, Entrepreneurship and Innovation 2014-2020 (EPAnEK) (MIS 5002739), cofunded by the European Union through the European Regional Development Fund and the Hellenic State through ESPA 2014-2020.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.