Black Sea thermohaline properties: Long-term trends and variations.

The current knowledge about spatial and temporal dynamics of the Black Sea's thermohaline structure is incomplete because of missing data and sparse distribution of existing measurements in space and time. This study presents 56 year continuous simulations of the Black Sea's hydrodynamics using the 3D General Estuarine Transport Model (GETM), without incorporating any relaxation toward climatological or observational data fields. This property of the model allows us to estimate independent temporal trends, in addition to resolving the spatial structure. The simulations suggest that the intermediate layer temperature is characterized by a weak positive trend (warming), whereas the surface temperature does not show a clear linear trend. Different salinity trends have been established at the surface (negative), upper (weaker negative) and main halocline (positive). Three distinct dynamic periods are identified (1960-1970, 1970-1995, 1995-2015), which exhibit pronounced changes in the Black Sea's thermohaline properties and basin circulation. Strengthening of the main cyclonic circulation, accompanied by intensification of the mesoscale anticyclonic eddy formation is found. Both events strongly affect the sea surface salinity but contribute in opposing directions. Specifically, strong composite large-scale circulation leads to an increase in sea surface salinity, while enhanced formation of mesoscale anticyclones decreases it. Salinity evolution with time is thus the result of the competition of these two opposing yet interdependent processes.

Introduction

The Black Sea has unique natural conditions like a positive net freshwater balance, strong basin wide cyclonic circulation, and the so called Rim Current, with well exhibited western and eastern gyres. The Cold Intermediate Layer (CIL) is defined by temperatures less than 8°C in the sub‐surface Black Sea's waters and over 90% of the basin deep water volume consists of anoxic water. Some of Europe's longest and largest rivers flow into it, including the Danube and the Dniepr. Fresh water from the rivers is spread along the shelf before entering the Rim Current and eventually is mixed with the waters from the basin interior. The interaction between the oxygen rich surface waters and the deep zones is limited by density stratification, which is amplified only under particular meteorological conditions, leading to a layered structure that affects the diversity of the organisms within the Black Sea [Murray et al., 1991; Oguz et al., 1992; Özsoy and Ünlüata, 1997; Ivanov et al., 1997]. Horizontally, two distinct regions can be outlined: one being the wide and shallow North‐Western Shelf (< 200 m), the second being the deep sea, which is bounded by the 1500 m isobath (Figure 1). The latter is mostly isolated from the direct riverine inflow, which is known to be a key driver in the shelf. However, the mesoscale eddies evolve along the periphery of the basin, as a part of the Rim Current dynamic structure. These effectively link coastal hydrodynamic and biogeochemical processes to those in the deep sea and thus provide a mechanism for two‐way transport between nearshore and offshore regions [Zatsepin et al., 2003]. These two regions show pronounced physical and biological differences [McQuatters‐Gollop et al., 2008]. Climate affects the Black Sea via atmospheric forcing and riverine inflow, which has been demonstrated to be a significant factor for the overall water balance and basin‐scale circulation [Oguz et al., 1995]. The Black Sea's hydrodynamic characteristics change irregularly in time, with estimated frequency from several years to several decades [Ginzburg et al., 2004; Oguz et al., 2003; Oguz et al., 2006; Kazmin et al., 2010; Zatsepin, 2007]. This leads to differences in the estimates of long‐term trends and, in most cases, these differences are fundamental as the identified trends can even be an order of magnitude apart. It is also crucial to study not only the trends that characterize the hydrodynamics of the surface layers, but also hydrological characteristics in the most stratified part of the pycnocline (at depths 50–200 m). Indeed, the presence of temperature or salinity trends in the intermediate layers of the Black Sea could be a regional consequence of global climate change. No positive long‐term temperature trend in the surface waters of the Black Sea has not been found [Polonsky and Lovenkova, 2004; Shapiro et al., 2010]. However, this does not mean that profound changes caused by regional variations of climate are absent [Mikaelyan, 1997; Oguz et al., 2003]. Isolation of trends in the deep layers of the sea is difficult because of the scarce and low quality of data for characterizing those layers [Ivanov et al., 2000 ; Knysh et al., 2011]. Simulating reasonably the temperature and salinity (T/S) evolution on seasonal, annual, inter‐annual and decadal scales is a key factor for determining possible trends in the Black Sea thermohaline structure.

Figure 1

