Coupled regional Earth system modelling in the Baltic Sea region

Non-linear 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 traditional standalone atmosphere and ocean models on which many of today's regional climate assessments rely on (e.g. EuroCordex, NOSCCA, BACC II). This promotes 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 into the region by using prescribed boundary conditions, and thus, permit a higher degree of freedom. In the past, several coupled model systems have been developed for Europe and the Baltic Sea region. This article reviews recent progress of model systems that allow two way communication between atmosphere and ocean models, models for the land surface including the terrestrial biosphere, as well as wave models at the air sea interface and hydrology models for water cycle closure. However, several processes that have so far mostly been realized by one way coupling such as marine biogeochemistry, nutrient cycling and atmospheric chemistry (e.g. aerosols) are not considered here. Compared to uncoupled standalone models, coupled earth system models models can modify mean near surface air temperatures locally up to several degrees compared to their standalone atmospheric counterparts using prescribed surface boundary conditions. Over open ocean areas, the representation of small scale oceanic processes such as vertical mixing, and sea ice dynamics appear essential to accurately resolve the air sea heat exchange in the Baltic Sea region and can only be provided by online coupled high resolution ocean models. In addition, the coupling of wave models at the oceanatmosphere interface allows 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 like e.g. heavy precipitation, storms, or tropical nights, and appear to be sensitive to coupling.

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 like e.g. heavy precipitation, storms, or tropical nights, and appear to be sensitive to coupling.
For the present-day climate, many coupled atmosphere-ocean and atmosphere-land surface models demonstrate added value with respect to single climate variables in particular when low quality boundary data were used in the respective standalone 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 for closing the water cycle remains problematic as the accuracy of precipitation provided by the atmosphere models is in most cases insufficient to realistically simulate the runoff to the Baltic Sea without bias adjustments.
Many regional standalone ocean and atmosphere models are tuned to well represent present day climatologies rather than accurately simulate climate change. More research is necessary about 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 to uncoupled ones.

