1 Introduction

The Baltic Sea exhibits outstanding SST changes over the last decades with for instance an increase of 1.35 °C in 1982–2006, corresponding to seven times the global rate (Belkin et al. 2009). This surface warming has multiple consequences as for instance a thermal stratification enhancement also reduces vertical mixing, or an increased risk of climate extremes such as marine heat waves. These physical changes alter also the biogeochemical conditions by, for instance, limiting the supply of nutrient to the euphotic zone. Ultimately, these changes can have an impact to various economic sectors such as fishing or tourism with the summer cyanobacteria blooms (e.g. Neumann et al. 2012; Meier et al. 2019a, b). Therefore it is critical to characterize and understand the drivers of these historical SST changes, disentangle the impact of climate change from the others changes (e.g. changes in aerosols concentrations), and evaluate accurately the future changes. Many studies have analyzed these changes using satellite and in situ SST datasets or with climate models for different periods (e.g. Belkin 2009; Gustafsson et al. 2012; Kniebusch et al. 2019; Lehmann et al. 2011; Siegel et al. 2006). All studies showed pronounced changes, but to our knowledge only Kniebusch et al (2019) analyzed detailed spatial and seasonal heterogeneities. Their results demonstrated that seasonal SST changes are very heterogeneous among the sub-basins of the Baltic Sea with for instance twice as much warming in the Baltic proper than in the Bothnian Sea in winter, and the opposite in summer (Kniebusch et al. 2019).

These differences are partly related to the heterogenous climatic and geomorphological characteristics of the Baltic Sea. The exchange of water masses with the open ocean is constrained by the narrow and shallow Danish straits. In addition, the Baltic Sea has a low average depth (~ 54 m, Fig. 1) and a strongly variable bathymetry that further limits the exchange between the Baltic Sea sub-basins (Fig. 1). This is reflected in long surface water renewal rates of around 30 years in the southern part and 40 years in the northern part (Meier 2007; Snoeijs-Leijonmalm and Andrén, 2017). In addition, the Baltic Sea is the second largest brakisch sea after the Black Sea with an average salinity of ~ 7.4 g/kg (Meier and Kauker 2003). Its water masses can be understood as a mixture of saline water from the North Sea and freshwater due to the presence of a very large number of rivers on its shores (Mohrholz et al. 2015). Thus the Baltic Sea salinity exhibits a strong north/south gradient, which is an additional factor favoring the formation of sea ice in the northern part. Finally, the vertical structure is characterized by a strong seasonal thermohaline stratification during summer and a perennial halocline in the Baltic Proper and in the Gulf of Finland which greatly reduces vertical mixing between the surface and the deeper layers.

Fig. 1
figure 1

Bathymetry (in m) of the Baltic Sea and the locations of its sub-basins: the Bothnian Bay and Bothnian Sea, the Gulf of Finland and Gulf of Riga, the Gotland, Bornholm, and Arkona basins, and the Danish straits. The Baltic Proper gathers the Gotland, Bornholm and Arkona basins

Full size image

These characteristics of the Baltic Sea (i.e. quasi-enclosed, shallow, strongly stratified) explain in major part the strong control of the atmosphere on its physical characteristics (Stigebrandt and Gustafsson 2003). Some studies (e.g. MacKenzie and Schiedek 2007; Meier et al. 2019a, b; Kniebusch et al. 2019) have shown a strong correlation between air temperature and SST. For instance, using sensitivity experiments, Meier et al. (2019a, b) showed that interannual variability in air temperature is crucial in representing the long-term trend in SST. Kniebusch et al. (2019) also showed in climate simulations that daily air temperature variability explains the majority of the daily variance in SST. Nevertheless, the latter study also revealed that the explained variability differs by 40–95% depending on location. The minimum values were found along the Swedish coast and in the northern part (Bothnian Bay and Bothnian Sea and Gulf of Finland) due to upwelling and the presence of sea ice in winter and spring isolating the surface waters from the atmosphere, respectively. Their results thus showed that parameters other than air temperature explain the daily variability of SST in the Baltic Sea, and suggested but did not clearly prove the same processes impacting long-term trends. In line with these results Barkhordarian et al. (2016) highlighted the effect of aerosols on recent SST trends. Based on observations and simulations, they showed that the recent decrease in regional aerosol concentration (which alter cloud cover and cloud albedo and thus insolation) due to the reduced industrial emissions from ~ 1980, is responsible for the inability of climate models to simulate the correct amplitude of observed warming over the period 1983–2014. To summarize, there is evidence that the SST trends are first-order driven by the air temperature but with a strong imprint of local characteristics such as changes in ice cover, cloud cover or surface winds.

All of these climate drivers are expected to continue in the future in response to anthropogenic forcings (BACC II Author Team 2015). However, their changes are not a linear function of greenhouse gas concentration, as there are many feedback loops that may alter past changes in the future. Therefore it is impossible to simply extrapolate past changes to predict future changes, and it is necessary to assess how and why past trends will change. Furthermore, the climate projections are affected by various sources of uncertainty (e.g. internal variability, atmospherics forcings) that need to be assessed by huge ensemble simulations.

In addition, to efficiently support marine spatial management it is necessary to objectively characterize spatial and seasonal differences in SST patterns. Clustering methods can be used to define well-separated homogeneous climate zones in, for example, minimizing intra-class variance and maximizing inter-class variance (Bador et al. 2015; Bernard et al. 2013). These methods have been widely used to define recurrent spatial patterns over time (e.g. Cassou et al. 2004; Dutheil et al. 2020; Grams et al. 2017), but are not yet widely used to spatially separate temporal trends.

