Sensitivity of the Baltic Sea overturning circulation to long-term atmospheric and hydrological changes

The sensitivity of the overturning circulation in the Baltic Sea is analyzed with respect to long-term changes in atmospheric and hydrological conditions by using two state-of-the-art ocean circulation models: RCO (Rossby Centre Ocean model) with a reference simulation and various sensitivity experiments as well as MOM (Modular Ocean Model). Historical reconstructions since 1850 lasting for more than 150 years are considered in order to identify coherences between the overturning stream function and surface wind, river runoff, major Baltic inflows, salinity and water temperature. Long-term time series are evaluated statistically for several sub-basins of the Baltic Sea for the two different models and for the sensitivity experiments concerning the inter-annual, multi-decadal, and centennial variability. We found that the simulated overturning circulation has a response time scale with respect to wind or runoff of around 30 years. The overturning circulation will decrease basin-wide under anomalous wind conditions which hamper the deep water flow within the Baltic Sea, under a warmer climate or when river runoff increases. However, a global sea level rise would reinforce the overturning circulation. Overall, multi-decadal variations of the overturning circulation are anti-correlated to and caused by the wind parallel to the cross-section of the overturning circulation. They are not caused by river runoff and the coherence with major Baltic inflows is small. Hence, the overturning circulation does not act as a good proxy for major Baltic inflows. river runoff, and sea level. Changes in wind or runoff affect the circulation after ~30 years in both models. The circulation reinforces when global sea level rises, but it decreases when the wind hampers the water mass exchange with the North Sea, when river runoff increases or under warmer climate. These counteracting responses allow no clear projection of the future overturning circulation of the Baltic Sea. Overall we found that pronounced 30-years variations of this circulation are caused by wind hampering surface water outflow of the Baltic Sea but not by river runoff and only little by saltwater inflows.


Introduction
The ocean's overturning circulation is a large-scale circulation system of surface and deep currents which transports large amounts of water, heat, salt, and a variety of constituents around the globe and impacts the Earth's climate crucially by ocean-atmosphere interactions (Schmittner et al., 2013). One of the most prominent examples is the Atlantic Meridional Overturning Circulation (AMOC; e.g. Cunningham et al., 2007;Msadek and Frankignoul, 2009;Kanzow et al., 2010) which is a basin-wide circulation in the latitude-depth-plane and is typically quantified by a meridional transport stream function. Maxima of this stream function at any latitude and depth specify the total amount of water moving meridionally above and below this depth in

Accepted Article
This article is protected by copyright. All rights reserved.
opposite directions. Instead of taking depth as a vertical coordinate there are also approaches to calculate the stream function as a function of density, salinity, or temperature (e.g., Döös et al., 2004;Lozier et al., 2019). However, drivers of the AMOC are unknown. Recently, Lozier et al. (2019) investigated forcings which control the variability of the meridional overturning circulation in the North Atlantic and found that besides buoyancy (density anomalies) the wind variability may be important on decadal time scales. Mechanisms are very likely different but a similar question about the causes of the variability of the overturning circulation might be asked for coastal seas such as the Baltic Sea with an estuarine-like circulation.
The Baltic Sea is a shallow, semi-enclosed sea in northern Europe with limited water exchange with the adjacent North Sea and with a freshwater supply of a surrounding catchment area that is about four times as large as the surface area of the Baltic Sea itself (Figure 1). Due to the large river runoff into the Baltic Sea (Bergström and Carlson, 1994) and the relatively small mixing (because tidal amplitudes are small), the Baltic Sea consists of a two-layer system with a brackish surface layer and a more saline deep water layer separated by a perennial halocline (Knudsen, 1900). The saltwater is supplied intermittently by small, medium or large saltwater inflows from the North Sea through the Danish straits. Only the so-called major Baltic inflow (MBI) events are able to ventilate the deep water in the central Baltic Sea and hence are of great importance for the environmental conditions below the permanent halocline (Mohrholz, 2018). Definitions of MBIs have developed over the past decades. First quantitative definitions were provided by Wolf (1972). In following studies, the criteria for MBIs were further developed (e.g., Franck et al., 1987;Matthäus and Franck, 1992;Fischer and Matthäus, 1996;Lehmann and Post, 2015;Mohrholz, 2018). The classical criteria for identifying MBIs were introduced by Matthäus and Franck (1992) who took a stratification coefficient and the bottom salinity at the light vessel Gedser Rev in the southwestern Baltic Sea into account. The strong stratification dominating mainly in the southwestern and central sub-basins of the Baltic Sea is connected with an estuarine-like overturning circulation which is presumably driven by the baroclinic pressure gradient, wind induced mixing, up-and downwelling, river runoff and sea level variability in Kattegat (Meier, 2005). It is noteworthy, that the total river runoff, the volume averaged salinity and the barotropic saltwater inflows show a considerable multi-decadal variability with a period of about 30 years but no statistically significant trend (e.g., Winsor et al., 2001;Meier and Kauker, 2003;Fonselius and Valderrama, 2003;BACC Author Team, 2008;BACC II Author Team, 2015;Mohrholz, 2018;Kniebusch et al., 2019).
In previous studies, mean in-and outflows of the overturning circulation were calculated from ocean circulation models (e.g., Meier and Kauker, 2003;Meier, 2005Meier, , 2007Placke et al., 2018) including their sensitivity to the external forcing such as wind, river runoff and the sea level variability in the Kattegat (Meier, 2005). Meier (2005) found that in steady state increased freshwater inflow and wind speed result in decreased salinity whereas increased amplitude of the Kattegat sea level results in increased salinity. Defining the age of water as the time elapsed since a water particle left the sea surface, the impact of changing freshwater or sea level in the Kattegat on the average age of water is comparatively small. This sugggests an approximate invariance of stability and deep water transport in steady state (Meier, 2005).
The overturning stream function (OSF) of the Baltic Sea was calculated for the first time on transects by Döös et al. (2004) using depth, density, salinity, and temperature coordinates. They found a closed overturning cell from the Kattegat up to the northern Baltic proper (i.e., the northern Gotland Basin; for the locations see Figure 1). With the help of passive tracers, Meier

