Warming of Baltic Sea water masses since 1850

The Baltic Sea is among the fastest-warming seas globally in recent decades affecting biogeochemical conditions such as euxinic areas but also pelagic and benthic marine ecosystems. It is therefore crucial to understand how this heat gain is distributed vertically. We used reconstructed atmospheric forcing fields for 1850–2008 to perform an ocean climate simulation that adequately captures climatogical temperature and salinity profiles. Then, a water mass classification distinguishes three water masses corresponding to the classical view, warm and fresh surface waters, cold and fresh intermediate waters, and cold and salty bottom waters, and two transition water masses. The temperature trends show a similar three layers pattern with fast warming at the surface (~ 0.06 K decade− 1) and bottom (> 0.04 K decade− 1) and slow in the intermediate layers (< 0.04 K decade− 1). The slow warming in the intermediate layer is explained by both weakly warmed water winter convection and the summer surface thermocline isolating the intermediate layers. The warming in the deep layers is related to warm surface inflow from the North Sea and Baltic proper in the southern and northern Baltic Sea respectively. Furthermore, sensitivity experiments show that the warming magnitude is controlled by rising air temperature while the vertical distribution of heat gain is related to surface wind conditions. Finally, the North Atlantic Oscillation and Atlantic Multidecadal Oscillation are well correlated with the temperature minimum and thus modulate the magnitude of warming in the intermediate layers on shorter time scales. This study provides a new picture of the Baltic Sea’s warming and suggests that this complexity is essential for understanding the influence of climate change on marine ecosystems.


Introduction
The Baltic Sea's temperature response to anthropogenic climate change (henceforth climate change) has significant implications for its overall functioning: hydrodynamics, biogeochemical cycles, and marine ecosystem. For instance, the vertical heat distribution influences the temperatures, the overturning circulation Placke et al. 2021), the stratification (Hordoir and Meier 2012), the nutrient supply to the surface , the surface dissolved oxygen and oxygen consumption at the bottom Krapf et al. 2022), the cyanobacteria blooms (Kahru et al. 2020;Wasmund 1997) or the body condition of cod (Receveur et al. 2022). In addition, the Baltic Sea has strong seasonal characteristics, such as the summer thermocline, spring blooms, winter and summer North Sea inflows, and winter ice cover (Feistel et al. 2008). With observed warming of 1.35 °C between 1982 and 2006 that exceeds the global average rate by a factor of seven (Belkin 2009), it is crucial to understand how this heat gain is distributed vertically and seasonally. Dutheil et al. (2022) have already analyzed seasonal variations in long-term trends in sea surface temperature (SST). This study goes beyond by analyzing the long-term trends in temperature as a function of depth in order to find the fingerprint of climate warming in different water masses of the Baltic Sea.
The Baltic Sea is a semi-enclosed marginal sea with a large latitudinal extent [54°N-66°N], resulting in a strong horizontal temperature gradient (~ 10 °C difference between the southern and northern Baltic Sea). In addition, large annual freshwater inflows from a large catchment area (4 times the size of the Baltic Bergström and Carlsson 1994) combined with a long turnover time of the total freshwater supply (35 years; Meier and Kauker 2003) are responsible for a strong horizontal salinity gradient (~ 10 g kg − 1 difference between the southern and northern Baltic Sea). The Baltic Sea comprises several basins with limited exchanges between them due to a complex bathymetry with several straits ( Fig. 1; Jakobsson et al. 2019). These different characteristics are the source of the significant horizontal differences in salinity and temperature between these basins. The vertical structure is qualitatively more homogeneous between the basins and can be decomposed into three layers: fresh and warm water at the surface, cold and fresh water in the intermediate layer, and cold and salty water at the bottom. This picture is actually more complex on a seasonal scale. The water column is well mixed up to the perennial halocline in winter. Then in spring, the surface thermocline develops, and the cold intermediate layer (CIL) appears. The CIL is characterized as a minimum temperature between the thermocline and the perennial halocline with a thickness estimated at 20-50 m (Liblik and Lips 2017). During spring, the temperature in the CIL is often lower than the temperature of maximum density (TMD). Thus, its depth can be explained by lateral advection of higher saline dense water (Chubarenko et al. 2017). The summer thermocline isolates the CIL and thus persists until fall/winter, when wind-induced mixing and winter convection erodes the surface thermocline and mixes the surface and subsurface waters, leading to the gradual disappearance of the CIL (Stepanova 2017).
The dynamics of the deep waters in the southern and central Baltic Sea are primarily driven by the inflows of salty and oxygenated water from the North Sea. The Major Baltic Inflows (MBIs) are generally barotropic in origin and occur in winter when weather conditions are favorable (e.g., Schimanke et al. 2014), but smaller inflows can also occur in summer driven by baroclinic pressure gradients (Mohrholz et al. 2006). The strong MBIs can propagate to the northern Gotland Basin and the Gulf of Finland, but their frequency is sporadic (one to a few per decade; e.g., Mohrholz 2018). Hence, a pronounced seasonal cycle of the deep water is not found except in the southwestern Baltic Sea (i.e., Arkona Basin and Bornholm Basin, see Fig. 1; Krapf et al. 2022) which is influenced by weaker and more frequent inflows (one to a few per year; Mohrholz et al. 2006). The northern Baltic Sea (i.e., Bothnian Sea and Bothnian Bay; Fig. 1) is separated from the central Baltic Sea by the Åland Sea, closed by a few sills that limit the water exchange (Jakobsson et al. 2019). Thus, the water exchange between the northern and southern Baltic Sea is considered an estuarine-like type, with surface water inflows from the northern Gotland Basin leading to a deep-water renewal in the Bothnian Sea (e.g., Marmefelt and Omstedt 1993). Giesse et al. (2020) showed in numerical experiments that winter convection has a small impact on deep water renewal in the northern Baltic Sea despite a lower vertical density gradient than in the south. These recent results thus support the classical view of deepwater ventilation by Baltic proper inflow (Marmefelt and Omstedt 1993).
Two of the key climatic oscillations for the Northern Europe climate are the North Atlantic Oscillation (NAO; Hurrell 1995) and the Atlantic Multidecadal Oscillation (AMO; Knight et al. 2005). These two oscillations modify the weather patterns over Northern Europe and thus alter precipitation, river runoff, sea level pressure, or wind fields (e.g., Cassou et al. 2004;Knight et al. 2006;Lehmann et al. 2011). These changes have a direct impact on the water characteristics of the Baltic Sea, with, for example, a positive correlation between NAO and SST or temperature in the CIL (e.g. Mohrholz et al. 2006), a negative correlation between NAO and ice extent (Börgel et al. 2020) and a negative correlation between AMO and salinity (Börgel et al. 2018).
These climatic oscillations modulate trends in the Baltic Sea's water characteristics on short to medium time scales. Thus, long time series are needed to identify the footprint of climate change on them. Numerous studies, using either in-situ or modeled data, have shown over different periods a significant increase in sea surface temperatures (Belkin 2009;Dutheil et al. 2022;Gustafsson et al. 2012;Kniebusch et al. 2019a;Lehmann et al. 2011;Siegel et al. 2006). This average Baltic Sea surface warming has been attributed to changes in air-sea heat fluxes in response to increasing air temperature (e.g., Kniebusch et al. 2019a), while the spatial and seasonal pattern of this surface warming can be explained by winter ice cover and surface winds ). In addition, Mohrholz et al. (2006) related the recent temperature increase in the Bornholm Basin to summer baroclinic inflows, and Kniebusch et al. (2019a) showed significant bottom warming of the entire Baltic Sea without exploring the cause. Most of these studies were focused on surface or bottom trends, and only a recent study (Liblik and Lips 2019) explored the vertical temperature trends from insitu data and showed significant and similar warming in surface and deep waters (0.4-0.6 °C decade − 1 ). Surprisingly no changes were detected in the CIL over 1982-2016, but these trends calculated over a 35 years period mixed the effects of climate change and the natural variability. Therefore, there is a need to explore the influence of climate change on Baltic Sea temperature more widely.
The present study aims to (1) characterize the vertical water masses in the Baltic Sea to (2) examine their seasonal variation in temperature and salinity. Then (3) we will use this framework to investigate the long-term trends in temperature, vertically and by water mass. Finally, we will explore (4) how the external (i.e., air temperature, surface winds, river runoff and sea level) and remoted (i.e., NAO and AMO) forcings modify these long-term trends. To achieve these objectives, we used a well-established general circulation model for the Baltic Sea and performed a reference historical climate simulation and sensitivity experiments for the 1850-2008 period associated with statistical methods of water mass classification. The paper is organized as follows: Sect. 2 describes the model configuration, the performed simulations, and the statistical analysis. Section 3 describes the climatological results, and Sect. 4 explores the long-term trends. Finally, Sect. 5 discusses the uncertainties associated with our results, compares them to the literature, and provides a summary with concluding remarks.