Thus, the aim of the present study is to (1) build homogeneous areas of seasonal SST trends in the Baltic Sea and (2) understand the processes explaining their characteristics. These two objectives will be achieved by using a reference historical climate simulation and sensitivity experiments for the 1850–2008 period associated with statistical methods of classification. In addition, we performed a large set of climate simulations under future conditions to (3) evaluate how global warming modifies historical changes and (4) assess the uncertainties related to the atmospheric forcings and the internal variability. These future simulations are identical to Meier et al. (2021) and were performed for two RCP scenarios, three sea-level rise scenarios, two nutrient load scenarios and four CMIP5 (Coupled Model Intercomparison Project) models, i.e. 48 climate change simulations. The paper is organized as follows: Sect. 2 describes the model configuration, the performed simulations, and the statistical analysis. Section 3 describes the results of the historical period and Sect. 4 of the future period. Finally Sect. 5 discusses the uncertainties associated with our results, compares them to the literature and provides a summary with concluding remarks.

2 Methods

2.1 RCO-SCOBI model

In this study the Rossby Centre regional Ocean climate model (RCO; Meier et al. 1999, 2003) coupled to a Hibler-type sea ice model has been used in the same configuration that in Meier et al. (2019a, b) and Kniebusch 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, and in the case of outflow, a modified Orlanski radiation condition is used (Meier et al. 2003). 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 (Gustafsson et al. 2012). The SCOBI model comprises the dynamics of nitrate, ammonium, phosphate, oxygen, hydrogen sulfide, three phytoplankton species (including nitrogen-fixing cyanobacteria), zooplankton, and detritus (Eilola et al. 2009). 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 2001). Inputs to the bulk formulae are state variables of the atmospheric planetary boundary layer, including 2-m air temperature, 2-m specific humidity, 10-m wind, cloudiness, and mean sea level pressure, and ocean variables such as SST, sea surface salinity, sea ice concentration, albedo, and water and sea ice velocities.

2.2 Historical simulation

The RCO historical simulation (1850–2008) has been forced at the surface by the High RESolution Atmospheric Forcing Fields (HiResAFF) data set developed by Schenk and Zorita (2012). This data set was built by applying the analogue 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 Meier et al. (2012) and Gustafsson et al. (2012), respectively. This data set was already evaluated in Gustafsson et al. (2012) and Meier et al. (2019a, b) and used in Kniebusch et al. (2019), Radtke et al. (2020) and Placke et al. (2021). This simulation will hereafter be referred to as REF.

2.3 Sensitivity simulations

In addition to REF, three sensitivity simulations were performed in order to investigate the influence of long-term trends in air temperature, surface winds and cloud cover on SST trends. These simulations are similar to REF, but with the following modified forcing data:

In the TAIR and WIND simulations, the air temperature and surface winds from 1904 are used and repeated for all 159-years of simulation respectively. The year 1904 was chosen because it corresponds to a cold period, i.e. without global warming of air temperature.

In the CLOUD simulation, the cloud cover from 1940 is used and repeated for all 159-years of simulation. The year 1940 was chosen because it corresponds to the year closest to the multiple year (1850–2008) climatology. All simulations performed over the 1850–2008 period are summarized in Table 1.

Table 1 List of simulations over 1850–2008
Full size table

2.4 Climate change simulations

To assess how the historical SST trends will be modified in the future two sets of climate change simulations were carried out.

First, Dieterich et al. (2019) produced an ensemble of scenario simulations with a regional coupled ocean–atmosphere model, called RCA4-NEMO. The atmospheric model RCA4 was run at a  horizontal resolution of 0.22 degrees and 40 vertical levels in the EURO-CORDEX domain (Jacob et al. 2014). Coupled to it is the North Sea-Baltic Sea model NEMO-Nordic at a horizontal resolution of two nautical miles and 56 vertical levels. Atmosphere and ocean are coupled at a 3 h frequency. This coupled model has been applied to downscale eight different CMIP5 models driven by three Representative Concentration Pathways (RCPs) each. In this study, four CMIP5 models (MPI-ESM-LR, EC-Earth, IPSL-CM5A-MR, HadGEM2-ES) and the greenhouse gas concentration scenarios RCP 4.5 and RCP 8.5 were used.

Then, the atmospheric field data at hourly to 6-hourly frequency are extracted from these regional climate simulations to force the RCO-SCOBI model over the 1976–2099 period. This regional model cascading is fully described in Meier et al. (2021). All climate change simulations analyzed in this study with their different forcings are listed in Table 2.

Table 2 List of climate change simulations
Full size table

2.5 Statistical methods

2.5.1 Trend calculation

First, the monthly average of fields (i.e. air temperature at 2 m, zonal and meridional components of surface winds, cloudiness and sea surface temperature) is computed from the 2-daily snapshot outputs model. Then at each grid point the linear trend is computed 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, thus the extreme values will have less influence on the trend calculation than with the least squares method. 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%.

2.5.2 Hierarchical clustering algorithm

To identify homogeneous areas of SST trends in the Baltic Sea, a statistical classification method is used. To that end, the seasonal SST trends are first computed over two periods: 1850–2008 and 2006–2099. Then all the longitude, latitude pairs of the seasonal SST trends 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, then it calculates for each class 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. The dendogram shows how is organized the classification and inform on the distance between each class. One of the advantages of hierarchical clustering is the inclusion of subdivisions in the upper divisions, so increasing the number of classes considered increases the accuracy of the clustering without changing the results.