Accepted Article
This article is protected by copyright. All rights reserved.
(2007) identified the three-dimensional pathways of inflowing high-saline and freshwater parcels from the Danish straits' entrance area and from the rivers, respectively. On average, the vertical overturning circulation of the inflowing saltwater is closed by rather patchy upwelling in the area of the Landsort Deep in the northwestern Gotland Basin according to model simulations (Meier, 2007). Seasonal changes of the outflowing surface branch of the overturning cell from the Baltic Sea into the Kattegat, which contains brackish water that originates mainly from the rivers in the northern Baltic Sea, are controlled by the zonal wind, although the annual mean outflow is given by the total runoff (Hordoir and Meier, 2010). Thus, strong westward directed wind can accelerate the surface outflow while strong eastward directed wind hinders it. Hence, the interannual variability of the seasonal freshwater outflow maximum is highly correlated with the North Atlantic Oscillation (NAO, e.g., Hurrell, 1995;Visbeck et al., 2001). In idealized scenario simulations based on two regionalized climate models under the RCP 4.5 and 8.5 greenhouse gas concentration scenarios, which neglected changes in the regional water cycle and consequently in runoff, Hordoir et al. (2018) showed that the overturning circulation in the Gotland Basin will decrease at the end of this century due to the

Accepted Article
This article is protected by copyright. All rights reserved. increased thermal stratification. A similar study was performed for the Gulf of Bothnia consisting of the Bothnian Sea and Bothnian Bay (Hordoir et al., 2019). From an ensemble of regional coupled atmosphere-ocean simulations, Gröger et al. (2019) analyzed projections including both changes in atmospheric and hydrological forcing. They found that the outflows from the southwestern Baltic Sea increased and that the inflows into the deeper layers were reduced by ~25 %. In the Baltic proper, their ensemble of scenario simulations did not show a robust change in the maxima of the OSFs.
In the present study, the OSF variability of the Baltic Sea is further analyzed based upon results from two ocean circulation models in order to consider model uncertainties. The two state-of-the-art ocean circulation models are the Rossby Centre Ocean model (RCO; Meier, 2007) and the Modular Ocean Model (MOM; Neumann et al., 2002;Griffies, 2004). For both models, historical simulations are performed, starting in January 1850 and ending in 2008 for RCO (Meier et al., 2012) and in 2009 for MOM, respectively. The study addresses the following key questions: 1) How different are the OSFs in the two models?
2) Is the inter-annual variability of the OSF modified by salt transports from the North Sea, in particular during MBIs?
3) What causes the pronounced multi-decadal variability of the OSF found in the historical simulation? 4) How sensitive is the OSF to atmospheric and hydrological changes on centennial time scales and what are the controlling processes?
The manuscript is organized as follows. In Section 2 the model setups, their forcing data, the determination of the OSF, of the salt transport associated with MBIs, and of the entrainment velocity are described as well as the experimental strategy. The results from the two different ocean circulation models and the sensitivity experiments performed with RCO are presented in Section 3. There the long-term mean structure of the OSF from the southern to the northern Baltic Sea is discussed and its inter-annual, multi-decadal, and centennial variability is analyzed in relation to atmospheric and hydrological conditions within different sub-basins. In Section 4 we discuss our results. The main findings of this study are summarized and concluded in Section 5 where also prospective ideas are outlined.

Model setups
For the present investigations two ocean circulation models are used. The Rossby Centre Ocean model (RCO) is a Bryan-Cox-Semtner primitive equation circulation model with a free surface (Killworth et al., 1991) and a lateral open boundary in the northern Kattegat at 57.7°N . The subgrid-scale mixing is parameterized using a k-ε turbulence closure scheme with flux boundary conditions (Meier, 2001). In horizontal direction no explicit diffusion is applied, but a harmonic parameterization of viscosity is chosen (Meier, 2007). The simulations are based on an Arakawa B-grid with no-slip lateral boundary conditions. In vertical direction level coordinates are used. The spatial

