Introduction

Among the world’s coastal seas, the Baltic Sea (Fig. 1) has the largest annual maximum hypoxic area (70,000 km2), followed by the northern Gulf of Mexico (23,000 km2), the north-western Black Sea (20,000 km2) and the East China Sea (12,000 km2)1. Hypoxia occurs when biogeochemical oxygen consumption exceeds ventilation by water transports and the dissolved oxygen concentration at the bottom falls below 2 mL O2 L−12. The spread of hypoxia has considerable consequences for ecosystem functioning and structure3,4 and is today perhaps one of the most severe environmental threats to coastal oceans worldwide1,3,5,6,7. As in many coastal seas, nutrient inputs into the Baltic Sea have been high since the 1950s8, the result of a growing human population, intensified agriculture and changing land use8,9. In the 1980s, nutrient loads reached their highest level, including peaks in nitrogen (N) and phosphorus (P) loads that were roughly three and seven times larger than preindustrial values, respectively8. Although the loads have decreased since the 1980s, present-day N and P loads are still about two and three times larger than preindustrial values, respectively, and the Baltic Sea has yet to return to its pre-1950s status8. Under anoxic conditions, phosphorus fluxes from the sediments into the water column are considerably larger than those occurring under oxic conditions10,11. As a result, nitrogen fixation has increased and earlier and more frequent cyanobacterial (blue–green algal) blooms have reinforced eutrophication12. This internal nutrient source counteracts the reduction in external loads that have been achieved since the 1980s8. Despite an improvement in water quality, indicated by ecosystem indicators in the shallow Kattegat (for the location, see Fig. 1), the offshore oxygen conditions in the Baltic Sea have yet to change substantially13. Thus, the environmental conditions in the Baltic Sea today are still characterised by (1) deoxygenation and hypoxia both in the deep offshore waters2,14 and in coastal zones15,16, (2) intensified cyanobacterial blooms17 and (3) reduced water transparency18.

Fig. 1: Baltic Sea bathymetry.
figure 1