2.5.3 Classification tree

A classification tree (using rpart package in R environment; Therneau et al. 1997) is applied to the seasonal trends of the following variables (explicative variables): air temperature at 2 m, zonal and meridional components of surface winds, cloudiness, and absolute (not trend) ice concentration. This method identifies the hierarchical thresholds discriminating the SST trends clusters identified previously (predictive variables). These explicative variables are chosen because they modify the heat fluxes at air-sea interface. An auto-correlation analysis was performed for removing the auto-correlated variables (r > 0.8), the variables retained are summarized in Table 3. The algorithm constructs a classification tree from the explicative variables that minimizes the classification error (based on Gini impurity index) of the predictive variables. Each split is defined by a simple rule based on a single explicative variable. For simplicity, the classification tree is cut once each SST trend class is defined. Thus at the end, a part of longitude, latitude pairs are well classified (all longitude, latitude pairs belonging to the correct SST trends clusters) and another part is bad classified. Finally, for estimating the errors done by the classification tree it is possible from all the longitude, latitude pairs well classified to reconstruct a map of SST trends clusters. Given the diversity of futures simulated by the four atmospheric forcings, it is difficult to extract a consensus tree for the future period. Therefore we limited this analysis to the historical period in REF.

Table 3 List of explicated variables used in classification tree. The auto-correlated variables (r > 0.8) were removed
Full size table

2.5.4 Upwelling frequency

The upwelling frequency has been calculated using the same method as proposed in 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 a visual analysis and has shown very good skills (Lehmann et al. 2012).

3 Results: past period

3.1 SST trends

Figure 2 displays the annual and relative seasonal SST trends over the Baltic Sea computed over the 1850–2008 period. The SST trends are characterized by a spatial average of 0.047 K.decade-1 and a weak spatial standard deviation of 0.008 K.decade−1 (Fig. 2a). This spatial distribution of annual SST trends exhibit an east/west gradient of ~  + 0.02 K.decade−1 with stronger warming in the eastern Bothnian Sea and Bothnian Bay and stronger warming in the western parts of southern Baltic. This rather weak spatial gradients of annual trend hide a strong heterogeneity at the seasonal scale between the Baltic Sea sub-basins. The relative SST trend is highest in the northern Baltic Sea in summer (~ 0.04 K.decade−1; Fig. 2d) and minimum in winter (~ − 0.03 K.decade−1; Fig. 2b) while the amplitude of seasonal changes are weaker in the southern part of the Baltic Sea and are centered on spring and fall (~ + 0.02 K.decade−1 and ~ − 0.01 K.decade−1).

Fig. 2
figure 2

Annual SST trend and relative seasonal SST trends (in K.decade−1) computed over the 1850–2008 period (calculated as the seasonal SST trend minus the spatially averaged annual SST trend). a Annual, b DJF, c MAM, d JJA, e SON. The hatched areas represent the regions where the trends are significant at a threshold of 95% from a Mann–Kendall non-parametric test

Full size image

3.2 Spatial classification of SST trends

As illustrated the seasonal SST trends display clear differences between Baltic Sea sub-basins. To separate these sub-basins into distinct spatial clusters of SST trends a hierarchical clustering algorithm is applied. The dendogram has been cut at two thresholds (three clusters and five clusters, Fig. 3a) in order to vary the spatial scale and thus the level of precision.

Fig. 3
figure 3

a Dendogram associated to b, c spatial clustering of seasonal SST trends over 1850–2008 computed with a hierarchical classification algorithm. Two thresholds are considered: b 3 classes and c 5 classes

Full size image

The first threshold considered separates the Baltic Sea into 3 sub-regions of SST trends (Fig. 3b). Cluster 1 (C1 in green) encompasses the Bothnian Sea and the Gulf of Finland, cluster 2 (C2 in blue) is located in the Bothnian Bay and cluster 3 (C3 in red) is located in the southern half of the Baltic Sea including the Gotland Sea, the Gulf of Riga, the Arkona and Bornholm basins and the Danish straits (see Fig. 1 for the location of each sub-region). The second threshold (5 clusters; Fig. 3c) divides C3 in three new sub-regions of SST trends: the cluster 3a (C3a in red; Fig. 3c) includes the Arkona Basin and the Danish straits, cluster 3b (C3b in yellow) is the western part of Gotland Sea and cluster 3c (C3c in purple) is the eastern part and the Gulf of Riga, while C1 and C2 remain unchanged (Fig. 3b).

C1 and C2 are characterized by a maximum of SST trend in summer and a minimum in winter, while maximum and minimum of SST trend occur respectively in spring and summer for C3 (Fig. 4a). The strongest seasonal amplitude occur for C2 with a difference of 0.09 K.decade−1 between the winter and the summer, while C1 and C3 display a weaker difference with 0.035 K.decade−1 and 0.021 K.decade−1 of seasonal amplitude, respectively.

Fig. 4
figure 4

Boxplots representing the spatial distribution for each cluster of SST trends (in K.decade−1). a 3-classes and b 5-classes. The whiskers represent the 10th and the 90th percentiles and the outliers are not represented

Full size image

C3b and C3c exhibit a SST trend very similar in winter and spring and are mainly separated by their SST trends in summer and fall with a larger value for C3b than for C3c (0.058 vs 0.036 K.decade−1 and 0.044 vs 0.033 K.decade−1 respectively; Fig. 4b). Finally the SST trends for C3a are larger than for C3b and C3c at all seasons (except in summer for C3b) with a maximum difference of 0.19 K.decade−1 in winter (Fig. 4b).