Accepted Article
This article is protected by copyright. All rights reserved.
resolution is 2 nautical miles in the horizontal and 3 meters in the vertical. The temporal resolution of the output data is two-daily. The model setup is the same as in Meier et al. (2019) driven by reconstructed three-hourly High RESolution Atmospheric Forcing Fields (HiResAFF; Schenk and Zorita, 2012) as atmospheric forcing and a combination of various observational and model datasets as hydrological forcing. The analyses of Meier et al. (2019) comprise a reference simulation and 12 sensitivity experiments, which all have the same experimental setup but modified forcing data. For the present study the reference simulation is utilized together with a subset of 6 sensitivity experiments which take certain conditions in air temperature, wind, runoff, and sea level rise into account. These experiments are explained in detail in Section 2.4.
The Modular Ocean Model (MOM version 5.1) is a global numerical ocean model based on the hydrostatic primitive equations (Pacanowksi and Griffies, 2000;Griffies, 2004) which has been adapted to the Baltic Sea with an explicit free surface, an open boundary condition with respect to the North Sea in the western Skagerrak at 8.3°E, and freshwater riverine input (e.g., Neumann et al., 2017). For vertical and horizontal mixing the K profile (KPP; Large et al., 1994) and the Smagorinsky (Smagorinsky, 1963) schemes are used, respectively. No explicit lateral diffusion is applied. The simulations are based on vertical z*-coordinates and on an Arakawa B-grid in the horizontal plane with no-slip lateral boundary conditions. The horizontal resolution is 3 nautical miles and the vertical layer thickness is 2 meters. The temporal resolution of 2-D and 3-D data is daily and monthly, respectively. Similar to the RCO simulations three-hourly HiResAFF data are taken as atmospheric forcing for MOM. The wind forcing data for the MOM simulation are interpolated to three-hourly data from daily data with a gustiness which is up to 1.3fold of the mean wind speed. I.e., linearly interpolated three-hourly data have been multiplied by a random factor between 1.0 and 1.3 as suggested by Brasseur (2001, their Table 6) who used this simple method in weather forecasting to determine wind gusts in regions where information about wind gusts is not available. This approach improves the representation of mixing in ocean circulation models. The river runoff is based on HELCOM (2015) and has been increased by about 20 % as a result of the calibration procedure. The simulations cover the years 1850-2009. Shortcomings due to the horizontal model resolutions are compensated such that both used models are adequate for the current climate simulations. The present version of RCO with an improved flux-corrected advection scheme was in detail evaluated by Meier (2007). Although the horizontal resolution of the model is too coarse to resolve the details of the topography in the Danish straits, MBIs and the temporal evolution of salinity in the various sub-basins are well simulated compared to monitoring data. In particular, the MBIs in the years 1993  and 2014  are simulated realistically. Even smaller inflows such as the summer inflow in August/September 2002 can be tracked from the sills in the Danish straits to the Gotland Deep in accordance with available observations . In the present study, the reconstructed atmospheric forcing for 1850-2008 (Schenk and Zorita, 2012) is of course less reliable than reanalysis data such as UERRA (http://www.uerra.eu/). However, a detailed evaluation using sea level, lightship, sea ice extent, environmental monitoring and hypoxic area data nevertheless confirmed that the model results are still satisfactory . This article is protected by copyright. All rights reserved.
For the MOM setup, the same resolution-related restrictions as for the RCO setup are valid. However, the model is able to reasonable reproduce salt water inflows in long-term simulations (e.g., Neumann and Schernewski, 2008) and on an event basis (e.g. Neumann et al., 2017).

Hydrological and atmospheric forcing data
The temporal evolutions of the hydrological and atmospheric forcing data of the two models RCO and MOM are shown in Figure 2. Runoff and wind data for the reference simulation RCO-REF are equally used for all RCO sensitivity experiments except for RUNOFF (which uses a climatological monthly mean river runoff), FRESH (which has increased runoff) and for W1904 (which uses wind from a year without an MBI). The total annual river runoff of the whole Baltic Sea excluding Kattegat and Skagerrak has in MOM by definition a mean value of approximately 16900 m 3 s -1 which is 20 % higher than in the reference simulation of RCO (with approximately 14100 m 3 s -1 on average) and hence coincides closely with the RCO sensitivity experiment FRESH ( Figure 2). All other sensitivity experiments have a similar long-term mean river runoff for the Baltic Sea area. The illustrated low-frequency winds of the two models are well correlated with a correlation coefficient of 0.95 because the wind forcing in both models is based on HiResAFF but the high-frequency components of the wind are differently determined. Moreover, the underlying wind values in RCO lie on the same grid as all other parameters (like, e.g., water temperature or salinity) and are only available on wet grid points while for MOM the wind forcing is calculated from an atmosphere grid that differs from the used parameter grid in size and resolution and has wind values on each grid point. The different parameterizations and data interpolation in the models as well as the imposed gustiness in MOM result in 10 % lower zonal and 30 % higher meridional wind speeds for MOM in contrast to RCO ( Figure 2).

Accepted Article
This article is protected by copyright. All rights reserved. The wind during the year 1904 which is used (namely regularly repeated plus a minor random shift between years) in the sensitivity experiment W1904 differs from the climatological mean wind during 1850-2008 used in REF. Figure 3 illustrates that the zonal and meridional climatological monthly mean winds are predominantly eastward and northward directed with a higher standard deviation in winter than in summer. However, the experiment W1904 on the long-term average has westward directed wind at the beginning of each year and southward directed wind in summer with very small standard deviations between years.

Accepted Article
This article is protected by copyright. All rights reserved.

Result analysis
For both models the overturning circulation in the Baltic Sea, which is expressed by the OSF, is calculated from the volume transport in zonal and meridional directions, respectively. Following Döös et al. (2004) the zonal OSF is determined by with x, y, z, t being the zonal, meridional, vertical, and time coordinates, H(x,y) being the depth level as well as u and v being the zonal and meridional velocity. Hence, the volume transport in zonal direction (u dy dz) is integrated over the meridional coordinate (latitude) at each longitude and depth from the southern coasts (yS) to the northern coasts (yN) and the meridional volume transport (v dx dz) is integrated over the zonal coordinate (longitude) at each latitude and depth from the eastern coasts (xE) to the western coasts (xW). The resulting volume transport is subsequently integrated over depth from the bottom towards the depth z and then averaged between the start time t0 and the end time t1. For the present study annual averages are used. The obtained OSF (in m 3 s -1 ) is a cross-section in the longitude-depth-plane (for zonal OSF) or latitude-depth-plane (for meridional OSF). For water bodies like the Baltic Sea where flows take place partly in zonal and partly in meridional direction the zonal and meridional OSF can be calculated

Accepted Article
This article is protected by copyright. All rights reserved.
separately in the appropriate basins and can then be combined to one OSF that describes the overturning circulation through the entire sea from the entrance region to the end basin.
The regions in which the OSF is determined and evaluated are illustrated in Figure 1 ). Within these areas the meridional OSF is calculated. The maximum of OSF(x0,z0) (maxOSF(x0,z0)) or OSF(y0,z0) (maxOSF(y0,z0)) is equal to the maximum zonal or meridional volume transport at longitude x0 or latitude y0 between the bottom and the depth z0.
For investigating MBIs into the Baltic Sea we follow Matthäus and Franck (1992) who defined an MBI as an inflow event with (1) a stratification coefficient G = 1 -So/Sb being ≤ 0.2 for at least 5 consecutive days at the station Gedser Rev (with So being the surface salinity and Sb being the bottom salinity) and (2) Sb being ≥ 17 g kg -1 . In the present study we consider the salt transport associated with MBIs through the Drogden Sill (12.7°E -13.0°E at 55.6°N) and Darss Sill (54.4°N -54.9°N at 12.6°E; locations see Figure 1) by summing up all salt transports through both transects that supply water with a salinity S ≥ 17 g kg -1 .
The magnitudes of the OSF determined for the different RCO sensitivity experiments differ from those of the reference simulation. This is attributed to changes in the entrainment velocity which is differently pronounced under the certain hydrological and atmospheric forcing conditions. Entrainment is a special form of mixing which leads in a two-layer-system (consisting of an active and a passive layer) directly to an enhanced volume flux of the active layer. One of the first attempts calculating entrainment into the Baltic Sea deep water from observational data was done by Kõuts and Omstedt (1993). In the present study the entrainment velocity we of the entire Baltic Sea is calculated for the different RCO sensitivity experiments following Stigebrandt (1985, his Eq. 8): with the factor ε=0.05 for the buoyancy flux B<0 and ε=1 for B≥0. Here m0=0.6 is an empirical constant, g is gravity, Δρ/ρ is the relative density difference between the upper and lower layers, and h is the mixed layer thickness. The friction current velocity u* is determined from Eq. 4

