Climate Signals in the Black Sea From a Multidecadal Eddy-Resolving Reanalysis

Ocean reanalyses are becoming increasingly important to reconstruct and provide an overview of the ocean state from the past to the present-day. In this article, we present a Black Sea reanalysis covering the whole satellite altimetry era. In the scope of the Copernicus Marine Environment Monitoring Service, the Black Sea reanalysis system is produced using an advanced variational data assimilation method to combine the best available observations with a state-of-the-art ocean general circulation model. The hydrodynamical model is based on Nucleus for European Modeling of the Ocean, implemented for the Black Sea domain with a horizontal resolution of 1/27°× 1/36°, and 31 unevenly distributed vertical levels. The model is forced by the ECMWF ERA5 atmospheric reanalysis and climatological precipitation, whereas the sea surface temperature is relaxed to daily objective analysis fields. The model is online coupled to OceanVar, a 3D-Var ocean data assimilation scheme, to assimilate sea level anomaly along-track observations and in situ vertical profiles of temperature and salinity. Temperature fields present a continuous warming in the layer between 25 and 150 m, where the Black Sea Cold Intermediate Layer resides. This is an important signal of the Black Sea response to climate change. Sea surface temperature shows a basin-wide positive bias and the root mean square difference can reach 0.75°C along the Turkish coast in summer. The overall surface dynamic topography is well reproduced as well as the reanalysis can represent the main Black Sea circulation such as the Rim Current and the quasi-permanent anticyclonic Sevastopol and Batumi eddies. The system produces very accurate estimates of temperature, salinity and sea level which makes it suitable for understanding the Black Sea physical state in the last decades. Nevertheless, in order to improve the quality of the Black Sea reanalysis, new developments in ocean modeling and data assimilation are still important, and sustaining the Black Sea ocean observing system is crucial.


INTRODUCTION
The Black Sea is the largest land-locked basin in the world with an area of 4.2 × 10 5 km 2 , a volume of 5.3 × 10 5 km 3 and a maximum depth of 2200 m (Özsoy and Ünlüata, 1997). It is connected to the Marmara Sea and Azov Sea through the straits of Bosphorus and Kerch, respectively. It is an estuarine basin, characterized by a positive net freshwater balance, mainly due to the outflow of some of the largest European rivers such as the Danube and Dniepr, and a high-rate of precipitation which in total exceeds the total evaporation most of the time over the basin (Kara et al., 2008;Volkov and Landerer, 2015). The resulting salinity of about 18 psu in the upper layer forms a strong stratification all over the basin where a saltier water of Mediterranean origin, crossing the Marmara Sea and the Bosphorus Strait, becomes the major source of ventilation for the anoxic lower layer (Ünlülata et al., 1990;Stanev and Beckers, 1999;Stanev et al., 2001). Another main characteristic of the Black Sea is the Cold Intermediate Layer (CIL) formed at the depth of the winter convection (Özsoy and Ünlüata, 1997). The upper layer circulation of the Black Sea is dominated by the Rim Current, a quasi-permanent cyclonic jet following the bottom topography which interacts with several anti-cyclonic eddies (e.g., Batumi and Sevastopol) along its pathway in the basin (Oguz et al., 1993;Korotaev et al., 2003).
The evolution of remote sensing has been crucial to understand some of the above-mentioned Black Sea processes, since it provides high temporal and spatial resolution observations (Korotaev et al., 2001). Kubryakov and Stanichny (2015) investigated the seasonal and interannual variability of the Black Sea eddies and found a relationship between the eddy properties and the intensity of the Rim Current using altimeter observations. In addition, sea surface temperature observations have helped to detect recent warming in the Black Sea as a response of climate change (Ginzburg et al., 2004;Shapiro et al., 2010;Mulet et al., 2018). However, the major challenge for studying the ocean dynamics in the Black Sea is the historical scarcity of sub-surface observations. Although this situation has been improved in the recent years with the first deployment of Argo floats after 2002 (Grayek et al., 2015), the number of in-situ observations significantly increased only after 2010. For instance, profiling floats contributed to the detection of a recent warming in the Black Sea and the reduction of the cold-water content in the CIL (Akpinar et al., 2017;Stanev et al., 2019).
Numerical ocean models represent a powerful complementary tool to investigate the three-dimensional state of the Black Sea circulation in time in absence of dense observations. Kara et al. (2005) used an eddy-resolving model to investigate the effects of ocean turbidity on upper-ocean circulation features including sea surface height and mixed layer depth. From a 56-year model simulation, Miladinova et al. (2017) revealed that temperature has a seasonal cycle at the surface, decreasing with depth down to the CIL. Next, the same simulation was used to investigate the formation and changes of the CIL and revealed that the cooling capacity of the CIL is highly variable and decreased drastically in the last decade of the simulation (Miladinova et al., 2018). Gunduz et al. (2020) related the reduced events of CIL formation in recent years to the amplified response to climate change of the Black Sea. Although current ocean models are highly sophisticated, including improvements in parameterization of physical processes of unresolved scale and incorporating numerical techniques that are optimal for ocean regions dynamically different, they still have some limitations and incorporate uncertainties from several sources (Lima et al., 2019). Therefore, they are not completely appropriate for providing accurate ocean monitoring indicators when used alone, nor to fully study the oceanic dynamics in the Black Sea.
Ocean reanalyses reconstruct the ocean state with a long integration of an ocean model constrained by atmospheric surface forcing and observations via data assimilation (Haines, 2018;Storto et al., 2019a). They provide a four-dimensional time series of the ocean state to study ocean dynamics and unravel sources and impacts of ocean variability. Ocean reanalyses can also provide initial and boundary conditions to other models as in downscaling simulations (de Souza et al., 2020) and uncoupled seasonal forecast initializations (Balmaseda, 2017). In the Black Sea, Knysh et al. (2011) conducted a pioneering investigation utilizing an ocean reanalysis. They applied a simple data assimilation scheme to ingest available in-situ observations from 1971 to 1993.
In this work we present a Black Sea reanalysis (BS-REA) that covers the altimeter era starting from 1993 until 2018. This reanalysis system has been continuously developed in the scope of the Copernicus Marine Environment Monitoring Service (CMEMS, Le Traon et al., 2019) since 2016 (Lima et al., 2020b). It is based on an eddy-resolving ocean model coupled with an advanced data assimilation scheme, which is very innovative for the Black Sea region. Here, we present a recently upgraded version in both model and data assimilation components, which we believe will help the community for a better understanding of the physical properties and dynamics of the Black Sea. Our objective is to ensure the best representation of the sea circulation and its thermohaline structure, as well as to provide more accurate ocean monitoring indicators that can help to understand the Black Sea response to climate change.
This article was organized as follows: in section 2 we outline the BS-REA configuration in detail; in section 3 we present the main characteristics of the BS-REA and discuss the results, and finally in section 4 the conclusions are drawn.