3.3 Drivers of SST trends classes

As shown by Meier et al. (2019a, b) there is a significant correlation between long-term air and sea water temperature trends. However there are some discrepancies between the trends of air temperature and SST especially according to the season and sub-basins as illustrated by the spatial pattern correlation calculated here and varying between − 0.66 (winter) and 0.55 (fall). The other identified main drivers of SST trends that modify the heat and radiation fluxes at the surface are cloud cover, sea ice concentration and u–v wind components. For determining the importance of these variables in the SST trend clustering, a classification tree was generated based on their seasonal linear trend.

The results of these classification trees are displayed in Fig. 5. The best predictor separating clusters C1, C2 and C3 is the ice concentration in spring and summer. Then, the best predictor to separate C3a and C3b clusters is the air temperature in winter and the meridional wind trends in fall to separate C3b and C3c.

Fig. 5
figure 5

Decision tree explaining the SST trend classes from the seasonal trends of explicative variables. a 3 classes and b 5 classes. The colored numbers represent the most frequent class at each node, the small numbers below the colored numbers are the number of (longiutde, latitude) pairs for each class with from left to right the classes 1 to a 3 and b 5, and the bold text at each node represents the threshold of the variable separating two branches

Full size image

It is possible to evaluate the error made by this statistical method to make a first assessment of the robustness of our results. Using all longitude, latitude pairs well classified by the classification tree, we reconstructed the SST trend classes to visualize the error made (Fig. 6). The left panel show an error rates of 8% achieved by the classification tree to reconstruct the SST trend clustering. In other words, this error rate indicates that 92% of all longitude, latitude points in the SST clustering can be determined with only 2 thresholds based on summer and spring ice concentration. These statistical results can be explained by two physical processes involving a warming regulation by ice cover. First ice cover isolates the sea and the air, thus limiting the heat fluxes; secondly the albedo of the sea ice is higher than the albedo of the sea, thus reflecting more of the incident solar radiation. During winter and spring the average sea ice cover extends to the southern Bothnian Sea isolating the underlying sea water. Conversely, during the melting season (in spring in the Bothnian Sea and in summer in Bothnian Bay), the surface albedo decreases, the solar radiation increases and the air–sea coupling increases with an increase of sensible and latent ocean heat gain, yielding a strong enhancement of the SST warming. Thus, the results of classification tree seem physically consistent and implies that ice concentration is the key parameter explaining the seasonal warming differences between the three identified sub-regions during the historical reconstructed period.

Fig. 6
figure 6

Reconstruction of SST trends clusters from the well classified (longitude, latitude) pairs by the regression tree. The both thresholds are considered: a 3 classes and b 5 classes

Full size image

The error rate of the classification tree for 5-classes is higher (23%) and deserves further analysis. The classification tree shows that intense warming in Arkona basin and Danish straits (C3a) in winter and the larger warming in the western (C3b) than in the eastern (C3c) Baltic proper in autumn is explained by an increase of air temperature and a decrease in the v-component of the wind at these respective seasons. To test the findings by the classification tree and for assessing the potential role of trends in surface winds, air temperature and cloud cover on SST trends, the WIND, TAIR and CLOUD sensitivity simulations are analyzed and compared to REF (see Table 1). The top and bottom panels in Fig. 7 display the winter and fall SST trends difference between REF and WIND and between REF and TAIR respectively. These results reveal that the air temperature trends yields an enhanced warming in C3a compared to C3b and C3c in winter (0.015 K.decade−1) representing 81% of the warming contrast between C3a and C3b, C3c in REF (Fig. 7c). However, the surface wind trends yields a homogeneous decrease of SST trends in C3a, C3b, C3c in winter characterized by a very weak warming contrast of 0.004 K.decade−1 between C3a and C3b, C3c (Fig. 7a). Therefore these results highlight the importance of air temperature trends in winter in agreement with the results of the classification tree. The air temperature effect can be explained by a lower increase in sensible heat fluxes (oriented towards the atmosphere in winter) in C3a than in C3b, C3c (0.26 vs 0.31 W.m−2.decade−1; black lines separates C3a and C3b, C3c in Fig. 8a) due to a reduced air-sea temperature anomalies in C3a.

Fig. 7
figure 7

SST trend difference (in K.decade−1) between a, b REF and WIND and c, d REF and TAIR in a, c winter and b, d fall computed over the 1850–2008 period

Full size image
Fig. 8
figure 8

a DJF trends in sensible heat fluxes (in W.m−2 decade−1), b SON upwelling frequency trends (in % decade−1) and c SON trends in sensible and latent heat fluxes (in W.m−2 decade−1). All these trends are computed over the 1850–2008 period in REF. The black lines delimit a the boundaries of C3* and c the boundaries between C4 and C5

Full size image

In addition, Fig. 7b reveals that the wind trends in Baltic proper yields an east/west warming contrast in fall. This warming contrast (calculated between C3b and C3c) in fall is of 0.008 K.decade−1 representing 73% of warming contrast between C3b and C3c in REF (0.011 K.decade−1). Conversely, Fig. 7d shows that air temperature trends induces a weaker warming contrast between C3b and C3c of 0.004 K.decade−1 (i.e. 36% of warming contrast in REF). These results are in agreement with those of the classification tree separating C3b and C3c. The wind effect on SST trends can be explained by two processes. First, Fig. 8b shows significant negative trends in upwelling frequencies in fall along the Swedish coast exceeding 0.5% decade−1. The reduced offshore transport of cold anomalies could explain the higher warming exhibits in western part than in eastern part of the Baltic proper. Another possible mechanism is a change in wind speed leading to changes in heat fluxes and horizontal transports. The Fig. 8c shows an increase of heat fluxes (latent and sensible) in the C3b and negative changes in C3c (− 0.08 vs 0.09 W.m−2 decade−1; black lines separates C3b and C3c in Fig. 8c). In fall, the heat fluxes are directed towards the atmosphere, thus an increase of these fluxes means an enhancement of loss heat for the sea and a reduced warming. Finally, the analysis of the CLOUD simulation strongly implies that cloud cover trends have no significant effect on SST trends (not shown).