Accepted Article
This article is protected by copyright. All rights reserved.
with the surface stress τ and the in-situ seawater density for the mixed layer ρ (i.e., at the surface) which are both given from the model. The buoyancy flux B through the sea surface is calculated following Stigebrandt (1985, his Eq. 6): with α being the thermal expansion coefficient for seawater in the mixed layer, β being the respective haline expansion/contraction coefficient, both are functions of temperature and salinity, Cp being heat capacity, and S being surface salinity. Furthermore, Qin is the incoming heat flux and F is the sum of precipitation minus evaporation and river runoff per basin area from the model output.
The thickness of the mixed layer h is determined for each wet grid cell of the model. For grid cells with a water depth ≥35 m the mixed layer depth is defined at the depth where a maximum of the vertical density gradient occurs in the thermocline range. This range is defined between the surface down to 27 m (42 m) depth for grid cells with a total water depth of 35-50 m (greater than 50 m). The maximum of the vertical density gradient needs to be greater than a certain threshold which is defined to be the depth-average over the annual mean vertical density gradient per grid cell. If no maximum exists in the thermocline range, an equal maximum is searched below, i.e., in the halocline range. For shallow water grid cells (with a total water depth <35 m) the equivalent maximum is searched within the whole vertical density gradient profile (independent of a thermocline or halocline range). If no maximum is found the mixed layer thickness is undefined such that no entrainment velocity can be calculated. This is often the case in shallow regions during wintertime.
The long-term mean entrainment velocity for RCO is then determined at each individual grid cell as follows: The annual time series of the calculated entrainment velocity of all 159 years are averaged to one mean annual time series. This 159-years mean annual time series has only valid values at time steps for which the entrainment velocity is calculated in each of the 159 years. I.e., at most grid cells in the deep regions (depths greater than 100 m where the entrainment velocity can be calculated at almost every time step throughout the year) the 159-years mean annual time series consists of almost wholeyear-data. In shallower regions the entrainment velocity can often not be determined during wintertime such that the 159-years mean annual time series has predominantly valid data only during summertime. For most of the very shallow coastal regions (0-35 m depth) even no entrainment velocity determination is possible throughout the whole year. Averaging the 159-years mean annual time series gives the 159-years mean entrainment velocity value at each grid cell.