Ocean Model
The present BS-REA hydrodynamic model is configured for the Black Sea region (the Azov Sea is not included) and it is based on NEMO v3.6 implicit free-surface implementation (Madec and The Nemo team, 2016), with a horizontal resolution of 1/27 • in the zonal direction and 1/36 • in the meridional direction, and 31 unevenly spaced vertical z-levels. This horizontal spatial resolution is chosen in order to have the same cartesian resolution in latitudinal and longitudinal directions, around 3 km at the model domain latitudes, which is conformed to an eddy-resolving scale; the Rossby radius of deformation in the Black Sea is approximately 20 km (Hallberg, 2013). The BS-REA horizontal spatial domain is shown in Figure 1.
The model is forced by the ECMWF ERA-5 atmospheric reanalysis (Hersbach et al., 2020) at the surface with a 0.25 • of spatial resolution and 1-hour time frequency. The atmospheric forcing variables are: the zonal and meridional components of 10 m wind (in m s −1 ), total cloud cover (in %), 2 m air temperature (in K), 2 m dew point temperature (in K) and mean sea level pressure (in Pa). Precipitation (in kg/m 2 s) over the basin is obtained from GPCP rainfall monthly database (Adler et al., 2003;Huffman et al., 2009), from which monthly climatological means are estimated considering the period 1979-2019. The momentum, heat and water fluxes are computed at the airsea interface based on the bulk formulae originally developed for the Mediterranean Sea (Castellari et al., 1998;Pettenuzzo et al., 2010) and applied as in the Black Sea forecasting system (Ciliberti et al., 2020).
The model bathymetry is based on the GEBCO gridded dataset at 30" resolution 1 in the Black Sea basin. The bathymetry is improved around the Bosphorus Strait with a high-resolution dataset, extensively described in Gürses (2016). Once acquired the high-resolution dataset, an optimal barycentric interpolation method is used to interpolate scattered bathymetric data on the regular spatial grid. The coastline is revised to account for and properly represents the coastal peculiarities and structures in the basin by using the NOAA shoreline dataset 2 . The river locations are remapped considering the new bathymetry.
For the river runoff, we use a monthly climatological mean estimate for the period 1960-1984 and provided by the SESAME project (Ludwig et al., 2009). The total number of rivers is 72, including the major ones such as Danube, Dnieper, Rioni, Dniester, Sakarya and Kizilirmak. The Danube runoff is distributed over five grid points to better represent its major branches, i.e., Chilia, Sulina, St. George. This special treatment accounts that the Chilia is the greatest one with three sub-branches. One is located in the south, in the Romanian territory, while the other two are in Ukraine. Sulina and St. George are located in the larger Danube floodplain, which occupies around 3500 km 2 . Thus, the distribution of the Danube discharge over its three main branches follows Panin et al. (2016); the Chilia spread 52% of the total discharge, while the remaining 48% is distributed in the Sulina (20%) and St. George (28%) branches, respectively. The salinity of the river waters is assumed to be zero.
Since the current model configuration of the BS-REA has closed lateral boundaries, the Bosphorus Strait net transport is parameterised as a river by means of surface boundary conditions while temperature and salinity are relaxed to a previous estimate. The net transport is computed iteratively from a simulation and a series of assimilation runs. A first iteration, which is a simulation, adopts a monthly climatology (Kara et al., 2008) and integrates for the whole reanalysis period. Then, every following iteration imposes the net outflow corrected by E-P-R estimates from the previous one in order to balance the water budget; evaporation (E) is model-derived and depends on each integration whereas precipitation (P) and river runoff (R) are monthly climatology as described above. In the BS-REA, a final correction, estimated from the CMEMS SSALTO/DUACS Delayed-Time Level-4 sea level anomalies measured by multi-satellite altimetry observations (Taburet et al., 2019), is applied to the freshwater balance at a single grid point adjacent to the Bosphorus Strait in order to impose Frontiers in Marine Science | www.frontiersin.org the observed trend and variability in the mean sea surface height. T and S are relaxed toward a monthly climatological profile computed from a high-resolution multi-year simulation (Aydogdu et al., 2018), to properly represent the water mass properties exchanged between the Mediterranean and Black Seas via the Bosphorus Strait. This relaxation is applied at five grid points surrounding the location of the Bosphorus Strait with a time frequency of 1 hour.
Finally, we restore the SST over the basin to the gridded CMEMS SST product (Buongiorno Nardelli et al., 2013). The restoring is done by added a damping term to the surface heat flux with a constant coefficient of dQ/dT = −200 W/m 2 /K.