In summary, our analyses clearly show that (1) the ice cover, by modifying albedo and air-sea heat fluxes, explains the north/south gradient of warming on a seasonal scale, (2) the enhanced warming in the western Baltic proper compared to the eastern part in fall is mainly due to changes in surface winds supported by a subsequent decrease in the upwellings frequency along the Swedish coast and an increase of heat fluxes in the eastern Baltic proper, (3) the enhanced warming in the southwestern Baltic Sea in winter is mainly explained by the trends in air-sea temperature anomalies leading to a lower increase in sensible heat fluxes (oriented to the atmosphere at this season, i.e. heat loss for the sea) compared to the surrounding areas.

4 Results: future period

Considering that climate change has a non-linear effect on SST trends due to several climate feedback loops, climate simulations were performed over the 1976–2099 period to assess how our previous results are modified for the future period. These simulations were carried out for two RCP scenarios. However, it was found that the spatial pattern of SST trends between both RCP scenarios are extremely similar (not shown) as already mentioned in Dieterich et al. (2019) and Gröger et al. (2019). Therefore, we will only describe and discuss the changes for the RCP8.5 scenario in the following.

4.1 Relative SST trends

Figure 9 top panels show the relative annual and seasonal trends in SST over the future period for the RCP8.5 scenario and their differences with the past period. The spatial pattern of relative annual SST trends corresponds to a tripolar warming with maximum SST trends in the Bothnian Sea and lower warming in the north (Bothnian Bay) and south of this sub-region (Fig. 9a). In addition, noteworthy is the lower warming of the shallow coastal zone. At the scale of the sub-basins and seasons, other spatial patterns emerge. In the Bothnian Bay, the warming is lower than the annual spatial mean in winter and spring, and higher in summer and fall, while in the Bothnian Sea the warming is always higher than the annual spatial mean for all seasons (Fig. 9b–e). In the Baltic proper the warming is higher than the annual spatial average in spring and summer, and the warming in the Danish straits is always lower than the annual spatial mean as in the other shallow coastal regions (Fig. 9b–e).

Fig. 9
figure 9

a–e Relative annual and seasonal SST trend (in K.decade−1), calculated as the SST trend minus the spatial averaged annual SST trend, over 2006–2099 period in RCP8.5. f–j Difference (future–past) of relative SST trend normalized by the annual, spatially averaged trends for 1850–2008 and 2006–2099. From left to right: Annual, DJF, MAM, JJA, SON

Full size image

These relative seasonal trends differ from those calculated over the past period (Fig. 3b–e), and to analyze these differences we show in Fig. 9 bottom panels the difference (future minus past) of the relative SST trends normalized by their annual spatially averaged trends. The normalization enables these relative changes to be compared independently of the mean warming associated with each period. During winter and spring the future period exhibits a stronger warming in the northern part and less intense in the southern part of the Baltic Sea, reducing the north/south warming gradient simulated in the past period (Fig. 9g, h). In summer, the warming is less intense in the northern part and the western part of the Baltic proper (Fig. 9i). Finally in fall, the warming is everywhere more intense in future than in past except in the Danish straits (Fig. 9j). These changes in summer and fall damp the east/west warming gradient found in the past period in the Baltic proper. To conclude, these results show that the spatial and seasonal trends in SST are not stationary with an acceleration of warming in the North Baltic Sea at all seasons except in summer and a deceleration in the South Baltic Sea in winter and spring. These changes imply that the warming in the Baltic Sea is not a linear response to increasing greenhouse gas concentration due to non-linear climate feedback processes.

4.2 Spatial classification of SST trends

In order to analyze how SST trends found in the past are in the future climate, we performed a spatial clustering of seasonal SST trends using the future simulations. This analysis is made for the four CMIP5 models forcing and for the multi-model mean. Figure 10 displays the results of this classification for the two thresholds (3 and 5-classes) and reveals significant changes compared to the past period. Furthermore, these changes are consistent among CMIP5 model forcings as indicated by the hatched areas that represent the region where 3 out of 4 models agree.

Fig. 10
figure 10

Spatial classification of SST trends computed with a hierarchical classification algorithm over 2006–2099 period for a–e 3-classes threshold and f–j 5-classes threshold under RCP85 scenario. This classification is computed from the a, f multi-model mean (MMM), b, g MPI-ESM-LR (MPI), c, h EC-EARTH, d, i IPSL-CM5A-MR (IPSL), e, j HadGEM2-ES (HadGEM2) models. The hatched areas in panels a and f represent regions where 3 out of 4 models agree

Full size image