Introduction
Climatic and environmental changes on regional scales are traditionally investigated using standalone models resolving processes specific for only one single environmental compartment (e.g. the terrestrial and marine biospheres, the hydrosphere, the atmosphere, the ocean, and the cryosphere, etc.). 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 (BACC II Author Team, 2015, Bøssing-Christensen et al., 2015 or the Euro-Cordex region  fall into the category of stand-alone ocean models or stand-alone atmosphere models whereas fewer assessments are based on coupled systems Meier et al., 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.
Climate change on regional scales can be much stronger than one would expect from external forcing such as greenhouse gases, solar radiation, etc. 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 run at typically coarse grid resolution. Recognizing these issues, there is an emerging demand for regional ESMs suitable for application at higher grid resolutions over a regional domain, accounting for biophysical coupling between the atmosphere and surface. National -e.g. MERGE: http://www.merge.lu.se/about-merge -and international strategic programmes are now stepping up their efforts to meet this demand. The Baltic Earth program (www.baltic.earth), 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).

Study area
Compared to 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 (Figure 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 while 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. Altogether, this makes Europe's climate fairly variable on small spatial scales and sensitive to perturbations in the large-scale atmospheric circulation patterns as e.g. reflected in 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 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 temperate humid climate in the southwest and 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) or 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.

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 degree of freedom. Consequently, coupled models can drift from observed conditions which makes their tuning more difficult but allows a more realistic interaction between models. Therefore, in their hindcast modes, standalone 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 millions of model time steps. Furthermore, in regard to future climate simulations the model boundary conditions like e.g. temperature are basically unknown and have to be derived in the standalone case from the output of available global climate models. Although advanced methods exist to make this data usable for high resolution models (e.g. Hay et al., 2007;Chen et al., 2016) there is evidence that in uncoupled ocean simulations the solution of the ocean models is too tightly controlled by the global model (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 for 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 a given time step to the next rather than predicting the variable at a given time. From this follows that the capability of a model to reproduce the transient behaviour 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 space and time scales 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 taken. A larger extension to the open Atlantic or Arctic Ocean substantially increases the memory of the system which has consequences for the model spinup and the economic operation of the model system.
In the past 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 e.g. the east coast of the U.S. (COAWST -Coastal Ocean-atmospherewave-sediment transport model, Warner et al., 2008;. 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 stepstones on the roadmap 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 focus on the specific demands of the Baltic Sea region within Europe. Thus, the focus is mainly set on the main physical feedbacks as they emerge in 2-way coupled earth system components. Therefore, this article does not aim at being overarching as some important components are not considered such as biogeochemical nutrient cycling on land and in the Baltic Sea or atmospheric chemistry (e.g. effect of aerosols).
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 & 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. 28 studies were included in their review, including three 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, and others replaced forest with bare soil. Modelled deforestation in boreal regions resulted in local cooling consistent with observations, but with a less consistent, slight cooling modelled in temperate regions in contrast to observations that indicate a slight warming in those zones. Goa et al. (2014) used the REMO RCM to investigate the biophysical effects of extensive peatland drainage and afforestation in Finland during the 20 th century. Simulations were made for a model domain centred 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 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 2 m air temperature in April, with the highest values being found over the most intensive forestation areas. In contrast, there was a slight cooling (< 0.1 K) over the growing season (May to 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, too, 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 applied lateral boundary forcing (pressure, humidity, temperature, and wind) every 6 h from ERA-Interim reanalyses. Sea surface temperature and sea ice extent were prescribed according to observations. Three simulations were performed, differing only in the land cover map used, for the 30-yr period 1981-2010. The control simulation used RCA4's standard, present-day land cover map defined in the ECOCLIMAP (Champeaux et al. 2005) product. This map reflects the considerable agricultural activity in Europe, with large areas with low forest cover in central, western and southern Europe. To effect maximum afforestation, Strandberg and Kjellström 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 focused on winter (DJF) and summer (JJA) season 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 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 attribute this to a land-ocean coupling whereby warmer and drier air from the surroundings (due to reduced ET over land) favours 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 cannot be explained by albedo changes, since this is reduced in these months. Since 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 in general afforestation 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, though DJF mean temperatures were reduced in 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 for 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, and their effects on extremes of precipitation. 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 to the deforestation and control simulations, with differences of 20-80% near the Baltic Sea region 10% in regions near the west European coast and increasing towards the east to reach 80%. This resulted in a reduction in winter precipitation extremes of up to 25% across the domain. and the biophysical influence (lower part in blue) of re-/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 re-and afforestation or enhanced forest productivity on near-surface temperatures in the Baltic Sea region. According to this, 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. As for the biogeochemical effects, these are increased carbon storage, weakening the radiative forcing, and more aerosols, which reduce the solar radiation reaching the land surface by additional scattering at the particles and more clouds. And as for the biophysical effects, these are increased roughness length, enhancing 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 re-and afforestation or enhanced forest productivity depends on the significance of the warming effect compared to 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 into account the biogeochemical and biophysical effects of land use change and management) 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 demonstrate clearly 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)  surfacetemperatures 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 identified important differences between developed and underdeveloped forests, the latter having lower evaporation rates as a result of their lower LAI and canopy height.

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 understanding 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 has 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 competitive interactions of growing vegetation stands comprising individuals or cohorts of different plant functional types. A number of DGVMs and land surface models (LSM) 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; 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(Haverd et al. , 2014(Haverd et al. , 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 (though see Kumkar et al (2020) for an application using CLM4.5) apart from the LPJ-GUESS (Smith et al. 2001 DGVM. 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 (Prentice et al. 2007; e.g. boreal or temperate, evergreen or deciduous, and broadleaved or needle leaved trees in the Baltic Sea region, and herbaceous species) or species (Hickler et al. 2012) that compete in natural and managed stands (forestry, crops and pasture), in response to climate, atmospheric CO2, and nitrogen (N) availability. As the stand structure evolves in response to environmental change and impacts the availability of key resources, then 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 (forest, natural, cropland etc.) 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 fertilisation (Olin et al. 2015a, b) and clearcut and continuous-cover forest management (Lagergren & 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.  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, burnt area), variables relating to the state and functioning of the soil (water content, C and N content, temperature, runoff, N leaching, 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, CH4 and N trace gases, BVOCs, and net ecosystem carbon exchange.

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 as well as historic and pre-industrial 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 the climate scenarios considered, the authors observed an increase in plant productivity and leaf area LAI, and a northward and upward advance of the boreal forest treeline 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) re-parameterised the most common European tree species in LPJ-GUESS and forced the model with AOGCM climate scenario output, downscaled to a spatial resolution of 10 × 10′. Climate change and CO2 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 broad-leaved forest replaces boreal conifer forest in the Baltic Sea region.

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 modelled natural and managed ecosystems to respond to climate and environmental change, and 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 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 CO2 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 (broadleaved or needleaved) 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 associated LAI. The specific forest PFT types simulated by LPJ-GUESS are aggregated into needleleaved and broadleaved trees before By changing the relative fractions and types in RCA, the LPJ-GUESS fields determine and update dynamically 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 broadleaved 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 CO2 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 CO2 concentrations are not updated over the limited domain covered by RCA. Thus, though 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-feedbackwere 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 in the full simulation period from the long-term means from LPJ-GUESS output in the coupled simulation for . 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 GCM forcing. Wramneby et al. (2010) showed that the snow-masking effect of forest expansion and greater LAI in the Scandinavian mountains, and consequent reductions in albedo, enhanced the winter warming trend. In central Europe, the stimulation of photosynthesis and plant growth caused by the increased CO2 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; 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). Zhang et al. found that warming and CO2 increases promote productivity increases, LAI increases and treeline 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.

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 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 using a bulk formula. Standalone atmosphere models usually read prescribed fields of sea ice and SST from reanalysis data sets or model data to calculate fluxes.

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 to their stand alone atmosphere and ocean modules. This is investigated in a number of studies (e.g. van Pham et al., 2014;Tian et al., 2013Gröger et al., 2015Primo et al., 2019;Kelemen et al., 2019;Akhtar et al., 2019;Cabos et al., 2020;Gröger et al., 2021). 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 4°W and 59°N. As forcing as well as at the lateral model boundaries the authors applied ERA-interim reanalysis data. Noteworthy, the coupled system showed systematically lower T2m air temperatures in the Baltic Sea and North Sea in the long term mean compared to the uncoupled atmosphere model. Consequently, interactive coupling reduced the models mean bias in T2m compared to the E-OBS observational data set 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. far remote from the coupled areas. Analysis of the modeled wind field indicated these areas were situated downwind from the coupled domain North Sea and Baltic Sea implying a more complex pattern of atmospheric advection of temperature anomalies which were not further investigated.
However, similar experiments of 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 bondaries prescribed by ORAS4 (Balsameda et al., 2013) the coupled system showed generally warmer near surface air temperatures over the Baltic Sea compared to 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 to 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 building up atmosphere-ocean feedbacks during summer and winter that control the SST in the Baltic Sea. 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 to observed SSTs. That was mainly caused by a cold bias of the ERA40 data set over Europe. In the coupled mode, driven by the same data set 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 to the uncoupled simulation). This is shown in Figure 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 reduces in the Baltic Sea compared to the ocean standalone model which uses prescribed atmospheric boundary conditions that can not respond to SST anomalies .

Summer thermal short circuit
During summer the above feedback loop is bypassed by the inverse thermal air sea contrast. Since 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 negative effect on wind strength ( Figure 4b). This was demonstrated at stations in stratified areas of the Baltic Sea and in the North Sea using lead correlation analysis between 10 m wind velocity and SST after a short term event of strengthened winds . During the first 70 hours after the event wind and SST were negatively correlated (with a peak r=0.7 at around 30 hours) implying decreasing SST with stronger wind mixing with colder deep waters. After 70 hours, correlation turned to positive values with a peak (r=0.7) around 130 hours as the colder water surface stabilizes the atmospheric boundary layer and thus promotes lower wind speeds. In the following wind mixing ceases again giving rise to heat gain from the warmer atmosphere again Gröger et al. (2015).
These results highlight the importance of thermal air sea coupling in mid-latitude marginal seas and are supported by a number of different studies. Tian et al. (2013) drove a coupled ocean atmosphere model by ERA-Interim reanalysis data (ERA-I). Similar to Gröger et al. (2015) the authors found too 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 by Tian et al. (2013). Moreover, Keleman et al. (2019) applied a model which 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 of winter precipitation patterns over eastern Europe through the altered representation of SSTs. However, SSTs themself were not validated so that a general conclusion about the performance of the whole system can not be drawn. Primo et al. (2019) used the same model showing that the air sea coupling can improve the representation of heat and cold waves but concluded that a general judgement of the whole system is difficult to draw but 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 thus are sensitive to the coupling. In particular, we here focus on three indices with importance for human health, agriculture, and ecosystem of the Baltic Sea: (1) number of tropical nights defined as nights where the daily minimum temperature does not fall below 20°C (e.g. Fischer and Schär, 2010;Teichmann et al., 2018;Meier et al., 2019;Gröger et al., 2021), (2) number of frost days defined as days where the daily minimum temperature falls below 0°C, and (3) number of warm periods defined as at least 3 consecutive days where 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 reference period 1970-1999 from both the coupled and the uncoupled simulation. 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 coupling we refer to Gröger et al. (2015) and . Figure 5a displays the simulated number of tropical nights simulated with the coupled ocean atmosphere model RCA4-NEMO (Wang et al. 2015;Gröger et al., 2019;Gröger et al., 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 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 limits also the exchange with colder waters from deeper layers. These two processes support a strong warming of the Baltic Sea during summer. Consequently, that 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 like north of the Black Sea or in parts of the Atlantic west off northern Iberia (Figure 5a, left). Figure 5a (right) demonstrates that the above described thermal effect over the Baltic Sea is much stronger pronounced in the coupled model including the dynamic ocean model. As a result, over the southern Baltic Sea the number of detected tropical nights within the reference period increases between 50% and 100% compared to the atmosphere standalone 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 while they are completely absent over the southern Mediterranean and southern part of the Atlantic. Again, the different thermal behavior between the North Sea and the Baltic Sea is obvious. The Baltic Sea behaves thermally like the continents supporting 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 underduring winter. By contrast the adjacent North Sea effectively damps the occurrences of frost days. However, a pronounced east west gradient is visible with lower 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 a similar thermal behavior as discussed for the Baltic Sea. The southern part of the eastern North Sea, namely the German bight is only shallow and supports rapid cooling. Altogether, this builds up the east west gradient seen in  Figure 5b (left). Generally, these thermally forced processes are represented in both the coupled and uncoupled model version but in the coupled model the number of frost days is significantly reduced nearly everywhere over the Baltic Sea (Figure 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 positive feedback on wind speed and 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 Figure 5c. In the Baltic Sea such periods are an important precondition for the occurrence of cyanobacteria blooms during summer. Also here, the southern Baltic Sea turns out to be a hot spot for the effect of interactive air sea thermal coupling as the number of such periods in the coupled model exceeds the corresponding number in the standalone atmosphere model by several orders of magnitude. Here the effect of the aforementioned summer thermal shortcut is seen. Enhanced mixing by winds brings cooler waters from depth to the surface. The cooler surface water likewise cool the air temperatures which imposes a stabilizing effect on the atmospheric boundary layer over sea thereby damping wind strength again (see Gröger et al., 2015 for details). This facilitates the development of longer lasting warm periods under summer.
So far, we discussed thermal effects of damping (in case of frost days) and amplifying (warm days, tropical nights). The associated feedback loops are more realistically represented in the fully coupled model as previously explained. However, a significant effect apart from the active coupling area over land seems to be rather weak. On the other hand we note that the interactive coupling area must be considered small compared to the whole domain. Thus extending the area of interactive coupling by e.g. including also the Mediterranean or larger parts of the North Atlantic may result in more intense effects of using coupled models Keleman et al., 2019;Akhtar et al., 2019;Cabos et al., 2020).

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, ERA40 SSTs often used in uncoupled atmosphere simulations are based on weekly or even monthly frequency (Fiorino, 2004;Uppala et al., 2005). Consequently, many studies compared coupled and uncoupled models with respect to the representation of extreme events in hindcast simulation. More recent products like e.g. ERA5 may improve the situation by providing a higher spatial and temporal resolution.
Jeworrek et al. (2017) diagnosed a better representation of atmospheric conditions favorable for the occurrence of convective snowbands. The authors attributed this improvement to a more accurate simulation of SSTs and subsequent air sea heat and moisture fluxes in case 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 snowbands has been recognized early by Gustafsson et al. (1998). Comparing two regional ocean atmosphere models in coupled and uncoupled mode Ho- Hagemann (2015 found that interactive air sea coupling can alter extreme precipitation over land. Ho-Hagemann et al. (2015 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 to the stand-alone atmospheric model CCLM. This was demonstrated for several flood events in Central Europe. A diagram displaying cause and effects as well as feedbacks between earth system components is shown in Figure 6 to explain the physical mechanism behind the better representation of heavy rainfall. The main effect is an altered SST in the coupled model further influencing wind speed and evaporation (so called Wind-Evaporation-SST (WES) feedback, Xie & Philander (1994)). When wind speed increases over an area, evaporation increases and so the latent heat flux is enhanced often leading to ocean cooling over this area. 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 article about the influence of atmospheric ocean interactions on heavy rainfall over Europe is available by Ho-Hagemann and  The Impact of air-sea coupling on simulations of mid-latitude cyclones was recently investigated using an ensemble approach by Ho-Hagemann et al. (2020). The coupled model GCOAST-AHOI reduces the large spread of wind speed, mean sea level pressure, surface temperature, cloud cover as well as radiation fluxes amongst ensemble members of the atmospheric model during the storm Christian occurring from 27 to 29 October 2013 in northern Europe.

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 2 Ho- Hagemann et al. (2015Hagemann et al. ( , 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).

Models
Atmosphere and others  A recent assessment of regional coupled modeling (Schrum, 2017) emphasized that both the location and extension of the coupled region (Sein et al., 2014), the coupling frequency (Fang et al., 2009), as well as the quality of initialization and boundary forcing (Wei et al., 2014) is critical. In the following, we focus on the size of the coupling area and consider its potential impact on climate simulations. We will elaborate options about suitable sizes for the coupling area as well. In studies where the coupling domain is comparatively large to the total domain often comprising multiple seas with a large heat inventory, the air sea coupling effect is often found to exten dt far inland (e.g., Somot et al. 2008, Ratnam et al. 2009, Ho-Hagemann et al. 2015. 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 Europe-Atlantic region.
A known problem for many atmospheric models (Vidale et al. 2003) is a dry bias over large areas of mid-latitude 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- . 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 of the atmosphere models, the realistic simulation of air sea feedbacks enhances the large scale wind and evaporation via altering the sea surface temperature and the land sea heat contrast, and therefore, reduces the dry bias.

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 reflects even 90 % of incoming solar radiation instead of high absorption of 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 as well as surface and bottom roughness since those control radiation, heat, moisture and momentum fluxes at the atmosphere-ice-ocean interface. The growth rate and temperature of sea ice impact on salt flux between ice and ocean. However, due to low salinity of the northern Baltic which experiences annual ice cover, that mechanism doesn't have a significant effect.
Theoretical frameworks of presently used sea ice models in climate applications were established already in the 1970's. The Tthermodynamical evolution of ice and snow is based on the classical heat conduction law, resolved numerically firstly by Maykut and Untersteiner (1971). The model resolves vertical temperature and salinity structure and surface temperature iteratively. The present community sea-ice models LIM-3/SI 3 (Rousset et al., 2015;Vancoppenolle et al., 2009) and CICE (Hunke and Dukowitc, 1997) apply this classical framework but include detailed parametrizations of snow, flooding, snow-ice formations, melt ponds, albedo and brines. For climate applications, the vertical structure is usually resolved for 1 to 3 layers. The oneOne layer model assumes a linear temperature profile on ice. This is a valid approximation for Baltic Sea where thermodynamically grown sea ice rarely exceeds one meter.
The Mmomentum 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 island. In present coupled models, two rheological solutions are commonly used. The Hibler (1979) viscous-plastic rheology (VP) implies that under very low strain rates the bulk and shear viscosities are constant and the model produces linear viscous behavior; otherwise the viscosities are calculated according to the plastic flow rule. VP rheology resolves non-linear behavior of Baltic sea ice dynamics well (Leppäranta et al. 1997). LIM-3/SI 3 model apply VP-rheology and it's a rheological choice of the NEMO-Nordic model (Pemberton et al. 2017;Hordoir et al., 2019). However, since VP-rheology is computationally demanding, numerically more feasible elastoviscous-plastic rheology (EVP; Hunke and Dukovich, 1997) is used in the CICE model. It also widely used in Baltic Sea applications Jakacki and Meler, 2018;Janecki et al, 2018).
The Tthird element of the sea ice models is a resolving 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. Choice of the minimum ice thickness has an impact on the modeled ice edge. Ridging of ice is taken into account since ice thickness can freely increase during the convergent ice motion, although 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 real behavior of the pack ice in a continuum scale.
These three governing equations are strongly coupled. Firstly, sea ice mobility is non-linearly related to the ice thickness and concentration. Thin or low concentration pack ice drifting close to free drift speed but even 0.5 meters thick solid ice cover can be stationary under action of strong winds in the Bay of Bothnia. In turn, ice motion generates fractures and leads on ice pack which enhance mobility, but more importantly increase sea ice mean thickness by new ice growth in leads and formation of pressure ridges in compression. Due to ice dynamics, mean ice thickness in the coastal boundary zone is thicker than in landfast ice regions (Oikkonen et al. 2016;Ronkainen et al, 2018). Mobility of pack ice has large consequences on formation of coastal leads which are local sources of heat and moisture.
In the Baltic Sea, correct modeling of landfast ice, which can extend several kilometers from shore, is essential. In the CICE, 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, 1966). This approach has been used in several applications and implemented on the Nemo-Nordic (Pemberton et al., 2017 ).
A study of Baltic Sea climate variability based on the Hibler model type coupled with Bryan-Semtner-Cox, z-coordinate baroclinic ocean model (the Rossby Centre Ocean model, RCO), developed by M. Meier (Meier 2002a, 2002b, Meier and Dorsher 2002, was performed in the early 2000's. RCO and the University of Helsinki sea-ice model (HIM), has been used for analysis of future ice conditions of the Baltic Sea region, (Haapala et al. 2001). Based on these models the authors made two 10-year simulations representing preindustrial and future scenarios of Baltic Sea ice conditions. Both models, on a global scale, delivered similar results, however on a regional scale there was a large variation between model results. Studies show dramatic decreases of ice extent, also the calculated ice thickness was lower in the scenario simulation and it is expected, based on this simulation, that Archipelago Sea and Quark will not be covered by ice in the future. The influence of greenhouse gas emissions (A2 and B2 IPCC scenarios that represent the climate of the late 21st century -2071-2100) on ice conditions was performed based on 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, a reduction of annual maximum ice thickness up to 60% and a decrease of ice days over 90%. All of the results are locally dependent. They also concluded that total ice area depends on air temperature with less influence from other physical factors.
Recently two present day coupled ice-ocean models have been developed for the Baltic Sea regionice part (LIM3.6) was evaluated in the NEMO-NORDIC model (Pemberton et al., 2017), that covers Baltic and North Seas and B-CESM (Jakacki and Meler, 2018), only for the Baltic Sea area. 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.

Coupling strategies and pitfalls 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 include the Earth System modeling Framework (ESMF;Hill et al., 2004), CPL7 (Craig et al., 2012) 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., 2012). On the other side, there is the concept of the OASIS and YAC coupler (Valcke, 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, 2013) has been relaxed with the inclusion of the Model Coupling Toolkit (MCT; Larson et al., 2005; in the OASIS3-MCT coupler (Valcke etl 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 ESMF to build a regional ESM with a focus on coastal processes in the North Sea and Baltic Sea. The Modular Earth Submodel System (MESSy; Jöckel et al., 2010) originated in the need for atmospheric chemistry being 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 (Tab. 2.2.3.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, since 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 e.g. 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., 2016).
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 with the complex coastlines and so far a relatively low resolution of 0.11 to 0.22 degrees in the atmosphere models and 1/60 to 1/20 degree in the ocean models. With 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 ambiguous assignment of grid cells is the remapping be carried out on an exchange grid that is the joint set of grid cells between two participating grids (Balaji et al., 2005, Bauer et al., 2020. The exchange grid also addresses the vertical flux between different model components but has so far not yet been used in RCMs of the Baltic Sea region. 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., 2019b). 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 aimed at demonstrating improvements in simulating aspects of the present day climate. However, a general statement about model performance of coupled versus standalone models can not be drawn as many studies showed added value of coupling especially in those cases when the uncoupled model was driven by low quality forcing data (Tian et al., 2013;Gröger et al., 2015;. Furthermore, andand the uncoupled models often had a lower resolution at the air sea boundary. In turn, this implies the application of coupled models in particular in future climate projections as global climate models can have large biases and constraints of the governing model that drive the uncoupled model may be too tight (Mathis et al., 2018).

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 designated waves) are mainly generated by the wind; the wave field is thus 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 (wave height, wave age, wave steepness, etc). The transport of momentum is the key property, but waves also influences the exchange of heat and mass, as well as the turbulence of the lower atmosphere and the upper ocean. Indirectly, the 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 (ECMWF, 2017), 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,τoc/τ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. It leads to the level off of the drag coefficient 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 as well as the mixing Wu et al., 2017a). Wave breaking can transfer a significant amount of energy flux into the ocean (Melville et al., 2002), which mainly affects a 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), α CB u ❑ 3 , in the surface boundary (Craig and Banner, 1994), where u ❑ is the ocean-side friction velocity and α CB is a wave-related parameter which is usually treated as a constant (e.g., 100) without considering the impact of wave status.
Due to the periodic motion, waves induce a net drift in the wave propagation direction (u s =u L −u, the difference between the Eulerian velocity, u, and Lagrangian velocity, u L ), which is defined as Stokes drift (Stokes, 1880). The Stokes drift can impact on 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, K-Profile Parameterization turbulence scheme (McWilliams and Sullivan, 2000), k −ϵ turbulence model (Axell, 2002), etc. 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).
Nonbreaking-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 the 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, more detailed description of the wave-ice interaction as well as their parameterization in models refer to Zhao et al. (2015) and Squire (2018).

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 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 model, such as NEMO-WAM (Alari et al., 2016;Wu et al., 2019b;Staneva et al., 2017), two-way coupled model GETM-WAM . 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 under young waves since they extract more momentum for their growth (Janssen and Viterbo, 1996). Accordingly, wave coupling is more important under extreme conditions (young wave dominating). Based on the simulation of 23-year storms using WRF-SWAN model with a domain covering North and Baltic Sea, Larsén et al. (2019) found that the wave influence improves the model performance compared with the uncoupled model Sea spray under high winds can significantly impact the momentum and heat flux inat 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);Zhao et al. (2006). Together 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 have shown 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 include radiation stress, Stokes drift, bottom friction modification, and turbulent kinetic energy due to wave friction . They found that the coupled system does not have a significant influence over the open North Sea, however, it significantly improves 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 indicates that the wave-related processes improve the storm surge simulation in extreme conditions. Compared to the other waverelated processes investigated in their study, the sea-state-dependent momentum flux plays a dominating role which agrees with the results of Wu et al. (2019b).

Mesoscale features
Coastal upwelling occurs frequently in the summer months ofin the Baltic Sea, which is mainly induced by the divergence (convergence) of the wind stress (Lehmann and Myrberg, 2008). Ocean waves affect the coastal upwelling through altering the ocean-side stress (Eq. 1) and upper-ocean mixing (mainly through CSF and breaking wave-induced energy flux). Based on an offline waveocean 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 in the eastern coast of the Baltic Proper. In their study, three wave-related processes were implemented in the ocean model, i.e., CSF, sea-state dependent energy and momentum flux. With a similar ocean model set up in the Baltic 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 during 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 cancelled the influence of 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 up to 10% which varies with location (Fig. 7). And the combined effect of those wave-related processes increases the weak upwelling frequency (−4 < ΔT < 2.5 C) but decreases the strong upwelling frequency (ΔT < − 6 C), here (ΔT) ∘ ∘ is the SST difference from the zonal mean temperature (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 those 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 through adding an extra wave-age dependent coefficient to the mixing length . Simulations show that the swell impact on the atmospheric mixing can alter the mean surface wind up to 0.3 m/s in the Baltic Sea and the North Sea Wu et al., 2016). The magnitude and direction of the change depend on the wave field and environmental conditions Wu et al., 2016).
Based on a three summer month simulations, 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 up to 1 degree in the Baltic Sea, in which, the sea-state dependent momentum and TKE flux play anda warming role but 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 coupled model with an improved representation of the air-wave-sea interaction processes. Based on two 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 onin the summer months than in winter months (Wu et al., 2019a) which is because the wind speed is higher in winter and sea-statedependent momentum/energy flux and the Stokes drift related processes are more important for upper-ocean mixing.

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 one half of the precipitation over land areas (Lind and Kjellström, 2009), and about twice as much as 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., 2020). Decadal variations in Baltic Sea salinity are caused largely by the accumulated runoff to the Baltic Sea (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 North-Western European Shelf .
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 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 Baltic Sea is almost decoupled from the open ocean so that land-borne nutrients significantly contribute to determining ecosystem productivity (Thurow, 1997;Österblom et al., 2007).
Consequently, river runoff is also an important component for the coupled system modeling over the Baltic Sea region. Hagemann et al. (2020) provide an overview on the current state of high resolution modeling of river runoff (or discharge) within the framework of regional coupled system models (RCSM). Note that the present section 2.4 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 been recently 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, e.g. 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.

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 either 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 . To reflect the projected increase in precipitation in the northern part of the Baltic Sea (e.g., Donnelly et al., 2014), in the Bothnian Sea and Bothnian Bay a linearly increasing discharge to +10% in 2100 is used. Thus, only a few RCSM setups exist where a discharge model is included, but its resolution is rarely higher than 0.5°. Examples comprise 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 andSud, 1998) at 0.5° (Dell'Aquila et al., 2012;Sevault et al., 2014) and LARSIM (Bremicker, 2000) at 1/6° over Northern Europe (Lorenz and Jacob, 2014).
The latter 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. (2015Sein et al. ( , 2020a 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 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.

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, currently these studies cover only 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 landsurface models over the Danish Skjern River catchment (area: 2500 km 2 ).
Over Korea, the discharge model TRIP has been coupled to a RCM at 0.5°, 0.25°, and 0.125° in preparation of future RCMS studies (Lee et al., 2015). Nguyen-Quang et al. (2018) used a highresolution (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°a nd 0.25° resolution. To our knowledge, none of the hydrological models used in the studies listed above and by Hagemann et al. (2020) has been used in a fully coupled RCSM setup that can be applied for climate time scales and large-scale areas.
Very recently, Hagemann et al. (2020) developed a high-resolution version of the HD model that is globally applicable at 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 the Helmholtz-Zentrum 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.

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 on the one hand the water cycle in coupled general circulation models (GCMs) is often not closed (Liepert and Lo, 2013). On the other hand, GCMs and RCMs suffer from substantial biases, especially with regard to precipitation and the hydrological cycle (Flato et al., 2013;Kotlarski et al., 2014), which affect the simulated discharge. An example comes from 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 increasing discharge in the Bothnian Bay and the Bothnian Sea by 10% in 2100  or use discharge data from uncoupled hydrology models Meier et al., 2018;Meier et al., 2019aMeier et al., , 2019b. 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, e.g., to an extraordinary strong freshening in future scenarios (Sein et al. 2020a).

LAND-ATMOSPHERE
Previous works have shown the potential for substantial land-atmosphere coupling in the Baltic Sea region, through land use change, natural vegetation dynamics and land management. RESMs that downscale CMIP6 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 BVOC (Kulmala et al. 2014;Hantson et al. 2017) and aerosols precursor from natural and managed vegetation and wildfires, CH4 and N2O emissions from wetlands and agriculture.

ATMOSPHERE-OCEAN-SEA ICE
With respect to coupled atmosphere-ocean modeling the choice of the size of the coupled domain has always to be carefully considered and in most cases the economical costs to drive the model, sets an upper limit. When the atmosphere domain covers several seas as it does in Euro-CORDEX domain it might be feasiblefieasible to coupled only a few of them to save computationalconputational 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 as there important small scale processes take place determining the shelf-open ocean exchange as demonstrated for the North Sea (Gröger et al., 2013) that 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 and 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 emphasized in simulations for the Arctic Ocean that inclusion of the North Pacific in coupled domain destroys 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 .
Along the coastal zone small scale oceanic processes like upwelling can build up 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 for coastal zone dynamics can be expected from higher resolution of atmosphere models up to a few kilometers together with the transition to convection permitting models instead of hydrostatic models. With respect to future scenario climate simulation more research is necessary to assess whether or not the coupled and uncoupled models reveal different climate changes signals in scenario simulations in both the hydrographic and atmospheric properties of the Baltic Sea region. First attempts to address this question have been made (e.g. Gröger et al., 2021Meier et al., 2021).

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 off-shore wind energy applications in the coastal zone). It has also been shown that wave information improves the description of the 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 introduction of as accurate information of the earth system as possible is essential. For this, coupling to other high resolution models is a prerequisite. This includes primarily atmosphere, ocean and waves, but can be expected to also be beneficial with submodels for ice, aerosols and etc. Improvements by including wave models into regional earth system models can thus be expected in particular for features like convective precipitation, cold-air outbreaks, lake-snowfall, polar lows and mid-latitude storms.

HYDROLOGY
Socio-economic 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. . 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 which has not yet been implemented in hydrological models. This is especially important in the context of future climate change when 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 such as water abstractions, e.g. for irrigation, and regulation, e.g. by 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 as well as modules to simulate water withdrawals, e.g. for irrigation. Apart from prescribed scenarios, the future development of these impacts might be realized with simple economic modules.

INTERNAL VARIABILITY
Hydrodynamical earth system components, i.e. ocean and atmosphere primitive equation models are characterized by generating 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 usually long term averages and statistics are considered this is no problem provided that the averaging period largely exceeds the frequencies of internal variability of the system. It is neither problematic if short term events are analyzed, provided that only the statistics of these events over an appropriate long period is interpreted. Contrary to this, i.e. when single events are analyzed like e.g. the timing of certain storms (e.g. Ho- Hagemann et al., , 2020 or wind induced coastal upwelling events in the Baltic Sea. In these cases the model solution can be substantially modulated by internal variability and so the comparison of coupled vs uncoupled systems may be misleading. In turn, this requires a profound estimation of the robustness of the model solution with respect to the initialization, boundary conditions (Weisse et al., 2020;Ho-Hagemann et al., 2020) and even the computing platform (Geyer et al., 2021). For this, different methods have been developed mainly based upon the generation of large ensembles taking into account small perturbations in the initialization (e.g. Giorgi et al., 2000;Weisse et al., 2000;Weisse et al., 2003;Sieck and Jacob, 2016).
Comparing ensembles of coupled and uncoupled regional ocean atmosphere models Ho-Hagemann et al. (2020) has been recently demonstrated that interactive air sea coupling can substantially reduce internal model variability compared to uncoupled atmosphere models. In case of the storm Christian occurring from 27 to 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 parameterization of cloud-radiation interaction in the atmospheric model and (2) lack of an active two-way air-sea interaction.
In a similar approach of Ho-Hagemann et al. (2020), Wiese et al. (2020) also found the reduction of internal model variability in ensembles of an interactive atmosphere-wave coupled model compared to those of a stand-alone atmospheric model. The role of internal variability is still not sufficiently investigated in more complex coupled systems involving more components like e.g. fully coupled ocean-wave-atmosphere models as pointed out by Wiese et al. (2020) or coupled ocean-atmosphereland vegetation models.

OTHER COMPONENTS
For the current generation of RCMs a substantially lower projected 21st century warming compared to their driving global models has been demonstrated and the neglection of scenarios for time varying aerosolsaerosoles 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 aerosoles is a major challenge in RCMs and has been recognized within the CORDEX Flagship Pilot Study (FPS-Aerosol https://www.hymex.org/cordexfp s-aerosol/ wiki/doku.php?id=start).
Despite the advancements in including the biophysical feedback from the land vegetation as outlined in section 2.1, analogous feedbacks from the marine biogeochemical cycles on the physics are so far not implemented in RCSMs. The main feedback is due to the altered penetration of solar insolation by marine biota that further influence 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., 2020b) 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 of cyanobacteria blooms in a future warmer climate (e.g. Saraiva et al., 2019a;Saraiva et al., 2019b;Meier et al., 2019a;Meier et al., 2019b).

Conclusions and key messages
For the Baltic Sea region previous research identified a number of important feedback loops between earth system compartments that alter both the mean climate as well as 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 North and central Europe while an amplification of warming up to 1 K is obtained for southern Europe . Likewise, implementation of scenarios for land management in the Baltic Sea region can alter temperature locally to up to a few Kelvin (Strandberg and Kjellström, 2019). The alterations emerge from a number of coupled landatmosphere processes involving biophysical feedbacks such as changes in albedo and roughness length and evapotranspiration and related changes in LAI and vegetation-type as well as biochemical feedbacks from high carbon dioxide concentration.
Unlike the land surface where evaporation depends also 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 mainly identified thermal air sea coupling as most significant (e.g. Kjellström et al., 2005;Gröger et al., 2015;Gröger et al., 2021) with the predominant impact on simulated SSTs. For the Baltic Sea deviations between prescribed SST and sea ice fields and modelled quantities can add up to a few Kelvin (Tian et al., 2013;Gröger etal., 2015). In turn, good quality SST fields appear a prerequisite for the representation of extreme events over land such as convective snowbands, 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 demonstrate the importance of interactive wave coupling in mediating the atmosphere-ocean exchange of heat, momentum and mass which enables 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;Hagemann et al., 2020) and are an important step towards the closure of the water cycle. However, they likewise 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 in particular problematic for the Baltic Sea where runoff constitutes a major contribution to the halocline structure and the freshwater budget. On the road to Many studies aim at demonstrating added value of interactive coupling by direct comparisons of single variables between coupled and uncoupled models. However, this may be misleading as in many cases, these improvements simply reflect bad quality of prescribed boundaries used in the stand alone model (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 bad 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 parent biased global model by generating its 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 economic standalone models are a good tool provided that forcing and boundary data is of good quality. For future climate simulations where this is not the case, rather coupled models will be the first choice. The higher costs of coupled climate simulations will require to reduce the size of individual ensemble simulations which hampers a robust assessment uncertainties due to global models and scenarios. However, advanced techniques for reducing the ensemble size by by conserving the model ensemble spread are currently under development (e.g. Wilcke and Bärring, 2016).
Author contributions. WH and PM wrote the section about land biosphere coupling, HTMHH, HEMM, CD, MG wrote the section about ocean atmosphere coupling, AR and LW wrote the section about atmosphere ocean wave coupling, SH and CD wrote the section about hydrology coupling, and JH and JJ wrote the section about atmosphere ice ocean coupling. The sections of conclusions, uncertainties and gaps represent the joint effort of all authors. The final manuscript was compiled by MG.
Competing interests. The authors declare that they have no conflict of interest.
Code/Data availability. The review paper uses exclusively published data. We refer to the referenced literature and corresponding authors for reasonable requests.