Accepted Article
This article is protected by copyright. All rights reserved.
For the investigation of the Baltic Sea overturning circulation in combination with hydrographic and atmospheric conditions, the relevant data are processed on an annual basis. Annual means are centered on wintertime, i.e., the data with its original temporal resolution from each model are averaged from the beginning of July of one year until the end of June of the subsequent year. Hence, winter conditions which are more variable than those in summer and which drive important processes like for instance MBIs are covered in their whole temporal extent. In order to focus on the general long-term (lowfrequency) variation of parameters over the more than 150 years long simulation period, the data are low-pass filtered using 10-year running averages which are shifted by one year each. Spatial means of parameters, which serve as representative means for the whole Baltic Sea, as well as the total river runoff are calculated for the Baltic Sea area excluding the Kattegat and Skagerrak.
For finding correlations between two time series wavelet coherences are computed according to Grinsted et al. (2004) and magnitudes for the individual periods are integrated over time in order to find long-term mean results. Furthermore, it must be noted that water temperature in this study refers to potential temperature for RCO simulations, but to conservative temperature for MOM. For water temperatures and salinities typical in the Baltic Sea the difference between both temperatures is on average about 0.5°C, with larger deviations (but typically smaller than 1°C) for warm water and lower deviations for cold water (cf. IOC et al., 2010).
The sensitivity experiments performed with the model RCO take modified forcing conditions in air temperature, wind, runoff, and sea level rise into account. An overview of these used experiments performed for the years 1850-2008 is shown in Table 1.

Long-term mean overturning circulation
The long-term mean OSF along the cross-section from the southern to the northern Baltic Sea as simulated by the two models RCO and MOM is presented in Figure 4.

Accepted Article
This article is protected by copyright. All rights reserved.

REF and lower in all other sub-basins of the central and northern
Baltic. The differences in the overturning circulation between the two models (Figure 4a-b) can partly be explained by the different atmospheric and hydrological forcing ( Figure 2) with higher runoff and higher mean wind speed in MOM compared to RCO (see Section Methods and below). Due to the higher runoff, the outflow in MOM is also higher than in RCO-REF.
It must be noted that the Bornholm Basin is characterized by two vertical circulation cells which are on top of each other. The surface cell is wind-driven with an upper branch of inflowing water due to the Ekman dynamics and an outflow at medium depths due to mass conservation. The second vertical cell is driven by saltwater inflows at the bottom and has an outflow at medium depths together with the wind-driven upper cell. The resulting OSF shows a minimum at about 10 m depth where the surface inflow takes place and it reverses at about 20 m depth where the outflow happens.
The characteristic of the long-term mean overturning circulation simulated by the individual sensitivity experiments of RCO is shown in Figure 4c-h. The OSF magnitudes simulated in T1904, T2008, and RUNOFF differ very little from the reference run. For the experiments FRESH and W1904 the overturning circulation weakens when compared to RCO-REF over a wide area with positive magnitudes and strengthens in few areas where negative magnitudes occur. In contrast, the experiment GSLR shows a large-area intensification of the overturning circulation. These results imply that in experiments with warmer or colder climate compared to REF and when runoff is kept constant instead of being inter-annually variable the overturning circulation is less affected. The overturning circulation decreases basin-wide when river runoff increases or when winds prevent MBIs which hence leads to lower salinity and smaller density gradients. However, in case of a sea level rise by 1.5 m the overturning circulation intensifies.

Accepted Article
This article is protected by copyright. All rights reserved.

Inter-annual variability
The relation between maxOSF records and the salt transport associated with MBIs is investigated from Figure 5. During the analyzed period the salt transport is little correlated with the maxOSF in the Bornholm Basin and Gotland Basin with correlation coefficients of about 0.5 and 0.2, respectively. Of course, a time lag must be expected between MBIs entering the Bornholm Basin and the maxOSF in the Gotland Basin. During the period with the highest maxOSF (in the year 1961 in the Bornholm Basin and in the year 1944 in the Gotland Basin) no maxima in salt transport associated with MBI events were observed. Otherwise, the strongest salt transport in the year 1977 is only related to a minor peak in the maxOSF in the Bornholm Basin but to no peak in the Gotland Basin. Hence, MBIs and annual maxOSF maxima are not well correlated. Not all MBIs replace the deep water in the Gotland Basin and the impact of the salt transports into the Baltic Sea on the maxOSF is overshadowed by other drivers such as the lowfrequency variability of the wind (see the discussion below).

Accepted Article
This article is protected by copyright. All rights reserved.