C2 is now confined to the northern part of the Bothnian Bay, and thus C1 has expanded northward. This difference is probably explained by the reduction in spring ice cover, with an ice concentration which exceeds 2% in only one CMIP5 model forcings and restricted to the extreme northern part of Bothnian Bay (hatched area bounded by the red lines; Fig. 11b–e). Thereafter, C3 retracted southward to about 57° N while C1 expanded as much. During the past period the distinction between C1 and C3 was related to lower warming in C1 than in C3 in winter and spring due to the presence of sea ice reaching a concentration of ~ 7% in the Åland Sea in spring (hatched area bounded by the red lines; Fig. 11a). As the sea ice disappeared in this season in the Bothnian Sea, the albedo is reduced and the air–sea coupling via heat fluxes is enhanced, leading to an intensification of the warming in this region compared to the past period. Thus the warming contrasts between the Baltic proper and the Bothnian Sea are reduced, which explains the southward extension of C1.

Fig. 11
figure 11

Ice concentration (in %) in spring for a past and b–e future period. For the future period the four CMIP5 model forcings are considered b MPI, c EC-EARTH, d IPSL, and e HadGEM2.1. The hatched areas bounded by the red lines show the regions where the ice concentration is higher than 6% in a and higher than 1% in b–d

Full size image

For 5-classes classification and the multi-model mean, C3b and C3c are now mainly separated along a zonal axis instead of a meridional axis for the past period (Fig. 3b), with C3b in the northern part of the Baltic proper and C3c in the southern part (Fig. 10f). This new cutting for the multi-model mean is coherent with the results of Fig. 9i, j, that shows a more intense warming in the eastern part of Baltic proper in future than under past conditions, reducing the east/west warming contrast during these seasons. The previous separation between C3b and C3c was explained by a decrease in the frequency of upwelling and a reduced heat loss in the western part of Baltic proper in fall (Fig. 12a, f). According to future simulations, the decrease in the upwelling frequency in fall appears only for one CMIP5 model forcing, while one of the models projects an increase (Fig. 12b–e). Furthermore, in two CMIP5 model forcings the trends in heat fluxes are lower in western than in eastern Baltic proper (Fig. 12h, j) and the opposite for the two others CMIP5 models forcings (Fig. 12g, i). These changes in upwelling frequencies and heat fluxes tend to attenuate the warming contrasts between the western and eastern Baltic proper and may therefore explain this new cutting for C3b and C3c.

Fig. 12
figure 12

Upwelling frequencies trends (%.decade−1) in fall for a past and b–e future period. Trends in sensible and latent heat fluxes (in W.m−2.decade−1) in fall for f past and g–j future period. Trends in sensible heat fluxes (in W.m−2.decade−1) in winter for k past and l–o future period. For the future period the four CMIP5 model forcings are considered, from left to right: MPI-ESM-LR (MPI), EC-Earth, IPSL-CM5A-MR (IPSL), HadGEM2-ES (HadGEM2)

Full size image

Finally, C3a is now confined to the Kattegat in the multi-model mean (Fig. 10f) whereas it extended to the Arkona Basin in the past. This new cutting is also coherent with the results of Fig. 9g, h, that show a less intense warming in Arkona and Kattegat regions in future than in past conditions, reducing the east/west warming contrast during winter and spring. Nevertheless this new cutting is not a common picture shared by individual clustering (Fig. 10g–j), suggesting that these results could be an artifact related to the multi-model mean. For the past period, the sensitivity experiments have revealed that the air temperature trends was the main factor explaining the enhanced warming in C3a, related to lower changes in sensible heat fluxes for C3a than for C3b, C3c. This analysis of trends in sensible heat fluxes for the future period shows the same kind of results, with a lower increase in sensible heat fluxes in southwestern Baltic Sea than in surrounding areas (Fig. 12k–o). These results are consistent with the individual clustering for the future period that show no significant changes for C3a.

In summary, the climate change will alter likely the spatial pattern of seasonal SST trends, with a more (less) intense warming in the northern Baltic Sea in winter and spring (summer) in response to the vanishing of sea ice, thus damping the north/south warming gradient observed over the past period. Furthermore, the east/west warming gradient in Baltic proper should be reduced in the future due to reduced changes in upwelling frequencies and in heat fluxes. Finally, the east/west warming gradient in the southwestern Baltic Sea should not change due to the persistence of relative decrease in sensible heat fluxes in the western part. Overall, contrasts in seasonal SST trends between Baltic Sea sub-basins should be attenuated compared to the period, in response to the reduced influence of sea ice-related climate feedback with the vanishing of ice cover.

5 Discussion, summary and conclusions

5.1 Discussion

Our study is a complement to the work by Kniebusch et al. (2019) on SST trends over past period and extends it for future scenarios. From statistical analyses, we have identified five sub-regions of SST trends, and related them to common drivers explaining these spatial and seasonal differences. In agreement with the Kniebusch et al. (2019)’s assumption, our results show that differences in SST trends between the northern and southern Baltic Sea are explained by summer and spring ice cover. Furthermore, our study reveals that the stronger SST trends in the western part than in the eastern part of the Baltic proper are supported by both reduced upwelling frequencies along the Swedish coast and an increase of heat fluxes in the eastern Baltic proper in response to surface wind changes. Contrary to the analyses of Kniebusch et al. (2019) who suggested also an effect of cloud cover, the analysis of a sensitivity experiment with a fixed cloud cover (CLOUD) suggests only negligible effects of long-term cloud cover changes on SST trends. However we can questioned the reliability of solar radiation trends in the atmospheric forcings. Indeed, a recent study suggested that the underestimation of SST trends in climate models over the recent period (1983–2014) would be due to the strong decrease in regional industrial aerosol emissions over the same period, changes that are not accounted for in these models (Barkhordarian et al. 2016). The reduction in aerosol concentrations enhances the radiative forcing through several mechanisms, such as changes in cloud albedo, cloud lifetime or cloud cover (e.g. Myhre et al. 2013). The solar radiation trend averaged on the Baltic Sea and calculated over 1980–2008 in our historical simulation (+ 1.7 W.m−2.decade−1) is consistent with observations over 1983–2014 period (~ 2 W.m−2.decade−1; Barkhordarian et al. 2016). Thus, despite the non-inclusion of these processes in the HiResAFF atmospheric forcings, the variability of the solar radiation is correct, probably due to the large-scale variability at the boundaries of the domain which is well represented. However, future changes in regional industrial aerosol emissions are not accounted for in the climate change scenarios, and could imply uncertainties in solar radiation trends and in SST trends, thus modifying partly the results of this study.