Observations
The BS-REA assimilates sea level anomaly (SLA), temperature and salinity observations. The specific products assimilated are: (i) in-situ T/S profiles from both SeaDataNet 3 (Pecci et al., 2020) and CMEMS NRT in-situ product (von Schuckmann et al., 2016) and (ii) along-track sea level anomaly from all available missions, pre-processed and distributed by the CMEMS Sea Level TAC (Taburet et al., 2019). For SLA assimilation, the choice of the mean dynamic topography (MDT) is a key point and can impact the quality of results (Yan et al., 2015). In BS-REA, a modelbased MDT is computed using a 20-year (1993-2012) mean of sea surface height derived from a model integration with the assimilation of T and S profiles only.
The in-situ instrumental errors assume different values for T and S and vary in the vertical dimension based on statistics derived from Ingleby and Huddleston (2007), whereas the in-situ representation errors vary horizontally on the model grid according to previous model statistics with respect to observations and adopt same values for T and S. Both components of in situ errors are constant over time. For SLA observations, the instrumental error is set to 4 cm, and the representation errors monthly and spatially vary following Oke and Sakov (2008).

Data Assimilation Scheme
The data assimilation scheme is the OceanVar (Dobricic and Pinardi, 2008;Storto et al., 2011), a three-dimensional variational (3D-Var) assimilation algorithm. The 3D-Var scheme aims to iteratively find an optimal analysis field, x a , that minimizes a cost function (Eq. 1).
, where x is the unknown ocean state, equal to the analysis x a at the minimum of J, x b is the background state, d = y − H (x b ) is the misfit between an observation y and its modeled correspondent (in the observation space) where H, the observation operator, maps the model fields at the observation location. The method accounts for the background and observation uncertainties through the error covariance matrices B and R, respectively. The observational error covariance matrix R is diagonal in the observation space and includes the sum of instrumental and representation errors and an error component according to the time of each observation with respect to the analysis time, i.e., the observation error is multiplied by a weight depending on the absolute temporal distance between observation and analysis.
OceanVar was originally developed for the Mediterranean Sea (Dobricic and Pinardi, 2008) and later extended to global ocean applications (Storto et al., 2011(Storto et al., , 2014. In OceanVar, in order to avoid the inversion of the B matrix and to precondition the minimization of the cost function, the B matrix is defined as B = VV T , in which V is decomposed in a sequence of linear operators: Hence the V operator is used to model the background error covariance matrix and includes correlations among variables and each of its linear operators are described below. In addition, a new control variable v is used for the minimization step by considering the transformation v = V + δx and thereby δx = Vv; the superscript "+" indicates a generalized inverse. The inclusion of the control variable in Eq. 1 results in a rearranged cost function, as follows: Thus, the variational cost function is solved with the incremental formulation (Courtier, 1997) and the preconditioning of the cost function minimization is achieved through a change-of-variable transformation from the physical (Eq. 1) to the control space (Eq. 2).
OceanVar is a multivariate scheme, i.e., the state vector, x, can contain the following model state variables: T, S, SLA u and v. However, only the first three variables are employed in the present BS-REA implementation; each control vector element is a linear combination of SLA, T, S. The assimilation of in situ profiles includes a background quality-check according to Eq. 3, which rejects observations in the case the square departure from the background (d 2 ) exceeds the sum of the background (σ 2 b ) and observation (σ 2 o ) error variances by a threshold value (α). This threshold is currently set to 11 for both S and T.
For the minimization of J, the balance of the two terms in Eq. 2 defines the shape and magnitude of the analysis increments. The V v operator consists of background-error T and S vertical covariances that are extracted empirically from a model integration with the assimilation of T and S profiles using the full model resolution; the same above-mentioned integration that is used to compute the model-based MDT. The daily temperature and salinity anomalies with respect to the monthly mean are calculated to generate a set of monthly EOFs (Empirical Orthogonal Functions, only the first 15 modes are retained). V h represents horizontal correlations that are modeled through a first-order recursive filter (Farina et al., 2015), with a fixed correlation length-scale of 20 km. Determined by V η , the SLA is covaried with T and S through a balance model (dynamic height) that imposes local hydrostatic and geostrophic balance among T, S, and SLA increments (Storto et al., 2011), according to the equation: where δη and δρ are, respectively, the sea level anomaly and density increments, so that δρ is integrated in the vertical from the bottom depth h b to the surface. The ρ(T, S) relation is calculated with the 1980 United Nations Educational, Scientific and Cultural Organization (UNESCO) International Equation of State (IES 80; Fofonoff and Millard, 1985). We assume a "level of no motion" at 1000 m, which corresponds to the depth h * where horizontal velocities are considered practically zero. This implies, through geostrophy, that the corresponding pressure increment δp * h vanishes too, which results in the equation: Once the analysis increments are computed with OceanVar, the method of incremental analysis update (IAU) is used to spread the analysis increments in the first time-steps during the model initialization (Bloom et al., 1996). As a further reading on the data assimilation scheme, we refer to Dobricic and Pinardi (2008) and Storto et al. (2011).

Bias Correction
All data assimilation systems are affected by biases due to imperfect numerical models, inaccurate observations and limitations in the assimilation scheme itself (Dee, 2005). From a previous experiment with the assimilation of T, S, and SLA, we detected the evolution of systematic biases in T and S over time periods with very sparse in situ observations. For example, we have noticed drifts in temperature below 300 m starting in 1996 when the number of in situ profiles drastically reduces while altimeter observations are available. Since such drift was not present in another experiment with the assimilation of only in situ profiles, we conclude that it was generated by the SLA assimilation conducted alone.
In order to prevent those drifts, BS-REA employs a large-scale bias correction (LSBC) below 300 m throughout integration. The LSBC is formulated as: We define the estimated bias b = x − x clim as the difference between the instantaneous temperature and salinity fields with respect to T and S climatologies, which is computed for the period 1993-2018 from the above mentioned experiment that assimilates only in situ profiles; dx dt denotes the T and S tendencies, whereas M (x) represents all dynamical and thermodynamical processes and boundary conditions involving T and S during the NEMO integration. The operator L is the estimator of the model bias. It consists of a low-pass spatial filter, configured to filter out spatial scales shorter than 20 km, and is formulated as a firstorder Shapiro filter (Shapiro, 1970) that uses 250 iterations. The final bias is subtracted from the tendency, as in the incremental algorithm (Bloom et al., 1996) with a relaxation coefficient of 1200 days in order to not deplete the seasonal variability.

Numerical Experiments: Strategy and Setup
Following a spin-up of 5 years (1988)(1989)(1990)(1991)(1992) with T and S assimilation, the BS-REA starts from 1993, as soon as the altimeter observations are available, with an assimilation cycle of 2-days. That is, if the model initializes at time t, the next analysis is performed at time t + 2. The observation window is 4 days centered at the analysis time, i.e., each cycle includes observations from 2 days before and after the analysis time. Table 1 summarizes the main aspects of the BS-REA configuration, which are also described in Lima et al. (2020a). For comparison, we also present a control experiment, covering the same period of BS-REA, with exactly the same set up for air-sea interaction, such as the same atmospheric forcing and heat flux correction using the analyzed SST, but without data assimilation and LSBC.

RESULTS AND DISCUSSION
In this section we present the assessment of the BS-REA. Estimated Accuracy Numbers (EAN), which include bias and root mean square difference (RMSD), are computed using the daily outputs of the reanalysis and compared to observations using a quasi-independent approach since the validation is done by comparing the daily-averaged BS-REA fields with respect to both assimilated and rejected observations. Moreover, we provide ocean monitoring indicators such as the temperature, salinity, and ocean

BS-REA Evaluation
In Figure 2, the seasonal maps of the SST bias and RMSD are shown. There is a predominance of positive SST bias all over the basin while a negative bias manifests in limited zones such as the western Anatolian coast in summer and autumn, in river influenced areas in the northwestern shelf during the whole year and in the vicinity of the Azov Sea except in spring. The BS-REA exhibits the lowest RMSD in spring, whereas the highest RMSDs are reached in summer and autumn. For instance, the RMSD exceeds 0.75 • C along the upwelling region centered at 33 • E (Sur et al., 1994;Özsoy and Ünlüata, 1997) in the Turkish coast, where we believe that overestimated surface winds from  the atmospheric dataset may intensify the upwelling events in summer and autumn. In general, BS-REA performs better than the control experiment in terms of bias and RMSD ( Table 2). The only variable where it is similar or even slightly lower is the SST which is strongly controlled by the atmospheric forcing and the SST relaxation, both of which are the same in the two experiments. For temperature, the highest RMSDs for the layer 10-100 m are 0.63 • C and 1.50 • C for BS-REA and the control, respectively.
While the control experiment has a negative bias of −0.25 • C in the upper layers, BS-REA retains a quite reduced positive bias of 0.01 • C.
In Figure 3, we show the vertical temperature error and bias profiles for different subregions in the Black Sea. The RMSD is relatively higher in the northwestern region (dark blue) that is under the influence of the Danube River where a maximum RMSD close to 2.25 • C arises around the thermocline. The other two regions with relatively large errors are the northeastern (light blue) and southwestern ones (green), which, respectively, may be related to the absences of the Azov Sea and an open Bosphorus Strait in the BS-REA configuration. Bias profiles manifest the largest discrepancies with the observations above 40 m where there is a predominance of a negative bias (above −0.5 • C), except in the northwestern region affected by a positive bias. The Hovmöller diagrams of the temperature bias and RMSD (Figure 4, upper panels) reveal a clear seasonal pattern such that the values are low (high) in winter (summer). The highest errors are in the thermocline, where the prevalence of negative biases is evident each summer from 10 m down to 60 m. There is an evident lack of in-situ observations between 1997 and 2003 in the Black Sea which limits the reanalysis system to be constrained only by altimeter observations for a long time period.
BS-REA shows significant improved skills also for salinity with respect to the control experiment ( Table 2). The RMSD is reduced in the entire water column from 0.66 PSU to 0.41 PSU (0.77 PSU to 0.16 PSU) in 0-10 m (10-100 m). The bias is also decreased from surface down to 500 m, mainly in the layer 10-100 m where the BS-REA bias is 0.001 PSU. BS-REA presents a slightly negative bias of −0.02 PSU in 0-10 m.
The vertical error profiles for different regions show that BS-REA represents with less quality the salinity above 20 m, in particular in the northwestern region (Figure 5). In this region, the RMSD overcomes 1.2 PSU at surface and the bias reaches −0.2 PSU at around 10 m depth. This is probably due to the limitation of imposing monthly climatological runoff such as  for the Danube River, which may cause a poor representation of salinity close to the river mouth. Unfortunately, we do not have a long and uninterrupted time series of river discharges for the Black Sea to be used for a more accurate parameterization. In deeper levels, the salinity RMSD is relatively higher only in the southwestern region in the layer 60-100 m where the RMSD exceeds 0.4 PSU. This may be due to the parameterization of the Bosphorus Strait in the current model configuration in which we relax the model towards a climatology from a model simulation. Nevertheless, bias is low in the southwestern Black Sea and comparable to other regions. Hovmöller diagrams show that both salinity bias and RMSD remain low over time (Figure 4, bottom panels). However, we note large RMSD that may exceed 1.5 PSU near the surface, mainly during some temporal intervals before 2008. The mean RMSD of sea level anomaly is 2.25 cm for BS-REA which corresponds to a reduction of ∼39% with respect to the control ( Table 2). Time series of SLA RMSD present a continuous reduction of the values during the first years of BS-REA integration. Error values fluctuate around 2 cm since 2005, whereas the control error ranges between 3 and 4 cm, sometimes exceeding 4 cm (Figure 6). Horizontal maps of RMSD reveal minor seasonal differences with the largest values close to the shelf areas (not shown), where there is a dominance of the mesoscale activities along the Rim Current. For example, relatively high errors can be found in the Crimean Peninsula, where there is a regular activity of the Sevastopol eddy, and in the southeastern region, which is related to the presence of the Batumi eddy. These eddies are quasi-stationary anticyclonic features that have been examined in the Black Sea (Kubryakov and Stanichny, 2015;Kubryakov et al., 2018).
Initially, the basin-averaged temperature time series shown in Figure 7 demonstrate a high seasonal variability in 0-25 m, with values above (below) 20 • C (7.5 • C) in most of the summers (winters). The seasonal signal weakens between 25 and 150 m and disappears below 150 m. The temperature trends are estimated in two different periods (1993-2018 and 2005-2018) and indicate an overall warming of the basin especially in the period 2005-2018 with a decreasing trend in deeper layers. The values are 0.083 • C year −1 (0.12 • C year −1 ) in 0-25 m and reduce to 0.0041 • C year −1 (0.0092 • C year −1 ) in 150-300 m for the period 1993-2018 (2005-2018). For comparison, Ginzburg et al. (2004) Degtyarev (2000) also noted a positive temperature trend of 0.016 • C years −1 in a less thick layer (50-100 m) from 1985 to 1997. After 2008, the formation and presence of the CIL is observed only in 2012 and, to a lesser extent, in 2017. Figure 8 (bottom left panel) shows the basin mean temperature anomaly with respect to a reference climatology evaluated from the same    (Mulet et al., 2018;Lima et al., 2020c) which will be discussed in section 3.3. In Figure 9, similar analyses for salinity reveal that salinity trends reduce in depth and are higher during the most recent period (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)

Ocean Heat Content
The investigation of the Black Sea ocean heat content follows Lima et al. (2020c), who defined the ocean heat content like anomalies with respect to the reference period of 2005-2014, following the equation below: with ρ 0 equal to 1020 kg m −3 and c p equal to 4181.3 J kg −1• C −1 are, respectively, the density and specific heat capacity; and dz indicate the a certain ocean layer limited by the depths z 1 and z 2 ; T m corresponds to the monthly average temperature and T clim is the climatological temperature of the corresponding month. In this study the ocean heat content is estimated as the deviation from the reference period of 1993-2014. its seasonal components. In winter and spring, the negative values of SSH extend to the easternmost coast. In summer and autumn, the western basin presents similar SSH properties while, in the eastern basin, the negative values are more restricted to the inner basin. The upper layer Black Sea circulation structures are consistent with the SSH gradients showing a seasonal variability with the only exception of the permanent Rim Current encircling the entire basin and forming a large-scale cyclonic gyre. The mean upper layer circulation develops around the Rim current together with the Batumi gyre in the easternmost basin and smaller scale eddies along the Anatolian coast. The Rim current bifurcates into two branches after the Crimean Peninsula with a smaller one recirculating in the northwestern shelf and merges back to the main branch around 30.5 • E. The Rim current accelerates along the Turkish coast around 32 • E, then detaches from the shelf and penetrates into the deep basin before going again close to the coast around 35.5 • E. In winter, the eastern and western gyres are less defined. Following the Rim Current, the Batumi anticyclonic eddy is well defined in summer, seems to be more confined near the Georgia coast in autumn, but it is less apparent in winter and spring. Next, the presence of Sevastopol anticyclonic eddy is very clear near the southwest of the Crimean Peninsula in spring and summer, whereas it is less distinguishable in winter and autumn. All these circulation patterns are consistent with the previous estimates described in Oguz et al. (1993); Özsoy and Ünlüata (1997), Korotaev et al. (2003) and Gunduz et al. (2020).

CONCLUSION
The BS-REA system shows very satisfactory skills compared to the model simulation, which highlights the importance of using data assimilation to improve the model representation. BS-REA also has the ability to represent the main Black Sea circulation, the Rim Current, as well as the mesoscale features in the Black Sea, such as the quasi-stationary anticyclones Sevastopol and Batumi eddies, respectively, near the Crimean Peninsula and southeastern region. Notwithstanding, the BS-REA has shown a reduced ability to represent the impact of Danube waters on the sea, which is possibly due to the current model configuration such as the application of monthly climatological runoff. Furthermore, the absence of the Bosphorus and Kerch straits negatively impacts the BS-REA representation in regions adjacent to the Azov and Marmara Seas.
The system is very suitable for understanding the physical state of the Black Sea in recent years and allows to obtain more accurate ocean monitoring indicators for the sea, which are important to understand its response to climate change. The temperature analyses have indicated a recent faster warming of the Black Sea that has impacted its CIL formation. Since 2009, the disappearance of CIL is evident, although some weaker CIL sporadic events are detected in 2012 and 2017. Additional investigations show a relative reduction in the ocean heat content in these years, which coincides with the reemergence of the CIL.
Trends in temperature, salinity and ocean heat content reveal a warming and salinification of the Black Sea, especially in the past few years. However, since trends based on short records are very sensitive to the beginning and end values of the time series and cannot in general reflect long-term climate trends, longer time series are needed to confirm these tendencies. This requires a continuous improvement of the BS-REA system through new developments in ocean modeling and data assimilation together with the maintenance of the Black Sea ocean observation system. In addition, for future work, we consider comparing our results with global models such as those from the Ocean Model Intercomparison Project (Lin et al., 2020;Chassignet et al., 2020) and global ocean reanalyses (Storto et al., 2019a,b), which can also allow us to quantify uncertainties through an ensemble of model results in the Black Sea.
In order to further improve the reanalysis, the next generation of the Black Sea systems will include a revision of the hydrodynamical core and new capacities from the data assimilation scheme. Regarding the core model, the new version will use higher resolution in vertical (e.g., from 31 to 121 z-levels with partial steps) and upgrade to NEMO v4.0. The Bosphorus Strait is going to be represented as an open boundary thanks to the inclusion of the Marmara Sea box in the numerical grid: it will ingest the high-resolution model solutions provided by the Unstructured Turkish Straits System (U-TSS, Ilicak et al., 2021) -T, S, SSH, U, V -with the scope to optimally interface the Black Sea with the Mediterranean Sea. Such new developments, together with the revision of the land forcing and data assimilation scheme to account for high resolution EOF, will be part of the new Black Sea forecasting system  that entered in service in May 2021 and will be uptaken by the BS-REA in the near future.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://doi.org/10. 25423/CMCC/BLKSEA_MULTIYEAR_PHY_007_004.

AUTHOR CONTRIBUTIONS
LL led the work, developed the reanalysis system, and contributed to the elaboration of this work at all levels. SCi led the hydrodynamical model development and provided a valuable scientific contribution to improve this work at different levels. DA, SCa, and MI contributed to the hydrodynamical model development. SM contributed to the scientific scope, providing important insights to improve this work at different levels. AA contributed to improving the model setup, particularly for the correction of the model freshwater budget and provided a valuable scientific contribution to improve this work at different levels. RE contributed to define the validation strategy. AC contributed to set up the bias correction scheme. AA, RE, AC, and EJ contributed to define the data assimilation strategies. RL provided the access to all available observational dataset. SCr, LS, and FP contributed to set up the operational procedures and interfaces. EP, GC, and EC provided useful comments that helped to guide this work. All authors contributed to the article and approved the submitted version.