Multi-decadal variability
In the following, we consider the temporal evolution of the maxOSF within different subbasins, which are depicted in Figure 4. Figure 6a-c shows timeseries of the maxOSF for the RCO reference simulation and sensitivity experiments as well as for the MOM simulation. In all considered sub-basins maxima occur around the year 1940, 1960 and in the 1990s. In the Bornholm Basin and the Bothnian Sea further maxima occur around 1880 and 1920, respectively. The Gotland Basin, Bothnian Sea, and Bothnian Bay show a common additional maximum in the 1970s. Note that these patterns hold for most of the simulations but not for all of them. In general, there is a time shift between the simulated maxima: those in the north occur later in time than those in the south. The low-pass filtered timeseries of the maxOSF with a cutoff period of 10 years for the different model runs per sub-basin correlate predominantly very well with the respective reference run of RCO and have correlation coefficients higher than 0.9 except for the experiment W1904 (around 0.5) and MOM (around 0.6).  Figure  6d) and the multi-decadal variations in maxOSF in the Bornholm Basin and Gotland Basin for REF are low with correlation coefficients of about 0.5 and 0.3, respectively. Moreover, the sensitivity experiment W1904 shows a low-frequency variability in salt transport by MBIs (Figure 6d) which is probably caused by the sea level variations in Kattegat, i.e., by the zonal wind field over the North Sea.

Accepted Article
This article is protected by copyright. All rights reserved.

Accepted Article
This article is protected by copyright. All rights reserved. Hence, the multi-decadal maxOSF variations in both the Bornholm Basin and the Gotland Basin are caused by the wind parallel to the cross-section of the OSF and not by the river runoff or MBIs. More eastward directed wind would imply a lower maxOSF in the Bornholm Basin because it results in a blocking of the outflow in the surface layer which consequently results in smaller inflows due to mass conservation. This is in accordance to the results by Meier and Kauker (2003) and Hordoir and Meier (2010). On the other hand more northward directed wind would imply a lower maxOSF in the subbasins in the central and northern Baltic Sea.
Within the Bornholm Basin and the Gotland Basin, the maxOSF and the mean salinity are correlated with an approximate time lag of 5 to 10 years (Figure 7) because higher wind parallel to the cross-section of the OSF and opposite to the outflow direction implies lower maxOSF and after a typical response time lower salinity. About 50 % of the mean salinity variations is explained by runoff in accordance to Meier and Kauker (2003), Börgel et al. (2018), and Kniebusch et al. (2019). Although mean salinity and maxOSF in the Bornholm Basin and the Gotland Basin are higher in MOM than in RCO the variability between the two models is similar and likely explained by the same mechanisms. The maxOSF also explains parts of the multi-decadal variations in the mean water temperature within each sub-basin because a larger transport of water masses into the deep layer would change the deep water temperature according to the temperature of the inflowing water. However, for the northern sub-basins multi-decadal variations in water temperature are low. In all sub-basins, a considerable temperature trend over the entire record caused by global warming is found.
The sensitivity experiments support our findings that the multi-decadal variability of the maxOSF is mainly driven by the wind parallel to the cross-section of the OSF. The multidecadal variability of the maxOSF is very similar for REF and the sensitivity experiments except for W1904 which differs significantly from REF both in magnitude and amplitude. This holds in particular for the Gotland Basin and the Bothnian Sea ( Figure 6). In the Bornholm Basin, the amplitudes of the maxOSF-variation are similar for all experiments because the sea level in Kattegat is the same. Hence, volume flows through the Danish straits affect the maxOSF in the Bornholm Basin but not in the Gotland Basin.
In the experiments GSLR and W1904, the maxOSF is systematically higher and lower,

Accepted Article
This article is protected by copyright. All rights reserved.
respectively. In all other experiments, the deviations to REF are relatively small and will be explained in detail in Section 3.4.
Also the time-integrated wavelet coherences between the multi-decadal variations in maxOSF and the parameters runoff, mean salinity, mean water temperature, and mean meridional wind speed at 10 m height in the different sub-basins for both models (not illustrated) have mainly low values (0.3-0.4 on average) in the short-period range and higher values (0.4-0.5 on average) on time scales greater than about 30 years.

Accepted Article
This article is protected by copyright. All rights reserved.

Centennial changes
Furthermore, the long-term mean maxOSF, average salinity, average water temperature as well as the average zonal and meridional wind at 10 m height for the individual subbasins are analyzed. The mean values with their standard deviations for all investigated model simulations from RCO and MOM for the different sub-basins are shown in Figure  8. Note, that the wind field for RCO is the same in all experiments within each sub-basin except for the sensitivity experiment W1904. When compared to REF the experiment W1904 reveals stronger meridional and weaker zonal wind in the northern Baltic Sea and slightly enhanced wind in the south. Changes in the forcing cause in some of the sensitivity experiments considerably different water temperatures and salinities. In T2008 the water temperature is noticeable higher than in REF, namely by about 2°C in the southern sub-basins and about 1°C in the north. In W1904 the water temperature is slightly increased compared to REF because the surface layer (with higher temperatures than at greater depths) is thicker in this experiment. Mean salinity in the experiments FRESH and W1904 is decreased compared to REF in most of the sub-basins, by about 1 and 1.5-2 g kg -1 , respectively. Further, in GSLR the salinity is about 1.5-2 g kg -1 higher than in REF, which is in agreement with Meier et al. (2018). For MOM the temperature is lower and salinity is higher than for RCO-REF in the southern and central Baltic Sea.