In addition, several other sources of variability (e.g. internal, natural, atmospheric forcings) may modify our results for the future period. To assess their effect on SST trends we used a pool of climate change simulations performed for four atmospheric forcings in three sea level rise scenarios and two nutrient load scenarios, i.e. twenty four climate change simulations for each RCP scenarios. Considering that sea level rise and nutrient load scenarios have a small effect on SST trends (Löptien and Meier 2011) we used these six simulations for each atmospheric forcings to assess the internal variability. These simulations are identical to those in Meier et al. (2021). The top panels in Fig. 13 display the annual SST trends uncertainties related to the atmospheric forcings and internal variability. Despite the robustness of the SST trends spatial pattern (p-value < 0.05 everywhere) found in Fig. 9a–e, this figure reveals an important dependency of SST trends to atmospheric forcings with a spread (spatially averaged) of ± 0.13 K.decade−1 from the multi-model mean and a spatial deviation standard of 0.014 K.decade−1. However the effect of internal variability is half the size with a spread of ± 0.06 K.decade−1 from the multi-model mean while the spatial deviation standard is similar (0.015 K.decade−1). These results highlights that the uncertainties in the mean intensity of warming are mainly due to atmospherics forcings but that the uncertainties in the spatial pattern of warming are due to both internal variability and atmospheric forcings. In addition, a calculation of the SST trends by 30-year slice period every one year over the entire period (2006–2099) shows that annual SST trends are variable over time. In RCP4.5 scenario, the natural variability appears to modulate these trends with successive periods of increasing and decreasing SST trends with a period of about 30 years. However, in the RCP 8.5 scenario, SST trends gradually increase over the first 40 years of the period reaching a maximum of 0.5 K.decade−1 over the period 2045–2075, before declining slightly from 2060 onwards, as in the RCP 4.5 scenario, a result of the pronounced natural variability in this scenario as well.

Fig. 13
figure 13

SST trends deviation from the multi-model mean calculated from the a four atmospheric forcings and b the three SLR scenarios and two nutrient load scenarios. The spatially averaged deviation is noted in the bottom right-hand box. c SST trends calculated by 30-year slice period every one year over the entire period (2006–2099) in RCP4.5 (blue line) and in RCP8.5 (red lines) scenarios. The lines represent the multi-model mean and shading around the multi-model mean represent the deviation due to the atmospheric forcings

Full size image

In line with recent studies (e.g. Meier et al. 2019a, b; Kniebusch et al. 2019) that have shown that short- and long-term SST variability is dominated by air temperature variability, our study focused primarily on atmospheric processes to explain spatial and seasonal SST trends. Nevertheless, we showed that changes in meridional winds were consistent with changes in upwelling frequency and could explain differences in warming between sub-basins, thus highlighting the importance of oceanic processes. Indeed, the modulation of several oceanic processes not considered here (e.g. vertical mixing, Ekman pumping, horizontal and vertical convection) in response to changes in atmospheric forcing can explain seasonal and local differences in warming. For instance, changes in stratification induced by changes in winds and melting sea ice alter the depth of the mixed layer (Meier et al. 2021) and thus modify the turbulent mixing between surface and sub-surface water and finally influencing the local warming values. Similarly, circulation changes due to changes in surface winds or vertical thermohaline structure could modify the distribution of heat between sub-basins. The study of these processes is outside the scope of this paper, but will be evaluated in a future study focusing on vertical changes in temperature trends.

Another source of uncertainty is the lack of air–sea coupling in our simulations that can alter the local warming and modify our results. First, Gröger et al. (2015, 2021) showed that there is a thermal response to air–sea coupling that influence the summer and fall air-sea heat fluxes. At these seasons, wind events draw cooler waters from depth to the surface. Upwelling cools the atmospheric boundary layer and so increases the static stability of the atmosphere with a negative effect on wind speed (summer short circuit). In addition, Ho-Hagemann et al. (2017) showed in CCLM model that the air–sea coupling enhances the moisture convergence over the Baltic Sea reducing the summer dry bias, which likely modify the shortwave radiation changes and thus the warming. Likewise, the changes in upwelling frequencies and the reduced SST trends associated, modify the land/sea thermal contrast and thus the low-level atmospheric circulation that likely feedbacks on upwelling efficiency. Regional coupled ocean–atmosphere configurations have been developed recently (e.g. Dieterich et al. 2019) and could be used to study the effect of air–sea coupling on climatological and extreme temperature changes. Finally, given the importance of ice concentration on SST trends, the reliability of its modeling is crucial. Although the skill of the RCO model on this point is correct (e.g. Meier 2001; Kneibusch et al. 2019), the uncertainty associated with the evolution of this variable in the future could alter our results.

