Coupled regional Earth system modeling in the Baltic Sea region

Nonlinear responses to externally forced climate change are known to dampen or amplify the local climate impact due to complex cross-compartmental feedback loops in the Earth system. These feedbacks are less well represented in the traditional stand-alone atmosphere and ocean models on which many of today's regional climate assessments rely (e.g., EURO-CORDEX, NOSCCA and BACC II). This has promoted the development of regional climate models for the Baltic Sea region by coupling different compartments of the Earth system into more comprehensive models. Coupled models more realistically represent feedback loops than the information imposed on the region by prescribed boundary conditions and, thus, permit more degrees of freedom. In the past, several coupled model systems have been developed for Europe and the Baltic Sea region. This article reviews recent progress on model systems that allow two-way communication between atmosphere and ocean models; models for the land surface, including the terrestrial biosphere; and wave models at the air–sea interface and hydrology models for water cycle closure. However, several processes that have mostly been realized by one-way coupling to date, such as marine biogeochemistry, nutrient cycling and atmospheric chemistry (e.g., aerosols), are not considered here.

In contrast to uncoupled stand-alone models, coupled Earth system models can modify mean near-surface air temperatures locally by up to several degrees compared with their stand-alone atmospheric counterparts using prescribed surface boundary conditions. The representation of small-scale oceanic processes, such as vertical mixing and sea-ice dynamics, appears essential to accurately resolve the air–sea heat exchange over the Baltic Sea, and these parameters can only be provided by online coupled high-resolution ocean models. In addition, the coupling of wave models at the ocean–atmosphere interface allows for a more explicit formulation of small-scale to microphysical processes with local feedbacks to water temperature and large-scale processes such as oceanic upwelling. Over land, important climate feedbacks arise from dynamical terrestrial vegetation changes as well as the implementation of land-use scenarios and afforestation/deforestation that further alter surface albedo, roughness length and evapotranspiration. Furthermore, a good representation of surface temperatures and roughness length over open sea and land areas is critical for the representation of climatic extremes such as heavy precipitation, storms, or tropical nights (defined as nights where the daily minimum temperature does not fall below 20 ^∘C), and these parameters appear to be sensitive to coupling.

For the present-day climate, many coupled atmosphere–ocean and atmosphere–land surface models have demonstrated the added value of single climate variables, in particular when low-quality boundary data were used in the respective stand-alone model. This makes coupled models a prospective tool for downscaling climate change scenarios from global climate models because these models often have large biases on the regional scale. However, the coupling of hydrology models to close the water cycle remains problematic, as the accuracy of precipitation provided by atmosphere models is, in most cases, insufficient to realistically simulate the runoff to the Baltic Sea without bias adjustments.