Accepted Article
This article is protected by copyright. All rights reserved. The different OSF magnitudes determined for the different RCO sensitivity experiments might be explained by changes in entrainment velocity under the certain hydrological and atmospheric forcing conditions of each experiment following Stigebrandt (1985). Therefore we calculated the long-term mean entrainment velocity for the whole Baltic Sea for each sensitivity experiment and compared it to that of the reference simulation ( Figure 9). The entrainment velocity for RCO-REF is predominantly positive with magnitudes of around 1*10 -5 m s -1 over wide areas of the Baltic proper and the Gulf of Finland (erosive deepening). The magnitudes are increased in the Landsort sub-basin (approximately 2*10 -5 m s -1 ) as well as in parts of the Bothnian Sea and Bothnian Bay (partly up to almost 5*10 -5 m s -1 ). Strongly varying entrainment velocity values in the Kattegat are due to the shallow bathymetry which makes the determination of the mixed layer depth and hence the entrainment velocity more variable than at deeper locations.

Accepted Article
This article is protected by copyright. All rights reserved.
Differences between the entrainment velocity determined from the sensitivity experiments and REF are smallest (predominantly around zero for most regions of the Baltic Sea) for T1904 and RUNOFF (which both have a similar maxOSF as REF), but also for the experiment FRESH. Although the increased runoff in the experiment FRESH does not lead to a noticeable entrainment velocity change when compared to REF it leads to the mentioned slight decrease in maxOSF because of a weaker perennial stratification (Meier, 2005). In the experiment T2008 the entrainment velocity is smaller than in REF almost throughout the whole Baltic Sea indicating that increased air temperature leads to weaker entrainment from below into the surface layer and consequently to a weaker maxOSF. An intensification and shallowing of the thermocline reduces the erosion of the halocline due to higher shielding (Hordoir et al., 2018;Gröger et al., 2019). For W1904 and GSLR the entrainment velocity in the Baltic proper is predominantly smaller and higher than in REF, respectively. Due to the long response time scale of the Baltic Sea of about 30 years (Meier and Kauker, 2003), changes in stratification do not impact multidecadal variations of the OSF but affect time scales longer than the response time scale.

Accepted Article
This article is protected by copyright. All rights reserved. Contrary to the claim by Hordoir et al. (2018), we found that maxOSF is not a good proxy for MBIs. The determined low correlation coefficients of 0.3 to 0.5 between salt transports associated with MBIs and the maxOSF indicate that less than 25 % of the variability of maxOSF is explained by MBIs. This is due to the fact that MBIs which are defined as periods with intensive saltwater inflow through the Danish straits do not necessarily propagate from the entrance region of the Baltic Sea to its central regions or, if they do they might have a travel time of several months. Moreover, MBIs are limited to the bottom layer while the OSF is a vertical averaging concept. Indeed, annual mean salt transports associated with MBIs and annual mean OSF maxima are uncorrelated. It is rather the wind parallel to the cross-section of the OSF which influences the maxOSF in the various sub-basins on decadal (Bothnian Bay) to multi-decadal time scales (Bothnian Sea and Gotland Basin) rather than river runoff or MBIs. Due to simple Ekman dynamics applied to a semi-enclosed basin in conjunction with mass conservation, wind opposite to the outflow direction of the sub-basin slows both mean out-and inflows down (Krauss and Brügge, 1991; see also Meier and Kauker, 2003).
Our results may indicate that a warmer future climate in the 21st century, for instance as recently projected using ensemble scenario simulations by Saraiva et al. (2019), may lead to a slight decrease of the maxOSF in most sub-basins (except in Bornholm Basin). Besides their projected increase in air and hence water temperature the projected changes in climate will be associated with increased river runoff (which would lead to decreased salinity), but particularly with a global mean sea level rise which in turn causes an increase in saltwater inflows as well as salinity and hence leads to stronger vertical stratification (Meier et al., 2017). The latter would hinder vertical oxygen fluxes (notwithstanding that stronger inflows transport more oxygen into the Baltic Sea) and thus may lead to an expansion of hypoxic areas. The sensitivity experiments in this study show that increased river runoff into the Baltic Sea slightly weakens the maxOSF as well. But a global sea level rise in turn would lead to an increase in the maxOSF which could counteract the effects by the local temperature and runoff increase. Hence, making a statement concerning the change of the OSF in the future climate is difficult as the magnitudes of the projected changes in the multiple drivers are rather uncertain. Uncertainties in, for instance, sea level projections are large (IPCC, 2013). Hordoir et al. (2018) investigated four scenario simulations driven by two different GCMs and by two greenhouse gas concentration scenarios, RCP 4.5 and 8.5. However, they neglected changes in river runoff and global sea level. In their two RCP 4.5 simulations, maxOSFs are projected to decline by 5 and 6.5 % at the end of the century which is somewhat larger than in our sensitivity experiment T2008. However, the warming in the scenario simulations by Hordoir et al. (2018) might also be larger than 2°C as in our experiment. Gröger et al. (2019) investigated an ensemble of five scenario simulations with a regional, fully coupled atmosphere-sea-ice-ocean model. They assumed a linear increase in runoff of 10% by the end of the twenty-first century. All ensemble members by Gröger et al. (2019) indicated a strengthening of the zonal, wind driven near surface overturning circulation in the southwestern Baltic Sea towards the end of the century whereas at the same time the overturning cell of the deep water is declining. In the Baltic proper, the strength of the meridional overturning shows no robust change, although three out of five ensemble members