Water depths (in m) and the locations of international monitoring stations listed by the Helsinki Commission (HELCOM, http://helcom.fi) at which oxygen concentrations are regularly recorded. The current annual sampling frequency is illustrated by coloured circles. The domain of the Baltic Sea model is limited, with open boundaries in the northern Kattegat (black line). The so-called Baltic proper comprises the Arkona Basin, the Bornholm Basin and the Gotland Basin.

Full size image

In addition to eutrophication, the Baltic Sea has also been affected by the changing climate (see also Supplementary Note 1). During the period 1982–2006, the sea surface temperature (SST) in the Baltic Sea increased by 1.35 °C, the largest warming among all coastal seas worldwide19. For the 21st century, projections of the Baltic Sea suggest higher water temperatures, less sea-ice cover, higher sea levels and increasing acidification compared with the present climate20.

Further spreading of hypoxia will presumably depend on the intensity of both human activities on land (e.g., fertiliser use, wastewater treatment, burning of fossil fuels) and climate change20. Scenario simulations for the Baltic Sea have suggested an expansion of the hypoxic area in a future climate because of increased nutrient loads owing to enhanced river runoff, reduced oxygen fluxes from the atmosphere to the ocean and the accelerated recycling of organic matter in response to higher water temperatures21. Furthermore, more frequent and longer-lasting events were projected in areas with intermittent hypoxia22. The expanding areas of anoxic bottom water will likely lead to larger phosphorus fluxes from the sediment into the water column23, whereas the feedback mechanisms inherent in biogeochemical cycling may prevent the system from returning to its pre-1950s condition24,25.

To counteract nutrient enrichment in the water, in 2007 the countries bordering the Baltic Sea agreed upon a plan to reduce nutrient loads, the so-called Baltic Sea Action Plan (BSAP)26. However, the country-wise-allocated, maximum allowable loads set by the BSAP to restore the good ecological status of the Baltic marine environment were calculated without taking into account the impact of a changing climate. Restoration plans similar to the BSAP have been implemented for Chesapeake Bay (the Chesapeake Bay Programme, since 1983; see https://www.chesapeakebay.net/) and the Gulf of Mexico (the Gulf Hypoxia Action Plan, in 2008; see https://www.epa.gov/ms-htf/gulf-hypoxia-action-plan-2008). However, these nutrient load abatement strategies similarly failed to consider a changing climate, despite reports that global warming may accelerate the expansion of hypoxia21,27,28. Although numerical models have been developed to investigate the combined impact of nutrient load and climate change29,30,31, the spread of the projected changes in multi-model ensembles, i.e., the “uncertainty” of regional climate projections, is considerable30,32. The uncertainties are related to (1) natural variability, (2) differences between the climate models (henceforth, model uncertainty), (3) the unknown future greenhouse gas concentration pathway, expressed by various plausible representative concentration pathways33 (RCPs) and (4) the unknown future nutrient load pathway, expressed by various plausible shared socioeconomic pathways34,35 (SSPs). A general drawback of previous ensembles of regional scenario simulations was that they were too small (<10 members), because of limited computational resources29. Consequently, studies on the uncertainties of future projections for coastal seas remain rare32. However, as decision-makers in environmental protection agencies and similar organisations increasingly seek quantitatively reliable climate information, robust uncertainty estimates of projections are urgently needed. To account for these needs, we developed a method based upon an ensemble of 48 scenario simulations to analyse and quantify the relative importance of uncertainties and their sources. We then used the ensemble spread of the changes for subgroups of the ensemble and an approximation of the internal variability to quantify the various uncertainties.

Our results demonstrate that the importance of natural variability in the detection of trends has so far been underestimated. For variables, such as SST, bottom oxygen concentration and Secchi depth (SD) as a measure of water transparency, for example, this natural variability is the largest source of uncertainty during the coming decades. An exception is salinity, in which case the uncertainties in the changes of the regional water cycle and global sea-level rise (SLR) prevent a prediction of future salinity changes. Nevertheless, the results of regional climate modelling can be used in the planning of observational programmes to monitor the success of mitigation actions, not only for the Baltic Sea but also for other coastal seas with multiple stressors and complex oceanography. An additional finding of our study was that, despite climate change, the BSAP may soon lead to significantly improved bottom oxygen concentrations in the north-western Gotland Basin, whereas in shallow coastal zones any signal indicative of climate change or nutrient load reductions may be difficult to detect. Thus, assuming the scenarios of this study, the detection of systematic changes requires that, in the worst case, records of bottom oxygen concentration from coastal observatories must be longer than 70 years.

Results

Future projections

Using an ecosystem model for the Baltic Sea11,36 and a dynamical downscaling approach29 (Fig. 2), we regionalised scenario simulations of four Earth System Models (ESMs) driven by two greenhouse gas scenarios, medium (RCP 4.5) and high (RCP 8.5) (Table 1). Two nutrient load scenarios were applied: the BSAP (following SSP 1) and a business-as-usual or reference scenario (REF) that assumes the same average nutrient loads for the future (following SSP 2) as during a reference period (2010–2012)29,37. Nutrient inputs during the latter were modified by taking into account the projected increase in river discharge (Fig. 3). Finally, three global SLR scenarios, low, medium and high, complemented our ensemble of 48 projections. The natural or internal variability of the ecosystem was measured by the interannual fluctuations in all simulations on time scales up to 30 years. Approximately, such fluctuations arise in the absence of any radiative or biogeochemical forcing because the response of the ecosystem to these drivers would be on longer time scales. The different dimensions of the ensemble were designed to answer the following questions: (1) How does the ecosystem of the Baltic Sea react to changes in different external drivers? (2) How do natural variability and global climate model uncertainties affect scenario simulations of the Baltic Sea ecosystem?

Fig. 2: Dynamic downscaling approach for the Baltic Sea region.
figure 2

The various components of the modelling approach are explained in the Methods section.

Full size image
Table 1 List of 48 scenario simulations.
Full size table
Fig. 3: Projections of river discharge and nutrient loads from land and atmosphere into the entire Baltic Sea.
figure 3

Upper panel: The low-pass filtered runoff data (in m3 s−1), using a cutoff period of 30 years, of four regionalised ESMs (illustrated by different line types) under RCP 4.5 (green) and RCP 8.5 (red) scenarios. Lower panels: bioavailable phosphorus (in 106 kg P yr−1, left panels) and nitrogen loads (in 109 kg N yr−1, right panels) from land (upper panels) and atmosphere (lower panels) under RCP 4.5, BSAP (blue), RCP 4.5, REF (green), RCP 8.5, BSAP (orange, not visible) and RCP 8.5, REF (red) scenarios. Nutrient loads during the historical period are depicted in black. Past and future nutrient loads were compiled by Saraiva et al.29. Note the different scales of the y axes. Atmospheric phosphorus deposition is small compared with the phosphorus loads from land.

Full size image

River runoff and nutrient loads

In RCP 4.5 and RCP 8.5, river runoff is projected to increase 2–22% and 7–22%, respectively (Fig. 3). The ranges correspond to the ensemble spreads. The percentages were calculated from the differences between historical (1976–2005) and future (2069–2098) periods and reflect the enhanced hydrological cycle, especially in the northern Baltic Sea region during winter. Note that the river runoff shows a pronounced multidecadal variability that affects the temporal evolution of salinity in the Baltic Sea38.

Since the 1950s, both nutrient loads from the Baltic catchment area and atmospheric deposition have increased while the human population has grown such that untreated wastewater discharges have increased and agricultural practices have intensified8 (Fig. 3).

Water temperature

The ensemble spread in Fig. 4 indicates that volume-averaged water temperatures at the end of the century will be 1.3–2.2 °C (RCP 4.5) and 3.0–4.2 °C (RCP 8.5) higher than the ensemble mean water temperature for the historical period (Table 2). The ranges correspond to the ensemble spreads (maxima minus minima) shown in Fig. 4.

Fig. 4: Records of selected climate indices and ecosystem indicators and their uncertainties.
figure 4

Left panels: from top to bottom, annual mean volume-averaged temperature (in °C), volume-averaged salinity (in g kg−1), hypoxic area (in 104 km2), volume-integrated primary production (in 109 kg C yr−1) and volume-integrated nitrogen fixation (in 108 kg N yr−1) for the entire Baltic Sea, including the Kattegat (see Fig. 1) in historical (≤2005, black lines) and scenario simulations (>2005, coloured lines) driven by four regionalised ESMs (illustrated by different line types) under RCP 4.5, BSAP (blue), RCP 4.5, REF (green), RCP 8.5, BSAP (orange) and RCP 8.5, REF (red) scenarios. A spin-up simulation since 1850 was performed as illustrated by the evolution of hypoxia. Right panels: related uncertainties (or ensemble spread expressed in standard deviations of ensemble members compared with the ensemble mean): (1) natural variability with periods <30 years (black line) (Eq. 2), (2) climate model uncertainty estimated from four ESMs (blue line) (Eq. 3), (3) the unknown greenhouse gas concentration pathway estimated from RCP 4.5 and RCP 8.5 (red line) (Eq. 5), (4) the unknown nutrient load emission pathway estimated from BSAP and REF (green line) (Eq. 6) and (5) the unknown mean SLR estimated from low, medium and high scenarios (yellow line) (Eq. 4). The dotted curve denotes the square root of the sum of the five variances that express the various sources of uncertainty (Eq. 8, right-hand side). Owing to the filtering approach, natural variability can only be calculated until 2069.

Full size image
Table 2 Changes in sea surface temperature.
Full size table

The spatial pattern of SST changes between present and future climates was found to depend on the season (Supplementary Fig. 1). In the northern Baltic Sea, maximum warming is delayed by several months compared with the atmosphere because of the persisting sea-ice cover39. According to the simulations, maximum warming of the surface water will occur during summer (June to August) instead of during winter (December to February), as in the atmosphere20. From winter to autumn, maximum warming is projected to shift from the Gulf of Finland (for the location see Fig. 1), via the central Bothnian Sea and the central Bothnian Bay to the northernmost coastal zone of the Bothnian Bay because of the insulation effect of the seasonal sea-ice cover. The amount of warming in the northern Baltic Sea would be larger in summer than during all the other seasons.

Salinity

The projected ensemble mean changes in volume-averaged salinity were smaller than the spread of the ensemble members (Fig. 4). Contrary to the results of previous studies that projected a decreasing salinity because of the increasing river runoff caused by the intensified regional water cycle21, our ensemble does not indicate a change in salinity (Fig. 4). The main reason was that previous scenario simulations did not account for the SLR, which according to our results would more or less completely compensate for the increasing runoff. A higher mean sea level in the world’s oceans would lead to more intense saltwater influxes from the adjacent North Sea into the Baltic Sea, causing larger vertical salinity gradients and a stronger stratification within the Baltic Sea40. Note that the regionally generated SLR components, i.e., the thermal and haline expansions and wind-driven changes in the Baltic Sea, are small compared to the remotely generated SLR.

Oxygen concentration

The simulations showed that the fate of bottom oxygen concentrations in the Baltic Sea will differ between shallow and deep areas. In shallow waters without perennial vertical stratification, sea bottoms are well ventilated from the sea surface and the residence time is short. In deep areas, a perennial stratification caused by vertical salinity gradients (the halocline) hinders ventilation of the deep water, and predominantly lateral water transports supply oxygen. Hence, in areas below the halocline, changes in the bottom oxygen concentration mainly depend on the nutrient load scenario that controls oxygen consumption (Fig. 5). However, the RCP and SLR scenarios impact the fate of oxygen concentrations as well. For example, under the combined BSAP and low SLR scenario, oxygen conditions below the halocline are projected to improve everywhere according to either RCP scenario. However, under a combined BSAP, high SLR and RCP 8.5 scenario, oxygen concentrations below the halocline may even decrease in certain regions such as the Bothnian Sea. A higher sea level in the Kattegat will cause a stronger stratification and, hence, smaller vertical oxygen transports from the sea surface to the deep water, with the increased inflow of saline and oxygen-rich water from the North Sea ventilating only the deeper parts of the Baltic Sea, a mechanism well known from observations10,40.

Fig. 5: Changes in the summer mean bottom oxygen concentration.
figure 5

Changes in the summer (June to August) mean bottom oxygen concentration (in mL O2 L−1) between the periods 1976–2005 and 2069–2098 under RCP 4.5 (upper panels) and RCP 8.5 (lower panels) driven by BSAP (first and third columns) and REF (second and fourth columns), and low (first and second columns) and high (third and fourth columns) SLR scenarios.

Full size image

Under the REF scenario, oxygen conditions in the Baltic proper below the halocline are projected to improve as well, reflecting the delayed response to the historical decline in nutrient loads since the 1980s owing to the long response time of the Baltic Sea8. However, these changes will be much smaller than under the BSAP scenario (Fig. 5).

In shallower areas without a halocline, all projections indicated a relatively small (<1 mL O2 L−1) but systematic decrease in bottom oxygen concentrations (Fig. 5) because the oxygen solubility of warm water is less than that of cold water. The exceptions are limited to river-dominated coastal areas close to river mouths with strong stratification, where the oxygen concentration would increase owing to the historical reductions in nutrient inputs. Hence, in shallow areas of the Baltic Sea the degree of oxygen decline is mainly controlled by the RCP scenario, with a minor impact of the SLR scenario.

Depending on the nutrient load scenario, the size of the hypoxic area at the end of the century is projected to change between −78 and −22% (BSAP) or between −39 and +21% (REF) (Fig. 4). The RCP or SLR scenario had less impact on the size of hypoxic area than did the nutrient load scenario. However, the differences in the projected size of the hypoxic area according to the simulations driven by the four ESMs were considerable (Table 3). Results similar to those for the hypoxic area were obtained for primary production and nitrogen fixation by cyanobacteria, but with somewhat different response time scales for the BSAP scenario (Fig. 4). At the end of the century, nitrogen fixation, accompanied by harmful cyanobacterial blooms, will almost vanish under the BSAP scenario, regardless of the climate scenario, because surface nutrient concentrations, in particular phosphate, will be much smaller than during the historical period (not shown).

Table 3 Changes in hypoxic area.
Full size table

Water transparency

In addition to oxygen concentrations, an important ecosystem indicator is the SD, defined as the depth at which a black-and-white or white disk lowered into the water is no longer visible. For the design of the BSAP, only these two ecosystem indicators were used by the Helsinki Commission (HELCOM) and long records of oxygen concentration and SD observations were used to evaluate our ecosystem model during the historical period31. With the implementation of the BSAP, by the end of the century considerable increases in SD of up to 2 m, particularly in the coastal zone of the Baltic proper, are projected (Fig. 6). These changes are not significantly different between the various SLR or RCP scenarios. Under the REF scenario, SD changes are much smaller and slightly differ between RCP scenarios. On average, under the BSAP and REF scenarios the SD will be 0.6 m and 0.1–0.2 m deeper, respectively.

Fig. 6: Changes in the summer mean Secchi depth.
figure 6

Summer (June to August) mean Secchi depth changes (in m) between the periods 1976–2005 and 2069–2098 under RCP 4.5 (upper panels) and RCP 8.5 (lower panels) and under BSAP (left panels) and REF (right panels) scenarios. Ensemble mean results of all ESM and sea-level rise scenarios are depicted.

Full size image

Uncertainty

Figure 4 shows that the spread in the projections can be attributed to the sources of uncertainty, i.e. the natural, interannual variability (here with periods <30 years), model uncertainty, and the choices of the RCP, SSP and SLR. Uncertainties in temperature before ~2065 are dominated by natural variability and model uncertainty. After 2065, the largest source of uncertainty is related to RCP. Uncertainties caused by either the SLR or SSP are negligible. For salinity, the most important uncertainty during almost the entire 21st century is caused by the differences between the ESMs derived, inter alia, from the large spread in the river runoff projections (Fig. 3). However, the uncertainty in the SLR is also important and becomes dominant after 2070.

Uncertainties in the size of the hypoxic area before 2025 are dominated by the natural variability (Fig. 4), whereas after 2025 the largest source of uncertainty is related to the SSP. However, all other factors contribute to the uncertainty. Hence, the predictability of the size of the hypoxic area over long time scales requires accurate knowledge of ESM sensitivities and the various drivers of the projections. Similar results were found for primary production and nitrogen fixation (Fig. 4). For both, the model system predicts that, already in 2024 and 2029, respectively, the uncertainty due to SSP will be larger than the natural variability. Interestingly, the results showed a significant decrease in the natural variability of nitrogen fixation over time and that, at the end of the century, the uncertainty contribution due to an unknown SLR will be as large as the uncertainty due to the RCP.

The natural variations of SST, bottom oxygen concentration and SD are significantly larger than the corresponding model uncertainty but the opposite is true for bottom salinity (Supplementary Figs. 2–5). For salinity, model differences, above all with respect to precipitation, rather than the natural variability prevented the detection of signals.

For SST, the natural variability is most pronounced in upwelling regions along the coasts (Supplementary Fig. 2). Owing to the prevailing south-westerly winds, upwelling is a common feature along the Swedish east coast. The largest model uncertainty occurs in the northern Baltic Sea, very likely related to the ice breakup in seasonally ice-covered regions41.

For bottom salinity, the natural variability is largest along the pathway of inflowing saltwater and at the depth of the halocline at the slopes of the bathymetry, indicating the vertical motion of the halocline accompanied by saltwater inflows (Supplementary Fig. 3). The model uncertainty is largest below the halocline because the wind fields and river runoff of the regionalised ESMs differ significantly, resulting in different saltwater inflows and bottom salinities.

For bottom oxygen concentration, the spatial pattern of the natural variability is similar to the bottom salinity pattern because saltwater inflows below the halocline usually transport higher concentrations of dissolved oxygen than those of the surrounding water (Supplementary Fig. 4). In addition, vertical salinity gradients in the Baltic Sea change over time, affecting vertical oxygen fluxes along the slopes of the Gotland Basin and in the Gulf of Finland. The model uncertainty shows a similar pattern but its magnitude is much smaller than the natural variability. In the coastal zone, both the natural variability and the model uncertainty are smaller than in the deep water.

For the SD, the natural variability is largest in the Baltic proper, in particular, in the coastal zone but horizontal gradients between the coastal zone and the open sea are smaller than those of the other variables discussed herein (Supplementary Fig. 5). As in the case of the SST or bottom oxygen concentration, model uncertainty is smaller than the natural variability.

In summary, patterns of natural variability and model uncertainty can be mainly explained by regional processes such as coastal upwelling events, saltwater inflows and the vertical motion of the halocline. With the exception of salinity, the model uncertainty was found to be smaller than the natural variability.

Time of emergence

The time at which the signal of climate or environmental changes can be distinguished from the noise of the natural variability is referred to as the time of emergence (TE), and it is a key variable in climate projections and risk assessments42. Figure 7 shows TE relative to the reference period 1991–2020 for the winter SST, the winter bottom oxygen concentration and the summer SD based on 12 regionalised climate simulations (four ESMs combined with three SLR scenarios) driven by RCP 8.5 and the BSAP. Areas with the shortest TE will be the first in which climate and environmental changes relative to the reference period 1991–2020 will be identifiable. During the historical period (1976–2005), nutrient loads decreased considerably (Fig. 3), but as our aim was to calculate TE with potential relevance for stakeholders, we shifted the reference period from the historical period to 1991–2020 to allow a comparison with the currently low nutrient loads vs the future. If the historical period is applied as a reference, then past changes in nutrient loads primarily dominate the response of the ecosystem model, because the impact of a changing climate will be smaller than that of changing nutrient loads.

Fig. 7: Time of emergence of sea surface temperature, bottom oxygen concentration and Secchi depth changes.
figure 7

Time of emergence (after 2000) of the winter (December to February) sea surface temperature (upper row), winter bottom oxygen concentrations (middle row), and summer (June to August) Secchi depth (lower row) changes under the RCP 8.5 and BSAP scenario relative to the reference period 1991–2020. From left to right, first quartiles, medians and third quartiles are depicted. Only areas with robust signals are shown, i.e., >50% of the models have a time of emergence within the period 2006–2068. The seasons with the earliest TE in large areas are presented. The Secchi depth is derived from simulated underwater light fields.

Full size image

For the three variables shown in Fig. 7, the seasons were selected according to the TE. For the SST and bottom oxygen concentration, areas with statistically significant changes are larger during winter than during summer whereas for the SD the opposite is true.

The first signals indicative of SST changes are projected to occur in winter sometime during 2040–2050, in the centres of the northern Baltic Sea and the Gotland Basin. The range was estimated from the first and third TE quartiles of the ensemble members. In coastal upwelling areas in the central Baltic Sea, SST changes will occur later, during 2050–2070.

For bottom oxygen concentration changes, the shortest TE was determined in the north-western Gotland Basin, occurring roughly during the period 2020–2030 and in all seasons (in Fig. 7 winter values are shown). Hence, in this area, the impact of the BSAP and previous efforts of nutrient load abatements will be detected relatively soon. In the eastern Gotland Basin, however, the bottom oxygen concentration below the halocline is also controlled by saltwater inflows from the North Sea, thus extending the TE to sometime between 2030 and 2050. In shallow areas without halocline, the pronounced natural variability is much larger than any signal, preventing the detection of changes in bottom oxygen concentrations until at least 2070. These areas along the southern and eastern coasts of the Baltic proper are characterised by water depths shallower than ~40 m. Note that the horizontal resolution of the model was too coarse to properly resolve river mouth estuaries and many of the fjords along the Swedish coast. Some of these estuaries are vertically stratified, resulting in prolonged residence times. Hence, for these coastal areas our model might fail to calculate realistic TE values.

In the Gotland Basin, the TE of the summer SD changes is between 2030 and 2060, with a pronounced east–west gradient (Fig. 7). Thus, the increased SD occurring in response to the implementation of the BSAP will be detected first in the western and then in the eastern Gotland Basin. In the Bothnian Sea, a robust TE cannot be determined before 2070.

Discussion

Underestimated uncertainties

Our method probably underestimates the uncertainties inherent in the different sources because the ensemble used is still too small. In particular, model uncertainty is likely to be underestimated. To test this hypothesis, further ESMs could be regionalised, following the dynamical downscaling method (Fig. 2), and added to the ensemble. We applied a cluster analysis of eight regionalised ESMs, including the four of this study, to investigate model uncertainties. For this analysis, selected climate indices, i.e., temperature, precipitation, total cloud cover, two near-surface wind components, ocean mixed layer depth, surface and bottom water temperature and surface and bottom salinity, were used. The results suggested that the selected sub-group of ESMs represents the overall uncertainties of atmospheric forcing relatively well (Supplementary Fig. 6).

However, model uncertainties are caused not only by global but also by regional models. Marine ecosystem models, in particular, are subject to large uncertainties because, unlike the physical models, they are not completely based on first principles such as the conservation of energy and momentum. The sources and sinks of the prognostic variables of biogeochemical cycling are highly parameterised because the underlying processes are not well understood. Although scenario simulations with various ecosystem models for the Baltic Sea have been assessed before30,32, coordinated multi-model experiments are currently not available. The disadvantage of uncoordinated experiments, however, is that the various sources of uncertainty cannot be distinguished. Nonetheless, the method presented in this study can serve as a prototype to estimate uncertainties and identify their sources in high-resolution projections for coastal seas, as discussed below.

Long response time

Another difficulty in analysing the various drivers of the marine ecosystem is associated with the long response time of coastal seas such as the Baltic Sea, which in the latter case is on the scale of several decades43,44. Given the delayed response in bottom oxygen concentrations to the nutrient load reductions implemented since the 1980s, the additional effect of the BSAP, represented in the model by a linear nutrient load reduction between 2012 and 2020 (Fig. 3), can approximately be addressed by comparing the projections for the BSAP and REF scenarios. For this comparison, the slight increase in nutrient loads in the REF scenario owing to the increase in river runoff (Fig. 3) can be neglected. According to the results of our model, the improved bottom oxygen conditions below the halocline are at least in part owing to historical changes in the nutrient supply before 2012 (Fig. 5). In REF, a robust TE of bottom oxygen concentration changes was not found for the eastern Gotland Basin, but in the north-western Gotland Basin the mean signal was larger than the natural variability (not shown). Hence, the early emergence of improved oxygen conditions in the north-western Gotland Basin reflects not only the BSAP but also the historical nutrient load reductions between 1980 and 2012. However, an improvement in the oxygen conditions of the entire basin requires full implementation of the BSAP.

CMIP6 vs CMIP5

Additional global scenarios consider the new greenhouse gas pathways developed under CMIP645,46. They include a larger number of greenhouse gas scenarios in order to explore vigorous mitigation actions aimed at keeping global warming below 1.5 or 2.0°C but also overshoot scenarios of high emissions during the first half and negative emissions during the second half of this century45. Accordingly, slightly altered total responses can be expected. It should be noted, however, that the methodology developed in this study is also fully applicable to CMIP6 (and any other) scenarios, which will be investigated once regionalisations of CMIP6 simulations are available.

Apart from their inclusion of new scenarios, in the CMIP6 models, the equilibrium climate sensitivity to a doubling of greenhouse gas concentrations is higher than in the predecessor CMIP5 models (1.5–4.5 °C for CMIP5 vs 1.8–5.6 °C for CMIP6 models47), at least as determined in initial assessments. However, the CMIP6 models in which warming is highest at the end of the century do not well simulate the recorded warming trend during the historical period and should therefore be excluded from policy-relevant assessments45,48,49,50,51.

Finally, it should be noted that in CMIP5 the historical period (the period for which observation-based greenhouse gas concentrations are available) ends in 2005, whereas in CMIP6 it ends in 2014. However, since the differences in greenhouse gas concentrations between CMIP5 scenarios and historical simulations are very similar during the first two decades after the scenario simulation starts, a large impact on our results was unlikely. In addition, internal variability is the largest uncertainty source during the early stage of the scenario periods while differences between radiative forcing are quite low.

Natural variability

Changes in all variables studied were affected by a pronounced natural variability. Hence, projected changes such as increasing the water temperature of the Baltic Sea cannot be attributed to global warming alone. This finding has large consequences for the detection of future climate change signals with the aid of environmental monitoring programmes (see below) and for future projections. The ensemble sizes of scenario simulations in most climate change studies of coastal seas in general and of the Baltic Sea in special were too small30,32. Further research on natural variability is needed to disentangle climate change and natural variability, especially on decadal to multidecadal time scales.

Novel scenario simulations

Despite underestimated uncertainties and the usage of older forcing scenarios, the presented scenario simulations for the Baltic Sea are unique. The results differed from those of earlier studies based on CMIP321 because in this study (1) nutrient load scenarios were revised37,52, (2) initial conditions were improved using a long spin-up simulation with a start in 185053,54 and (3) various SLR scenarios were considered, such that vertical stratification was higher. In particular, the revised, lower nutrient loads of this study led to a significant improvement in the ecological status of the Baltic Sea under the BSAP scenario. Hence, our ensemble of scenario simulations suggested that under the combined effects of projected temperature increase, SLR and nutrient load reduction the long-term response of the oxygen budget below the halocline would be dominated by the anticipated nutrient scenario. Owing to the relatively large number of projections, a better understanding of the signal-to-noise ratio was gained, in which the latter is controlled by considerable climate model differences and natural variability limiting the predictive capacity of any future projections.

Implications for management

Our results have major implications for the management of coastal seas that currently suffer from eutrophication, oxygen depletion and other environmental threats such as overfishing, contaminants, marine litter, marine traffic, invasive species, offshore activities and habitat degradation. In this context, the Baltic Sea well represents coastal seas that are threatened by multiple stressors. Reusch et al.55 argued that “the Baltic Sea can serve as a time machine to study consequences and mitigation of future coastal perturbations”.

Coastal seas are particularly vulnerable to a changing climate because of the potential impacts related to SLR, temperature increases and deoxygenation56. However, there have been few investigations of eutrophied coastal seas under a changing climate; examples are studies in Chesapeake Bay57,58,59, the northern Gulf of Mexico60 and the Baltic Sea21. The impacts of climate change on coastal seas will be manifold. For instance, SLR may cause flooding, sediment transport and erosion and even higher salinities40,61. A recent study showed that the potential contributions of melting of ice sheets to the SLR are larger than previously estimated, and a global total SLR exceeding 2 m by 2100 lies within the 95% uncertainty ranges for a high emission scenario62. We also considered these high SLR scenario simulations and found a considerable impact of such an SLR on the size of the hypoxic area of the Baltic Sea.

Further, our results suggest that mitigation actions such as the BSAP will be successful also in future climates, regardless of the climate scenario. The nutrient load reductions will result in an improved ecological status in which cyanobacterial blooms have vanished, the hypoxic area has become smaller and water transparency has improved. The improved water transparency will allow light to penetrate the thermocline more frequently and may eventually cause subsurface phytoplankton blooms63. However, as both warming and SLR will counteract the efforts of the BSAP, larger nutrient load reductions might be needed to reach the same environmental targets defined without considering a changing climate. In principle, our approach will allow the calculation of revised maximum allowable loads, including their uncertainty ranges, under future climate scenarios.

Our modelling approach can serve as a prototype for similar studies of the Baltic Sea and other coastal seas because it has the following advantages: (1) Calculations of the impact of nutrient load reductions on marine biogeochemistry consider climate change. This may be particularly important for assessments of coastal seas in which the response times of marine ecosystems to changes in external nutrient loads are long. (2) The dynamical downscaling approach is performed using a regional climate model with high resolution. Compared with global models, a coupled atmosphere–ice–ocean model better resolves regional orography, the regional land-sea mask, processes such as deep cyclones and regional feedback mechanisms such as the ice-albedo and surface wind-mixed layer feedbacks64,65,66. (3) Future nutrient loads are calculated using a terrestrial biogeochemical model, considering the impact of climate change on the land surface. (4) An ensemble of scenario simulations enables the estimation of uncertainties in the projections and is able to reveal the sometimes limited predictive capacity of the model system attributable to model differences and natural variability.

Implications for marine monitoring

Further, our TE results have implications for the design of cost-intensive, long-term environmental monitoring programmes. We suggest focusing on regular measurements of bottom oxygen concentration below the halocline at locations less affected by saltwater inflows, such as the north-western Gotland Basin. In general, regular observations in the coastal zone for the detection of bottom oxygen concentration changes might only be of value if the records are homogeneous and long enough compared to the pronounced natural variability. Currently, such observational records do not exist for the Baltic Sea. Figure 1 shows the locations of international monitoring stations with regular sampling, as listed by HELCOM. Many of them are located in the shallow coastal zone. Our method may help to improve the current observation system and serve as a prototype in the design of cost-effective monitoring programmes also in other coastal seas prone to eutrophication.

Methods

Modelling approach

The regional climate model RCA4-NEMO67,68 was used to regionalise four ESMs: MPI-ESM-LR, EC-Earth, IPSL-CM5A-MR, HadGEM2-ES (Gröger et al.; with references for the ESMs therein69). The selection of the ESMs from available models was motivated with the help of a cluster analysis (see Supplementary Methods, Supplementary Fig. 6). The scenario simulations of the ESMs are part of the latest Intergovernmental Panel on Climate Change (IPCC) assessment report70. The atmospheric surface fields resulting from dynamical downscaling were used to force the marine ecosystem model RCO-SCOBI for the Baltic Sea, including the Kattegat11,36 (Fig. 2). With its aid, both ocean physics44 and marine nutrient cycles for nitrogen and phosphorus, including concentrations of phyto- and zooplankton as well as dead particulate organic matter (detritus), were calculated31. This so-called NPZD type model considers nutrient cycling both in the water column and in the sediments71,72.

In the RCO-SCOBI model11, SD as a measure of water transparency is calculated as SD = 1.7/k(PAR), where k(PAR) is the coefficient of the underwater attenuation of the photosynthetically available radiation73,74. In the model, the concentrations of phytoplankton and detritus control the attenuation of light in the water column. Consequently, changes in the SD are expressed by changes in the phytoplankton and detritus concentrations, whereas the concentration of coloured dissolved organic matter is assumed to be constant over both time and sub-basin areas, because more-detailed information, for instance from observations, is not available. This simplified approach follows that of earlier studies of the Baltic Sea75.

Surface air temperature and precipitation from the same four regionalised ESMs were also used to force the hydrological model E-HYPE (Hydrological Predictions for the Environment, http://hypeweb.smhi.se), a process-based, high-resolution multi-basin model established for Europe (Fig. 2). In this study, E-HYPE was used to calculate river runoff and nutrient loads from the entire Baltic Sea catchment area76,77.

The applied dynamical downscaling strategy with several one- and two-way coupled models was chosen because the regional climate model RCA4-NEMO does not contain a marine or terrestrial biogeochemical component (Fig. 2). The horizontal and vertical grid resolutions of the coupled physical-biogeochemical ocean model are 3.7 km and 3 m, respectively. Compared with other Baltic Sea models, the skill of RCO-SCOBI relative to observations during historical periods driven by either downscaled atmospheric reanalysis data78,79 or regionalised climate forcing30 is good. For instance, Meier et al.30 showed that during the historical period the group of climate simulations included in our study, using RCO-SCOBI, belonged to the three best out of eight groups of Baltic Sea models. The cost function used as the criterion for the model intercomparison was calculated from the mean biases divided by the standard deviations of the observations at the most frequently monitored stations. These locations are also representative of the different sub-basins of the Baltic Sea. For the model intercomparison, temperature, salinity and dissolved oxygen, inorganic nitrogen and inorganic phosphorus concentrations during 1980–2005 were assessed30. Furthermore, the important saltwater inflow events and vertical stratification are well represented in RCO-SCOBI25,36,43,80.

Scenarios

In this study, 48 scenario simulations for the Baltic Sea were carried out based upon a matrix of (1) four driving ESMs68,69; (2) two greenhouse gas concentration scenarios following RCP 4.5 (medium) and RCP 8.5 (high)33; (3) two nutrient load scenarios following regionalised SSP 1 (a low case scenario, called BSAP26) and SSP 2 (a business-as-usual or REF scenario that takes into account the impact of climate change and the current socioeconomic conditions29,37,52) and (4) three linearly increasing SLR scenarios from the year 2005. The projected mean sea-level changes in the year 2100 are: (scenario 1) 0 m, (scenario 2) the ensemble mean of RCP 4.5 (0.54 m) and RCP 8.5 (0.90 m) IPCC projections56,81 and (scenario 3) the 95th percentiles of the low- (1.26 m) and high-case (2.34 m) scenarios following Bamber et al.62 (Table 1). Hence, we added another dimension to the ensemble to consider the uncertainties in SLR on global scales. Previous sensitivity studies carried out with a step-wise increase in the mean sea level indicated the importance of the SLR for the Baltic Sea ecosystem40. Here, we applied a linearly increasing mean sea level by deepening the water depth at all grid points every 10 years.

Under the low scenario (scenario 1), SLR is neglected, thus following previous scenario simulations. For water exchange between the Baltic Sea and the North Sea, the sea level relative to the seabed of the sills in the entrance area limits transport40. The sills located in the two connecting shallow and narrow straits (Great Belt and Öresund) are only 8 and 18 m deep (Fig. 1). Therefore, during wind-driven inflow events, the water column is mixed vertically and the amount of inflowing saltwater is proportional to the cross-sections of the straits at the sills approximately40. Based on the eustatic SLR and the glacial isostatic adaptation, the relative SLR during the period 1915–2014 was estimated to be 0–1 mm year−1 82. Hence, a very optimistic SLR scenario would be that the future SLR simply follows the lower estimate of the past SLR. The mean sea-level conditions (scenario 2) were compiled according to the IPCC Special Report on the Ocean and Cryosphere in a Changing Climate56 without considering the glacial isostatic adaptation and the gravitationally induced adjustments for the Baltic Sea region that were made due to preferential melting of the Antarctic ice sheet81,83. For the high scenario (scenario 3), experiments with an SLR according to the 95th percentile, as suggested by Bamber et al.62, were carried out. Low- and high-SLR scenarios were used to approximately mirror the entire range of potentially possible sea-level scenarios adapted to the Baltic Sea. SLR scenarios are not independent of RCP scenarios because larger greenhouse gas concentrations cause accelerated SLRs (Table 1). However, we treated SLR as an additional dimension of the uncertainty space because the interactions between the atmosphere-ocean system and the ice sheets of Greenland and Antarctica are not properly considered by the ESMs used in CMIP562.

As mentioned above, the river runoff and nutrient load data for nitrogen and phosphorus were calculated with the hydrological model E-HYPE77 and driven by air temperature and precipitation fields regionalised from the ESMs. However, the hydrological forcing data of IPSL-CM5A-MR under the RCP 8.5 scenario were not available. Instead, we used the forcing data of MPI-ESM-LR under RCP 8.5 because a cluster analysis following Wilcke and Bärring84 suggested that under RCP 8.5 air temperature and precipitation changes for the Baltic Sea region in MPI-ESM-LR and IPSL-CM5A-MR are similar.

Nutrient load scenarios defined plausible pathways of nutrient inputs from rivers, point sources and atmospheric deposition following SSP 1 and SSP 229,37,52. In the BSAP scenario, nutrient loads decreased linearly from the values in 2012 (average for 2010–201285) to the maximum allowable input defined in the mitigation plan in 2020. Thereafter, the nutrient loads remained constant until 2098 (Fig. 3).

In our opinion, the low and high RCP, SSP and SLR scenarios represent optimistic and pessimistic scenarios, respectively, and follow plausible future pathways. However, RCP 2.6 scenario simulations are missing and the ranges of the SLR scenarios may be too large. Nevertheless, the chosen ranges were not likely to impact our results considerably.

Model evaluation

The standard deviations of SST, bottom salinity, bottom oxygen concentration and SD from four regionalised ESMs (MPI-ESM-LR, EC-Earth, IPSL-CM5A-MR, HadGEM2-ES) for the period 1977–2005 were compared with the standard deviations of ocean reanalysis data available for the period 1971–199986 (Supplementary Figs. 2–5). Despite minor regional differences between the various models and the observations caused by differences in the atmospheric and hydrological forcings, the large-scale patterns and magnitudes of the standard deviations in the historical climate simulations and reanalysis data were similar. Thus, we argue that the quality of the simulated variability in our scenario simulations was sufficient for the purpose of this study. The natural variability \(\sqrt {\bar V(r,s)}\) was smaller than the total standard deviation because periods >30 years were excluded.

Uncertainty analysis

The different dimensions of the ensemble studied here (ESMs, SLRs, RCPs and SSPs; see Table 1) allowed us to examine how the marine ecosystem responds to different external drivers and how uncertainties in global climate projections affect projections of the Baltic Sea ecosystem. The simulated ensemble spread of the projected changes was used to estimate the uncertainties87. Each scenario simulation x(m, l, r, s, t) as a function of time t and driven by ESM m (m = 1, …, Nm), SLR scenario l (l = 1, …, Nl), RCP scenario r (r = 1, …, Nr) and SSP scenario s (s = 1, …, Ns) with Nm = 4, Nl = 3, Nr = 2 and Ns = 2 was separated into an average for the period 1991–2020 \(\bar x_0\left( {m,l,r,s} \right) = \bar x\left( {m,l,r,s,t_0} \right)\) with the centred time t0 = 2005, a low-pass filtered component with a cutoff period of 30 years \(\bar x\) and a residual ε:

$$x\left( {m,l,r,s,t} \right) = \bar x_0\left( {m,l,r,s} \right) + \bar x\left( {m,l,r,s,t} \right) + \varepsilon \left( {m,l,r,s,t} \right)$$
(1)

Hence, the differences in initial conditions owing to model differences and natural variability during the historical period were removed and not considered in the analysis. Because no more than one climate realisation was available for each of the four downscaled ESMs, we estimated the uncertainty of natural variations V(t) from the high-frequency part of each scenario simulation and calculated the ensemble mean of the variances \(\sigma _{t|t^\prime }^2\) of 30-year running time windows t’ є {1, …, 30} at time t:

$$V\left( t \right) = \frac{1}{N}\mathop {\sum }\limits_{m,l,r,s} \sigma _{t|t^\prime }^2\left( {\varepsilon \left( {m,l,r,s,t^\prime } \right)} \right)$$
(2)

with N = Nm Nl Nr Ns = 48. This approach assumes that the low-pass filtered time series mainly contains the signal of anthropogenic climate or nutrient load changes of each individual ESM and neglects the low-frequency natural variability on time scales larger than 30 years. The model uncertainty M(t) was estimated from the ensemble mean variances in the low-pass filtered projections driven by the different ESMs averaged for all RCPs, SSPs, SLRs:

$$M\left( t \right) = \frac{1}{{N_lN_rN_s}}\mathop {\sum }\limits_{l,r,s} \sigma _m^2\left( {\bar x\left( {m,l,r,s,t} \right)} \right)$$
(3)

The scenario uncertainties related to the future SLR (L(t)), RCP (R(t)) and SSP (S(t)) are the corresponding ensemble mean variances of the multi-model means x(l,r,s,t):

$$L\left( t \right) = \frac{1}{{N_rN_s}}\mathop {\sum }\limits_{r,s} \sigma _l^2\left( {{\boldsymbol{x}}\left( {l,r,s,t} \right)} \right)$$
(4)
$$R\left( t \right) = \frac{1}{{N_lN_s}}\mathop {\sum }\limits_{l,s} \sigma _r^2\left( {{\boldsymbol{x}}\left( {l,r,s,t} \right)} \right)$$
(5)
$$S\left( t \right) = \frac{1}{{N_lN_r}}\mathop {\sum }\limits_{l,r} \sigma _s^2\left( {{\boldsymbol{x}}\left( {l,r,s,t} \right)} \right)$$
(6)

with the multi-model mean

$${\boldsymbol{x}}(l,r,s,t) = \frac{1}{{N_m}}\mathop {\sum }\limits_m \bar x\left( {m,l,r,s,t} \right)$$
(7)

The low-pass filtered total variance \(\bar U(t)\) of the anomalies of all scenario simulations is approximately equal to the low-pass filtered sum of the five uncertainties V(t), M(t), L(t), R(t) and S(t):

$$\bar U\left( t \right) = \overline {\sigma _{m,l,r,s}^2(x\left( {m,l,r,s,t} \right) - \bar x_0\left( {m,l,r,s} \right))} \cong \overline {V\left( t \right) + M\left( t \right) + L\left( t \right) + R\left( t \right) + S\left( t \right)}$$
(8)

From the requirement that the signal-to-noise ratio be equal or larger than two, the time of emergence TE(m,l,r,s) of each scenario simulation was calculated:42

$$\frac{{Signal}}{{Noise}} = \frac{{\bar x\left( {m,\,l,\,r,\,s,T_E} \right) - \bar x_0\left( {m,l,r,s} \right)}}{{\sqrt {\frac{1}{2}\left( {\sigma _{T_E|t^\prime }^2(\varepsilon \left( {m,l,r,s,t^\prime } \right)) + \sigma _{t_0|t^\prime }^2(\varepsilon \left( {m,l,r,s,t^\prime } \right))} \right)} }} = 2$$
(9)

Figure 7 shows the median, first and third quartiles of the TE(m,l,r,s) of 12 projections driven by four ESMs combined with three SLR scenarios. By definition, robust signals required that >50% of the models had a time of emergence within the period 2006–2068. Our method expands previous approaches42,87,88 by adding additional dimensions of SSPs and SLRs to the ensemble of regional climate simulations, taking multiple stressors into account, and by separating noise from anthropogenic signals by filtering. Hence, in this study, the natural variability was defined as variations with periods smaller than 30 years. However, the approach does not consider the natural variability on longer time scales such as variations caused by the Atlantic Multidecadal Oscillation39,89. For further information, the reader is referred to the Supplementary Methods.