RCO model
In this study, the Rossby Centre regional Ocean climate model (RCO; Meier et al. 1999) coupled with a Hibler-type sea ice model (Mårtensson et al. 2012) has been used in the same configuration as in Kniebusch et al. (2019a) and Meier et al. (2019). The horizontal and vertical resolutions are of 3.7 km and 3 m, respectively. The subgrid-scale mixing in the ocean is parameterized using a k-ε turbulence closure scheme with flux boundary conditions (Meier 2001). A flux-corrected, monotonicity-preserving transport scheme is embedded without explicit horizontal diffusion (Meier 2007). The model domain comprises the Baltic Sea and Kattegat with lateral open boundaries in the northern Kattegat. In the case of inflow temperature, salinity, nutrients (phosphate, nitrate, and ammonium), and detritus, the values are nudged toward observed climatological profiles. In the case of outflow, a modified Orlanski radiation condition is used . Daily sea level variations in the Kattegat at the open boundary of the model domain were calculated from the meridional sea level pressure gradient across the North Sea using a statistical model . Fluxes between the atmosphere, ocean, and sea ice (heat, radiation, momentum, and matter) are parameterized using bulk formulae adapted to the Baltic Sea region (Meier 2002). Inputs to the bulk formulae are atmospheric state variables, including 2 m air temperature, 2 m specific humidity, 10 m wind, cloudiness and sea level pressure, as well as ocean variables such as SST, sea surface salinity, sea ice concentration, albedo, and water and sea ice velocities. The RCO model is a state-of-the-art ocean circulation model for the Baltic Sea and was recently applied in many climate studies, for instance, by Dutheil et al. (2022), Kniebusch et al. (2019a, b), Meier et al. (2018Meier et al. ( , 2021Meier et al. ( , 2022, and Saraiva et al. (2019a, b).

Reference simulation (REF)
The RCO historical simulation  has been forced at the surface by the HiResAFF data set developed by Schenk and Zorita (2012). This data set was built by applying the analog method to assign regionalized reanalysis data to the few available observational stations in the early periods. Thus, consistent multivariate forcing fields are obtained without artificial interpolation. River runoff and riverine nutrient loads were reconstructed following Gustafsson et al. (2012) and Meier et al. (2012). This data set was already evaluated in Gustafsson et al. (2012) and Meier et al. (2019) and used in Kniebusch et al. (2019a), Radtke et al. (2020), Placke et al. (2021) and Dutheil et al. (2022). This simulation will hereafter be referred to as REF.

Sensitivity simulations
In addition to REF, five sensitivity simulations were run to investigate the influence of river runoff, sea level, surface winds and long-term trends in air temperature on seawater temperature trends. These simulations are like REF but with the following modified forcing data: • In the TAIR simulation, the air temperature from 1904 is used and repeated for all 159 years of simulation. The year 1904 was chosen because it reflects a cold period. • In the WIND simulation, the surface winds from 1904 are used and repeated for all 159 years of simulation. • In the RUNOFF simulation, the river runoff is fixed to the monthly climatology computed over 1850-2008, while in the FRESH simulation, the river runoff is increased by 20% compared to 1850-2008. • In the SLR simulation, the sea level is 24 cm lower than in 1960.
All simulations performed are summarized in Table 1.

Trend calculation
First, the monthly average of sea temperature is computed from the 2-daily snapshot outputs of the model. Then at each grid point, the linear trend is calculated with the Theil-Sen estimator (Sen 1968;Theil 1950). The trend computed with this method is the median of the slopes determined by all pairs of sample points. The advantage of this expensive method is that it is much less sensitive to outliers. The trends are computed seasonally and annually. In the last case, the annual cycle is removed before computing the linear trend. The significance of trends is evaluated from a Mann-Kendall non-parametric test with a threshold of 95%.

Water mass classification
The water mass classification partitions the ocean into distinct pieces whose origin, evolution and fate can be measured, studied and modeled (Groeskamp et al. 2019). This approach was first developed to overcome the scarcity of direct measurements of the ocean for tracking ocean circulation (Broecker and Peng 1982). Subsequently, several theoretical developments were done to study the transformation mechanisms of water masses in order to understand their life cycle during transport (Groeskamp et al. 2019). Here, we use this approach to vertically discretize the water masses in each Baltic Sea basin and then track their monthly evolution, long-term temperature trend and understand the mechanisms responsible for their changes. This approach is particularly interesting in our case to follow the evolution of the CIL whose life cycle depends on interactions at the air-sea-ice interface and internal processes. A statistical classification method is used to identify the watermass areas in the Baltic Sea. To that end, monthly climatologies of temperature and salinity are first computed. Before classification, the temperature and salinity values are normalized and standardized by basin. This first step is justified by the strong horizontal variations of temperature and salinity in the Baltic Sea and the fact that we focus on vertical variations of water masses. Indeed without this first normalization step, the water masses would not only be separated along the vertical but also along the horizontal, which would give a similar number of water masses along the vertical but multiplied by the horizontal diversity of the basins, which would simply lead to more confusion in our results (Fig. S2). The boundaries of each basin are shown in The yellow dashed lines from A to G represent the North to South vertical sections across the basins of the Baltic Sea, and from H to I and from J to K the west to east vertical sections across the Bothnian Sea and Bothnian Bay, respectively. Four stations (BO3, US5b, SR5, BY15) are represented by red dots Fig. 1. Finally, all the longitude-latitude pairs of the monthly climatologies of temperature and salinity are classified using a hierarchical clustering algorithm (euclidean distance and ward aggregation, using the 'Cluster' package in the R programming environment; Kaufman and Rousseeuw 1990). This classification algorithm first considers all pairs of longitude-latitude as an individual cluster. For each class, it calculates the distance to the other classes and merges the two closest classes. This process is repeated until a single class is obtained, and thus a so-called dendrogram is constructed. This dendrogram is cut to get five classes (not shown). One of the advantages of hierarchical clustering is the inclusion of subdivisions in the upper divisions, thus increasing the number of classes considered increases the accuracy of the clustering without changing the results.

North Sea inflows
The North Sea inflows were detected in RCO following Meier and Kauker (2003) using a salinity threshold of 17 g kg − 1 . This study focuses on the influence of North Sea inflows on bottom temperature trends. Thus, based on the results of Mohrholz et al. (2006) we modified the inflow volume threshold to 10 km 3 in the Arkona Basin for at least six consecutive days to consider all North Sea inflows (winter major and minor barotropic and summer minor baroclinic). Then, we calculated the annual frequency, the annual averaged volume, the mean duration, and the mean temperature of these North Sea inflows.

Mixed-layer depth and thermocline strength
The mixed-layer depth (MLD) was calculated following de Boyer Montégut et al. (2004), using a threshold value of 0.2 °C for the difference between the near-surface water temperature at 10 m depth and the temperature at the MLD. The thermocline strength was calculated as a finite-difference of the temperature at each depth level (from 0 to 100 m), then the maximum value at each longitude-latitude was extracted.

Upwelling
The upwelling frequency has been calculated using the same method as proposed by Lehmann et al. (2012). This method is based on the temperature difference between the coastal SST and the surrounding water. Thus, to detect an upwelling event, we calculate, from the 2-daily snapshots, the temperature difference of each pixel with the zonal mean corresponding to this pixel. An upwelling is detected if this difference is lower than − 2 °C. Finally, a mask is applied for removing all the points located beyond 28 km from the coast. As this method is based on a difference with the zonal mean, it is limited for regions where the coastline is mainly oriented along an east/west axis, as in the Gulf of Finland. Nevertheless, this automatic method has been compared to visual analysis and has shown excellent skills (Lehmann et al. 2012).

Validation of REF
Over the recent period  the observed vertical profiles of annual temperature (black lines in Fig. 2) are characterized by high surface temperatures varying between 6 and 9 °C, depending on the location, and minimal temperatures (~ 3 °C) between 30 and 60 m that correspond to the CIL. Below this intermediate layer, a progressive increase in temperature with depth is observed. The warming below the CIL is less pronounced when moving northward, with finally no temperature increase below the minimum temperature depth at station BO3 in the Bothnian Bay. The shape of the vertical profiles and their variation with location are well represented by REF (blue lines in Fig. 2) with values within a one standard deviation envelope, but some biases are noticeable. The REF simulation underestimates the surface temperature from 0.5 to 1 °C in the Bothnian Bay and the Bothnian Sea and overestimates it by 0.5 °C at the BY15 station in the Gotland Basin. Furthermore, the REF simulation overestimates the temperature below the CIL at SR5 and BY15 but underestimates it at the BO3 and BY5 stations.
The vertical profiles of annual salinity show a different pattern between the northern (BO3 and SR5) and southern (BY15 and BY5) stations. At the northern stations, observations show a gradual increase in salinity with depth, while at the southern stations, there are two salinity jumps, the first near the surface (~ 30-40 m) and the second in the deeper layers (~ 80-90 m). These patterns are well represented in REF. However, the simulation underestimates salinity (~ 0.5-1 g kg − 1 ) in the northern Baltic Sea and below 50 m in the Bornholm Basin. Figure 3 extends this validation with temperature/salinity diagrams at these four stations. The observations (black dots) show a characteristic shape with a decrease of temperature until a minimum and associated with a weak increase of salinity. This minimum corresponds to the core of the CIL. Then, the temperature and salinity increase in the deeper layers. This shape is very well represented in REF (blue dots in Fig. 3), but the salinity is underestimated (~ 0.5-1 g kg − 1 ) in the northern Baltic Sea. To extend the validation period, Fig. 2 and Fig. 3 show also a comparison between observed and modeled values for 1900-1959 at the BY5, BY15 and SR5 stations. At these three stations the REF simulation shows similar skills to represent the observed annual temperature and salinity profiles and the T/S diagrams than for 1 3 . We note that about 40 years, 30 years and 20 years of data are available at SR5, BY15 and BY5 stations respectively for 1900-1959. The reader may also refer to Meier et al. (2019) for a validation of seasonal variations in salinity and temperature values. In addition, Placke et al. (2018) showed that the RCO model is within the state of the art of regional models over the Baltic Sea. These validations highlight the ability of the RCO model to represent the vertical structure in temperature and salinity of the Baltic Sea and thus indicate that it is appropriate for our study.

Climatology of vertical temperature
A cross-section from North to South of the Baltic Sea crossing most basins (yellow dashed lines in Fig. 1) is analyzed to characterize the vertical structure in temperature of the Baltic Sea. Figure 4 shows the annual and relative seasonal temperature along this cross-section. As described in Sect. 3.1, the vertical annual temperatures (Fig. 4a) are characterized by maximal temperatures in surface separated by a thermocline from minimal temperatures between 30 and 60 m, and a slight warming below this depth. A negative North-South temperature gradient is also visible. The relative seasonal temperature reveals that during the winter, the surface and sub-surface water are well mixed with low vertical gradients of temperature until 60 m (Fig. 4b). In this season, the positive temperature anomalies in deep layers of Bornholm Basin are probably related to the winter inflow from the North Sea. The CIL appears in spring with the core located around 40 m (Fig. 4c). At the same time, a surface thermocline starts to develop in the southern Baltic Sea. During the summer, the surface thermocline (~ 20-30 m) is developed, and the core of CIL deepens around 50 m (Fig. 4d). Finally in fall, the surface thermocline is progressively eroded by wind-induced mixing and the core of the CIL reaches 60 m before to completely disappear in winter (Fig. 4e).

Classification of water masses
A water mass classification based on the monthly climatology of temperature and salinity has been performed to better characterize the evolution of the vertical structure in temperature and salinity of the Baltic Sea and especially to follow the life cycle of the CIL. The classification algorithm has determined five classes (Fig. 5). Class C1 (in green) occurs from the surface to 60 m during winter and early spring and is characterized by fresh (from 2 to 8 g kg − 1 ) and cold water (0-5 °C; Fig. 6). Then during spring, C1 is progressively replaced in the surface layers by the class C4 but persists in the intermediate layers (30-60 m) in summer, almost disappears in fall except in the southern Baltic Sea before re-emerging in December (Fig. 5). C1  ods 1960-2008 and 1900-1959 respectively. The shaded envelopes around the vertical profiles represent one standard deviation calculated as the interannual variability over the period considered. The red dots in Fig. 1 show the locations of these stations presents exactly the water characteristics and dynamics to the CIL. Below the class C1 between 50 and 80 m, class C2 (in blue) is present all year (Fig. 5). C2 is characterized by cold (> 6 °C) and saline water (3-10 g kg − 1 ), explaining its deeper origin compared to C1 and is present all year (Figs. 5 and 6). C3 (in red) defines the summer surface water with high temperature (10-18 °C) and low salinity (> 7 g kg − 1 ). C4 (in yellow) is a transition water mass that occurs primarily in spring and fall but persists in summer below C3 near the thermocline (Fig. 5). C4 is characterized by intermediate temperature (5-13 °C) and low salinity (> 7 g kg − 1 ). Finally, class C5 (in purple) is in the deepest layers of the Baltic Sea between the permanent halocline and the bottom and is present all year (Fig. 5). This class is thus characterized by cold (> 8 °C) and saline water (3-14 g kg − 1 ). Figure 7 displays the monthly variations of temperature, salinity and volume by basin in order to provide a deeper understanding of the dynamics of these water masses. First, C1 shows strong seasonal variations in the Bothnian Bay and the Bothnian Sea (Fig. 7). From November to March, the The vertical and horizontal lines represent one standard deviation calculated as the interannual variability over the period considered for temperature and salinity respectively. Red dots in Fig. 1 show the locations of these stations temperature decrease is accompanied by a salinity increase (~ 0.2 g kg − 1 ), probably due to the brine rejection during the sea ice formation and vertical mixing between the surface and sub-surface (Fig. 7a). The volume of C1 increases from November to December then it is maximal and stable until March (Fig. 7b). In contrast, the temperature increase from April to October is associated with a weaker change in salinity (0.1 g kg − 1 ) and a strong volume decrease of C1, suggesting a heat gain mainly by diffusivity and radiation. In the southern Baltic Sea, the absence of sea ice formation induced no salinity change in winter. The salinity increase from May to October is probably due to the deepening of C1 associated to the decrease of its volume. In fall, a strong decrease of salinity is modeled for C1 in all basins, probably related to high precipitation and runoff rates in this season (e.g., Meier et al. 2006;Saraiva et al. 2019) and to the windinduced mixing with fresher surface waters as indicated by the strong volume increase at this period.
For C2, in the Bothnian Bay and the Bothnian Sea, the decrease in temperature during winter is much lower than C1 and is associated with a weak change in salinity. This class is between 60 and 80 m below the well-mixed, cold winter convection zone (i.e., C1), suggesting that this weak winter temperature decrease in C2 is related to diffusive heat loss as supported by the stable volume of C2 at this season. From April to October, temperature changes in C2 are associated with salinity changes and volume increase, suggesting either inter-layers mixing or convection along the slopes. Longitudinal transect reveals that C2 is connected in May and June to fresher coastal surface water (Fig. S3a-f), explaining partly the salinity decrease during these months. The analysis of the frequency of upwelling reveals that this connection is due to more frequent upwelling events during May and June than in April (Fig. S3g-i). In the Gotland Basin, monthly variations in temperature and salinity are minimal for C2, suggesting that this class is isolated from others with little exchange of matter and heat as supported by the very stable volume (Fig. 7b).
For C3 and C4, a strong increase in salinity occurs in spring. These two classes are in the coastal regions in spring and expand to the centers of the basins in summer ( Fig.  S3a-f). The salinity of the coastal regions is very low in spring due to increased runoff from precipitation and snowmelt. Thus this expansion from the coast to the center of the   Fig. 1 indicate the transect sections basins is accompanied by an increase in volume and explains the strong salinity increase in spring. Finally, C5 is completely isolated by the permanent halocline and thus presents no seasonal variation except in Bornholm Basin. In this basin, from December to May the temperature decreases and is associated with a salinity increase and the opposite the rest of the year. Similar behavior is modeled for C2 in the Bornholm basin. These salinity and temperature variations for C2 and C5 in the Bornholm Basin are likely explained by a strong influence of weak and frequent inflows from the North Sea.

Results: long-term temperature trends
In this section, first, the long-term trend (over 1850-2008) in temperature as a function of depth and water masses from the REF simulation will be analyzed. Then, the sensitivity simulations to external forcing (air temperature, surface winds, and runoff) and the effect of the natural variability are also analyzed to further improve our understanding of the mechanisms responsible for these trends. Figure 8 displays the annual and relative seasonal temperature trends along the North-South transect (yellow dashed lines in Fig. 1). The annual temperature trends exhibit a vertical three-layers pattern similar to climatological temperatures (Fig. 4a vs. Fig. 8a) with high temperature trends in surface layers (~ 0.06 K decade − 1 ), minimal temperature trends in intermediate layers (< 0.04 K decade − 1 ), and higher temperature trends in deep layers (> 0.04 K decade − 1 ). The trends and absolute values of relative temperature in winter and spring also show a similar spatial pattern (spatial pattern correlation > 0.6; Fig. 4b-c vs. Fig. 8b-c). In winter, negative relative values are found everywhere except in the deep layers of the southern Baltic Sea. In spring, the same spatial pattern is observed as in winter. The only notable difference is the positive relative surface trends in the southern Baltic Sea. In summer, strong positive values of relative warming (up to 0.06 K decade − 1 ) are observed at the surface in the northern Baltic Sea, associated with strong negative values in the sub-surface (up to -0.06 K decade − 1 ). In autumn, the strong positive surface values disappear, and the negative subsurface values are less pronounced (up to -0.04 K decade − 1 ) than in summer but extend to the southern Baltic Sea.

Mechanisms of formation
This section focuses mainly on the intermediate and deep layers because these areas have been little studied in terms of temperature trends, unlike the surface. Seasonal relative temperatures and their trends exhibit strong similarities, especially in the intermediate layers ( Sect. 4.1). Thus, we test whether the slow warming of the intermediate waters is due to the translation of cold and weakly warmed surface waters into the intermediate layers with winter convection, then isolated from the surface by the summer thermocline and thus conserved until autumn. To evaluate this potential mechanism, we first show in Fig. 9 the monthly evolution of vertical temperature trends averaged by basin. For intermediate layers (40-60 m), the temperature trends follow this seasonal cycle: they increase until spring, then decrease to reach their minimum value in August-September, and finally increase in the second half of the fall until the following spring. The basins differ in that the maximum warming is reached later when moving northwards, from March for Bornholm Basin to June for Bothnian Bay. This result suggests that winter convection fuels the CIL in the southern Baltic Sea while in the northern it might be formed by thermal convection during warming in spring as the temperature in the CIL is lower than the TMD in winter. This explanation is supported by the timing when the sea temperature matches the TMD (Fig. S4) and the strong vertical gradient of the winter temperature trends in the northern Baltic Sea. These seasonal variations of temperature trends in the intermediate layer with maximum values in winter/spring and minimum values in late summer show that the slow warming modeled in the intermediate layer do not originate solely from winter convection and that other processes are also at work. Fig. 7 Monthly variations of (a) Temperature/Salinity (T/S) diagram and (b) volume of water masses (in 10 3 km 2 ) averaged by basin (rows) crossed by the North-South transect (yellow dashed lines in Fig. 1). Each cluster of water mass is represented by a marker: C1 is a square, C2 is a solid circle, C3 is a triangle, C4 is a diamond and C5 is a ring. The color of the markers represents in (a) the month of the year, from dark blue (January) to dark red (December) and in (b) distinguish the different clusters (C1 in green, C2 in blue, C3 in red, C4 in yellow and C5 in purple)  This analysis horizontally averages coastal and offshore waters, thus mixing different water masses and their properties (e.g., long-term trends) as different water masses coexist at the same depth (Fig. S3). Therefore, in Fig. 10 we analyzed the monthly evolution of temperature trends by water masses to disentangle this effect. In the Bothnian Bay and the Bothnian Sea, the monthly variations in temperature trends for C1 (green) and C2 (blue) are slightly more variable than described previously, with minimum values in winter and summer and maximum values in fall and spring. The variations in spring and summer temperature trends for C1 and C2 are in phase with C3 and C4; only the amplitude is damped. This result suggests a significant influence of Theses monthly variations in northern Baltic Sea can be explained by various mechanisms. In winter, two processes can explain the weak sub-surface warming (C1 and C2).
First, wind-induced mixing of surface and subsurface waters prior to sea ice formation. As freezing is delayed by climate change this process can take place over a longer time period. The second process is winter convection (due to brine rejection) of surface waters that are weakly warmed due to the isolation of the ice cover . Then in May and June, the surface warming is maximal due to the positive albedo feedback in response to the melting of the sea ice. As the surface temperature (C4) reaches the TMD at this time (Fig. S4), thermal convection translates the surface heat gain in the water column up to 60 m ( Fig. 9a-b). In summer, the surface warming gradually decreases with the reduction of shortwave radiation. The formation of a surface thermocline weaker than in the southern Baltic Sea only partially isolates the subsurface (C1 and C2) but enough to dampen the decrease. Finally in fall, the windy condition mixes the warm surface water (C3 and C4) with the subsurface, increasing temperature trends in C1 and C2.
In contrast, in the Gotland Basin, the decrease in temperature trends modeled between spring and summer in C3 (-0.05 K decade − 1 ) and C4 (-0.09 K decade − 1 ) is associated with no change in C2 and a weak decrease in C1 (-0.01 K decade − 1 ), highlighting a lower influence of surface warming on the sub-surface than in the northern basin. The absence of thermal spring convection as the water temperature is already larger than TMD at this season (Fig. S4) and a stronger summer thermocline can explain this decoupling between the surface and the sub-surface.
Finally, the seasonal variations in temperature trends in the Bornholm basin are significantly different from those in the other basins. First, C1 and C2 are phase-shifted with a decrease in temperature trends in C1 from April to October and an increase after and the opposite in C2. In contrast, C1 is partially (from May to December) in phase with C4 and C2 with C5. This shows that the warming in C1 is rather surface controlled, while the warming in C2 is directly related to the bottom conditions. In addition, the Bornholm Basin is the only basin to show seasonal variations in temperature trends in C5. These variations agree with seasonal variations in surface warming in the Danish straits and Kattegat (see Dutheil et al. 2022), suggesting that the bottom temperature trends in Bornholm Basin are controlled by North Sea inflows. Figure 11a shows the annual time series of North Sea inflows temperature (a 10 years running average is applied) and its long-term trend. C5 and North Sea inflows have similar temperature trends (> 0.1 K decade − 1 ) supporting the idea that North Sea inflows drive temperature trends in C5.
To summarize, the warming of the intermediate layer in the northern Baltic Sea is controlled by various seasonal processes such as wind-induced mixing, haline and thermal convection, and shortwave radiation. In the Gotland Basin, surface processes have a weaker influence on the temperature trends in the intermediate layer than in the northern Baltic Sea, probably due to a stronger stratification and the absence of sea ice. Finally, the bottom temperature trends in the Bornholm Basin can be related to the North Sea inflows.

External forcings
This section aims to understand which external drivers influence long-term temperature trends. Dutheil et al. (2022) highlighted the impact of surface winds and air temperature trends on surface warming. Hence, we performed sensitivity simulations for these two drivers. Figure 12 compares the annual and summer temperature trends along the northsouth transect in the REF, TAIR, and WIND simulations (summer is the season where the differences between the sensitivity experiments are most notable; see Fig. S5). No significant annual warming is modeled in TAIR except in the deep layers of the southern Baltic Sea, emphasizing the influence of the North Sea inflows on temperature trends in this area. As the open boundaries conditions at Kattegat are set to the monthly climatological conditions, we analyzed three indices (mean volume, duration, and annual number of North Sea inflows), characterizing potential changes in the intensity and frequency of North Sea inflows (Fig. 11b-d) and which could explain the bottom warming modeled in the southern Baltic Sea. This analysis reveals no change in duration, a slight decrease in mean volume (-0.1 km 3 decade − 1 ; Fig. 11c) and an increase in annual number of North Sea inflows (+ 0.6 inflows decade − 1 ; Fig. 11b). The latter can lead to these bottom temperature trends in the southern year running average is applied for all time series. The temporal correlation between the AMO or NAO and the temperature minimum for each season is indicated in Table 2 Baltic Sea observed in TAIR. An additional mechanism which is not evaluate here would be a shift in seasonality, with more warm inflows during summer and fall instead of cold inflows during winter. Despite no significant annual warming in TAIR (Fig. 12c), the spatial pattern of relative summer temperature trends are very similar to REF (spatial pattern correlation of 0.89) with a reduced magnitude (Fig. 12b vs. d). As only the air temperature trends are suppressed in TAIR, this summer warming pattern is driven by the other surface forcings, which remain identical to REF.
Thus, it is probably explained by the thermal isolation of the summer thermocline, which remains similar to REF since the wind-induced mixing is not modified. However, the absence of air temperature warming in TAIR reduces the sensible heat flux directed towards the sea during this season and may explain the smaller magnitude of this summer warming in TAIR than in REF.
Conversely, setting the surface wind conditions to 1904 induces a different spatial pattern of summer warming compared to REF (pattern correlation of 0.69), with weaker warming at the surface and much more vertically extensive negative relative trends than in REF. The summer of 1904 is windier than the REF summer climatological conditions (Fig. S1), resulting in a weakening of the thermocline and a deepening of the mixed layer depth in response to the increased wind-induced mixing (Fig. 13). Thus, the heat gain is more vertically extensive and reduces the warming contrast between the surface and the sub-surface. In addition, the windier conditions also enhance the latent heat flux, which further reduces the warming. In summary, changes in latent heat fluxes and wind-induced mixing in response to windier conditions alter the vertical distribution of heat gain and hence surface and subsurface warming. Finally, we also explored the influence of changes in runoff with the FRESH and RUNOFF simulations. The temperature trends in FRESH and RUNOFF show negligible differences relative to REF (Fig. S5), probably due to the long period considered in these simulations (159 years), as the stratification does not change after the spin up period.
In conclusion, the atmospheric conditions are responsible for both the magnitude and spatial pattern of vertical temperature trends on the surface and sub-surface. The air temperature trends dominate the magnitude of warming, while the climatological conditions of surface winds control the vertical distribution of heat and thus explain the vertical structure of surface and sub-surface warming.

Remote drivers
To further investigate the atmospheric control on temperature trends, we analyze the relationship between natural variability (i.e., AMO and NAO) and the temperature in CIL for the Bothnian Sea and Gotland Basin. Our classification of water masses is based on temperature and salinity monthly climatologies, thus it is not possible to follow the CIL over the period 1850-2008 with this method. Therefore, we consider the minimum temperature in the water column (Tmin) as a reference to track the temperature evolution in the CIL. Figure 14 compares the temporal evolution of normalized and standardized AMO and NAO indices with Tmin for each season. Overall, Tmin correlates well with the AMO and the NAO in winter, spring, and sometimes even summer (Table 2). This result is consistent with the high auto-correlation (from 0.6 to 0.8) between Tmin in winter and the other seasons. Thus, the warming in the intermediate layers is enhanced during positive phase of NAO and AMO, and damp during negative phase. Tmin is better correlated with the AMO in the Bothnian Sea (except in spring) and the NAO in the Gotland Basin (except in winter). The correlations decrease over the year, except for the NAO in spring, with a higher correlation than in winter. This result can be explained by the modification of the air-sea heat fluxes by the ice cover. Indeed, the persistence or melting of sea ice in spring is strongly related to the NAO (correlation of -0.65 between NAO and maximum sea ice cover; Fig.  S6). Furthermore, the monthly correlation between Tmin and NAO is maximum in March (0.67; not shown), which favors spring in this analysis and explains why this result also occurs in the Gotland basins where there is no sea ice.

Discussion
Overall, the REF simulation adequately represents the average salinity and temperature conditions, highlighting that our model configuration is suitable for analyzing temperature trends in the water column. Nevertheless, REF systematically underestimates salinity in the northern Baltic Sea (Figs. 2 and 3). Sensitivity simulations to runoff were performed and showed that increasing or decreasing the runoff does not change the temperature trends. Thus, the mean salinity has a negligible influence, probably because stratification remains similar without changing wind forcing, as shown by Meier (2005). This result implies that the salinity biases simulated in REF in the northern Baltic Sea probably have limited influence on our results.
In addition, the REF simulation tends to underestimate the sea ice cover (e.g., Kniebusch et al. 2019a). Thus air-sea heat fluxes and wind-induced mixing in winter and spring are probably overestimated in REF relative to observations. This bias may alter the temperature trends found on surface and sub-surface during these seasons in the northern Baltic Sea. However, Giesse et al. (2020) have shown that the haline convection due to brine rejection during the sea ice formation has non-significant impact on the deep-water renewal. Therefore, the bottom warming in the northern Baltic Sea might not be influenced by this bias and thus are likely related to surface water inflows from the Baltic proper. Likewise, bottom warming in the southern Baltic Sea are probably associated with North Sea inflows. The statistics of MBIs (event numbers, duration, volume) modeled are roughly well represented in the REF simulation compared to the observations (Mohrholz 2018). However, the monthly climatological conditions applied to the open Kattegat boundaries and land uplift are kept constant. Thus, potential changes (e.g., temperature or salinity) at this open boundary are not considered in this study. During the simulation period sea level has risen by about 25 cm relative to the sill bottoms causing an increase in salinity of about 0.35 g kg − 1 according to model simulations . This small salinity changes might have some impacts on stratification in the Baltic Sea due to increased saltwater inflows in response to the sea level rise. Thus, to assess the influence of sea level on Baltic Sea's warming we analyzed the SLR simulation. This simulation shows very similar temperature trends compared to REF, only the bottom warming in the Bornholm Basin during fall is slightly weaker than in REF (-0.02 K decade − 1 ).
A limitation of this study is the lack of air-sea coupling in our simulations. Indeed, Gröger et al. (2015Gröger et al. ( , 2021 showed a thermal response to air-sea coupling that influences the summer and fall air-sea heat fluxes. During these seasons, wind events draw colder waters from depth to the surface. Upwelling cools the atmospheric boundary layer and increases the atmosphere's static stability with a negative effect on wind speed (summer short circuit). However, the enhanced warming in the northern Baltic Sea due to the ice-albedo feedback causes higher atmospheric instability and larger wind speed ( ~ + 1 m.s − 1 ), as shown by Meier et al. (2011). These air-sea mechanisms could influence our results by changing the wind-induced mixing. Furthermore, Ho-Hagemann et al. (2017) showed with the help of the CCLM model that the air-sea coupling enhances the moisture convergence over the Baltic Sea, which likely modifies the shortwave radiation changes. Although Dutheil et al. (2022) found that the cloud cover trend has insignificant effects on the long-term SST trends, the reliability of the modeled cloud cover is questionable. Thus, this mechanism could alter our results in coupled air-sea simulations. Another limitation of our framework is the water masses classification. This classification is based on temperature and salinity monthly climatologies calculated over the 1850-2008 period. Thus, changes in the vertical or seasonal distribution of water masses during this period are not considered.
Overall, our results are qualitatively in agreement with the results of Liblik and Lips (2019), showing significant positive temperature trends in surface and bottom waters (0.4-0.6 °C per decade) and weaker changes in the CIL. The differences in the warming magnitude can be likely attributed to the different time periods considered in both studies. The main difference between these two studies concerns the bottom temperature trends in the northern Baltic Sea. Liblik and Lips (2019) showed no bottom temperature trend in the northern Baltic Sea, whereas our analysis revealed a significant temperature trend of 0.02-0.03 K decade − 1 . This difference can probably be attributed to an observation period of 35 years in Liblik and Lips (2019), analyzing thus the multi-decadal variability and not the long-term trend and to a lack of data in these basins coupled with sampling biased toward summer (their Fig. 2). To assess the influence of the period considered, we analyzed the temperature trends in the REF simulation during 1982-2008 (Fig. S7), the most like Liblik and Lips (2019). In agreement with Liblik and Lips (2019), Figure S7 shows non-significant warming of ~0.2°C. decade-1 in the deep layers of the northern Baltic Sea during 1982-2008. This highlights that more than 35 years of data collection are needed to extract significant trends in the bottom warming in the northern Baltic Sea.
Finally, our results have large implications because many Baltic Sea components are dependent of the temperature. For instance, Krapf et al. (2022) showed that the yearly variations of hypoxic and anoxic areas are strongly related to the bottom temperature. Likewise, Receveur et al. (2022) have also highlighted the role of bottom temperature warming in the recent decline of western Baltic Sea's cod stock. It is therefore clear that the assessment and understanding of the Baltic Sea's warming over the whole water column deserves more attention in order to accurately assess the impacts of climate change. An unresolved question here is how climate change will influence vertical temperature trends in the future considering the changes projected by climate models on the climatological conditions such as ice cover, stratification, and precipitation. As our runoff sensitivity simulations show, changes in salinity, such as those induced by increased precipitation over northern Europe in response to climate change (e.g. Christensen et al. 2022), will likely have a small effect on the vertical distribution of heat gain and temperature trends. However, changes in cloud cover, the decline in ice cover or an enhancement of stratification could alter long-term temperature trends.

Conclusion
This study investigates the influence of climate change on the Baltic Sea's warming under past conditions. A historical simulation (1850-2008) was performed with a regional oceanic model using HiResAFF atmospheric forcings that capture the main observed temperature and vertical salinity characteristics. This includes the annual mean vertical profiles of temperature and salinity but also their co-variations in temperature/salinity diagrams for the main basins. This accurate historical simulation of the temperature and salinity conditions justifies using our regional model configuration to assess the vertical temperature trends under past conditions. Then, a water mass classification based on a monthly climatology of temperature and salinity was performed to objectively separate the different water masses and characterize their vertical and seasonal variations and their long-term temperature trends. This analysis allowed us to distinguish five water masses, three corresponding to the classical three-layer view, surface, intermediate and deep, each characterized by different salinity and temperature values explaining their depth, and two being transition water masses. The calculation of the temperature trends has also revealed a picture in three layers with high temperature trends in the surface (0.06 K decade − 1 ) and at the bottom (> 0.04 K decade − 1 ) and low temperature trends in the intermediate layers (< 0.04 K decade − 1 ). Actually, the monthly variations of the warming reveal a more complex picture with several processes driving the vertical heat gain distribution. The slow warming of the intermediate layer has been related to: (1) the weakly warmed surface waters translating to the intermediate waters due to winter convection and wind-induced mixing, (2) the summer surface thermocline thermally isolating the intermediate layers. The faster warming of bottom layers can be related to the inflow of warm surface waters from the North Sea to the bottom of the Baltic Sea. In addition, sensitivity simulations to surface wind and trends in air temperature showed that the magnitude of surface and subsurface warming is primarily controlled by rising air temperature, while surface wind conditions drive the vertical distribution of this heat gain. However, variations in runoff do not impact long-term temperature trends. Finally, we analyzed how natural variability (i.e., NAO and AMO) influences the temperature in the intermediate layers.
Our analysis showed that Tmin correlates with AMO in the northern Baltic Sea and with NAO in the southern Baltic Sea. The correlation is largest in winter for AMO and in spring for NAO; this detail can probably be explained by the influence of NAO on the persistence of sea ice in spring which modulates the air-sea heat fluxes. Our results reveal for the first time the complex warming of the Baltic Sea and may have important implications for a better understanding of the climate change impacts on marine ecosystems.