Accepted Article
This article is protected by copyright. All rights reserved.
indicate at least a northward expansion of the main overturning cell. This result might be explained by differences in changes of the north-south wind component. In the Bothnian Sea, all ensemble members show a significant weakening of the meridional overturning (Gröger et al., 2019).
As our study shows that the wind may have a strong influence on the OSF further detailed analyses should be performed with focus on changes in future wind conditions such as the northward shift of the westerlies during summer reported by Gröger et al. (2019).
Our results support earlier findings by Meier (2005). He showed with the help of sensitivity experiments that changes in the deep water transport depend on the time scale of the applied forcing anomaly. Long-term changes of freshwater and saltwater inflows and of lowfrequency wind anomalies on time scales long compared to the turnover time of the freshwater content cause the estuary to drift into a new steady state with altered salinity. However, stability and deep water transport are almost constant because mixing does not significantly change (Meier, 2005;his Fig. 3). This conclusion is true as long as the depth of the halocline does not significantly change compared to the mixed layer depth, e.g. in FRESH and RUNOFF. In our experiments W1904 and GSLR the halocline depth substantially deepened and shoaled, respectively, causing a significant change in the erosion of the halocline by wind-induced and convective mixing.
In a simple two-layer model such changes might be explained by changes in the entrainment velocity across the bottom of the mixed layer depth following Stigebrandt (1985). In REF, the first term for the calculation of the entrainment velocity (2 m0 u* 3 /(g Δρ/ρ h), see Equation 3) is dominant everywhere, whereas the third term related to the haline buoyancy fluxes at the surface is negligible (not shown). In T2008 the thermocline is more intense compared to REF due to warming (increased Δρ during spring and summer, see Equation 1) causing a reduction of the entrainment velocity and a shielding of the halocline from the impact of surface turbulence (Figure 9). An additional effect is the increased buoyancy flux B at the sea surface during spring (not shown). However, the latter term is smaller compared to the first term even in the northern Baltic Sea indicating that in case of warming the resulting reduction of thermal convection is less important contrary to the claims by Hordoir et al. (2019). In summary, the entrainment velocity is smaller in T2008 compared to REF and, hence, the overturning circulation is smaller. To the contrary, in GSLR the entrainment velocity increases because of a shallower halocline depth. As a consequence mixed layer depth h is shallower during most of the year. The effect is partly counteracted by an increased density jump at the mixed layer base during spring and summer (not shown). In W1904 there are reversed results for the entrainment velocity (which is decreased) with shallower mixed layer depth during summer and autumn and decreased density jump at the mixed layer base during the entire year.

Conclusions
The conclusions of this study are summarized as follows:

Accepted Article
This article is protected by copyright. All rights reserved. 1) Differences between the two models RCO and MOM, e.g. in strength of the OSF in the various sub-basins are partly explained by different wind and runoff forcing. The models have a response time scale between maxOSF and wind as well as between maxOSF and runoff of about 30 years. Furthermore, the water temperature and maxOSF co-vary on intra-decadal and multi-decadal time scales for MOM and RCO, respectively.
2) The OSF is not a good proxy for MBIs because MBIs do not necessarily propagate from the entrance region to the Gotland Basin. On a multi-decadal time scale the maxOSF in the Gotland Basin depends mainly on the wind parallel to the cross-section and opposite to the upper branch of the OSF which hampers both out-and inflow under periods with anomalous strong wind. MBIs through the Danish straits only affect the maxOSF in the Bornholm Basin.
3) The response of the overturning circulation depends on the time scale T of the forcing: 3.1) For T >> 30 years (which is the approximate response time scale of the Baltic Sea) a global sea level rise by 1.5 m as simulated in the sensitivity experiment GSLR leads to a significantly stronger OSF in the whole Baltic Sea. With 20 % more river runoff and with an about 2 °C higher air temperature the overturning circulation is slightly weakened especially in the central and northern Baltic Sea. Lacking MBIs weakens the overturning circulation considerably. In all these cases, the vertical stratification in the Baltic Sea is weakened causing a reduced baroclinic pressure gradient.
3.2) For T ≤ 30 years the maxOSF is anti-correlated to wind especially in the southern Baltic Sea with partly phase shifts of 5 to 10 years between both parameters. A wind anomaly opposite to the outflow direction hampers the outflow and consequently the inflow of saline water into the sub-basin due to mass conservation. All other multi-decadal (or shorter) changes of runoff or MBIs are unimportant for the maxOSF because the time scale of the forcing anomalies is too short such that the overall stratification cannot adopt to the applied changes. Hence, the multi-decadal maxOSF variations with a period of about 30 years are caused by the wind parallel to the cross-section of the OSF and not by river runoff or MBIs.
4) The impact of projected warming and freshening on the maxOSF is relatively small.
Future global sea level rise may counteract the decline of the OSF due to warming. However, projections of future sea level are uncertain preventing conclusions on quantitative changes of the OSF in the future.