The results of this study allowed us to define areas of homogeneous SST trends, which could be useful in terms of marine ecosystems or public management of these sub-areas. The 5 sub-regions identified in this study show significant differences in SST trends by season, which can lead to very different pressures on marine ecosystems. For instance, a recent study shows that the temperature threshold triggering initial spawning of herring in the southwestern Baltic Sea is 3.5–4.5 °C (Polte et al. 2021). Their results revealed that the late seasonal onset of cold spells, the corresponding lengthening of the larval hatching period and early larval hatching peaks significantly reduced larval production and ultimately led to a reduced abundance of juveniles. In a general way, the critical time for most of marine species is the larvae period which is strongly affect by the regional and seasonal thermal conditions (e.g. Dodson et al. 2019; Hüssy et al. 2012), therefore regional and seasonal changes in SST may be a more relevant indicator to assess the thermal stress on the reproductive capacity of marine ecosystems. Our findings are also interesting in term of public management with the tourism which is an important economic sector in the region (HELCOM 2018). Indeed, there is many environmental indicators that favor tourism and are affected by the SST. For instance, the summer cyanobacteria bloom or the jellyfish reproduction (e.g. Holst 2012; Treible and Condon 2019) are related to the thermal conditions and have negative influence on the attractiveness of tourist areas. The information on SST trends can further be used for the spatial planning of the coastal Baltic Sea areas.

5.2 Conclusions

This study investigates the spatial and seasonal heterogeneities of SST trends in Baltic Sea under past and future conditions. A historical simulation (1850–2008) was performed with a regional oceanic model using HiResAFF atmospheric forcings. This simulation revealed a spatially averaged SST trend of ~ 0.05 K.decade−1 over this past period but with strong variability according to the seasons and the sub-basins of the same magnitude as the spatial average. A spatial clustering was then carried out for separating precisely the sub-regions of homogeneous SST trends across the seasons. This analysis was done for two thresholds, according to the first threshold the Baltic Sea is separated in 3 sub-regions: C1 encompasses the Bothnian Sea and Finland Gulf, C2 corresponds to Bothnian Bay, C3 encompasses all the southern part of Baltic Sea. According to the second threshold C1 and C2 remain identical but C3 is now divided in 3 new sub-regions: C3a encompasses Arkona Basin, the Danish straits and Kattegat, C3b corresponds to the western part of Baltic proper and C3c the eastern part. For analyzing the drivers of the seasonal and spatial SST trends differences identified by the clustering method, a classification tree was performed based on the seasonal trends of 5 explicative variables: air temperature, zonal and meridional components of surface winds, cloudiness, and ice concentration. This analysis revealed that the ice cover in spring and summer, isolating the sea to the atmosphere, explains the SST trends differences between C1, C2 and C3. Then, the more intense warming in C3b than in C3c during summer and fall are related to change in surface winds, yielding a decrease of the upwelling frequencies along the Swedish coast and an increase of heat fluxes in the eastern Baltic proper. Finally, the more intense warming in C3a than in C3b and C3c in winter and spring were related to the air temperature trends in winter associated with a relative decrease of sensible heat fluxes. A reconstruction of the initial SST trends cluster based on the classification tree showed a low rate of error and sensitivity experiments with fixed surface winds, air temperature and cloud cover confirmed these findings; reinforcing the robustness of our results. Climate change simulations over 1976–2099 were carried out with four CMIP5 forcings to assess possible changes in the spatial pattern of seasonal SST trends from past conditions. Spatial clustering based on these new simulations showed some consequent differences from past conditions. C1 extended in north and south, confining C2 in the extreme northern part of Bothnian Bay. C3b and C3c are separated along zonal axis instead of a meridional axis in past conditions, and C3a is confined to the Danish straits and Kattegat. The differences in C1, C2 and C3 between these two periods are associated with a reduction of north/south warming gradient in response to a strong decrease of ice concentration modifying the surface albedo and the air-sea coupling. The differences in C3b and C3c are associated with a reduction of east/west warming gradient in Baltic proper which can be explained by changes in upwelling frequencies and heat fluxes. These analyses have been done for the multi-model mean and the four CMIP5 models forcings. They show consistent changes through these different forcings, reinforcing the robustness of our results. However the changes in C3a for the multi-model mean are not shared by the clustering with the four CMIP5 models forcings, suggesting that these changes are only an artifact of the multi-model mean.

Finally, the study allows us to draw some conclusions for the SST trends in Baltic Sea over past and future periods. First, the annual SST trends hide important spatial heterogeneities at the seasonal scales that are important to consider for a precise assessment of their impacts on the marine ecosystems and the societies. Second, these seasonal trends are mainly driven by three atmospheric variables: the ice cover, the air temperature and the surface winds that modify the balance of heat fluxes at the air-sea interface. However the long-term trends in cloud cover and the changes in shortwave radiation associated have a secondary effect on SST trends. Third, climate change will alter the intensity of past SST trends but also their spatial pattern. Overall, contrasts in seasonal SST trends between Baltic Sea sub-basins should be attenuated in the future. RCO model project an acceleration of warming in northern part of Baltic Sea in winter and the opposite in summer due to the reduced effect of nonlinear climate feedbacks such as sea ice albedo. In Baltic proper, most of the simulations show a faster warming of the eastern part than the western part in summer and autumn which could be related to changes in upwelling frequencies and heat fluxes. Fourth, for a given RCP scenario the atmospheric forcings is the main source of uncertainty on the intensity of warming, while the uncertainties on the spatial pattern of warming are equally due to the atmospheric forcings and the internal variability. Finally, all these changes are strongly modulated by the natural variability.