Many regional stand-alone ocean and atmosphere models are tuned to suitably represent present-day climatologies rather than to accurately simulate climate change. Therefore, more research is required into how the regional climate sensitivity (e.g., the models' response to a given change in global mean temperature) is affected by coupling and how the spread is altered in multi-model and multi-scenario ensembles of coupled models compared with uncoupled ones.

1 Introduction

Climatic and environmental changes on regional scales are traditionally investigated using stand-alone models that resolve processes specific to only one single environmental compartment (e.g., the terrestrial and marine biospheres, the hydrosphere, the atmosphere, the ocean and the cryosphere). Many projections that served as a basis for the recent climate change assessments for the North Sea (NOSCCA; May et al., 2016; Schrum et al., 2016) and Baltic Sea (Bøssing-Christensen et al., 2015, 2021) or the EURO-CORDEX region (Jacob et al., 2014) fall into the category of stand-alone ocean models or stand-alone atmosphere models, whereas fewer assessments have been based on coupled systems (Meier et al., 2015, 2021). However, as environmental compartments interact with each other via mass, momentum and energy exchange, the interfaces and boundary layers between the compartments are of great importance.

The use of Earth system models (ESMs) that interactively couple atmospheric, marine and terrestrial energy, water and biogeochemical dynamics is becoming increasingly common practice in global climate assessments (e.g., IPCC 2013) and international coordinated protocols for climate simulations (e.g., the emission-driven CMIP6-C4MIP; Eyring et al., 2016; Jones et al., 2016). Their employment in regional climate change and impact studies (e.g., Döscher and Meier, 2004; Somot et al., 2008; Wramneby et al., 2010; Meier et al., 2011; Cabos et al., 2020; Sein et al., 2020, Zhang et al., 2018; Soto-Navarro et al., 2020; Dieterich et al., 2019a; Gröger et al., 2019, 2021) is, however, still rare.

Climate change on regional scales can be much stronger than one would expect from external forcing such as greenhouse gases or solar radiation alone (e.g., Vogel et al., 2017; Stuecker et al., 2018). This is because the direct forcing can be strongly modulated by feedback processes that act to amplify or dampen the change in climate variables. Biophysical (e.g., albedo- and evapotranspiration-mediated) feedbacks may significantly affect the interactions between the Earth system compartments at the regional scale. Moreover, these feedbacks are affected by small-scale physiographic features, such as mountain ranges and coastline features, which are poorly captured by global ESMs when they are run at typically coarse grid resolution. Recognizing these issues, there is an emerging demand for regional Earth system models (RESMs) that are suitable for application at higher grid resolutions over a regional domain and that account for biophysical coupling between the atmosphere and surface. National (e.g., MERGE; http://www.merge.lu.se/about-merge (10 September 2021) and international strategic program are now stepping up their efforts to meet this demand. The Baltic Earth program (https://www.baltic-earth.eu/, 10 September 2021), which brings together climate and environmental scientists around the Baltic Sea, emphasizes the importance of realistic coupled modeling to achieve an improved Earth system understanding of the Baltic Sea region, with “Multiple drivers for regional Earth system changes” being one of the “Grand Challenges” addressed by the program (Reckermann et al., 2021).

1.1 Study area

Compared with other continents, Europe and the adjacent European marine sectors are extremely diverse on especially small scales. Classifications based on precipitation, temperature distributions (Köppen, 1923) and other environmental factors (Metzger et al., 2005) distinguish between 5 and 18 different climate types ranging from polar climates in high mountainous areas and Iceland to temperate regions of humid or even dry nature (Fig. 1; e.g., Kottek et al., 2006; Beck et al., 2018). By contrast, vast areas of comparably uniform environmental conditions that can be found on the world's bigger cratonic continents are absent. Climatologically, the western part of Europe is influenced by the oceanic climate (i.e., linked to the large heat content of the North Atlantic), whereas the eastern part is of continental character with increased seasonality. Likewise, Europe is positioned between the fully polar climate in the north and subtropical climates to the south. On the whole, this makes Europe's climate fairly variable on small spatial scales and sensitive to perturbations in the large-scale atmospheric circulation patterns as reflected in, for example, the strong impact of the North Atlantic Oscillation (NAO; Hurrel, 1995; Scaife et al., 2008; Rousi et al., 2020). Thus, the demand for modeling Europe's climate includes high resolution as well as a comprehensive process description by the respective coupled model components. This makes this region a challenging test case for high-resolution Earth system modeling.

Figure 1 Climate classification based on E-OBS monthly mean temperature and precipitation (Cornes et al., 2018). Classes are defined after Köppen (1923).

The Baltic Sea climate is influenced by a temperate, humid climate in the southwest and a snow climate in the north and east (Fig. 1), with an enhanced seasonal cycle giving rise to highly variable meteorological conditions related to predominant weather regimes over the region (Hertig and Jacobeit, 2014) ranging from severe storms, summer heat waves and winter cold-air outbreaks (Smith and Sheridan, 2020) to prolonged dry periods in the southern part. Many of these phenomena are directly subject to local and regional thermal feedbacks between the atmosphere, the land and the ocean and, thus, require a realistic exchange of mass and energy as realized by interactively coupled regional Earth system models.

1.2 Towards Earth system modeling of the Baltic Sea region

From a theoretical point of view, the coupling of two or more interacting models to create a more comprehensive system implies that boundary processes formerly prescribed, parameterized or even neglected are explicitly simulated. This removes observational constraints and/or empirically derived relationships on the model solutions and, thus, increases the model's degrees of freedom. Consequently, coupled models can drift from observed conditions; this makes their tuning more difficult but allows for a more realistic interaction between models. Therefore, in their hindcast modes, stand-alone models can be expected to be closer to observations once prescribed boundary conditions are of good quality.

From the climate perspective, which envisages simulations over several decades or even centuries, the model should ideally be drift-free to integrate over several million model time steps. Furthermore, with respect to future climate simulations, the model boundary conditions, such as temperature, are basically unknown and have to be derived from the output of available global climate models in the stand-alone case. Although advanced methods exist to make this data usable for high-resolution models (e.g., Hay et al., 2000; Chen et al., 2013), there is evidence that the solution of the ocean models is too tightly controlled by the global model in uncoupled ocean simulations (Mathis et al., 2018).

Hence, Earth system models are the prime tool for simulating cross-compartment feedback loops (Claussen, 2001; Giorgi and Gao, 2018; Heinze et al., 2019). In turn, such feedbacks play an important role in mediating the response of the Earth climate to a given external forcing or perturbation. Consequently, the ability of ESMs to simulate such feedbacks is essential. However, the capability of coupled models to better represent cross-compartment feedbacks and more realistically model dynamical processes is often not adequately accounted for during model validation. By contrast, model validation usually aims to demonstrate the models' ability to represent mean climate. By their nature, climate models are tools to iteratively solve the change in a variable from one given time step to the next rather than to predict the variable at a given time. From this, it follows that the capability of a model to reproduce transient behavior, such as interannual variability or long-term trends, is essential to estimate if a climate model can yield reliable answers to how changes in climate forcing will likely impact on climate variables. However, regional models are often validated more by how well they reproduce a given climatology of the present-day climate rather than by how well they reproduce trends derived from the historical past (Kerr, 2013). As every compartment model has its own spectrum of internal spatial scales and timescales, the inertia of the system increases when including slower components. This can become important, especially when a decision about the size of the ocean model domain has to be made. A larger extension to the open Atlantic or Arctic Ocean substantially increases the memory of the system, which has consequences for the model spin-up and the economical operation of the model system.

Over the past few decades, several advancements in coupled modeling have led to a growing number of different regional climate models on the way to a fully comprehensive description of the Earth system for many regions of the world, such as the east coast of the US (COAWST – Coupled Ocean–Atmosphere–Wave–Sediment Transport model; Warner et al., 2008, 2020). Recent reviews of global and regional Earth system modeling have elaborated past and recent trends and summarized future challenges for further development (Schrum, 2017; Giorgi and Gao, 2018; Giorgi 2019; Heinze et al., 2019; Jacob et al., 2020). Realized stepping stones on the road map towards coupled Earth system modeling for the Baltic Sea region include coupled ocean–sea-ice–atmosphere models, coupled ocean–wave and atmosphere–wave models, coupled vegetation–atmosphere models, and coupled ocean–atmosphere–hydrology models. This article aims at reviewing the latest developments, and describes the problems and benefits of using coupled models, with a focus on the specific demands of the Baltic Sea region within Europe. Thus, the emphasis is mainly set on the main physical feedbacks as they emerge in two-way coupled Earth system components. Therefore, this article does not aim to be comprehensive, as some important components are not considered, including biogeochemical nutrient cycling on land and in the Baltic Sea or atmospheric chemistry (e.g., the effect of aerosols).

2 Current state of knowledge

2.1 Land surface–atmosphere coupling and modeling

2.1.1 Biophysical mechanisms

The Baltic Sea region's land surface and terrestrial ecosystems couple strongly to local climate, and biophysical interactions determine the exchanges of momentum, heat and water between the land and the atmosphere as well as the spatiotemporal dynamics of the atmospheric boundary layer (ABL), including the wind speed, surface and air temperature, humidity and precipitation, and the atmospheric radiative balance.

Changes in the properties of the land surface and ecosystems, whether natural or anthropogenic, will lead to climate change through a number of well-established biogeochemical and biophysical feedback mechanisms (IPCC 2019). The biogeochemical mechanisms include the release or uptake of greenhouse gases (primarily CO₂, although CH₄ and N₂O are of considerable importance in the Baltic Sea region; Gao et al., 2014) and the emissions of black carbon, aerosol precursors (e.g., biogenic volatile organic compounds – BVOCs) and organic carbon aerosols that can alter the atmospheric composition, including cloud condensation nuclei and the fraction of diffuse and global radiation (e.g., Kulmala et al., 2014).

Important properties of the land surface include its albedo, its roughness, the species composition, the properties and phenology of green vegetation (e.g., leaf area index – LAI) and plant physiology (e.g., leaf stomatal and canopy conductance). Local climate is altered as a result of changes to the shortwave and longwave radiation, the turbulent fluxes of sensible and latent heat (i.e., evapotranspiration – ET) and momentum (e.g., Bonan 2008; Anderson et al., 2011; Pielke et al., 2011; Mahmood et al., 2014; Ellison et al., 2017; IPCC, 2019).

Forests typically have a lower surface albedo than grassland, pastures and cropland. Thus, deforestation tends to increase the albedo, whereas reforestation and afforestation have the opposite effect. Furthermore, coniferous forest albedo (0.05–0.15) is lower than deciduous forest albedo (0.15–0.20) (Anderson et al., 2011; Bonan 2008). Hence, the species composition determines the net albedo in a given region. As snow albedo ranges from 0.45 to 0.95 depending on age, history and mechanical disturbance, the albedo of the land surface in the Baltic Sea region is particularly sensitive to the duration and extent of snow cover as well as to the underlying vegetation type (Anderson et al., 2011). Reforestation or afforestation leads to a lower albedo in periods of snow cover, with greater net radiation at the land surface, stronger sensible heat fluxes and a warming over the forested area as a result (Anderson et al., 2011).

Since trees are taller than grasses and crops, forested regions in the Baltic Sea region have greater roughness lengths and tend to couple more strongly to the atmosphere, creating more turbulence than grass-covered regions or cropland, with higher sensible and latent heat fluxes. An increase in the surface roughness in association with reforestation or afforestation, increased rates of tree growth or altered forest management practices can lead to stronger turbulent fluxes and strong turbulent mixing (Winckler et al., 2019).

Trees transpire more water than grass or crops as a result of their larger leaf area and deeper roots. Thus, forests have higher evapotranspiration rates and latent heat fluxes than grasslands, although irrigated croplands can also have high ET rates. Higher latent heat fluxes cool the surface and moisten the ABL. A decrease in the latent heat fluxes associated with deforestation tends to warm the surface and leads to higher sensible heat fluxes and a warming of the ABL. The ABL is also drier, which can lead to a reduction in precipitation.

2.1.2 Biophysical effects of land-use and land-cover changes on climate in the Baltic Sea region

A number of modeling studies have examined the influence of land-use and land-cover change (LULCC) on climate variables in northern European domains, including the Baltic Sea region. The standard approach (Gálos et al., 2012; Strandberg et al., 2014; Strandberg and Kjellström, 2019) is to alter the static land-cover input to the coupled model and to compare simulations with an unchanged, control simulation. This approach does not permit two-way coupling in which local climate changes resulting from the perturbation subsequently alter land surface properties and vegetation characteristics.

Perugini et al. (2017; see also IPCC, 2019) reviewed the published literature on the biophysical effects of anthropogenic land-cover change on temperature and precipitation in boreal, temperate and tropical regions. A total of 28 studies were included in their review, including 3 based on observations and 25 that were based on idealized regional and global climate model simulations designed to estimate the regional and global biophysical effects of complete deforestation or afforestation. To effect deforestation in their simulations, some authors replaced forest with grassland, whereas other authors replaced forest with bare soil. Modeled deforestation in boreal regions resulted in local cooling consistent with observations but with a less consistent, slight cooling modeled in temperate regions in contrast to observations that indicate a slight warming in those zones.

Goa et al. (2014) used the REMO regional climate model (RCM) to investigate the biophysical effects of extensive peatland drainage and afforestation in Finland during the 20th century. Simulations were made for a model domain centered on Finland but covering a large part of the Baltic Sea region. The model grid had a horizontal resolution of 18 km with 27 vertical levels up to 25 km in the atmosphere.

Maps from the Finnish National Forest Inventory (FNFI) were used to compute changes to the fractional coverage of REMO's 10 land-cover classes in Finland from the 1920s, when peatland drainage and forestation began, to the 2000s. Over this period, coniferous forest replaced large regions previously covered by peat bogs and mixed forest. Two 18-year (1979–1996) simulations were then made to compare the effect of the land-cover change, driven with 6-hourly lateral boundary conditions from the ERA-Interim reanalysis data and identical land cover outside of Finland.

Goa et al. (2014) found that the reduction in albedo associated with prescribed peatland forestation resulted in an increase of up to 0.43 K in the 2 m air temperature in April, with the highest values being found over the most intensively forested areas. In contrast, there was a slight cooling (<0.1 K) during the growing season (May–October), associated with greater ET from coniferous forests.

Strandberg and Kjellström (2019) used the Rossby Centre RCM RCA4 (Strandberg et al., 2015; Kjellström et al., 2016), a successor to RCA3 (Samuelsson et al., 2011), to investigate and attribute the climate impacts of maximum potential afforestation or deforestation in Europe (see also the study by Gálos et al., 2012, using the REMO RCM). Horizontal grid spacing in RCA4 is approximately 50 km over the EURO-CORDEX domain covering Europe, and the model has 24 vertical levels in the atmosphere. For their study, Strandberg and Kjellström (2019) applied lateral boundary forcing (pressure, humidity, temperature and wind) every 6 h from ERA-Interim reanalysis data. Sea surface temperature and sea-ice extent were prescribed according to observations.

Three simulations, which only differed with respect to the land-cover map used, were performed for the 30-year period from 1981 to 2010. The control simulation used the standard, present-day land-cover map from RCA4 defined in the ECOCLIMAP (Champeaux et al., 2005) product. This map reflects the considerable agricultural activity in Europe, showing large areas with low forest cover in central, western and southern Europe. To effect maximum afforestation, Strandberg and Kjellström (2019) used the LPJ-GUESS dynamic vegetation model (Smith et al., 2001; see below) to produce a map of potential natural forest cover for Europe in equilibrium with present-day climate. Finally, maximum deforestation was implemented by converting forest fractions according to the potential forest cover to grassland in the control simulation.

In their analysis, Strandberg and Kjellström (2019) focused on winter (December–January–February, DJF) and summer (June–July–August, JJA) seasonal mean temperatures and precipitation as well as daily minimum and maximum temperatures. Afforestation decreased albedo in both seasons, especially in regions in eastern Europe with long snow cover and little forest cover. In contrast, deforestation increased albedo throughout the region, especially in the northern Baltic Sea region in winter.

Afforestation in Europe generally resulted in increased ET, as trees have a larger leaf area and deeper roots. This leads to colder near-surface temperatures in JJA. Deforestation gave the opposite effects, with warmer near-surface temperatures due to decreased ET in summer. In the Baltic Sea region, deforestation-induced reductions in ET coincide with the largest differences in the fraction of forest (e.g., a reduction of 20 %–35 % in Scandinavia). Interestingly, ET increases by 20 % over the Baltic Sea in JJA. Strandberg and Kjellström (2019) attribute this to a land–ocean coupling, whereby warmer and drier air from the surroundings (due to reduced ET over land) favors increased ET over the Baltic Sea when it comes into contact with the sea.

Afforestation in winter resulted in a decrease in temperature in central and southern Europe that can not be explained by albedo changes, as this is reduced during these months. As ET is also low during winter, the changes were attributed to atmospheric circulation changes resulting from increased roughness.

The winter low-pressure systems simulated by RCA4 lose their energy earlier because of increased friction over afforested areas due to their greater roughness, and afforestation generally leads to a simulated winter climate with less cyclonic activity in central Europe. Associated with this, the mean geopotential height at 500 hPa increased by 100 m in the Baltic Sea region.

Strandberg and Kjellström (2019) also found that the biophysical effects of afforestation or deforestation in Europe on daily minimum and maximum temperatures were stronger than the impacts on mean near-surface temperatures. In the case of afforestation, although DJF mean temperatures were reduced throughout most of Europe, there was a particularly strong warming of daily minimum temperatures (up to 2–6 ^∘C in Germany) that could be attributed to increased cloud cover and reductions in outgoing longwave radiation. During summer, on the other hand, the marked changes in mean temperatures were mainly caused by respective changes in daily maximum temperatures (i.e., decreases in the case of afforestation and increases in the case of deforestation). Finally, by rerunning the simulations and confining the applied deforestation and afforestation changes to the western and/or eastern parts of the model domain, the authors also showed that the climatic effects of afforestation or deforestation in Europe were mainly local.

Belušić et al. (2019) followed up on the study of Strandberg and Kjellström (2019) and used a cyclone-tracking algorithm to study how the same idealized deforestation and afforestation scenarios affected the number and intensity of cyclones as well as how they affected precipitation extremes. Consistent with the results of Strandberg and Kjellström (2019), Belušić et al. (2019) found that the larger surface roughness after afforestation reduced the number of cyclones over Europe compared with the deforestation and control simulations: differences were 20 %–80 % near the Baltic Sea region, 10 % in regions near the western European coast and they increased towards the east to reach 80 %. This resulted in a reduction in winter precipitation extremes of up to 25 % across the domain.

Figure 2 Outline of the biogeochemical (upper part in brown) and the biophysical (lower part in blue) influence of reforestation/afforestation or enhanced forest productivity in the Baltic Sea region on near-surface temperatures. The overall effect on near-surface temperatures varies by season and region, depending, for instance, on snow cover and incoming solar radiation (adapted from May et al., 2020).

Figure 2 summarizes the various biophysical and biogeochemical influences of reforestation and afforestation or enhanced forest productivity on near-surface temperatures in the Baltic Sea region. According to the figure, the decrease in surface albedo (resulting in increased absorption of incoming solar radiation at the land surface) is the only effect that leads to a warming of near-surface temperatures, while all the other effects lead to a cooling. The biogeochemical effects shown in Fig. 2 are increased carbon (C) storage, which weakens the radiative forcing, and more aerosols, which reduce the solar radiation reaching the land surface due to additional scattering from particles and more clouds. The biophysical effects shown are increased roughness length, which enhances the turbulent fluxes of energy and momentum, and increased evapotranspiration in late spring and summer, which strengthens the fluxes of latent heat and weakens the fluxes of sensible heat. The magnitude of the overall cooling in the Baltic Sea region associated with reforestation and afforestation or enhanced forest productivity depends on the significance of the warming effect compared with the cooling from the other biophysical and biogeochemical effects.

Following up on their earlier, theoretical and observation-based studies in which changes in land management were shown to affect surface temperature to a degree similar to changes in land-cover type (Luyssaert et al., 2014), Luyssaert et al. (2018) used the ORCHIDEE-CAN land surface model (further developed to explicitly take the biogeochemical and biophysical effects of land-use change and management into account) coupled to the LMDZ atmospheric circulation model to investigate the trade-offs associated with using European forests to meet the climate objectives of the Paris Agreement. Their analyses clearly demonstrate that the biophysical effects of forest management must be taken into account in any assessment of climate mitigation strategies, with consequences for policy and forestry in the Baltic Sea region and beyond (e.g., in relation to the optimal balance of coniferous and deciduous forest in the region). Similarly, Kumkar et al. (2020) used offline simulations with the CLM4.5 land surface model to quantify the sensitivity of turbulent fluxes and land surface temperatures in Fennoscandia to scenarios of changes to forest composition and structure relative to their present-day values. They found that replacing conifers with deciduous forests could cool the surface by 0.16 K annually and by 0.3 K is summer months, mainly as a result of their higher albedo, and they identified important differences between developed and underdeveloped forests, with the latter having lower evaporation rates as a result of their lower LAI and canopy height.

2.1.3 Dynamic global vegetation models applied to the Baltic Sea region

Dynamic global vegetation models (DGVMs) are numerical models of terrestrial ecosystems that simulate the properties, dynamics and functioning of potential, natural and managed vegetation and their associated biogeochemical and hydrological cycles as a response to climate and environmental change. Prentice et al. (2007) summarize their historical development, design and construction principles as well as the processes typically included, their evaluation and examples of their application. DGVMs incorporate research and knowledge from different disciplines, including plant geography; plant physiography; biogeochemistry, including soil biogeochemistry; vegetation dynamics and demography; biophysics; agriculture; and forest management. DGVMs have been used to study the observed and expected impacts on terrestrial ecosystems in the Baltic Sea region resulting from climate and environmental change. An understanding of these impacts is a necessary first step to comprehending the dynamics in coupled RESMs, where ecosystem change is allowed to influence local and regional climate through the biophysical feedback mechanisms outlined above.

Recent works have shown that in order to realistically simulate ecosystem carbon balance, climate responses, and ecosystem recovery following disturbances due to land-use change, management interventions, and natural disturbance processes such as fires and storms (Fisher et al., 2018; Pugh et al., 2019), it is important to incorporate the size- and age-structure and demography of vegetation and ecosystems explicitly, and to account for the competitive interactions of growing vegetation stands comprising individuals or cohorts of different plant functional types. A number of DGVMs and land surface models (LSMs) are now moving in this direction, away from a traditional tiled or area-based (Smith et al., 2001) land surface representation, including the CLM4.5(ED) LSM (Moorcroft et al., 2001; Fisher et al., 2015, 2018) and the Functionally Assembled Terrestrial Ecosystem Simulator (FATES) vegetation demography submodel (Koven et al., 2020), the POP (Population Orders Physiology) module for woody demography in the CABLE LSM (Haverd et al., 2013, 2014, 2018), the RED module in the JULES LSM (Argles et al., 2020) and the SEIB-DGVM (Sato et al., 2007). To date, however, no demographic model has been applied to study the Baltic Sea region specifically (see Kumkar et al., 2020, for an application using CLM4.5) apart from the LPJ-GUESS DGVM (Smith et al., 2001, 2014). LPJ-GUESS explicitly represents the size, age structure, spatial heterogeneity and temporal dynamics of co-occurring cohorts of plant functional types (PFTs), i.e., classifications of plants according to their physical, phylogenetic and phenological characteristics (e.g., boreal or temperate, evergreen or deciduous, and broadleaf or needleleaf trees in the Baltic Sea region, and herbaceous species; Prentice et al., 2007) or species (Hickler et al., 2012) that compete in natural and managed stands (forestry, crops and pasture), in response to climate, atmospheric CO₂ and nitrogen (N) availability. As the stand structure evolves in response to environmental change and impacts the availability of key resources, the growth, survival and the outcome of competition are affected.

LPJ-GUESS represents different land use and management in separate stands (Lindeskog et al., 2013). The fraction of the grid cell covered by each stand (e.g., forest, natural and cropland) type can change in time, following external land-use datasets (e.g., Hurtt et al., 2020). LPJ-GUESS also allows for detailed management interventions for representative crops (represented as crop functional types – CFTs), grassland grazing, mowing and fertilization (Olin et al., 2015a, b) as well as clear-cut and continuous-cover forest management (Lagergren and Jönsson, 2017). Disturbances due to management actions such as forest clearing, prognostic wildfires and a stochastic generic disturbance regime induce biomass loss and reset vegetation succession (Smith et al., 2001). N-cycle-induced limitations on natural vegetation and crop growth, C–N dynamics in soil biogeochemistry and N trace gas emissions are included (e.g., Smith et al., 2014; Olin et al., 2015a, b) as well as BVOC (isoprene and monoterpene) emissions (Hantson et al., 2017).

LPJ-GUESS output variables describe the vegetation state (PFT/species composition, LAI, vegetation height, biomass, tree density, and burned area), variables relating to the state and functioning of the soil (water content, C and N content, temperature, runoff, N leaching, and loss of dissolved organic C and N), and climatically important fluxes to and from each simulated stand (evapotranspiration, gross and net primary productivity, autotrophic and heterotrophic respiration, fluxes from wildfires, CH₄ and N trace gases, BVOCs, and net ecosystem carbon exchange).

2.1.4 Modeling terrestrial ecosystems in the Baltic Sea region

LPJ-GUESS has been applied in many studies to simulate terrestrial ecosystems in the Baltic Sea region, under current, future, and historic and preindustrial climate conditions. Koca et al. (2006) simulated the impacts of climate change on natural ecosystems in Sweden in response to different regional climate change scenarios. In all of the climate scenarios considered, the authors observed an increase in plant productivity and LAI, and a northward and upward advance of the boreal forest tree line by the end of the 21st century. The current dominance of Norway spruce and to a lesser extent Scots pine was found to be reduced in favor of deciduous broadleaf tree species in future scenarios across the boreal and boreo-nemoral zones. These changes are consistent with earlier studies (Miller et al., 2008) of the effects of climate and biotic drivers on Holocene vegetation in Sweden and Finland, where observed changes to the northern distribution limits of temperate trees and species at the tree line were attributed to millennial variations in summer and winter temperatures.

Hickler et al. (2012) reparameterized the most common European tree species in LPJ-GUESS and forced the model with an atmosphere–ocean general circulation model (AOGCM) climate scenario output, downscaled to a spatial resolution of 10 arcmin×10 arcmin. Climate change and CO₂ increase resulted in large-scale successional shifts, with 31 %–42 % of the total area of Europe projected to be covered by a different vegetation type by the year 2085 depending on the scenario used. Consistent with the earlier results of Koca et al. (2006), trees replace tundra in arctic and alpine ecosystems, and temperate broadleaf forest replaces boreal conifer forest in the Baltic Sea region.

2.1.5 Regional Earth system modeling with interactive vegetation dynamics in the Baltic Sea region

Coupled regional Earth system models (RESMs) extend RCMs to include the terrestrial biosphere as an integral dynamic component interacting in a two-way coupling with the atmosphere, with representations of both vegetation dynamics and terrestrial biogeochemistry. Such a framework allows for modeled natural and managed ecosystems to respond to climate and environmental change and to influence local and regional climate through the biophysical feedback mechanisms outlined above.

RCA-GUESS was the first published and evaluated RESM (Smith et al., 2011) to include the terrestrial biosphere as an integral dynamic component, and it couples LPJ-GUESS to the RCA3 RCM (Samuelsson et al., 2011). In its RCA-GUESS configuration, LPJ-GUESS is driven by the daily mean temperature, soil water content, precipitation and downward shortwave radiation simulated by RCA3, and the CO₂ concentration is read from the same source used to force RCA.

In its uncoupled configuration, the land surface scheme of RCA3 uses ECOCLIMAP to specify the cover fractions of two vegetated land surface tiles, one representing the type of forest (broadleaf or needleleaf) and the other open land (including crops, pasture and grassland) for each grid cell in its domain.

In RCA-GUESS, LPJ-GUESS replaces the static ECOCLIMAP land-cover description and aggregates its vegetation fields to update the tile fractions, their type and their associated LAI. The specific forest PFTs simulated by LPJ-GUESS are aggregated into needleleaf and broadleaf trees before providing the information to RCA, while open land includes a varying coverage of herbaceous vegetation.

By changing the relative fractions and types in RCA, the LPJ-GUESS fields dynamically determine and update the surface albedo, LAI, surface roughness and conductance in RCA grid cells. For example, albedo is calculated in RCA using a weighted average of prescribed albedo values for needle leaved and broadleaf trees, open land vegetation, snow and bare soil. Similarly, the fluxes of sensible and latent heat are calculated as weighted averages of the individual tiles.

Terrestrial CO₂ exchange is simulated by LPJ-GUESS, enabling biogeochemical ecosystem responses to be assessed consistently with the biophysical land–atmosphere interactions (Zhang et al., 2014). However, the atmospheric CO₂ concentrations are not updated over the limited domain covered by RCA. Thus, although biophysical feedback loops are closed in RCA-GUESS, biogeochemical feedback loops are open.

Wramneby et al. (2010) used RCA-GUESS to identify hot spots of biophysical vegetation–climate feedbacks for future climate conditions in Europe. Two simulations – feedback and non-feedback – were run over Europe for 1961–2100 to isolate the effect of feedbacks from vegetation dynamics. In the feedback simulation, RCA and LPJ-GUESS were coupled throughout the entire simulation period. In the non-feedback simulation, land cover in RCA was prescribed and fixed over the full simulation period from the long-term means from LPJ-GUESS output in the coupled simulation for 1961–1990. The difference between the climate change signal (2071–2100 minus 1961–1990) from the feedback simulation and the corresponding signal from the non-feedback simulation was used to calculate the additional contribution of the vegetation–climate feedback to the background climate changes simulated by RCA3, driven by lateral AOGCM forcing.

Wramneby et al. (2010) showed that the snow-masking effect of forest expansion and greater LAI in the Scandinavian mountains as well as the consequent reductions in albedo enhanced the winter warming trend. In central Europe, the stimulation of photosynthesis and plant growth caused by the increased CO₂ concentration, longer growing seasons and warming mitigated the future warming through a negative feedback due to enhanced ET associated with the increased vegetation cover and LAI.

Zhang et al. (2014, 2018) applied RCA-GUESS over the Arctic domain of the Coordinated Regional Climate Downscaling Experiment (CORDEX-Arctic) – which includes the northern part of Baltic Sea region – to investigate the role of the biophysical feedbacks from vegetation in the Arctic region under three different climate scenarios (RCP2.6, RCP4.5 and RCP8.5, where RCP stands for Representative Concentration Pathway). Zhang et al. (2018) found that warming and CO₂ increases promote productivity increases, LAI increases and tree line advance into the Arctic tundra, with the consequence that two biophysical effects have the potential to alter the spatiotemporal signal of future climate change in the Arctic region: interactive vegetation results in albedo-mediated warming in early spring and ET-mediated cooling in summer, amplifying or modulating local warming and enhancing summer precipitation over land.

2.2 Ocean–atmosphere coupling

The treatment of ocean–atmosphere exchange of momentum, heat and mass substantially differs between coupled and uncoupled models. Basically, in the coupled mode, the ocean is driven by fluxes and sea level pressure calculated in the atmosphere model (e.g., Fig. 3) and is used to drive the coupled ocean model which, in turn, communicates simulated fields of sea ice and surface water temperature (SST) to the atmosphere model. By contrast, uncoupled ocean models use atmospheric forcing fields (passively coupled ocean, PCO) to calculate air–sea fluxes with a bulk formula. Stand-alone atmosphere models usually read prescribed fields of sea ice and SST from reanalysis datasets or model data to calculate fluxes.

Figure 3 Mass, momentum and heat exchange as realized in the atmosphere–ocean model RCA4-NEMO. The abbreviations used in the figure are as follows: ICA – interactively coupled atmosphere, PCA – passively coupled atmosphere, ICO – interactively coupled ocean and PCO – passively coupled ocean. Source: Gröger et al. (2015).

2.2.1 Impact on mean climate

One of the main questions addressed so far is if an interactively coupled atmosphere–ocean model would significantly change the long-term climate compared with their stand-alone atmosphere and ocean modules. This has been investigated in a number of studies (e.g., van Pham et al., 2014; Tian et al., 2013; Gröger et al., 2015, 2021; Primo et al., 2019; Kelemen et al., 2019; Akhtar et al., 2019; Cabos et al., 2020). Van Pham et al. (2014) developed a regional atmosphere–ocean general circulation model (RAOGCM) built upon the ocean model NEMO coupled to the atmosphere model COSMO-CLM (the COSMO model in Climate Mode, hereafter denoted as CCLM) for the EURO-CORDEX domain. The coupling domain encompassed the Baltic Sea and the North Sea until 59^∘ N and 4^∘ W. As forcing as well as at the lateral model boundaries, the authors applied ERA-Interim reanalysis data. It is noteworthy that the coupled system showed a systematically lower long-term mean 2 m air temperature (T2m) in the Baltic Sea and North Sea compared with the uncoupled atmosphere model. Consequently, interactive coupling reduced the model's mean bias in T2m compared with the E-OBS observational dataset as a reference. Interestingly, the authors found most significant changes between coupled and uncoupled runs over continental areas of central and eastern Europe (i.e., at a distance from the coupled areas). Analysis of the modeled wind field indicated that these areas were situated downwind from the coupled domain in the North Sea and Baltic Sea, implying a more complex pattern of atmospheric advection of temperature anomalies which was not further investigated.

However, similar experiments by Gröger et al. (2015), using nearly the same ocean model NEMO but coupled to the regional atmosphere model RCA4, were somewhat contradictory. In their setup forced by ERA40 reanalysis data with lateral boundaries prescribed by ORAS4 (Balsameda et al., 2013), the coupled system showed generally warmer near-surface air temperatures over the Baltic Sea compared with the uncoupled RCA4 model. Furthermore, no significant differences were found in air temperatures over land between coupled and uncoupled simulations. With respect to sea surface temperatures, Gröger et al. (2015) found the strongest differences in the Baltic Sea, where SSTs in both the coupled and uncoupled system were too cold in winter compared with satellite products. However, winter SSTs were significantly higher in the coupled model (thereby reducing the bias) due to seasonally varying feedback loops controlling the ocean–atmosphere heat exchange. Figure 4 sketches the main mechanisms comprising the atmosphere–ocean feedbacks during summer and winter that control the SST in the Baltic Sea.

Figure 4 Side (a) of the figure shows a positive winter feedback loop, and side (b) shows a negative summer short-circuit. Drawn after Gröger et al. (2015).

Winter thermal–wind mixing positive feedback

During winter, the Baltic Sea is usually warmer than the atmosphere, supporting a net heat flux out of the ocean.

During ocean offline simulations, in which the ocean model was forced with ERA40 atmospheric reanalysis data, simulated winter SSTs in the Baltic Sea showed a strong bias compared with observed SSTs. This was mainly caused by a cold bias of the ERA40 dataset over Europe. In the coupled mode and driven by the same dataset at the lateral boundaries, the bias nearly vanished due to a thermal feedback loop between the ocean and the atmosphere which resulted in stronger vertical mixing and increased transport of warmer deep waters to the surface (thereby reducing the cold bias at the surface compared with the uncoupled simulation). This is shown in Fig. 4a. In the coupled model, the atmospheric boundary layer is disturbed by warm anomalies generated in the Baltic Sea. This promotes stronger winds that, in turn, feed back to the ocean with stronger vertical mixing, thereby increasing heat exchange with warmer, deeper water layers. As a result, the ocean model's cold bias decreases in the Baltic Sea compared with the ocean stand-alone model which uses prescribed atmospheric boundary conditions that can not respond to SST anomalies (Gröger et al., 2015).

Summer thermal short-circuit

During summer, the above feedback loop is bypassed by the inverse thermal air–sea contrast. As the atmosphere is generally warmer than the ocean in summer, any wind-induced upward mixing of cold, deep water will tend to stabilize the atmospheric boundary layer with a negative effect on wind strength (Fig. 4b). This was demonstrated at stations in stratified areas of the Baltic Sea and in the North Sea using lead correlation analysis between the 10 m wind velocity and SST after a short-term event of strengthened winds (Gröger et. al., 2015). During the first 70 h after the event, the wind and SST were negatively correlated (with a peak r=0.7 at around 30 h), implying decreasing SST with stronger wind mixing with colder deep waters. After 70 h, the correlation turned to positive values with a peak (r=0.7) around 130 h as the colder water surface stabilized the atmospheric boundary layer and, thus, promoted lower wind speeds. Following this, wind mixing ceased again giving rise to heat gain from the warmer atmosphere Gröger et al. (2015).

These results highlight the importance of thermal air–sea coupling in midlatitude marginal seas and are supported by a number of different studies. Tian et al. (2013) drove a coupled ocean–atmosphere model with ERA-Interim reanalysis data (ERA-I). Similar to Gröger et al. (2015), the abovementioned authors also found overly cold SSTs in winter but with a substantial lower bias in the coupled model. No detailed feedback analysis was carried out, but it is likely that the positive winter feedback was also present in the model of Tian et al. (2013). Moreover, Keleman et al. (2019) applied a model that couples an atmosphere model for the EURO-CORDEX region to regional ocean models for the Mediterranean, the North Sea and the Baltic Sea. They found a significant improvement in winter precipitation patterns over eastern Europe through the altered representation of SSTs. However, the SSTs themselves were not validated; thus, a general conclusion on the performance of the whole system can not be drawn. Primo et al. (2019) used the same model and showed that air–sea coupling can improve the representation of heat and cold waves but concluded that a general judgment regarding the whole system is difficult to draw and depends on the variable considered.

Impact on climate indices

Apart from the thermal coupling effect on mean climate, we now consider temperature-related climate indices often used in regional climate assessments such as CORDEX (e.g., Jacob et al., 2014; Teichmann et al., 2018; Kjellström et al., 2018; Gröger et al., 2021). The indices are strongly related to the ocean heat content and are, therefore, sensitive to the coupling. In particular, we focus on three indices with importance for human health, agriculture and the Baltic Sea ecosystem: (1) the number of tropical nights, which are defined as nights when the daily minimum temperature does not fall below 20 ^∘C (e.g., Fischer and Schär, 2010; Teichmann et al., 2018; Meier et al., 2019a; Gröger et al., 2021); (2) the number of frost days, which are defined as days when the daily minimum temperature falls below 0 ^∘C; and (3) the number of warm periods, defined as at least 3 consecutive days on which the daily maximum temperature reaches 20 ^∘C. For our analysis, we take available data from hindcast runs, as described in Gröger et al. (2015). All indices are derived from the 1970–1999 reference period from both the coupled and the uncoupled simulations. Both runs are driven with ERA40 atmospheric reanalysis data at the lateral model domain boundaries. Lateral boundary conditions were set according to the ocean reanalysis ORAS4 (Balsameda et al., 2013). For details about the coupling, we refer the reader to Gröger et al. (2015) and Dieterich et al. (2019b).

Figure 5 The left panel of (a) shows the number of tropical nights diagnosed from the coupled regional ocean–atmosphere model during the 1970–1999 reference period. The right panel of (a) displays the difference (%) of the coupled minus the uncoupled model. The panels in (b) are the same as those in (a) but for the number of frost days. The left panel of (c) shows the number of periods of at least 3 consecutive warm days (days when the daily maximum temperature exceeds 20 ^∘C). The right panel of (c) displays the difference (%) of the coupled minus the uncoupled model.

Figure 5a displays the number of tropical nights simulated with the coupled ocean–atmosphere model RCA4-NEMO (Wang et al., 2015; Dieterich et al., 2019b; Gröger et al., 2019, 2021). A clear land–sea pattern is seen with only sporadic occurrences over land with the exception of the southern part of the Iberian Peninsula and the Pannonian Basin south of the Carpathians. The higher effective heat capacity of the ocean is responsible for the frequent occurrences over the Mediterranean and over the southern North Atlantic sector, which also includes the Bay of Biscay. The thermal effect of the Baltic Sea water body is obvious. Unlike the North Sea, which receives colder waters from the North Atlantic, the Baltic Sea has very limited exchange of waters with the adjacent North Sea. In addition, a strong seasonal thermocline during summer also limits the exchange with colder waters from deeper layers. These two processes support a strong warming of the Baltic Sea during summer. Consequently, this prevents the air temperature from falling below 20 ^∘C during warm periods in the second half of the summer. As a result, the Baltic Sea displays a range of tropical nights that matches the range found further to the south, such as north of the Black Sea or in parts of the western Atlantic off northern Iberia (Fig. 5a, left). Figure 5a (right) demonstrates that the abovementioned thermal effect over the Baltic Sea is much more pronounced in the coupled model including the dynamic ocean model. As a result, the number of detected tropical nights within the reference period increases by between 50 % and 100 % over the southern Baltic Sea in the coupled model compared with the stand-alone atmosphere model.

Figure 5b (left) shows the number of frost days during the reference period. A clear land–sea pattern is seen with strongly diminished occurrences over the open-ocean areas in the north, whereas they are completely absent over the southern Mediterranean and the southern part of the Atlantic. Again, the different thermal behavior between the North Sea and the Baltic Sea is obvious. With respect to its thermal behavior, the Baltic Sea is similar to the continents and supports a large number of frost days. This is related to the strong winter halocline that hampers wind-forced mixing and convective mixing with deep waters. Consequently, the upper-layer water body of the Baltic Sea can rapidly cool during winter. In contrast, the adjacent North Sea effectively damps the occurrence of frost days. However, a pronounced east–west gradient is visible with fewer occurrences in the west. Here, warmer waters from the Atlantic enter the North Sea and spread southward. The eastern part of the North Sea is influenced by low-salinity waters derived from the Baltic Sea. These waters flow northward along the Norwegian coast and impose a haline stratification there. This results in similar thermal behavior to that discussed for the Baltic Sea. The southern part of the eastern North Sea, namely the German Bight, is shallow and supports rapid cooling. Altogether, this creates the east–west gradient seen in Fig. 5b (left). Generally, these thermally forced processes are represented in both the coupled and uncoupled model versions, but in the coupled model, the number of frost days is significantly decreased over almost the entire Baltic Sea (Fig. 5b, right). Here, the aforementioned winter mixing feedback loop operates (i.e., stronger mixing in the coupled model increases the winter sea surface temperature with a positive feedback on wind speed, resulting in a significantly higher sea surface temperature and reducing the number of frost days). Over the North Sea, coupling generates a positive anomaly along a band between the 2–4^∘ E meridian. It is likely that this reflects shifts in the gradients caused by slightly altered flow paths of the water masses derived from the Baltic Sea and the North Sea.

Finally, the coupling effect of warm periods is displayed in Fig. 5c. In the Baltic Sea, such periods are an important precondition for the occurrence of cyanobacteria blooms during summer. Moreover, the southern Baltic Sea is found be a hot spot with respect to the effect of interactive air–sea thermal coupling, as the number of such periods in the coupled model exceeds the corresponding number in the stand-alone atmosphere model by several orders of magnitude. Here, the effect of the aforementioned summer thermal short-circuit is seen. Enhanced mixing by winds brings cooler waters from the depths to the surface. The cooler surface water also lowers the air temperatures; this imposes a stabilizing effect on the atmospheric boundary layer over sea, thereby further damping wind strength (see Gröger et al., 2015, for details). This facilitates the development of longer-lasting warm periods in summer.

So far, we have discussed the thermal effects of damping (in the case of frost days) and amplifying (warm days and tropical nights). The associated feedback loops are more realistically represented in the fully coupled model, as previously explained. However, effects outside of the active coupling area over land seem to be rather weak. On the other hand, we note that the interactive coupling area must be considered small compared with the whole domain. Thus, extending the area of interactive coupling by, for example, also including the Mediterranean or larger parts of the North Atlantic may result in more intense effects when using coupled models (Primo et al., 2019; Keleman et al., 2019; Akhtar et al., 2019; Cabos et al., 2020).

2.2.2 Impact on extreme events

Apart from the representation of mean climate, many air–sea coupling processes are important in generating hazardous events such as extreme precipitation, storm track paths or flooding (see also Rutgersson et al., 2021). Often these events are generated remotely over the open ocean and, thus, require a realistic representation of the ocean's surface. Therefore, when no high-resolution ocean model is coupled to the regional atmosphere model, the sea surface has to be represented by available products for sea surface temperatures (SSTs) and sea ice. These products are often limited in quality and frequency. For example, the ERA40 SSTs, which are often used in uncoupled atmosphere simulations, are based on a weekly or even monthly frequency (Fiorino, 2004; Uppala et al., 2005). Consequently, many studies have compared coupled and uncoupled models with respect to the representation of extreme events in hindcast simulation. More recent products such as ERA5 may improve the situation by providing higher spatial and temporal resolution data.

Jeworrek et al. (2017) diagnosed a better representation of atmospheric conditions favorable for the occurrence of convective snow bands. The authors attributed this improvement to a more accurate simulation of SSTs and subsequent air–sea heat and moisture fluxes in a case where the atmosphere model RCA4 was coupled to the ocean model NEMO setup for the North Sea and Baltic Sea. The importance of accurate SSTs for the representation of snow bands was recognized early by Gustafsson et al. (1998).

Comparing two regional ocean–atmosphere models in coupled and uncoupled mode, Ho-Hagemann (2015, 2017) found that interactive air–sea coupling can alter extreme precipitation over land. The abovementioned authors pointed out that the coupled model COSTRICE improved the low-level large-scale moisture convergence over the North Atlantic and the moisture transport towards central Europe. As a result, the simulated summer heavy rainfall improved compared with the stand-alone atmospheric model CCLM. This was demonstrated for several flood events in central Europe. A diagram displaying the cause and effects as well as feedbacks between the Earth system components is shown in Fig. 6 in order to explain the physical mechanism behind the improved representation of heavy rainfall. The main effect is an altered SST in the coupled model which further influences wind speed and evaporation – the so called wind–evaporation–SST (WES) feedback (Xie and Philander, 1994). When the wind speed increases over an area, evaporation increases and the latent heat flux is subsequently enhanced, which often leads to ocean cooling over the area in question. The lowered SST generates a horizontal SST gradient on the sea surface and also increases the land–sea heat contrast which, in turn, supports increasing wind speed. Stronger wind over the North Sea then generates a larger latent heat flux from the ocean to the atmosphere and intensifies the low over central Europe and the North Sea, which both support the large-scale moisture convergence from the North Sea to central Europe. A review of the influence of atmosphere–ocean interactions on heavy rainfall over Europe is available in Ho-Hagemann and Rockel (2018).

Figure 6 An air–sea feedback and interaction diagram. For each arrow, the initial status indicates that the source quantity is increasing, and the sign (− or +) indicates the changing tendency of the target quantity. Colors denote the group of changes or states over sea (blue), land (brown) and the land–sea interactions (green). Source: Ho-Hagemann et al. (2017).

An improved representation of extreme and mean temperatures in the CCLM atmosphere models is also reported in Primo et al. (2019) and Kelemen et al. (2019) when coupled to regional high-resolution ocean models for the North Sea, the Baltic Sea and the Mediterranean. Akhtar et al. (2014) analyzed 11 historical medicane events simulated by the atmosphere-only model CCLM and the coupled model CCLM-NEMOMED12 with different horizontal atmospheric grid spacings of 0.44, 0.22 and 0.08^∘. In this analysis, the coupled simulations improved significantly compared with atmosphere-only simulations at higher atmospheric grid resolution (0.08). The characteristic features of medicanes, such as warm cores and high wind speeds, are more intense in coupled simulations compared with atmosphere-only simulations. Akhtar et al. (2019) also demonstrated improved simulations of cyclones over the Mediterranean in the coupled system model of CCLM and two ocean models, NEMO-Nordic and NEMOMED12, compared with the atmosphere-only model with prescribed SSTs.

The impact of air–sea coupling on simulations of midlatitude cyclones was recently investigated by Ho-Hagemann et al. (2020) using an ensemble approach. The coupled model GCOAST-AHOI reduces the large spread of wind speed, mean sea level pressure, surface temperature, cloud cover and radiation fluxes amongst ensemble members of the atmospheric model during Cyclone Christian, which occurred in northern Europe between 27 and 29 October 2013.

2.2.3 Influence of the size of the coupling area

As outlined above, the size of the air–sea coupling area will influence how strong the coupling effect will be and how far it may propagate further over land. Table 1 lists several RAOGCMs applied for the Baltic region that have been developed for various areas within the recent decade. In these regional coupled system models (RCSMs), the atmospheric models cover different domains, such as the EURO-CORDEX domain, the European domain or the Arctic. The domain of the ocean model can be global with increased resolution over the North Sea and Baltic Sea as in MPIOM (Sein et al., 2015), or regional as in the other RCSM setups. For example, in a common NEMO setup used by several RAOGCMs (e.g., Van Pham et al., 2014; Gröger et al., 2015; Wang et al., 2015; Dieterich et al., 2019b; Akhtar et al., 2019; Primo et al., 2019), the Baltic Sea and North Sea regions are considered. In other coupled systems, such as Ho-Hagemann et al. (2015, 2017, 2020), the ocean model domain extends from the Baltic Sea and North Sea region into a part of the North Atlantic. Moreover, one atmospheric model can be coupled to more than one ocean model as in Akhtar et al. (2019) and Primo et al. (2019).

Table 1 Regional coupled climate models and their air–sea coupling domains.

A recent assessment of regional coupled modeling (Schrum, 2017) emphasized that the location and extension of the coupled region (Sein et al., 2014), the coupling frequency (Fang et al., 2009), and the quality of initialization and boundary forcing (Wei et al., 2014) are critical. In the following, we focus on the size of the coupling area and consider its potential impact on climate simulations. We will also elaborate on options with respect to suitable sizes for the coupling area.

Sein et al. (2014) pointed out that the choice of the coupled model domain based on simple geographical arguments is not sufficient, and decisions should be made based on the fundamental understanding of oceanographic and atmospheric processes and their feedbacks. Mikolajewicz et al. (2005) discussed the challenges regarding the global mass and energy balance arising from coupling between a global model and a regional model. Inconsistencies can occur when the inflow to the coupled region is calculated based on the global forcing while outflow is calculated based on the regionally coupled model solution. On the one hand, using a larger domain for the ocean model gives the ocean more degrees of freedom (Mikolajewicz et al., 2005) by putting the boundary conditions in the deep ocean in the North Atlantic and not in the North Sea; on the other hand, more air–sea coupling effects over the North Atlantic on the simulated climate over Europe can be taken into account.

In studies where the coupling domain covers a relatively large area inside the entire integration domain and commonly comprises multiple seas with a large heat inventory, the air–sea coupling effect is often found to extend far inland (e.g., Somot et al., 2008; Ratnam et al., 2009; Ho-Hagemann et al., 2015, 2017). Gröger et al. (2021) hypothesize that SST anomalies must have a critical extension and be linked to a sufficiently large heat content of the underlying ocean to impose a significant effect on large-scale atmospheric circulation. Otherwise, the fast transport (relative to the ocean) within moving air masses supports a rapid dispersion of temperature anomalies in the atmosphere. Li (2006) indicated that varying the SST over the Mediterranean Sea could initiate atmospheric teleconnections, which can influence precipitation over remote regions such as the Atlantic–European region.

A known problem for many atmospheric models (Vidale et al., 2003) is a dry bias over large areas of midlatitude continents. Sensitivity experiments with different regional models and different resolutions showed that interactive coupling can reduce this bias in those models that include parts of the North Atlantic in the coupled domain (Ho-Hagemann et al., 2017). A large part of precipitation over Europe is linked to moisture originating from the North Atlantic. Thus, a realistic moisture convergence over the North Atlantic–European region is essential to obtain good precipitation patterns. The authors concluded that, in the presence of precipitation biases in atmosphere models, the realistic simulation of air–sea feedbacks enhances the large-scale wind and evaporation via alteration of the sea surface temperature and the land–sea heat contrast and, therefore, reduces the dry bias.

2.2.4 Atmosphere–sea-ice–ocean modeling

The Baltic Sea is seasonally covered by sea ice, and the importance of this ice's influence on the general state of the Baltic Sea is unquestionable. Ice cover creates a barrier between the atmosphere and the sea that results in a direct impact on the exchange of mass, energy and momentum. Ice significantly modifies or even eliminates the interaction between the atmosphere and the sea. Furthermore, ice and snow can reflect up to 90 % of incoming solar radiation instead of the high absorption of this radiation by the sea surface. This albedo-related positive feedback effect is the main reason behind the amplification of climate change in the polar regions.

In a coupled modeling system, the most important prognostic sea-ice parameters are ice/snow surface temperature, albedo, ice concentration, growth rate of ice, and surface and bottom roughness, as these factors control radiation, heat, moisture and momentum fluxes at the atmosphere–ice–ocean interface. The growth rate and temperature of sea ice impact the salt flux between the ice and ocean. However, due to the low salinity of the northern Baltic, which experiences annual ice cover, this mechanism does not have a significant effect.

The theoretical frameworks of presently used sea-ice models in climate applications were established in the 1970s. The thermodynamical evolution of ice and snow is based on the classical heat conduction law, which was first numerically resolved by Maykut and Untersteiner (1971). The model resolves the vertical temperature and salinity structure and the surface temperature iteratively. The present community sea-ice models LIM-3/SI³ (Rousset et al., 2015; Vancoppenolle et al., 2009) and CICE (Hunke and Dukowitc, 1997) apply this classical framework but include detailed parametrization of snow, flooding, snow–ice formations, melt ponds, albedo and brines. For climate applications, the vertical structure is usually resolved for one to three layers. The one-layer model assumes a linear temperature profile on ice. This is a valid approximation for the Baltic Sea, where thermodynamically grown sea ice rarely exceeds 1 m.

The momentum balance equation of sea ice includes wind stress, bottom stress due to the ocean current, sea surface tilt, internal stress of the ice pack and the Coriolis force. The main uncertain term is the internal stress of sea ice. In the Baltic, this term can be dominant due to the large effect of the coastline and islands. In present coupled models, two rheological solutions are commonly used: the Hibler (1979) viscous-plastic rheology (VP) implies that the bulk and shear viscosities are constant and the model produces linear viscous behavior under very low strain rates; otherwise, the viscosities are calculated according to the plastic flow rule. VP rheology resolves the nonlinear behavior of Baltic sea-ice dynamics well (Leppäranta et al., 1998). The LIM-3/SI³ model applies VP rheology, and it is a rheological choice of the NEMO-Nordic model (Pemberton et al., 2017; Hordoir et al., 2019). However, as VP rheology is computationally demanding, numerically more feasible elasto-viscous-plastic rheology (EVP; Hunke and Dukovich, 1997) is used in the CICE model. It also widely used in Baltic Sea applications (Meier, 2002b; Jakacki and Meler, 2018; Janecki et al., 2018).

The third element of the sea-ice models is the resolution of ice thickness distribution g(h) (Thorndike et al., 1975). In the classical Hibler (1979) model, g(h) is approximated with two ice thickness categories: thin ice, interpreted as open water, and thick ice. The choice of the minimum ice thickness has an impact on the modeled ice edge. Ridging of ice is taken into account, as ice thickness can freely increase during the convergent ice motion, although the ice concentration is constrained to be 1.0 at maximum. To solve g(h) numerically, several ice categories are needed (Hibler, 1980; Flato and Hibler, 1995). An alternative approach is to solve the ice concentration and mass for each ice category or ice type in a Lagrangian ice thickness space (Bitz et al., 2001). Multi-category sea-ice models apply redistribution functions to describe an average evolution of the pack ice deformation processes. Several deformation processes, such as compacting, rafting and ridging, are possible during a single time step (Haapala et al., 2005). This mimics the real behavior of the pack ice on a continuum scale.

These three governing equations are strongly coupled. Firstly, sea-ice mobility is nonlinearly related to the ice thickness and concentration. Thin or low-concentration pack ice is drifting at approximately free drift speed; however, even if the ice is 0.5m thick, solid ice cover can be stationary under the action of strong winds in the Bay of Bothnia. In turn, ice motion generates fractures and leads on the ice pack which enhance mobility and, more importantly, increase the sea-ice mean thickness owing to new ice growth in leads and the formation of pressure ridges in compression. Due to ice dynamics, the mean ice thickness in the coastal boundary zone is thicker than in landfast ice regions (Oikkonen et al., 2017; Ronkainen et al., 2018). The mobility of pack ice has large consequences regarding the formation of coastal leads, which are local sources of heat and moisture.

In the Baltic Sea, the correct modeling of landfast ice, which can extend several kilometers from shore, is essential. In the CICE, the landfast ice regime is parameterized by introducing basal stress due to grounded ridges (Lemieux et al., 2016). A simplified approach is to assume that the landfast ice regime is dependent on sea depth (Palosuo, 1963). This approach has been used in several applications and was implemented in NEMO-Nordic (Pemberton et al., 2017 ).

A study of Baltic Sea climate variability based on the Hibler model type coupled with the Bryan–Semtner–Cox, z-coordinate baroclinic ocean model (the Rossby Centre Ocean, RCO, model), developed by Markus Meier (Meier 2002a, b; Döscher et al., 2002), was performed in the early 2000s. The RCO and the University of Helsinki sea-ice (HIM) models were later used for the analysis of the future ice conditions of the Baltic Sea region (Haapala et al., 2001). Based on these models, the authors carried out two 10-year simulations representing preindustrial and future scenarios of Baltic Sea ice conditions. On a global scale, both models delivered similar results; however, on a regional scale, there was large variation between model results. The abovementioned studies showed dramatic decreases in ice extent, and the calculated ice thickness was also lower in the scenario simulation; it is expected, based on this simulation, that the Archipelago Sea and Quark will not be covered by ice in the future. The influence of greenhouse gas emissions (the A2 and B2 IPCC scenarios that represent the climate of the late 21st century, 2071–2100) on ice conditions was examined using four RCOA (RCOA is a coupled RCO with the atmosphere model) model simulations. Each analyzed scenario was made using the same model but with different boundaries created by two different climate models. Results showed that the mean annual ice volume will decrease by about 80 % or more, which amounts to a reduction in the annual maximum ice thickness of up to 60 % and a decrease in ice days of over 90 %. All of the results are location dependent. These studies also concluded that total ice area depends on the air temperature, with less influence from other physical factors.

Recently two present-day coupled ice–ocean models have been developed for the Baltic Sea region: the ice part (LIM3.6) was evaluated in the NEMO-Nordic model (Pemberton et al., 2017), which covers the Baltic and North seas, and B-CESM (Jakacki and Meler, 2018), which only covers the Baltic Sea area, and the model is based on the Community Earth System Model, where sea ice is represented by CICE; and the oceanic part is the Parallel Ocean Program. Both are working well as present-day climate models.

2.2.5 Coupling strategies and pitfalls in comparing coupled and uncoupled models

In recent years, two different coupling architectures have been used and actively developed. On the one hand, there is a single executable concept that uses the coupler as a driver to call different Earth system components and handle the communication between them. Examples of this concept include the Earth System Modeling Framework (ESMF; Hill et al., 2004), CPL7 (Craig et al., 2017) or C-Coupler1 (Liu et al., 2014). This approach might require a substantial degree of adaptation of existing code to fit into the coupling framework. With respect to the performance of coupled systems, single executable design (CPL7) has been shown to be superior to multiple executable design (CPL6) for today's configuration of coupled GCMs (Craig et al., 2017). On the other hand, there is the concept of the OASIS and YAC coupler (Valcke et al., 2015; Hanke et al., 2016) that orchestrates the individual executables of the Earth system components via a communication library. This approach is less invasive for existing codes and requires only the insertion of communication calls without the need to fit into a common framework. The initial performance bottleneck in OASIS3 (Valcke et al., 2015) has been relaxed with the inclusion of the Model Coupling Toolkit (MCT; Larson et al., 2005; Jacob et al., 2005) in the OASIS3-MCT coupler (Valcke et al., 2015; Craig et al., 2017). MCT parallelizes the regridding between different Earth system components and the necessary communication.

The hierarchical approach with individual models as entities of a framework is typical for ESMs where the different components are developed within one institution. With MOSSCO (Lemmen et al., 2018), there is also an example of coupled model development in the Baltic Sea region that uses the ESMF to build a regional ESM with a focus on coastal processes in the North Sea and the Baltic Sea. The Modular Earth Submodel System (MESSy; Jöckel et al., 2010) originated from the need for atmospheric chemistry to be coupled to atmospheric dynamics and has evolved into a system with a coupled regional atmosphere component – COSMO (Kerkweg et al., 2018).

Most other efforts in coupled atmosphere–ocean modeling in the Baltic Sea region are based on the use of community models that are coupled with the community OASIS coupler (Table 1). The advantage is that specific development tasks can be distributed to different communities, and individual groups benefit from each other's expertise.

Traditionally, sequential coupling has been used and has the advantage that component models update the variables and fluxes with the most recent information from other coupled components. On the other hand, a sequential coupling is less efficient, as other components need to wait until the active component has updated its state. Nowadays, concurrent coupling is preferred where all coupled components run concurrently with state variables and fluxes that have been updated commonly during the last coupling time step. This implies that the coupling time step between coupled components needs to resolve the important processes that lead to feedbacks between, for example, atmosphere and ocean. A typical example would be the moisture flux in a RCM to investigate monsoon dynamics (Yang et al., 2019), the generation of medicanes (Akhtar et al., 2014) or coastal upwelling (Perlin et al., 2007).

A long known issue with coupled modeling is the mismatch of land–sea masks between atmosphere and ocean components of a coupled system (Jones, 1999). The Baltic Sea region is especially prone to this issue due to its complex coastlines and the currently relatively low resolution of 0.11–0.22^∘ in the atmosphere models and $\mathrm{1}/\mathrm{60}$–$\mathrm{1}/\mathrm{20}$^∘ in the ocean models. With the higher resolution of convection-permitting regional climate models (CPRCMs) approaching the resolution of Baltic Sea models (Belusic et al., 2019), the effects might become less prominent, but the conceptual issue remains. One way to deal with the ambiguous assignment of grid cells is remapping, which can be carried out on an exchange grid that is the joint set of grid cells between two participating grids (Balaji et al., 2006; Bauer et al., 2020). The exchange grid also addresses the vertical flux between different model components but has not been used in RCMs of the Baltic Sea region to date.

Comparing coupled and uncoupled ocean model runs usually involves simulations that use different bulk formulae to calculate air–sea fluxes in the atmospheric boundary layer in the coupled and uncoupled modes, respectively (e.g., Dieterich et al., 2019a). To separate the effect of coupling from other differences between coupled and uncoupled models, this needs to be addressed in more detail in future efforts. Nevertheless, many studies have aimed at demonstrating improvements in simulating aspects of the present-day climate. However, a general statement about model performance of coupled versus stand-alone models can not be drawn, as many studies have shown the added value of coupling especially in cases when the uncoupled model was driven by low-quality forcing data (Tian et al., 2013; Gröger et al., 2015; Ho-Hagemann, 2017). Furthermore, the uncoupled models often had a lower resolution at the air–sea boundary. In turn, this implies that the application of coupled models, in particular in future climate projections as global climate models, can have large biases and that the constraints of the governing model that drive the uncoupled model may be too tight (Mathis et al., 2018).

2.3 Ocean–wave–atmosphere coupling

One major characteristic of the atmosphere–ocean interface is the presence of surface gravity waves – the surface changes as a direct response to the atmospheric forcing. Surface gravity waves (hereafter referred to as waves) are mainly generated by the wind; thus, the wave field is strongly dependent on the wind field. For coastal areas, the over-water fetch is one additional aspect. Waves are characterized by a variety of properties (e.g., wave height, wave age and wave steepness). The transport of momentum is the key property of air–sea interactions affected by waves, but waves also influence the exchange of heat and mass as well as the turbulence of the lower atmosphere and the upper ocean. Indirectly, wave influences can affect the whole boundary layer/mixed layer in the atmosphere as well as the ocean.

As a buffer role at the air–sea interface, waves can be divided into growing waves (young waves) and decaying waves (swell) with very different impacts on the atmosphere and the ocean. Waves extract energy from the air-side stress when waves are growing. By contrast, they release momentum to the ocean in the presence of decaying waves. When considering the role of waves, the stress balance at the air–sea interface is expressed as follows (ECMWF, 2017):

$$\begin{array}{}\text{(1)}& {\mathit{\tau }}_{\mathrm{a}}={\mathit{\tau }}_{\text{oc}}+{\mathit{\tau }}_{\mathrm{w}}+{\mathit{\tau }}_{\text{ds}},\end{array}$$

where τ_a is the air-side stress, τ_oc is the ocean-side stress, τ_w is the wave-induced stress and τ_ds is the momentum flux from the wave field to mean currents. In traditional stand-alone or atmosphere–ocean coupled models, τ_a is identical to τ_oc without considering the role of waves. It has, however, been shown that the normalized ocean-side stress, ${\mathit{\tau }}_{\text{oc}}/{\mathit{\tau }}_{\mathrm{a}}$, can be larger (smaller) than 1.5 (0.85) in the Baltic and the North seas (Alari et al., 2016; Wu et al., 2019b; Qiao et al., 2021).

For high-wind-speed conditions, sea spray and airflow separation caused by breaking waves are important aspects. This causes the drag coefficient to level off under extreme winds (about 25–30 m s⁻¹ for the mean wind speed at 10 m). For swell conditions, the air–wave interaction is more complex, and waves influence the momentum flux, the turbulence of the atmosphere and the mixing (Nilsson et al., 2012; Wu et al., 2017a).

Wave breaking can transfer a significant amount of energy flux into the ocean (Melville et al., 2002), which mainly affects down to few significant wave height depths from the ocean surface (D'Asaro, 2014). The breaking-wave-induced energy flux is commonly parameterized as an additional input of turbulent kinetic energy (TKE), α_CBu³, in the surface boundary (Craig and Banner, 1994), where u is the ocean-side friction velocity, and α_CB is a wave-related parameter that is usually treated as a constant (e.g., 100) without considering the impact of wave status.

Due to their periodic motion, waves induce a net drift in the wave propagation direction: $u_{\mathrm{s}}=u_{\mathrm{L}}-u$, which is the difference between the Eulerian velocity, u, and the Lagrangian velocity, u_L, defined as Stokes drift (Stokes, 1880). The Stokes drift can impact the wave-averaged momentum equation through the Coriolis–Stokes force (CSF), the vortex force and a Stokes-corrected pressure (e.g., Suzuki and Fox-Kemper, 2016):

The CSF can alter the ocean Ekman transport, which is indicated by numerical simulations (Polton et al., 2005). The Langmuir turbulence (LT) induced by the vortex force is one of the most important wave-related processes, which can indirectly affect the whole mixed layer through the turbulent transport of the wave-induced turbulence (Belcher et al., 2012). Many turbulent closure schemes have been modified in order to implement the LT influences (Ali et al., 2019), such as the K-profile parameterization turbulence scheme (McWilliams and Sullivan, 2000) and the k−ϵ turbulence model (Axell, 2002). The Stokes-corrected pressure is small in the order of the Rossby number and can be neglected in coarse-resolution models. In addition, the Stokes drift can also affect the mass and tracer transport through the divergence of the sea surface height and tracer advection, respectively (McWilliams and Sullivan, 2000; Wu et al., 2019b).

Non-breaking-wave-induced ocean mixing proposed by Qiao et al. (2004) and Babanin and Haus (2009) is a direct stirring function in the ocean by waves. In coastal areas, the wave-induced bottom stress can also affect the circulation and sea water level (Davies and Lawrence, 1995).

In marginal ice zones, short waves are attenuated rapidly by ice, and long waves can propagate much further into the ice areas. Wave radiation stress on the ice is a term usually neglected in the ice model. The mechanism of the wave–ice interaction is complex; for a more detailed description of the wave–ice interaction as well as its parameterization in models, the reader is referred to Zhao et al. (2015) and Squire (2018).

2.3.1 Effects on the coupled system when introducing waves

To capture the wave–current and wave–atmosphere interaction, several models coupled with a wave model have been developed in the Baltic Sea and the North Sea region. The wave model has been coupled to regional climate models, i.e., RCA-WAM (Rutgersson et al., 2012; Wu et al., 2015) and COSMO-WAM (Wahle et al., 2017), and weather prediction models, i.e., WRF-WAM (Wu et al., 2017b) and WRF-SWAN (Larsén et al., 2019). On the ocean side, several wave–ocean coupled models have been tested: one-way coupled models, such as NEMO-WAM (Alari et al., 2016; Wu et al., 2019b; Staneva et al., 2017), and two-way coupled models, such as GETM-WAM (Staneva et al., 2016). Recently, a three-way coupled atmosphere–wave–ocean system (Uppsala University-Coupled Model, UU-CM) was developed and tested in the Baltic Sea region (Wu et al., 2019a).

Storm simulations

The wave feedback on the atmospheric simulation is significant for young waves, as they extract more momentum for their growth (Janssen and Viterbo, 1996). Accordingly, wave coupling is more important under extreme conditions (when young waves are dominant). Based on the simulation of 23-year storms using the WRF-SWAN model with a domain covering the North and Baltic seas, Larsén et al. (2019) found that the wave influence improves the model performance compared with the uncoupled model (mainly under a high-wind range).

Sea spray under high winds can significantly impact the momentum and heat flux at the air–sea interface. An effective roughness length parameterization considering sea spray influences was developed by Wu et al. (2015) based on the studies of Kudryavtsev and Makin (2011) and Zhao et al. (2006). Along with the heat flux parameterization considering sea spray (Andreas et al., 2015), the parameterizations were implemented into the RCA-WAM coupled model for storm simulations (Wu et al., 2015). The results showed that the sea spray impact on the wind stress can significantly improve the wind simulation during storms; in addition, the sea spray impact on the heat flux improves the temperature simulation.

The storm-induced surge has been investigated using ocean–wave coupled systems. Wave–current interaction processes were implemented in the GETM-WAM model, which included radiation stress, Stokes drift, bottom friction modification and turbulent kinetic energy due to wave friction (Staneva et al., 2016). Staneva et al. (2016) found that the coupled system did not have a significant influence over the open North Sea; however, it significantly improved the simulation in the coastal areas in terms of significant wave height, water level and current. In the recent study by Staneva et al. (2017), the CSF, sea-state-dependent momentum and energy flux were implemented in the NEMO model with external wave forcing data. The simulation of two storms in the North Sea indicated that wave-related processes improve the storm surge simulation under extreme conditions. Compared with the other wave-related processes investigated in their study, the sea-state-dependent momentum flux played a dominant role, which agrees with the results of Wu et al. (2019b).

Mesoscale features

Coastal upwelling occurs frequently in the Baltic Sea during the summer months and is mainly induced by the divergence (convergence) of the wind stress (Lehmann and Myrberg, 2008). Ocean waves affect the coastal upwelling by altering the ocean-side stress (Eq. 1) and upper-ocean mixing (mainly through CSF and the breaking-wave-induced energy flux). Based on an offline wave–ocean coupled model, Alari et al. (2016) found that the ocean model better captures the distribution of the sea surface cold water induced by a coastal upwelling event on the eastern coast of the Baltic Proper. In their study, three wave-related processes were implemented in the ocean model: CSF, sea-state-dependent energy and momentum flux. With a similar ocean model set up in the Baltic Sea and the North Sea area, Wu et al. (2019b) investigated the wave impact on the Baltic Sea coastal upwelling in terms of intensity and frequency from June to September in 2015. In addition to the three processes investigated in Alari et al. (2016), the Stokes drift impact on the mass and tracer advection was also included in their model. They found that the Stokes drift impact on the mass and tracer advection largely canceled the influence of the CSF on the coastal upwelling simulation. Compared with the control experiment (without wave-related processes), the simulation including the four wave-related processes changes the upwelling frequency by up to 10 %, although this varies with location (Fig. 7). Furthermore, the combined effect of these wave-related processes increases the weak upwelling frequency ($-\mathrm{4}<\mathrm{\Delta }T<\mathrm{2.5}$ ^∘C) but decreases the strong upwelling frequency ($\mathrm{\Delta }T<-\mathrm{6}$ ^∘C); here, ΔT is the SST difference from the zonal mean temperature (Wu et al., 2019b).

Figure 7 Panels (a)–(d) show the upwelling frequency during the period from June to September 2015, and panels (e)–(h) show the difference between the simulation considering wave influence and the control experiment. (This figure is adapted from Wu et al., 2019b.)

Two convective snow band cases in the Baltic Sea region were investigated using the RCA-WAM model by Jeworrek et al. (2017). Due to the wave feedback on the atmosphere through changing the sea surface roughness length, the RCA-WAM model processes a time shift of several hours in the maximum 10 m wind simulation. In the simulation of the two convective snow cases, the wave coupling has less influence than that from the ocean coupling (see above).

Climate simulations

The swell impact on the atmospheric mixing was introduced into the RCA-WAM regional climate model by adding an extra wave-age-dependent coefficient to the mixing length (Rutgersson et al., 2012). Simulations show that the swell impact on the atmospheric mixing can alter the mean surface wind by up to 0.3 m s⁻¹ in the Baltic Sea and the North Sea (Rutgersson et al., 2012; Wu et al., 2016). The magnitude and direction of the change depend on the wave field and environmental conditions (Rutgersson et al., 2012; Wu et al., 2016).

Based on a 3-summer-month simulation, Alari et al. (2016) found that the wave-related processes (sea-state-dependent momentum flux, TKE induced by wave breaking, and Stokes–Coriolis forcing) can change the sea surface temperature by up to 1^∘ in the Baltic Sea; here, the sea-state-dependent momentum and TKE flux play a warming role, but the Stokes–Coriolis forcing plays a cooling role along coastlines. Similar to the sea surface temperature, the Stokes–Coriolis forcing dominates the bottom temperature cooling along coastlines where the water depth is relatively shallow. However, in their study, the LT influence was not included in the model.

Recently, Wu et al. (2019a) developed a fully coupled atmosphere–wave–ocean model with an improved representation of the air–wave–sea interaction processes. Based on a 2-month-long (January and July 2015) simulation, they found that the coupling has a significant impact on coastal areas in the Baltic Sea. The wave–current interaction has a larger influence in the summer months than in the winter months (Wu et al., 2019a) because the wind speed is higher in winter, and the sea-state-dependent momentum/energy flux and the Stokes-drift-related processes are more important for upper-ocean mixing.

2.4 Hydrological coupling – closing the water cycle

In coupled system modeling, hydrological coupling usually refers to the closure of the water cycle between land and ocean via the interactive simulation of river runoff. River runoff is an important component of global and regional water cycles, especially in the Baltic Sea catchment where it comprises about half of the precipitation over land areas (Lind and Kjellström, 2009) and about 2 times the precipitation over sea areas (Jacob, 2001; Leppäranta and Myrberg, 2009). Freshwater inflows from runoff and precipitation affect the thermohaline circulation (Knudsen, 1900; Placke et al., 2021). Decadal variations in Baltic Sea salinity are largely caused by the accumulated runoff to the water body (Meier and Kauker, 2003; Väli et al., 2013; Radtke et al., 2020). On the one hand, the thermohaline circulation of the Baltic Sea is also influenced by inflows of highly saline water from the North Sea (that itself may be strongly impacted by precipitation and river runoff; Lehmann and Hinrichsen, 2000). On the other hand, river runoff into and net precipitation over the Baltic Sea mainly induce its outflow into the North Sea where it is an important source of stratification in the northwestern European shelf (Hordoir and Meier, 2010).

In addition, river runoff and the associated nutrient loads substantially influence the functioning of the marine ecosystem. These inflows from land, carrying fresh, nutrient-rich water, determine coastal physical conditions and the nutrient concentration. Therefore, they dominantly influence primary production and affect the variability of the whole ecosystem (e.g., Gustafsson et al., 2012; Daewel and Schrum, 2017). In the Baltic Sea, this becomes even more relevant, as the water body is almost decoupled from the open ocean so that terrestrial nutrients significantly contribute to determining ecosystem productivity (Thurow, 1997; Österblom et al., 2007).

Consequently, river runoff is also an important component for coupled system modeling over the Baltic Sea region. Hagemann et al. (2020) provide an overview of the current state of high-resolution modeling of river runoff (or discharge) within the framework of regional coupled system models (RCSMs). Note that Sect. 2.4 of this paper comprises excerpts and, hence, some overlap with the introduction section of Hagemann et al. (2020), as the associated information is relevant to the present review. In addition to traditional regional climate models (RCMs), RCSMs have recently been developed to conduct climate change studies at high spatial and temporal resolutions. In order to adequately represent biogeochemical cycles, a proper description of the transport of chemical species, such as nitrogen, phosphorus, carbon and silicon, into the ocean requires a very detailed representation of stream characteristics (such as flow paths, lakes, ponds, reservoirs, wetlands and floodplains) because the smallest water bodies may exhibit large parts of the retention on land (Bouwman et al., 2013). Therefore, RCSMs require a high-resolution discharge component to couple their atmosphere/land components to the ocean component and to adequately resolve smaller catchments and the day-to-day variability of discharge.

2.4.1 Coupling to coarse-resolution discharge models

The current discharge models applied in coupled (or Earth system) models for global or regional climate simulations usually do not fulfill this requirement. In global ESMs, discharge (or routing) models are frequently part of the coupled system (often as part of the land surface scheme), but their spatial resolution is usually 0.5^∘ (Roeckner et al., 2003; Guimberteau et al., 2012) or coarser (Lawrence et al., 2011; Milly et al., 2014; Best et al., 2011).

Instead of using hydrological coupling, many RCSMs use prescribed runoff, taken from climatology, observations or model data. For example, Gröger et al. (2015) coupled the Rossby Centre regional Atmospheric climate model RCA4 on a 24 km European domain with the NEMO ocean model over the Baltic Sea and North Sea, but river runoff was prescribed as a daily climatology of an E-HYPE ERA-Interim hindcast (Donnelly et al., 2016). To reflect the projected increase in precipitation in the northern part of the Baltic Sea (e.g., Donnelly et al., 2016), a linearly increasing discharge to +10 % in 2100 is used in the Bothnian Sea and Bothnian Bay. Thus, only a few RCSM setups exist where a discharge model is included, but its resolution is rarely higher than 0.5^∘; examples of such setups are the Hydrological Discharge (HD) model (Hagemann and Dümenil, 1998) at 0.5^∘ (Sein et al., 2015; Sitz et al., 2017; Elizalde, 2011), TRIP (Oki and Sud, 1998) at 0.5^∘ (Dell'Aquila et al., 2012; Sevault et al., 2014) and LARSIM (Bremicker, 2000) at $\mathrm{1}/\mathrm{6}$^∘ over northern Europe (Lorenz and Jacob, 2014).

The latter model is part of the fully coupled RCSM BALTIMOS for the Baltic Sea and links the atmospheric RCM REMO to the Baltic ocean/sea-ice model BSIOM. Lorenz and Jacob (2014) introduced the BALTIMOS system and showed first results from a simulation driven by analysis data from the European Centre for Medium Range Weather Forecast (ECMWF). Sein et al. (2015, 2020) introduced the RCSM ROM in which the global ocean–sea-ice–marine biogeochemistry model MPIOM/HAMOCC with regionally high horizontal resolution is coupled to the atmospheric RCM REMO and the global HD model (see above) via the OASIS coupler. They used forcing from the NCEP/NCAR reanalysis and ECHAM5/MPIOM (Roeckner et al., 2003; Jungclaus et al., 2006) historical simulations and evaluated their results for the North Atlantic and North European region, where they also specifically addressed the Baltic Sea catchment.

2.4.2 Coupling to high-resolution discharge models

Several studies exist where a RCM was coupled to a very high-resolution regional hydrology model that covers the full range of processes at the land surface. As such coupled systems are often limited by computational effort, these studies currently only cover short periods or relatively small catchments/areas (e.g., Mauser and Bach, 2009; Senatore et al., 2015; Shrestha et al., 2014). Larsen et al. (2014) presented results from a full two-way coupling of the HIRHAM RCM over a 4000 km×2800 km domain at 11 km resolution and the combined MIKE SHE-SWET hydrology and land surface models over the Danish Skjern River catchment (area of 2500 km²).

Over Korea, the discharge model TRIP has been coupled to a RCM at 0.5, 0.25 and 0.125^∘ in preparation for future RCSM studies (Lee et al., 2015). Nguyen-Quang et al. (2018) used a high-resolution (1 km) river network to define hydrological transport units within the grid boxes of the ORCHIDEE land surface model and applied this over the Mediterranean region using forcing data at 0.5 and 0.25^∘ resolution. To our knowledge, none of the hydrological models used in the studies listed above or in the study of Hagemann et al. (2020) have been used in a fully coupled RCSM setup that can be applied for climate timescales and large-scale areas.

Very recently, Hagemann et al. (2020) developed a high-resolution version of the HD model that is globally applicable at a 5 min resolution (HD5 model). In their study, offline HD5 model results were evaluated over Europe and the Baltic Sea catchment. In order to prepare high-resolution scenario simulations over Europe and the Baltic Sea, the HD5 model has already been coupled within GCOAST (Geesthacht Coupled cOAstal model SysTem), which is the RCSM developed at Helmholtz-Zentrum Hereon, Geesthacht. This coupling is necessary to ensure that the simulated discharge and other hydrological variables are consistent with the climate variables so that interactions between the different compartments of the regional Earth system can be considered. Ho-Hagemann et al. (2020) introduced a GCOAST subset in which OASIS3-MCT couples the atmospheric RCM CCLM over the 0.11^∘ EURO-CORDEX domain, the HD5 model over Europe, and the ocean–sea-ice model NEMO-LIM3 over the Baltic Sea, North Sea and parts of the North Atlantic. Using this GCOAST-AHOI subset, they investigated the effects of air–sea coupling on internal variability of the regional atmospheric model.

2.4.3 Hydrology models in future climate scenarios

The aforementioned sensitivity of the Baltic Sea to yearly freshwater input implies that any coupled hydrology model needs to be very accurate in simulating the total yearly discharge over the catchment area. Consequently, the use of regionally coupled hydrology models in future climate scenarios is challenging because the water cycle in coupled general circulation models (GCMs) is often not closed (Liepert and Lo, 2013). Moreover, GCMs and RCMs suffer from substantial biases, especially with respect to precipitation and the hydrological cycle (Flato et al., 2013; Kotlarski et al., 2014), which affect the simulated discharge. An example of this is the regionally coupled atmosphere–ocean–sea-ice–marine biogeochemistry model ROM (Sein et al., 2015) which overestimates catchment precipitation by 11 %–14 %, translating into a 23 %–34 % overestimation in discharge when driven with different reanalysis datasets. Thus, current studies using future scenarios for the Baltic Sea mostly rely on rough estimates of a 10 % increase in discharge in Bothnian Bay and the Bothnian Sea by 2100 (Dieterich et al., 2019b; Gröger et al., 2019), or they use discharge data from uncoupled hydrology models (Saraiva et al., 2019a, b; Meier et al., 2018, 2019a, b). By contrast, the use of the online-coupled HD component is physically more consistent, as the water cycle is closed but may transfer biases from the atmosphere/land components into the ocean, leading, for example, to an extraordinarily strong freshening in future scenarios (Sein et al., 2020).

3 Uncertainties and gaps

3.1 Land–atmosphere

Previous works have shown the potential for substantial land–atmosphere coupling in the Baltic Sea region via land-use change, natural vegetation dynamics and land management. RESMs that downscale CMIP6 (Coupled Model Intercomparison Project Phase 6) outputs to the Baltic Sea region should, as a first step, use land-use/management scenarios consistent with CMIP6 protocols, i.e., from the LUH2 dataset (Hurtt et al., 2020). Simulations with RCA-GUESS, RCAO-GUESS (Zhang et al., 2020) and other RESMs have shown the importance of and potential for land–atmosphere–ocean interactions in the Baltic Sea region. Coupling new Earth system components is expected to further advance our understanding of climate and environmental change in the Baltic Sea region. These include land–ocean/freshwater interactions, including the leaching of carbon and nutrients from land, the use of simulated BVOCs (Kulmala et al., 2014; Hantson et al., 2017) and aerosols precursors from natural and managed vegetation and wildfires, and CH₄ and N₂O emissions from wetlands and agriculture.

3.2 Atmosphere–ocean–sea ice

With respect to coupled atmosphere–ocean modeling, the choice of the size of the coupled domain always has to be carefully considered, and in most cases, the economical cost of driving the model sets an upper limit. When the atmosphere domain covers several seas, as it does in the EURO-CORDEX domain, it might be feasible to coupled only a few of them to save on computational costs (e.g., Tian et al., 2013; Gröger et al., 2015). However, with respect to regional ocean studies, it appears reasonable to avoid lateral boundaries too close to the shelf break: important small-scale processes take place in this region that determine the shelf–open-ocean exchange, as demonstrated for the North Sea (Gröger et al., 2013), and benefit from high resolution. On the other hand, including slower components with longer internal timescales like parts of the open North Atlantic will allow the regional model to generate its own long time internal variability, which is likely out of phase with the decadal variation of the driving global model. One example is provided by Sein et al. (2014), who used simulations for the Arctic Ocean to emphasize that the inclusion of the North Pacific in coupled domain destroys the interannual correlation with reanalysis data. Hence, special nudging techniques might be necessary to adapt the regional ocean model large-scale circulation to that of the driving global ocean model, as this is practiced for atmospheric regional models (von Storch et al., 2000).

Along the coastal zone, small-scale oceanic processes like upwelling can create strong SST gradients along the coasts, further influencing land–sea–atmosphere dynamics. However, SST patterns are still strongly smoothed while communicated to the atmosphere, as the SST field is interpolated onto the atmosphere grid. Thus, more advancements in coastal zone dynamics can be expected from higher-resolution atmosphere models (up to a few kilometers) and the transition to convection-permitting models instead of hydrostatic models. With respect to future scenario climate simulations, more research is necessary to assess whether or not the coupled and uncoupled models reveal different climate change signals in scenario simulations in both the hydrographic and atmospheric properties of the Baltic Sea region. First attempts to address this question have already been made (e.g., Gröger et al., 2021, Bøssing Christensen et al., 2021; Meier et al., 2021).

3.3 Waves

Adding surface waves to a coupled system is becoming more important with increasing resolution, in particular, when detailed information is required in complex areas (such as for offshore wind energy applications in the coastal zone). It has also been shown that wave information improves the description of ocean mixing. There is, however, still limited knowledge concerning the impact of waves on the exchange of heat and mass. The interaction between waves and ice is also not well described in present models, which induces uncertainties in the freezing of ice as well as the wave properties in the marginal ice zone. This is important in polar regions but also potentially in the northern Baltic Sea.

For extreme events (often linked to local–mesoscale systems), high-resolution models are crucial, and the introduction of Earth system information that is as accurate as possible is essential. For this, coupling to other high-resolution models is a prerequisite. This primarily includes atmosphere, ocean and waves, but it can be expected to also be beneficial with submodels for parameters such as ice and aerosols. Thus, improvements from including wave models into regional Earth system models can be expected, in particular for features like convective precipitation, cold-air outbreaks, lake snowfall, polar lows and midlatitude storms.

3.4 Hydrology

Socioeconomic development is a major driver for nutrient loads into the Baltic Sea (Arheimer et al., 2012; Bartosova et al., 2019) that contribute to eutrophication and hypoxic conditions (e.g., Saraiva et al., 2019a, b). So far, only a few hydrological models include nutrient cycling to provide explicit estimates of nutrient inputs to the sea (Hundecha et al., 2016). This needs to be addressed in more detail in future model development. Hence, a major challenge is to include/combine (1) land-use change, (2) terrestrial carbon and nutrient cycling from fertilizers, (3) dynamic vegetation modeling and (4) long-term storage of nutrients in the soils.

The potential increase in weather extremes will have an effect on soil erosion and nutrient loads, which is a topic that has not yet been implemented in hydrological models. This is especially important in the context of future climate change, as some extreme conditions may become more frequent (e.g., Jacob et al., 2014). For the adequate representation of biogeochemical cycles, the hydrological model must not only consider the respective lateral transports, but it must also include detailed biogeochemical process descriptions (e.g., Tang et al., 2018).

However, even with a coupled hydrological component, human impacts and their future developments are currently not regarded in coupled models of the Baltic Sea region. Here, many rivers are strongly affected by human impacts (e.g., water abstraction for irrigation) and regulation (e.g., dams). Consequently, for the modeling of rivers and watersheds that are highly influenced by human activities, related processes need to be implemented into the respective hydrological model component. As pointed out by Hagemann et al. (2020), this includes the implementation of existing and planned dams and reservoirs (based on available global databases), their management of river flow regulations and modules to simulate water withdrawals (e.g., for irrigation). Apart from prescribed scenarios, the future development of these impacts might be realized with simple economical modules.

3.5 Internal variability

Hydrodynamical Earth system components (i.e., ocean and atmosphere primitive equation models) are characterized by the fact that they generate noise as a result of turbulent dynamics (e.g., Weisse et al., 2000; Penduff et al., 2018; Wiese et al., 2020; Geyer et al., 2021). In climate applications, where long-term averages and statistics are usually considered, this is no problem, provided that the averaging period largely exceeds the frequencies of the internal variability of the system. It is also not problematic if short-term events are analyzed, provided that the statistics of these events are only interpreted over an appropriately long period. Contrary to this, i.e., when single events are analyzed, such as the timing of certain storms (e.g., Ho-Hagemann et al., 2017, 2020) or wind-induced coastal upwelling events in the Baltic Sea, the model solution can be substantially modulated by internal variability, and the comparison of coupled versus uncoupled systems may be misleading. Thus, this requires a profound estimation of the robustness of the model solution with respect to the initialization, boundary conditions (Wiese et al., 2020; Ho-Hagemann et al., 2020) and even the computing platform (Geyer et al., 2021). For this, different methods have been developed that are mainly based upon the generation of large ensembles with the consideration of small perturbations in the initialization (e.g., Giorgi et al., 2000; Weisse et al., 2000, 2003; Sieck and Jacob, 2016).

Comparing ensembles of coupled and uncoupled regional ocean–atmosphere models, Ho-Hagemann et al. (2020) recently demonstrated that interactive air–sea coupling can substantially reduce internal model variability compared with uncoupled atmosphere models. In case of Cyclone Christian, which occurred between 27 and 29 October 2013 in northern Europe, the authors found that the larger uncertainty in the atmosphere-only simulation was caused by a combination of two factors: (1) uncertainty in the parameterization of the cloud–radiation interaction in the atmospheric model and (2) the lack of an active two-way air–sea interaction.

Using a similar approach to Ho-Hagemann et al. (2020), Wiese et al. (2020) also found a reduction in the internal model variability in ensembles of an interactive atmosphere–wave coupled model compared with those of a stand-alone atmospheric model. The role of internal variability has still not been sufficiently investigated in more complex coupled systems involving more components such as fully coupled ocean–wave–atmosphere models, as pointed out by Wiese et al. (2020), or coupled ocean–atmosphere–land vegetation models.

3.6 Other components

For the current generation of RCMs, a substantially lower projected 21st century warming has been demonstrated compared with their driving global models, and the neglection of scenarios for time-varying aerosols has been identified as a key process (Boé et al., 2020). However, the use of explicit atmospheric chemistry and transport models as interactive parts within global Earth system models is still not common and is mostly parameterized in CMIP6 models (Stevens et al., 2017; Fiedler et al., 2019). Thus, the representation of aerosols is a major challenge in RCMs and has been recognized within the CORDEX Flagship Pilot Studies (FPS-Aerosol https://www.hymex.org/cordexfps-aerosol/wiki/doku.php?id=start, last access: 10 September 2021).

Despite the advancements in including the biophysical feedback from the land vegetation (as outlined in Sect. 2.1), analogous feedbacks from the marine biogeochemical cycles on the physics are currently not implemented in RCSMs. The main feedback is due to the altered penetration of solar insolation by marine biota that further influences heat absorption and, thus, the vertical distribution of heat (Lengaigne et al., 2009). Regional studies for the Arctic (Lengaigne et al., 2009) and for the Indian Ocean (Sein et al., 2021) suggest that intense phytoplankton blooms strongly influence vertical mixing and thermocline dynamics. The latter two processes are essential for the Baltic Sea, especially in the context of projected increases in cyanobacteria blooms under a future warmer climate (e.g., Saraiva et al., 2019a, b; Meier et al., 2019a, b).

4 Conclusions and key messages

For the Baltic Sea region, previous research has identified a number of important feedback loops between Earth system compartments that alter both the mean climate and extreme event statistics such as heavy precipitation, storms and flooding. Regional Earth system modeling is the only way to represent feedbacks between different Earth environmental compartments in an adequate way (Heinze et al., 2019).

The coupling of atmosphere models to dynamical components for the ocean and land constitutes a major step towards Earth system modeling of the Baltic Sea region. Climate projections for the Baltic Sea region up to the end of the century reveal that the feedback of vegetation changes on climate warming is mostly negative in northern and central Europe, whereas an amplification of warming of up to 1 K is obtained for southern Europe (Wramneby et al., 2010). Likewise, the implementation of scenarios for land management in the Baltic Sea region can alter the temperature locally by up to a few kelvin (Strandberg and Kjellström, 2019). The alterations emerge from a number of coupled land–atmosphere processes involving biophysical feedbacks such as changes in albedo and roughness length, and evapotranspiration and related changes in the LAI and vegetation-type as well as biochemical feedbacks from high carbon dioxide concentrations.

Unlike the land surface where evaporation also depends on precipitation with a further effect on heat exchange, freshwater supply to the ocean (via runoff and precipitation) has no direct influence on evaporation over sea. Accordingly, coupled atmosphere–ocean model studies have mainly identified thermal air–sea coupling as most significant (e.g., Kjellström et al., 2005; Gröger et al., 2015, 2021; Ho-Hagemann, 2017), with the predominant impact on simulated SSTs. For the Baltic Sea, deviations between prescribed SST and sea-ice fields and their modeled quantities can add up to a few kelvin (Tian et al., 2013; Gröger et al., 2015). In turn, high-quality SST fields appear to be a prerequisite for the representation of extreme events over land, such as convective snow bands, heavy precipitation, or flooding due to the influence of SST and the presence/absence of sea ice on the large-scale atmospheric circulation. Many studies have demonstrated the importance of interactive wave coupling in mediating the atmosphere–ocean exchange of heat, momentum and mass which enables one to more reliably resolve processes such as the effect of breaking waves otherwise parameterized in either ocean or atmosphere models.

The first fully coupled atmosphere–ocean–hydrology models for the Baltic Sea region have been accomplished (Sein et al., 2015, 2020; Hagemann et al., 2020) and are an important step towards the closure of the water cycle. However, they also elucidate common problems related to strong biases in input precipitation/evaporation. Without bias correction, this leads to unrealistic river runoff to the Baltic Sea (Sein et al., 2020). This is particularly problematic for the Baltic Sea, as runoff constitutes a major contribution to the halocline structure and the freshwater budget of this water body.

Many studies have aimed at demonstrating the added value of interactive coupling by direct comparisons of single variables between coupled and uncoupled models. However, this may be misleading, as these improvements simply reflect the poor quality of prescribed boundaries used in the stand-alone model in many cases (see, e.g., Gröger et al., 2015). This has implications for downscaling the effect of future climate changes because global models can be strongly biased on the regional scale and, thus, provide low-quality input data (the so called “rubbish in rubbish out problem”, e.g., Hall, 2014). Here, the more complex regional coupled models can develop more independently from the biased parent global model by generating their own climate, as demonstrated by Mathis et al. (2018).

Finally, coupled models can be computationally highly demanding, especially when multiple components are included. For hindcast simulations of the historical climate, the more economical stand-alone models are a good tool, provided that forcing and boundary data are of good quality. For future climate simulations where this is not the case, coupled models will be the first choice. The higher costs of coupled climate simulations will require a reduction in the size of individual ensemble simulations, which hampers robust assessment uncertainties due to global models and scenarios. However, advanced techniques for reducing the ensemble size by conserving the model ensemble spread are currently under development (e.g., Wilcke and Bärring, 2016).