Bathymetry and location map of the Black Sea, including the boundaries of studied zones. The 1500 m isobath is drawn in magenta, while the boundary of the shelf and deep sea compartments, separated by the 200 m isobath is shown in green. The locations of two points from 2001 R/V Knorr cruise are denoted with A1 and A2, while with W1 and W1 are shown the locations where the WOA13 (https://www.nodc.noaa.gov/) salinity data is extracted. The black ‘o’ indicates the location of the Rim current speed plotted in Figure 13.

Understanding the evolution of the Black Sea's physical characteristics is further hampered by the lack of winter‐time hydrographic observations [Ivanov et al., 2000 ; Belokopytov, 2011]. Moreover, the time and spatial distributions of in‐situ measurements are very irregular and incomplete. Many validated models have already been used to simulate the Black Sea hydrodynamics and in this way give a 3D look inside the basin's characteristics, filling the gaps in observations. Generally, they can be divided in two main groups: data assimilation models and models with strong or weak relaxation. A model without any relaxation is reported in Capet et al. [2012], however neither sea surface salinity (SSS) nor water column salinity is presented and discussed. Numerous models developed in 21st century are based on data assimilation [Besiktepe et al., 2001; Korotaev et al., 2003; Cannaby et al., 2015; Dorofeev et al., 2013; Dorofeev and Sukhikh, 2016], and they required broadly available satellite data for sea surface temperature (SST) and altimetry to constrain model runs. Due to the lack of regular data for SSS, typically climatological data sets are used instead. Usually, such type of models failed to reproduce the vertical Black Sea thermohaline structure (e.g., CIL properties) in a realistic way, but are appropriate for operational forecasting (CMEMS, http://marine.copernicus.eu/). The second group involves stand‐alone ocean models with different relaxation schemes to climatological or observational SST or/and SSS time series. A variety of numerical approaches exist, including the Bryan‐Cox primitive equation model [Stanev et al., 1995], the Geo Hydrodynamics and Environment Research (GHER 3D) model [Stanev and Beckers, 1999], the Dietrich Center for Air Sea Technology (DieCAST) model [Staneva et al., 2001], the Princeton Ocean Model (POM) [Oguz et al., 1995], the GFDL Modular Ocean Model (MOM) [Stanev et al., 2004], which is a later version of the Bryan‐Cox model and the Harvard Ocean Prediction System (HOPS) [Besiktepe et al., 2001]. In Kara et al. [2005] the Hybrid Coordinate Ocean Model (HYCOM) is applied, while in Enriquez et al. [2005] ‐ the Proudman Oceanographic Laboratory Coastal Ocean Modelling System (POLCOMS). Newer adaptations of these models still use SSS relaxation mainly because of the difficulty in accurately reproducing the surface salinity as well as its large vertical gradients. A key feature of all models is their need to use SSS time series to operate, while at the same time, Black Sea SSS observational data is scarce and inhomogeneous. Despite being limited, the existing observational evidence can be used to some extend for relaxation in hindcast simulations [Korotaev et al., 2003; Dorofeev et al., 2013; Dorofeev and Sukhikh, 2016]. However, we crucially need a model capable to perform scenarios for the future as well. For example, to create scenarios of potential future climate change one would need to design in advance the surface salinity responses to certain changes if relaxation is needed in the modeling tools. Furthermore, realistic simulation of T/S fields is also important for the nutrients budget in the Black Sea, which is in turn important for the phytoplankton bloom and hence the changes in regional ecological characteristics [e.g., Oguz et al., 2003, 2004, and 2006].

This motivated us to develop a model capable to simulate the mesoscale circulations and thermohaline structure in the Black Sea for a continuous multidecadal period without any relaxation toward external fields. Our 3D hydrodynamic model comprises 3D GETM (http://www.getm.eu/) and General Ocean Turbulence Model (GOTM). GETM applies primitive equations with hydrostatic and Boussinesq approximations to calculate the 3‐D flow field and the free surface. The quality of the forcing data affecting our simulations has been analyzed in Miladinova‐Marinova et al. [2016b] and Miladinova et al. [2016a] and the most appropriate forcing data capable to assess the potential changes in the Black Sea dynamics has been selected. Our main aim is to develop a model that is able to reproduce the Black Sea's main physical features in a long‐term simulation without depending on relaxation data. Fine spatial and long‐term temporal variations of the Black Sea thermohaline properties have not been examined previously. To fill in this gap, this study emphasizes the realistic representation of the Black Sea circulation and mean thermohaline structure, with a particular focus on the mean SSS variation. As the properties of salinity, temperature and density are approximately uniform in the deep water [Murray et al., 1991] we discuss herein the dynamics of the upper surface down to 250 m depth. Additionally, we evaluate the long‐term trends of T/S and summarize the evidence for local changes based on trends established for the Black Sea basin.

The paper is organized as follows. Section 2 contains a brief description of the study area, model approach and setup. A short description of data sources used for model validation is also given. Extensive model comparison against independent basin mean data is performed in section 3. Possible factors for SSS inter‐annual and long‐term variations are discussed in section 4. Long‐term changes of the surface, as well as, water column thermohaline properties are given in section 5. Some aspects of the Black Sea mesoscale variability are also discussed in section 5. The summary of the study is presented in section 6.

Materials and Methods

Study Area

The Bosphorus Strait connects the Black Sea with the Mediterranean Sea via the Marmara Sea and the Kerch Strait is the connection with the Azov Sea (Figure1). The Black Sea shelf edge slope is steep and the shelf is narrow except for the North‐Western Shelf region. We consider three distinct regions: the deep basin, which is bounded by the 1500 m isobath, the intermediate area between 200 and 1500 m and the shelf area of less than 200 m depth (Figure 1). This study is focused on the description and analysis of the properties and dynamics of the interior deep basin region (magenta line in Figure 1)

Numerical Model

The hydrodynamic model comprises the 3D GETM and 1D GOTM initialized on a high resolution 2 × 2 min latitude–longitude horizontal grid. The model bathymetry grid is produced from ETOPO1 global bathymetric grid with horizontal resolution of 1 min. Linear programming procedure was applied to smooth slightly the bathymetry [Sikirić et al., 2009]. The maximum depth of the model domain is 2200 m with a 70 levels general vertical grid, which is compressed toward the surface. A detailed description of the GETM equations can be found in Stips et al. [2004] and Peneva and Stips [2005]. Recently, this model has been successfully applied to simulate the Wadden Sea and the North Sea [Gräwe et al., 2016; van Leeuwen et al., 2015]. One way to minimize dissipation and dispersion is to use a numerical method which satisfies the Total Variation Diminishing (TVD) property [Harten, 1983]. Flux‐limiter methods satisfy the TVD property and can switch between a second‐order approximation (when the field is smooth) and a first‐order approximation (when it is near a discontinuity). Flux‐limiter methods have been applied to different numerical approaches so that oscillations present in the numerical solution can be minimized. GETM does not involve any internal procedure that relaxes the ocean temperature and salinity toward climatology. In the case of the Black Sea long‐term simulations, the most stable simulations were achieved when using a fine horizontal and vertical grid and applying a second‐order scheme combined with the second‐order monotone scheme with the Superbee limiter [Burchard and Bolding 2002].

The meteorological forcing from the European Centre for Medium Range Weather Forecast (ECMWF) available from http://www.ecmwf.int, based on 6 hourly records has been applied. It involves ERA‐40 project (1958–2001) and ERA‐Interim project (1979–2015). The effect of the meteorological forcing data on the model performance is studied in Miladinova‐Marinova et al., [2016b]. The data sets coming from ECMWF have been chosen as the best available meteorological forcing because of two main factors: it includes long‐term data records and the best circulation features are obtained when our model is using ECMWF data. It takes about 2 years for the model spin up. GETM is reading the precipitation data from the meteorological file, but calculating evaporation and condensation internally. This method has the advantage that the calculated humidity flux is consistent with the calculated latent heat flux. Two model runs have been chosen in order to identify possible artifacts arising from initialization and forcing uncertainties: main run with ERA‐40 (1958–1979) followed by forcing with ERA‐Interim (1980–2015) and run2 forced with ERA‐Interim (1979–2015). In order to achieve a consistent representation of both meteorological data sets, all atmospheric variables are downscaled using spline interpolation to a 0.5o×0.5° grid from the original ECMWF Gaussian grid. This method also slightly improves the conservation of the wind speed curl.

Freshwater input has been estimated using the values from the Global Runoff Data Centre (GRDC, http://www.bafg.de/GRDC) runoff. Being an estuarine basin, the Black sea is very sensitive to variations in the fresh water input. This is because the resulting buoyancy flow, induced by the river runoff, is essential for establishing the basin circulation. Comparison between runoff data sets from different data centers has revealed similar climatological mean annual cycles for all rivers considered herein [Miladinova‐Marinova et al., 2016b]. We have chosen to force our model with the GRDC data because it contains long‐term daily records of the Danube River. The long‐term variability of the thermohaline forcing is given in Figure 2. Note that the difference of precipitation and evaporation is negative due to excess of evaporation over precipitation. Water exchange in the Bosphorus and Kerch Straits is simulated as a river flow that contains surface outflow/inflow and bottom inflow/outflow. Starting from the long term steady state water and salt budgets in the Black Sea, the monthly averaged volume fluxes have been estimated [Miladinova‐Marinova et al., 2016b] and used as a further forcing condition. Both inflow and outflow have their own seasonal cycles, however we fixed the annual inflow flux to about 10×103 m3 s−1 [Ünlüata et al., 1990] and then estimated the corresponding outflow. The annual climatology of the Bosphorus net flow is also shown in Figure 2.

Figure 2

Annual climatology × 103 m3 s−1 of river runoff, precipitation‐evaporation and Bosphorus net flow.

The model is initialized by means of temperature and salinity 3D fields coming from the MEDAR/MEDATLAS II project (http://www.ifremer.fr/medar). The MEDAR data set for the Black Sea reflects the main features known from observations – the strong halocline at 70–150 m, the CIL at approximately 25–70 m and the doming of the isohalines due to the cyclonic Rim current.

For estimating structural breakpoints in the analyzed time series, we apply the method developed by Bai and Perron [2003].

Data Used for Validation and Comparison

The Pathfinder data (http://www.nodc.noaa.gov/sog/pathfinder4km/) is preferred for the SST model–data comparisons from 1982 to 2009 as it outperforms the commonly used satellite observations [Casey and Cornillon, 1999], and also has a much finer resolution than the others, which is needed for the fine resolution Black Sea model used in this paper. From 2005 to 2015 reanalyzed thermohaline fields extracted from the Copernicus Marine Environment Monitoring Service (CMEMS, http://marine.copernicus.eu) are used for comparing the model result with independent results. The CMEMS hydrodynamics are supplied by the Nucleus for European Modelling of the Ocean (NEMO, v3.4). The model solutions are corrected by the variational data assimilation. For the SST supplementary comparison, data redrawn from Oguz et al. [2006] is used. The winter‐mean SST data in Oguz et al. [2006] are constructed using the monthly mean data (http://badc.nerc.ac.uk/data/hadisst) from all the available in situ measurements from 1875 to 2000, performed within the interior part of the basin for depths greater than 1500 m, as well as the Advanced Very High Resolution Radiation (AVHRR) satellite observations. The annual T and S variations at the intermediate layers are compared with reconstructed seasonal or annual hydrodynamic data from reanalyzed model products [Knysh et al., 2011].

Model Validation

First we benchmark for the long‐term steady state salt budget in the Black Sea as presented in Figure 3. Monthly mean temperature exhibits much stronger seasonal than inter‐annual variation. Conversely, mean salinity varies less seasonally and more inter‐annually. The overall mean temperature is 9.41 ± 0.05°C without a clear linear trend, whereas the mean salinity is 21.162 ± 0.0045 with a small decreasing trend of 0.0012 ± 2.10−19 per year. We might conclude that the model assumption of a long‐term salt balance is satisfied.

Figure 3

Deep basin volume averaged monthly mean temperature and salinity.

Since our main goal is to develop a fully prognostic model we validate whether it reproduces the mean system behavior for a long time period within acceptable bounds. It is therefore less important if the model is capable to represent completely a particular event in a particular location or period. For this reason the simulated mean state of the system has been compared with mean observations or independent model results.

Surface Characteristics

General circulation of the Black Sea consists typically of the cyclonic Rim Current, meandering along the continental slope, and a number of cyclonic and anticyclonic mesoscale eddies located inside the main current or in‐between it and the coast [Oguz et al., 1993; Ginzburg et al., 2002; Korotaev et al., 2003 ; Zatsepin et al., 2003]. A usual way to show that metrological forcing capture adequately the Rim Current structure is to compare the simulated surface circulation with well‐known patterns already described in the literature. Figure 4a illustrates climatological average surface circulation in May (for 1962–2000 and at 5 m depth) to be compared with Capet et al. [2012] and Staneva et al. [2001]. The upper layer circulation in the deep basin consists of strong cyclonic eddies combined with smaller and weaker anticyclones. In agreement with observations, the Rim Current speed is estimated at 0.3–0.5 m s−1, indicating that the chosen bathymetry smoothing and numerical approach are adequate. Several anticyclonic eddies along the coastal side of the Rim Current are well visible in Figure 4, namely, the Batumi, the Sevastopol, the Crimea, the Synop, the Sakarya, the Kizilirmak, the Bosphorus and the Kaliakra. The largest and most persistent anticyclonic eddy, the Batumi, lies east of 38°E. The surface circulation in Figure 4b could be compared with the 5 dbar objectively mapped dynamic height anomaly from hydrographic observations in Oguz et al. [1993]. The large, sub‐basin and mesoscale surface circulation features observed through the SO90 campaign, which was an extensive, quasi‐synoptic three‐dimensional coverage of most of the Black Sea, are adequately simulated by the model. The detailed comparison with Figure 2 in Oguz et al. [1993] reveals a similar cyclonically dominated basin interior obtained by both observations and fully prognostic simulations. Note the same strong cyclonic eddy in the central part of the basin even though its exact location is different. Concerning the anticyclonic motion, the main bifurcation of the flow in both studies took place near the eastern end of the basin at about 39.5°E, leading to the formation of the Batumi anticyclonic eddy. Moreover, the estimated maximum Batumi surface current speed agrees with the observed speed of about 0.3 m s−1. Figure 4c is a snapshot of daily flow field in the southern part of the basin. It is given to be compared with the altimeter‐derived circulation in Korotaev et al. [2003] for the same date. Our model points out the same five anticyclonic eddies along the Anatolian coast as is shown in Korotaev et al. [2003]. Moreover, the meandering of the Rim Current in Figure 4c agrees with the meandering reproduced by the use of altimeter data. The agreement is particularly good in the eastern part of the basin, where both models, our fully prognostic and altimeter data based model in Korotaev et al. [2003], indicate almost no mesoscale activity except the Batumi eddy. In contrast, on 20 August 1997 the Kaliakra eddy is well seen in Figure 4c, while it cannot be seen in Korotaev et al. [2003]. Our model simulates enhanced formation and speed of the mesoscale cyclones inside the Rim Current [Kubryakov and Stanichny, 2015].

Figure 4

Model surface circulations (a) in May (climatological average for 1962–2000), (b) in September 1990 and (c) for 20 August 1997 are represented by arrows. Isohalines are plotted in colours in (a) and (b) and indicated on the colour scale, while the circulation speed is indicated by colours in (c). The approximate locations of several popular anticyclonic eddies are also shown.

Now we validate our model by comparing the calculated SST and SSS with satellite and in‐situ data. SST and SSS are extracted from the surface outermost 2 m temperature and salinity fields. The ability of GETM to simulate SST is first examined using climatological mean maps for 1982–2009 and the bias from Pathfinder SST climatology for the same period. Climatological annual mean SST distribution is shown for the entire basin (Figure 5a). Regarding the spatial distribution of SST in the Black Sea, the most noticeable characteristic is its increase in the direction from northwest to southeast for all seasons. It is worth to note the decreased SST in the center of the eastern gyre, alongside the Rim Current in the western part of the basin and central part of the Anatolian coast. The annual mean Pathfinder SST, interpolated on the GETM domain and averaged over 1982–2009, is subtracted from the simulated annual mean SST averaged over the same period (Figure 5b). When simulations are compared with the Pathfinder SST climatology, the model usually gives a bias of less than ±0.6°C. Here it is notable that the bias in the western and central basin is largely negative except for the south‐western shelf, while it is positive in the eastern basin. Differences along northern continental boundaries and the Anatolian coast are relatively larger (from −0.7 to −1.0°C). Obtaining accurate SST fields from satellite and in situ measurements is complicated for semi‐enclosed or enclosed basins such as the Black Sea, because this relatively small basin is dominated by coastline orientation, shape and topography, making the spatial resolution of a given data set an important factor [Kara et al., 2008]. The fine spatial resolution of our model allows us to describe the effect of the Rim current on the spatial SST variation, namely the lowering of the SST due to the flow of colder water from the North‐Western Shelf. Additionally, Dufois et al. [2012] reported the existence of a warm bias nearshore in Pathfinder SSTs, thus monthly SST values in some regions are not representative during specific periods (i.e., summer and peak upwelling season). Ignoring the nearshore regions, the modeled and satellite spatial patterns are well correlated indicating the model's capability to represent long‐term SST time series.

Figure 5

(a) Calculated climatological mean SST for 1982–2009. (b) Bias of model climatology compared to Pathfinder climatology.

SSS is not well observed nor modeled, and comparisons with available data and model results are presented for a few time periods only. In Figure 6a is given the basin‐wide climatological SSS for 2005–2015 based on CMEMS results. CMEMS's model data is the only systematic SSS data for the last decades. The only distinctive SSS feature in Figure 6a is the salinity difference between the North‐Western Shelf and the rest of the basin. Our climatology for the same period (Figure 6b) indicates the presence of an inner shelf front, associated with the fresh water inflow from large the rivers on the shelf. Salinity of about 18 ‰ establishes an interface between the interior waters and those mixed with the river waters. The latter have low density and are observed as a strip around the basin. This distinguishing feature of the surface salinity distribution has been recorded and simulated previously [Besiktepe et al., 2001; Dorofeev and Sukhikh, 2016]. It appears that for this period our model underestimates the salinity of the strip, however, due to the lack of data, we could not evaluate the model error precisely. In Figure 6c is plotted the mean SSS in July 1992, which may be compared with SSS in Figure 5 in Besiktepe et al. [2001]. The SSS field plotted in Figure 6c is consistent with SSS from data‐driven primitive equation simulation in Besiktepe et al. [2001]. Differences between both SSS in July 1992 can be seen in the deeper part of the shelf and in the center of the eastern gyre. The former is mainly due to the initialization field of the data‐driven model and the latter can be associated with the higher resolution of our model. The last subplot in Figure 6 represents the 20 year SSS climatology in February from 1993 to 2012. In February the river runoff is low and river water is confined near the coast. The river waters extend down to the Bosphorus entrance region, in the form of a narrow coastal belt of current with salinities less than 14 ‰. The highest SSS are calculated in the centers of the main gyres and in the south‐eastern part about 40°E and 42°N, which are also seen in Figure 6c and to some extend in Besiktepe et al. [2001]. Figure 6d is intended for comparison with the Figure 8 in Dorofeev and Sukhikh [2016]. The overall SSS distribution is similar: high salinity in the basin interior which is decreasing toward the coast. A difference exists in the SSS distribution in the south‐eastern part, where our model gives lower SSS values except for the patch about 40°E and 42°N.

Figure 6

(a) Climatological mean SSS for 2005–2015 based on CMEMS data assimilation model. (b) Calculated herein climatological mean SSS for 2005–2015. (c) SSS in July 1992. (d) Simulated climatological mean February SSS for 1993–2012.

Figure 8

Variations of annual SSS averaged over the interior basin with depth greater than 1500 m from two different simulations. The black vertical lines represent the duration of the periods I, II and III. The black dotted line represents SSS trend from the main run (Table 2). WOA13 10 year mean SSS at locations W1 and W2 are given with triangles. CMEMS data is also plotted for comparison.

In Figure 7 surface temperature fields in the deep basin (depth deeper than 1500 m) averaged over the basin are compared with independent sea surface temperature data from Pathfinder, CMEMS and from Oguz et al. [2006]. Monthly mean SSTs from the Pathfinder and CMEMS are first interpolated on our model grid and then the corresponding deep basin averaged values are calculated. The inter‐annual and inter‐decadal variations in the Black Sea hydrography can be identified through analysis of the SST averaged over the interior basin with depths greater than 1500 m during winter (December–March) [Oguz et al., 2003 and 2006; Kazmin et al., 2010; Belokopytov, 2011]. For this reason, the deep basin average of annual mean SST (ASST) and winter mean (December–March) SST (WSST) are shown in Figures 7a and 7b, respectively. We first analyze several statistical metrics to validate monthly mean SSTs, obtained from the simulations, such as Pearson's linear correlation coefficient, R, the mean difference (MD), which is the difference between the mean values, and root mean square difference (RMSD) (Table 1). When performing metrics significance evaluation, the confidence intervals are calculated following Nicholls [2014]. The correlations of model ASST/WSST and these from the Pathfinder are good and statistically significant as the 95% confidence intervals of R are narrow. The MD and RMSD values are less than 0.65°C, displaying a reasonable agreement between model results and Pathfinder data. Concerning WSST, note the negative difference of −0.62 ±0.1°C between model and Pathfinder WSST. Even though the Pathfinder images offer good spatial coverage, we should be aware that satellite patterns are subject to certain accuracy limitations, resulting from cloud cover, sensor capabilities, calibration procedures, etc. Comparisons with the WSST in Oguz et al. [2006] demonstrate lower correlation but a statistically insignificant MD (Table 1).

Figure 7

The deep basin average of (a) annual mean SST (°C, ASST), (b) winter mean (December–March) SST (°C, WSST). Results from simulation are denoted with green circles. SST data obtained by 4 km monthly‐mean gridded AVHRR Oceans Pathfinder data set is denoted by brown squares, while that from CMEMS with orange diamonds. WSST redrawn from Oguz et al. [2006] are given with grey triangles. The black vertical lines represent the duration of the periods I, II and III. The horizontal dotted lines show the mean values (Table 2).

Table 1

Statistics of Deep Basin Averaged ASST and WSSTa

	Model and Pathfinder, (1982–2009)	Model and WSST in Oguz et al. [2006], (1960–2001)
ASST  R	0.93 (0.85, 0.97)
MD (°C)	0.28 (0.18, 0.38)
RMSD (°C)	0.37 (0.35, 0.39)
WSST
R	0.93 (0.86, 0.97)	0. 65 (0.43, 0.8)
MD (°C)	−0.62 (–0.72, −0.52)	0.012 (–0.14, 0.17)
RMSD (°C)	0.67 (0.63, 0.7)	0.51 (0.49, 0.52)

The 95% Confidential Intervals are Given in Brackets.

We further examine calculated SSS from the two different model runs, the main run starting in 1958 is compared with run2 started in 1979. Run2 is initialized with MEDAR thermohaline climatology appropriate for the early 80s. The comparison between both runs is designed to study the sensitivity of model SSS to the starting year and to identify any model drift. Long‐term evolutions of annual SSS from both runs in the deep interior basin are shown in Figure 8. The SSS correlation between the runs is 0.86 with the 95% confidence interval (0.74, 0.93) and p‐value 0 but a quite significant and constant bias is observed. Since the forcing conditions of both runs are identical, it turns out that the difference is largely due to the salinity field used for initialization. Obviously, the MEDAR climatology used for the initialization of the second run does not match the simulated salinity by the main run at the run2 starting point. Further simulations do not overcome the initial difference and subsequently the mean difference between the main run and run2 is −0.36 (‐0.4,‐0.34). Since both runs display similar oscillations and notably oscillation amplitudes overlap at the end, we can conclude that there is no clear evidence of a model drift.

Due to the lack of consistent SSS data we could not validate precisely the SSS evolution (Figure 8). In support to our findings, Shokurova et al. [2011] reported that per decade the maximum salinity values were measured in the mid‐1960s and in 1990s. The SSS displays gradual decrease over the last 50 years [Oguz et al., 2006; Kazmin et al., 2010; Belokopytov, 2011; Shokurova et al., 2011]. WOA13 10 year mean SSSs at W1 and W2 (see Figure 1 for their locations) display decreasing trend, too. Note, the lower SSS in W2 than in W1 as simulated herein (Figure 6d). CMEMS annual mean SSS exhibits strong oscillation with substantially lower amplitudes from 2000 to 2005. For the last 5 years of study, the agreement between CMEMS data and our simulations is good (Figure 8).

Intermediate Layer Properties

In addition to validation of the surface mean fields, a further validation of the T/S vertical structure is presented in Figure 10. Simulated profiles are compared with data collected by 2001 R/V Knorr cruise (http://www.ocean.washington.edu/cruises/). Data for two locations in the deep basin exist: A1 is from the Rim Current area and A2 is located in proximity to the western gyre center (see Figure 1). Since the observational data is taken once in a particular location and there is no information about the errors of the measurements, we present only an illustration of the data‐model inter‐comparison without any statistical metrics. Modeled vertical T profiles agree well with the data, in both qualitative and quantitative comparisons. However, the simulated upper halocline position does not correspond perfectly to the observational evidence for the selected time and location. It is important to note that T and S vertical profiles in Figure 9 are attained from a 41 year fully prognostic simulation. The most challenging issue of any numerical model concerning the Black Sea is to calculate accurately the extremely large vertical differences in salinity, temperature and density and to resolve the halocline position precisely.

Figure 10

Time series from simulations and observations: (a) temperature averaged over the interior basin at levels of 50 m and 200 m depth; (b) salinity averaged over the interior basin at levels of 50, 100 and 200 m depth from the main run and run2. Reanalysis data redrawn from Knysh et al. [2011] is given with filled symbols.

Figure 9

Vertical profiles in 23 May 2001 at A1 and in 25 May at A2: (a) temperature and (b) salinity. Data profiles are from 2001 R/V Knorr cruise.

Additionally, T and S averaged over 50 m and 200 m depths as well as S averaged at 100 m depth have been analyzed (Figure 10). These layers contain the most stratified part of the pycnocline and the reanalysis data of Knysh et al. [2011] is available for comparison from 1971 to 1993. Salinity at 100 m and temperature at 200 m from both runs agree well with the reanalysis data of Knysh et al. [2011], while closer to the surface (50 m) the agreement is not as good. Firstly, we have already shown the high inter‐annual and spatial variability of S in the surface layer. Secondly, the reanalysis is based on data which is heterogeneously distributed in time and space. Thus we could not expect a perfect agreement. Additionally, the density of the data points used for reanalysis is low in the center of the cyclonic gyres where the higher salinity levels occur in the upper halocline. In contrast, the low difference at 100 and 200 m is an important indicator of the model's ability to correctly describe the Black Sea's thermohaline structure.

A detailed validation of the mean model states against independent observational or model mean data indicates

Using ERA‐40/ERA‐Interim atmospheric forcing, the large‐ and mesoscale surface circulations are well represented.

Mean SSTs are close to satellite and in‐situ estimations.

The model represents observed SSS variations to an acceptable standard, however the mean SSS validation is based on scarcer quantitative and qualitative comparisons.

Thermohaline initialization affects mainly SSS and less SST, T and S at the intermediate layers.

Mean vertical thermohaline structure is kept in acceptable bounds during the 56 year run. Specifically, the very large vertical differences in salinity, temperature and density in relatively thin layers are reasonably tracked by the model.

Key Factors for SSS Variation

The deep basin mean SSS has a significant long‐term decreasing trend (Figure 8). The tendency toward dilution of the upper layer is usually attributed to the increase of fresh water input. However, at annual scale, we fail to find a linear correlation between SSS and freshwater input, nor independently with river runoff, evaporation or precipitation. In the same time the freshwater input exhibits a weak increasing trend of 0.46 mm yr−1, so the SSS response to freshwater input might be delayed due to the influence of hydrological conditions. In several studies [Oguz et al., 2006; Kazmin et al., 2010] it has been suggested that a combination of stronger surface cooling and intensified wind stress in cold years usually lead to substantial SSS increase, whilst mild years/winters are correlated with SSS decrease. Oguz et al. [2006] found the temporal variations of the upper 200 m layer‐average salinity anomaly to reflect the climatic signature in terms of phases and duration, and are inversely related to the temperature variations shown in Figure 7b. On the other hand, the linear correlation between SSS and ASST is −0.457(0.084, 0.717) and p‐value 0.019. Even though the p‐value is less than 0.05, which is the upper limit for the 95% confidence interval, the calculated interval itself is wide and we cannot consider the correlation as statistically significant. The correlation between WSST and SSS is even poorer. Thus, we did not find a linear correlation between SSS and severity of annual (winter) weather, so below we will present an alternative explanation for the salinity evolution.

Figures 11 and 12 are designed to illustrate the annual variation of SSS in two particular years: 1993, which is the year with the lowest WSST (6.92°C) through the simulation period, and 2001 which exhibits one of the hottest winters in the last four decades (WSST = 9.16°C). These two particular years have been chosen to display spatially the annual SSS evolution as they are examples for SSS decrease/increase during the corresponding year. The freshwater input, which is the sum of river runoff and precipitation minus evaporation, is approximately the same in both years and it is slightly below the climatological mean (Figure 2).

Figure 11

Model monthly mean SSS distribution in 1993.

Figure 12

Model monthly mean SSS distribution in 2001.

Typically, in January–April the river runoff is almost completely entrained in the mixing area on the shelf and does not penetrate into the interior part of the sea where the salinity is higher. The cold and fresh waters are transported by the Rim Current periphery along the coast and reach the south‐eastern part of the basin in April–May where they are mixed with higher salinity surface waters in the vicinity of the Batumi anticyclonic eddy (Figure 4). Thus in January–April, a strong density gradient occurs between the near surface waters of the central and the shelf parts (see the upmost panel in Figure 11) causing an intensification of the cyclonic motion, and consequently an upwelling amplification in the inner and deeper part of the sea. This also causes a downwelling amplification over the shelf areas. In Figure 13 current speeds are shown at 28.6°E and 42.6°N within the Rim Current area (Figure 1) from the surface to 250 m depth for the years of interest. For example, in the beginning of 1993 SSS is higher than the long‐term mean value and the peripheral fresh water transport prolongs until June due to both strong density gradient and cyclonic circulation (Figures 11 and 13a) which weaken from July to November and collapse into smaller scale cyclonic structures. The weakness of the composite peripheral current system, combined with an activated peripheral anticyclonic eddies, give rise to mixing of the nearshore low salinity waters with high salinity waters from the interior and the SSS in the basin interior decreases. Thus in a year with a severe winter a decrease of SSS may occur.

Figure 13

Monthly mean contours of current speed (m s−1) at 28.6°E and 42.6°N from the surface to 250 m depth.

On the contrary, in the hot winter 2000–2001 both SSS and cyclonic circulation are low (Figures 12 and 13b). The circulation intensification onset happens in February and the current increases in March–June (to 0.3 ms−1 in the upper 50 m layers). In 2001 the Rim current was still intense even after the beginning of the water warming in summer preserving its composite structure (clearly visible in Figure 12). As a result, the density gradient between near surface waters of the deep and the peripheral areas is kept strong and the interior SSS rises. A possible explanation for the amplification of the cyclonic circulation is the fact that strong northerly winds (from the NE–NW sector) prevailed over the North‐Western Shelf in summer‐fall of 2001 (about 58% occurrence), whereas the northerly winds account for only 42% in 1993 with lower wind speed than in 2001 (Figure 14). The wind is known to be a key factor in the formation of near sea surface currents and, in particular, northerly winds over the Black Sea support basin cyclonic motion. Once more a linear correlation between annual mean values of SSS and wind speed/direction has not been established. In summary, different combinations of weather conditions can lead to substantial increase/decrease of the deep basin SSS.

Figure 14

Wind roses for eight directional and five wind speed magnitude bands, representing the probability of occurrence in percentage of winds with given direction and speed estimated for the North‐Western Shelf from August to December: (a) in 1993; (b) in 2001. Data are extracted from the atmospheric forcing wind speed and direction at 10 m height.

Long‐Term Changes

Surface Properties

In fact, the overall (56 years) trends in winter and annual SST are not significant due to high p‐values, while SSS exhibits a significant negative trend (Figures 7 and 8, and Table 2). Several authors (Ginzburg et al., 2004; Oguz et al., 2003; Oguz et al., 2006; Kazmin et al., 2010; Zatsepin, 2007] suggest to cut the long‐time period into several shorter intervals including both warm and cold sub‐intervals. The interval duration varies through the studies but it is approximately within 5–10 years. Applying nonlinear regression analyses to ASST, WSST and SSS, the following break points can be identified: 1972 and 1997 for ASST; 1971 and 1998 for WSST; 1970 and 1995 for SSS. Consequently, three specific periods 1960–1970 (I), 1970–1995 (II) and 1995–2015 (III) have been discerned. For each of the identified periods we discuss the variation in thermohaline properties separately.

Table 2

Mean Values and Trends of Deep Basin Averaged SSTs and SSS From the Main Run in Three Certain Periodsa

Periods	I	II	III	1960–2015
ASST (°C)	16.45 ± 0.33	15.69 ± 0.23	16.10 ± 0.21	15.98 ± 0.15
WSST (°C)	8.89 ± 0.37	8.01 ± 0.29	8.36 ± 0.25	8.3 ± 0.17
SSS (‰)	18.44 ± 0.13	17.84 ± 0.07	17.53 ± 0.1	17.83 ± 0.1
ASST (°C yr−1)	−0.051 (0.4)	−0.033 (0.004)	0.052 (0.0013)	−0.0024 (0.62)
WSST (°C yr−1)	−0.051 (0.4)	−0.039 (0.007)	0.047 (0.02)	−0.006 (0.26)
SSS (‰ yr−1)	−0.048 (0.009)	−0.007 (0.06)	−0.01 (0.2)	−0.02 (0.)

In Brackets After the Trend Values are Given the p‐values of the First Order Coefficient in the Linear Fit Considering 95% Confidence Limit

The first period is a relatively warm period, in which WSST varies in the range 7.46–9.61°C and ASST in 15.4 ‐ 17.6°C. The simulations show that the first period contains more years with positive SST anomalies above the long‐term mean (Figure 7a) and calculated negative trends are insignificant for both ASST and WSST (Table 2). SSS displays negative trend in accordance with the existing hypothesis that SSS is negatively related to the severity of winters/years [Oguz et al., 2006; Kazmin et al., 2010]. Below it is illustrated that SSS decreases as a consequence of main cyclonic system disintegration. The second period of SSTs can be considered to be the coldest period, characterized by a range of WSST variation between 7 and 9°C, (ASST between 14.7 and 16.2°C) and is including five extremely cold years [Titov, 2003], namely 1976, 1985, 1987, 1992 and 1993. Within this period lie the lowest temperatures in the 56 year period of the study: in 1987 the ASST reaches minimum of 14.7°C (simulations) and 14.2°C (Pathfinder); in 1993 the WSST reaches the minimum of 7.0°C (simulations) and 7.2°C (Pathfinder). In period II both SSTs display a significant drop (II, Table 2). The mean salinity in II has been reduced substantially in comparison with I, however SSS experiences a weaker negative trend in II. The third period is a warm one that comprises a strong warming phase from 1998/1999 to 2001, followed by a sharp drop in 2002, and again a warming phase since 2009. The highest WSST is calculated in 2001, from which we extrapolate that it was an extremely hot winter. The period is not only dominated by persistent positive SST anomalies but the mean ASST of the cooling phase 2002–2008 (16.04 ± 0.26°C) is close to the long‐term mean values. Contrary, the estimated mean WSST is 8.08 ±0.13°C of the cooling phase is below the long‐term mean. However, the low winter temperatures in several years could not balance higher summer temperatures and a steep positive ASST trend is calculated for the third period (III, Table 2). For SSS the last period is characterized by a growth of amplitude and frequency of SSS oscillations. It is worth to note the highest mean values of all the three variables appear in period I, followed by the lowest SSTs in II period. Finally in III SSTs turn to rise while SSS continues to decline. Established SSTs trends do not persist for the entire simulation period. Due to a significant inter‐decadal SST variation, there is a high variance in the trends reported [Ginzburg et al., 2004; Oguz et al., 2003; Oguz et al., 2006; Kazmin et al., 2010, Zatsepin, 2007; Shapiro et al., 2010]. Moreover, trends based on short records are very sensitive to the beginning and end dates of the recordings and cannot in general reflect long‐term climate trends. The SSS decreases in the individual periods and on the whole in entire period. For comparison, monthly average salinity observations (1954–2008) in the surface layer of the west half of the abyssal part of the Black Sea has been analyzed by Belokopytov [2011] and a decrease of about 0.04 ‰ per decade is established. However this analysis is based on scarce data, which is further compromised by a missing period, between 1993 and 2005. Even if our model might overestimate the magnitude of the SSS trend, the trend is clearly present in the last 5–6 decades.

Further we examine how surface circulation fluctuates and explore the possible parallels between the strength of the Rim Current and the variation of SSS on decadal or multi‐annual scale. In order to achieve this, the climatological mean surface circulations are calculated for the periods separated by the SSS break points (Figure 15). Seasonal variability of the Black Sea eddies has been studied using hydrological data [Oguz et al., 1993], altimetry data [Korotaev et al., 2003] and numerical simulations [Staneva et al., 2001]. There is no clear evidence backing seasonal and on inter‐annual variability of the Black Sea circulations, because of the use of different data sets and different methods of analysis. One general assumption, which is not confirmed by all studies, is the strengthening of the basin scale cyclonic circulation in winter and the intensification of anticyclone activity in the warm period. Typically in September lots of anticyclonic eddies are formed between the shelf edge and the Rim Current. Climatological surface circulations in February and September for the three SSS periods are shown in Figure 15. The colour bar is fixed in the interval 0–0.6 m s−1 for all figures in order to make more apparent the comparisons of the current speed. Additionally, the scale of 0.5 m s−1 arrow length is also given in the figures.

Figure 15

Climatological model surface circulations (a) in February and (b) in September for the period I; (c) in February and (d) in September for II; (e) in February and (f) in September III. The colour bar represents the surface speed, while arrows show both speed and direction.

The first period is characterized by the lowest circulation activity – both cyclonic and anticyclonic (Figure 15a and 15b). In February the maximum current speed is about 0.3 m s−1 and the main current is separated in several smaller size cyclonic eddies. The effect of low surface circulation on the SSS is twofold: it promotes the mixing of low salinity shelf waters and higher salinity surface waters of the basin interior and reduces the upwelling in the basin interior. Hence, SSS will decrease considerably (Table 2). Formation of anticyclonic structures is present in September. Noteworthy examples include the Sinop eddy, as well as the Batumi and Sevastopol eddies.

In the second period the peripheral cyclonic system is composite and strengthened in both cold and warm season (Figures 15c and 15d). In February the mean current is almost completely separated in two main gyres with several small and powerful central cyclonic eddies. Two anticyclonic eddies along the north‐eastern shelf (red arrows in Figure 15c) are present in February as well. The speed of the large‐scale circulation increases in September, accompanied by the movement and intensification of the anticyclonic eddies. The Sevastopol eddy is well visible in Figure 15d, while the Batumi eddy is flattened between the main current and a strong cyclonic eddy in the easternmost part of the basin. Moreover, there is an evidence of anticyclonic eddy formation inside the main current in the warm period. Consolidation of the main current prevents the horizontal mixing of the surface waters, while the formation of anticyclonic eddies amplifies it. SSS decrease slows in this period (Table 2).

The overall impression of the climatological surface circulations in the last period (Figures 15e and 15f) is that the Rim Current is the most strengthened accompanied by extremely strong anticyclonic circulations. Note the maximum speed of both circulation types exceed 0.6 m s−1. The composite movement of the Rim Current is the main reason for the distinct difference between salinity in the Rim Current interior and exterior (Figure 6c). Since SSS exhibits significant oscillations in the last decades (Figure 8), there is a need to identify the main reasons for SSS oscillation. In Figure 16, are given the climatological surface circulations in February for 1990–1995 and 1995–2000. In the former period SSS increases and in the latter it decreases. Comparing both surface plots we can distinguish two main differences: the Rim Current has lower speed in the latter period and higher anticyclonic activity due to weakening of the main circulation. A key factor for the SSS decrease in this period is the intensification of the Batumi anticyclonic eddy, where the fresh waters from the North‐Western Shelf are mixed with higher salinity surface waters from the south‐eastern part of the basin.

Figure 16

Climatological model surface circulations in February: (a) for 1990–1995 and (b) for 1995–2000. The colour bar represents the surface current speed, while arrows show both speed and direction.

In summary, multi‐decadal variability of the Black Sea surface circulation indicates consolidation and amplification of the main cyclonic motion in the recent decades. Additionally, a large intensification of the anticyclonic motion is simulated in the periods of main current weakening.

Water Column Properties

Despite recent increase in the observational data and research of the Black Sea thermohaline structure (Murray et al., 1991; Özsoy and Ünlüata, 1997; Ivanov et al., 2000; BSC, 2008; etc.), we believe that our long‐term fine‐resolution simulations will allow to refine and complement previously established concepts. The temporal behavior of temperature and salinity averaged over the deep sea basin is the typically established way of studying thermohaline and dynamical fields in the basin interior (depth > 1500 m). The evolution of seasonal, inter‐annual and multi‐decadal variations in temperature and salinity from 0 to 250 m for the 56 year period is shown in Figure 17. To facilitate reading, the temperature color bar is set in the range 6–9°C and the 8°C isotherm is shown with a red dashed line. The deeper water masses possess almost vertically uniform S and T of about 22‰ and 9°C, respectively, while important variations of T/S occur in the upper 250–300 m. Temperature is seasonally variable at the surface, decreasing with depth until reaches the coldest intermediate layer. Below the CIL the temperature increases to 9°C at about 300–400 m depth and it remains unchanged until the bottom. Seasonal variation of the near surface salinity is definitely visible, namely, salinity of the mixed layer (ML) is lower in summer‐fall than in winter‐spring (in contrary to the near surface temperature). Over the year, there can be seen the main processes forming the thermohaline structure: winter water cooling, a new CIL formation, ML, summer water warming, seasonal thermocline, the increase of CIL temperature by fall, and the broken continuity of the cold layer in several years (i.e., the CIL disappears in some deep regions). Figure 17 involves perceptibly expressed mean CIL location depths. The upper and lower bounds of CIL are identified by the location of the 8°С isotherm. It can be seen that in several winters the CIL depth decreases mainly due to stronger heat loss. It is deepened as a result of surface water warming, vertical advection, and the turbulent diffusion of heat. By fall, the temperature of the CIL increases on average. Thus, the CIL thickness decreases and its depth increases from April to October due to atmospheric heating. The mean ML depths variates approximately in the range from 5 to 35 m. We will focus now on the changes in the thermohaline structure as a consequence of changing climate conditions.

Figure 17

Deep basin model mean contours of temperature (°C) and salinity from the surface to 250 m depth for 1960–1970, 1971–1981, 1982–1992, 1993–2003, 2004–2015. The 8°C isotherm is shown with a red dashed line.

The multiyear variability of T and S in the most stratified segment of the pycnocline, which is located between the 50 and 200m depth, clearly indicates strengthening of the stratification in 1960–1970 (first panel in Figure 17). CIL is largely destroyed in 1962 and 1966, since the deep basin mean temperature of the intermediate layers are above 8°C. After increasing in the first decade, T and S decrease in the following one. The mean CIL depth is almost 100 m without significant seasonal and inter‐annual fluctuation. The only one exception took place in the beginning of the 1976 when the CIL shifted up owing to considerable heat loss. The third and fourth decades are associated with severe winters and consequent intensification of the vertical mixing leading to deepening of the pycnocline. The mean depths occupied by CIL display a clear seasonal cycle, namely, they are about 75 m in cold season and deepened in the warm periods. The warming of the Black Sea surface waters since the end of the 20th century went hand in hand with the warming of the CIL waters. During several years, the seasonal cycle of salinity in the ML is missing (e.g., in 2001–2002 and 2013–2015). For the same years a significant disintegration of the CIL is found. Not only the area covered by CIL shrinks but also the CIL thickness decreases. Several authors [Oguz et al., 2003; Stanev et al., 2014] suggested CIL disappears due to the atmosphere warming, for example in 2001 or 2011 because in both years the WSST is higher than 9°C. Our simulations (Figure 17) coincide with the experimental and previous modelling evidence that CIL temperature increases considerably and less area is covered with intermediate layers of temperature lower than 8°С in the last decades. On the other hand, in 2001–2002 the cold layers, with temperature less than about 8.5°С, remain confined below the seasonal thermocline and are refilled in the subsequent year. Note that in 2010–2011 and 2014–2015 the layering is almost destroyed. Particularly in fall, waters with temperature of about 9°C occupy approximately the total water column below the thermocline. There is an important difference between the thermohaline characteristics in the beginning and at the end of the simulation period mainly associated with the last SST strong warming trend (Table 2). Obviously, the upper layer thermohaline characteristics are modified as a result of the CIL temperature rise.

In addition to the SST and SSS trends, we analyze the multiyear variability of annual mean T/S at distinct depths to see if there is an alteration of the T/S trends at these depths (Table 3). Indeed, the temperature at the surface does not display an obvious linear trend (Table 2). At 50 m depth, it exhibits a strong positive linear trend (Figure 10a and Table 3), whereas at 200 m the trend is weaker and positive. Since the strongest positive trend is found at 50 m (approximately located at the core of the CIL), it is a signal that a warming of the intermediate layers is happening prior to the surface layer. Analysis of salinity evolution at intermediate depths indicates a negative trend at 50 m depth, which weakens with depth and further alters to a positive trend at 200 m.

Table 3

Trends of T (°C yr−1) and S (‰ yr−1) at Certain Depthsa

	T at 50 m	T at 200 m	S at 50 m	S at 200 m
Trend (1960–2015)	0.009 (0.005)	0.005 (0.)	−0.012 (0.)	0.012 (0.)

The p‐values of the First‐Order Coefficient in the Linear fit in the Case of 95% Confidence Limits are Given in Brackets.

Conclusions

A high‐resolution 3D hydrodynamic model of the Black Sea is used to calculate its physical characteristics for the 1960–2015 period. It is the first successful Black Sea model based on GETM. The Black Sea peculiarities, such as the large shallow North‐Western Shelf, characterized by a huge drainage basin, the very steep continental slope and the deep basin interior make it difficult to develop a successful 3D model without relaxation to observations or other data. By means of applying higher order advection schemes TVD and superbee limiter, we can maintain the realistic thermohaline structure (both vertically and horizontally) in long‐term runs. For the Black Sea the non‐diffusive properties of the scheme are especially important and resulted in significantly improved preservation of the CIL, including a permanent halocline and a better representation of water mass characteristics. Furthermore, the usage of the TVD scheme resulted in an improved representation of frontal structures and fresh water runoff along the western coastal and shelf region. Careful adjustment of fresh water balance and coupled air‐sea heat flux exchange with feedback are also important factors for fully prognostic model runs. Our model possesses fine horizontal (2'x2') and vertical resolution (70 layers), which together with downscaled realistic atmospheric forcing from ERA40 and ERA‐Interim allow us to resolve the mesoscale eddy field. The mean results of the calculations are successfully validated by comparison with other T/S data based on hydrological, drifting buoys and satellite measurements. Specifically, the comparison with the altimeter‐derived models indicates similar cyclonic and anticyclonic circulation features. Simulated characteristics are analyzed to study the evolution of the upper and intermediate layer dynamics (until 250 m) during the last 5–6 decades. They are also used to distinguish the driving mechanisms behind T/S variability and to separate long‐term trends from short‐term variations. One particularly significant result is that our model can successfully reproduce the adequate, small‐scale features of the general circulation in response to the meteorological and thermohaline forcing. Because of its extensive validation, the model results will contribute to estimate fundamental Black Sea features and their variation, since observations in the Black Sea are generally sparse, inhomogeneous and hence unable to give a systematic picture of reality.

No long‐term trend in surface temperature is identified, however a significant decreasing trend in surface salinity of 0.02 ‰ yr−1 is established. Three distinct periods in the thermohaline properties could be identified using structural breakpoint detection, namely, from 1960–1970, 1970–1995, and 1995–2015. These periods also exhibit very different current regimes, with a very weak and disintegrated cyclonic circulation in the first period, followed by periods with a strong main cyclonic circulation (the Rim Current) in both cold and warm periods of the year. A moderate, well‐pronounced Rim Current exists in the second period, which becomes even stronger in the last period (current speed >0.5 m s−1), accompanied by similar intensification of the anticyclonic eddy activity. Powerful and long‐lived anticyclones are generated during years with strong winter Rim Current.

An analysis of thermohaline properties at certain depths allows us to reveal the following irregularities. The surface layers do not display a trend in temperature yet follow a negative trend in salinity, which weakens with depth and further alters to a positive trend at 200 m. In the layers comprising the permanent pycnocline, temperature and salinity are found to be increasing over time. Detailed discussion of free salinity evolution and CIL disintegration and vanishing are believed to be the highlights of this study. The successful validation and long‐term application of the model will be incorporated in forecasts/scenarios and biogeochemical simulations.