Ocean‐Wave‐Atmosphere Interaction Processes in a Fully Coupled Modeling System

A high‐resolution coupled ocean‐wave‐atmosphere model (Uppsala University Coupled model, UU‐CM) of the Baltic Sea and the North Sea with improved representation of ocean‐wave‐atmosphere interaction processes is presented. In the UU‐CM model, the stress on the air‐sea interface is estimated as a balance of four stress terms, that is, the air‐side stress, ocean‐side stress, wave‐supported stress (absorption of momentum by the wave field), and the momentum flux from waves to currents (breaking waves). The vector differences between these four stress terms are considered in the coupled system. The turbulent kinetic energy flux induced by wave breaking, the Stokes‐Coriolis force and the Stokes drift material advection terms are added to the ocean circulation model component. Based on two‐month‐long (January and July) simulations, we find that the ocean‐wave‐atmosphere coupling has a significant influence on coastal areas. The coupled system captures the influence of surface currents and local systems such as coastal upwelling and their impact on the atmosphere. The wave‐current interaction enhances the upper ocean mixing and reduces the sea surface temperature in July significantly. However, the pattern of the wave‐current processes influences on the ocean current and waves are complex due to the stress differences in both magnitude and direction.


Introduction
Ocean-wave-atmosphere (OWA) interaction processes control the energy, momentum, and gas transfer between the ocean and the atmosphere. Traditionally, atmosphere, ocean, and wave models are run independently from each other, the OWA interface conditions are prescribed using bulk formulae, and the interaction between the components are ignored. It is common to tune the model parameters and, in general, this improves model performance. However, tuning cannot fully capture the dynamic processes across the air-sea interface (Wu, Sproson, et al., 2017). Thus, building seamless models is one of the ways to capture those interaction processes properly (Hov et al., 2017).
The first attempts of coupling atmosphere and ocean models were carried out in the late 1960s and early 1970s (Bryan et al., 1975;Manabe & Bryan, 1969). These models were designed mainly for global climate studies. The resolution and the complexity of the models have been dramatically increased in past decades. More and more components have been added to the coupled models, for example, sea-ice, atmospheric chemistry, and the carbon cycle (Barth et al., 2000;Hibler & Bryan, 1987;Qiao et al., 2013). Compared with coupled climate models, coupled numerical weather prediction models have developed much more slowly. Because the temporal scale of the ocean is larger than for the atmosphere, it has long been thought that the weather forecasts have little to benefit from coupled models. However, studies have shown that atmosphere-ocean coupled models can improve the model forecast skills, especially during severe weather (Bao et al., 2000;Loglisci et al., 2004;Mogensen et al., 2017;Sandery et al., 2010). A great number of coupled models (regional and global coupled models) have been developed for a variety of purposes during the past two decades (Bruneau & Toumi, 2016;Eyring et al., 2016;Giorgi & Gao, 2018;Lewis et al., 2019;Peng et al., 2012).
The OWA interaction processes can directly affect the heat and kinetic energy transfer between the ocean and atmosphere. Accordingly, they can alter the atmospheric and oceanic dynamic processes from the submesoscale to the global scale. A variety of regional-scale phenomena can be better captured in a coupled model, such as monsoon intraseasonal oscillations in the Asian-western Pacific region (Fu et al., 2007), El Niño-Southern Oscillation (Newman et al., 2009), and Madden-Julian oscillation (Newman et al., 2009). At the mesoscale and synoptic scale, air-sea coupling processes improve model performance concerning the intensity and influence regions of extreme events (i.e., Jeworrek et al., 2017;Lengaigne et al., 2018;Srinivas et al., 2016). A recent study by Byrne et al. (2016) found that the atmosphere-ocean coupling can enhance the transfer of wind energy into the ocean.
The influence of surface gravity waves on air-sea interaction processes has mostly been ignored in the development of coupled models in both climate and numerical weather prediction since the temporal and spatial scales of surface waves are much smaller than the atmospheric and oceanic dynamic scales (Hasselmann, 1991). A notable exception is the coupled atmosphere-wave forecast system developed by the European Centre for Medium-range Weather Forecasts (Janssen, 1989(Janssen, , 1991(Janssen, , 2004. The influence of surface waves on air-sea interaction processes has received considerable attention in recent years, and a wave model has been implemented in some coupled models Janssen & Bidlot, 2018;Liu et al., 2011;Qiao et al., 2013;Warner et al., 2010). Surface waves affect the marine atmospheric boundary layer by modifying the sea surface roughness (Janssen, 1989(Janssen, , 1991. Many sea state-dependent roughness length parameterizations have been proposed in order to capture the wave influence on the marine atmospheric boundary layer (e.g., Drennan et al., 2003;Liu et al., 2012). The sea surface roughness determines the momentum and heat fluxes between the atmosphere and the ocean interior and can thus affect the climate and individual weather events (Moon et al., 2004;Wahle et al., 2017;Wu et al., 2016). Under swell conditions, surface waves can have a significant influence on the atmospheric mixing for neutral and a slightly unstably stratified atmosphere (Rutgersson et al., 2012;Semedo et al., 2009;. Wave-current interaction processes affect both the upper ocean turbulence and the wave field. The wave-current interaction can modulate the sea surface roughness (e.g., Munk et al., 2000) and wave breaking (e.g., Melville et al., 2005) through shoaling effects (Holthuijsen, 2007). The nonlinear interaction between waves and currents can contribute to extreme wave events (e.g., Janssen & Herbers, 2009), which are a serious threat to offshore activities. The indirect influence of currents on waves can come from the modulation of wind and atmospheric stability caused by the atmosphere-ocean interaction. Wave breaking injects turbulence kinetic energy (TKE) into the near-surface ocean layer (Terray et al., 1996), which decays rapidly with depth (Janssen, 2012). The Stokes production related to the shear of the Stokes drift is thought to drive the generation of Langmuir turbulence (McWilliams et al., 1997;McWilliams & Sullivan, 2000) and to affect a much deeper layer down to the bottom of the mixed layer (e.g., Sullivan & McWilliams, 2010) and indirectly even deeper (e.g., Breivik et al., 2015). Some studies indicate that the nonbreaking wave can also affect the upper ocean mixing (e.g., Babanin, 2006;Qiao et al., 2004). Recently, Christensen et al. (2018) discussed the importance of wave processes on the short-term predictions of oceanic drift.
Previous studies have implemented many wave-related processes into traditional models using off-line coupled models (e.g., Alari et al., 2016;Law Chune & Aouf, 2018;Staneva et al., 2017;Wu et al., 2019). In these off-line experiments, the ocean models get the wave parameters from free-standing wave model integration. Indirect (feedback) effects between the atmosphere, the wave field, and the ocean can obviously not be captured in forced mode. In the previous study by Wu et al. (2019), the impact of wave-current interactions on the ocean simulation was investigated in an off-line ocean model. In this study, the impact of the OWA interaction processes on the simulations is investigated using a fully coupled atmosphere-wave-ocean-ice model in order to capture the direct and indirect influence of the interaction processes. Four sensitivity experiments are devised to investigate the role of OWA interaction processes in a high-resolution model domain covering the Baltic and North Seas. Based on these four experiments, we explore the impact of going from a realistic stand-alone atmospheric setup (Exp1) and atmospheric-wave model setup (Exp2) to a coupled system (Exp3 and Exp4). This paper is organized as follows: Section 2 presents the OWA interaction processes implemented in the coupled model. Section 3 describes the fully coupled model and the experimental setup. The data sets used in this study are introduced in section 4. The results are analyzed in section 5 and discussed and summarized in sections 6 and 7, respectively.

Air-Sea Wave Interaction Processes
The OWA interaction processes investigated in this study are summarized in this section. We mainly focus on the processes related to surface gravity waves.

Wind-Driven Waves
The wind blowing over the ocean surface generates waves with wavelengths ranging from capillary waves (with wavelengths less than 2 cm) to swell waves (where the wavelength can exceed a kilometer). The development and dissipation of waves in the ocean surfaces are determined by the energy balance in the air-sea interface. In deep water, the balance equation for the wave action density spectrum N = F∕ can be written (WW3DG, 2016), Here, F represents the wave spectral density, is the intrinsic circular frequency, k the wave number vector, c g the group velocity, c k represents the spectral advection velocity, u the ocean current, ∇ and ∇ k are the divergence operators in horizontal and spectral space, respectively, S nl is the nonlinear transfer source term, S in is the wind input source term, and S ds is the wave dissipation source term. In shallow water, more source terms will be added into the right-hand side of equation (1), such as wave-ice interaction, wave-bottom interaction and depth-induced breaking (WW3DG, 2016).
The wind input term, S in , is the source of the wave growth, which is usually modeled as positive definite (i.e., energy and momentum are not allowed to flow from the wave field to the atmosphere). The wind input has been studied for more than six decades, but the mechanism of wind-generated waves is still not well understood (see the review by Hristov, 2018). The two wind-wave generation mechanisms traditionally assumed are Phillips's resonance theory (Phillips, 1957) and Miles's instability theory (Miles, 1957). Based on wind-wave interaction mechanisms, several wind input source terms have been proposed and implemented in wave models, such as WAM cycle 3 (Group, 1988), Tolman and Chalikov (1996), WAM cycle 4 (Janssen, 2004), ST4 (Ardhuin et al., 2010), and ST6 (Zieger et al., 2015). The wind input term takes the general form in which k is the wave number, the direction, and (k, ) is the wind-wave growth rate. In Ardhuin et al. (2010), the wind-wave growth rate is parameterized as a function of wind, wave age, and roughness length.

Momentum Flux at the Air-Sea Interface
It is commonly assumed that the air-side momentum flux, a , is identical to the ocean-side momentum flux, oc . This is valid only when the ocean waves are in equilibrium with the local wind (Pierson & Moskowitz, 1964), that is, fully developed. Most of the time surface waves either release more momentum to the ocean interior than they absorb from the atmosphere or vice versa. Thus, considering the buffer role of surface waves at the air-sea interface, the momentum flux balance is (ECMWF, 2017) (3) where in is the wave-supported stress, and ds is the momentum flux from the wave field to mean currents, which is always negative (since the dissipation of wave energy and thus momentum must always be a sink to the wave field).
The airside momentum flux is usually estimated based on the bulk formula, where a is the density of air, u 0 is the surface current, and the drag coefficient C d is usually estimated from Monin-Obukhov similarity theory, where is von Kármán's constant, z 10 = 10 m, and z 0 is the airside surface roughness length. Over the ocean, z 0 is usually estimated from the Charnock relationship, that is, z 0 = c u * 2 ∕g, where c is the Charnock coefficient in the airside (Charnock, 1955) and u * is the air-side friction velocity. The Charnock coefficient is usually treated as a constant or a function of the mean wind speed. Based on these parameterizations, the airside stress is indifferent to the sea state. This is clearly not an accurate way to estimate the airside stress.
When the 2-D wave spectra are available, the Charnock coefficient c can be calculated following Janssen (1991) wherê= 0.0095 (Bidlot et al., 2007). The wave-supported stress is calculated from the 2-D wave spectrum, The total momentum flux to the ocean interior becomes (8)

Stokes Drift-Related Processes
The Stokes drift, u s = u L − u, is a second-order effect of linear waves proportional to the steepness of the waves. It is the difference between the total Lagrangian velocity, u L , and u, the Eulerian velocity (what is actually calculated by an ocean model). Water particles spend more time in the forward motion under the crest and higher in the water column, where the velocity is higher, than in the backward motion under the trough where it moves against the wave propagation direction at a deeper level (where the velocity is lower; Stokes, 1847;van den Bremer & Breivik, 2018). This prevents particle orbits from being completely closed, and the net effect is a slow forward drift. The Stokes drift affects the upper ocean flow through the interaction with the Coriolis force, the Coriolis-Stokes force, see Hasselmann (1970), the Stokes advection, and the Stokes shear force (Annex A1 of . After averaging over the fast time scale of surface waves, the governing equations of Eulerian ocean models can be expressed as ∇ · u = 0 (12) Here, c is a scalar quantity, such as temperature and salinity; is the sea surface height; D u and D c are the parameterizations of subgrid-scale physical processes for momentum and tracer equations, g is the gravitational acceleration withẑ the vertical unit vector. The Stokes drift impacts on the mass transport through the divergence of the sea surface height is expressed by Equation (14). For a detailed description of the ocean model set up see Wu et al. (2019). In this study, and in contrast with Wu et al. (2019), the Stokes advection in the momentum equation is switched on. The Stokes shear force, which is the source of the Langmuir turbulence (Suzuki & Fox-Kemper, 2016), is not discussed here.

Wave Energy Flux to Ocean
Wave breaking injects TKE into the upper ocean, typically to a depth of the order of the wave height (Terray et al., 1996). The TKE flux can be calculated from the spectral source terms in a wave model It is often assumed (Craig & Banner, 1994) that the energy flux varies only with the friction velocity and not with the sea state. This allows a simple cubic relation with the waterside friction velocity The wave energy parameter CB is then treated as a constant of about 100 (Craig & Banner, 1994).
The near-surface mixing length is usually estimated using the roughness length which is thought to be determined by the sea state (Stacey, 1999). When sea state information is lacking, the roughness length is usually parameterized as where w ≈ 70, 000 is a Charnock coefficient for the water-side roughness length. When sea state information is available, the roughness can be estimated by where H s is the significant wave height, and w is in the range 0.85-1.6 (Law Chune & Aouf, 2018); here, we use w = 1.

Models
A regional fully coupled ocean/ice-wave-atmosphere system has been developed at Uppsala University, called UU-CM. At this stage, UU-CM includes four component models, that is, atmosphere, wave, ocean, and ice. The OASIS3-MCT (Valcke et al., 2013) is used as the coupler between the different components.
The exchanged variables are shown in Table 1. Each component is described in the following subsections.
• Atmospheric model: The Weather Research and Forecasting (WRF) Model (version 3.8) is used as the atmospheric component in the coupled system. WRF is a three-dimensional nonhydrostatic atmospheric model (Skamarock et al., 2005). The following schemes are used in this study. The Thompson scheme (Thompson et al., 2008) is used for the microphysics; the Rapid Radiative Transfer Model scheme (Mlawer et al., 1997) for the longwave radiation; Dudhia's scheme (Dudhia, 1989) for the shortwave radiation; Nakanishi and Niino PBLs surface layer scheme (Nakanishi & Niino, 2006) for the surface layer physics, and finally Grell-3 for the cumulus parameterization. By default, the roughness length in the WRF model is calculated using Here, U 10 is the 10-m wind speed. Initial and boundary conditions are taken from the ERA-Interim reanalysis (Dee et al., 2011). The horizontal resolution of the WRF model is 6 km with 31 vertical levels extending to 50 hPa. The WRF domain is shown in Figure 1 as the outer box.  Figure 1). NEMO has 56 vertical levels ranging from 1.5 m at the surface to around 22 m thickness for the lowermost model level. The following 11 tidal constituents are included in the NEMO configuration: M 2 , S 2 , N 2 , K 2 , K 1 , O 1 , P 1 , Q 1 , M 4 , MS 4 , and MN 4 . The Oregon State University tidal constituent database (Egbert & Erofeeva, 2002) and the climatological data (Janssen et al., 1999) for the temperature and salinity are used as the boundary condition.

Experiments
Four experiments were devised to investigate the impact of the OWA interaction processes on the coupled model. The exchange variables in the four experiments are listed in Table 2. The physical interaction processes at the air-sea interface are described in the following text.
• Exp1: Here, only the atmosphere (WRF) and wave models (WW3) are activated and there is only one-way forcing of the wave field by atmosphere (see Table 2). The boundary layer scheme of WRF yields 10-m wind vectors to WW3, but there is no feedback to the atmosphere model. The water level ( ) and the surface currents in equations (1) and (4) are here set to zero. • Exp2: Here the wave model is coupled to the atmospheric model and the feedback from the wave field is done by returning the sea state-dependent Charnock coefficient (equation (6)) to WRF (see Table 2). • Exp3: Two-way coupling to the ocean (NEMO) and sea ice model is activated in Exp3 (see Table 2). WRF gets sea surface temperature (SST) and surface currents from NEMO and returns forcing fields to NEMO and the ice component, which include wind stress, short-wave radiation, long-wave radiation, and net water flux. Here, the ocean-side stress is the same as the air-side stress, that is, oc = a . The TKE flux from breaking waves is estimated using equation (16), and the water-side roughness length is calculated using equation (17). • Exp4: In this experiment, in addition to the processes included in Exp3, the wave-current interaction processes listed in section 2 are switched on. The variables exchanged between WW3 and NEMO are shown in Table 1. WW3 calculates the refraction by NEMO surface currents and adjusts the wave propagation speed in accordance with the water level, in which the impact of tide-induced currents and sea level height on waves is also included in the system. The ocean-side stress is calculated following equation (20). The Stokes drift-related processes listed in section 2.3 are implemented in NEMO. The breaking wave-induced TKE flux and the TKE roughness length are calculated following equations (15) and (18), respectively. Note. The mean wind speed are from stations A1 to A9, the significant wave height data are from station W1 to W12, and the sea surface temperature data are from station O1 to O11.
Exp1 is the control experiment against which the other experiments are compared to assess the impact of atmosphere-wave, atmosphere-ocean, and wave-ocean interaction on the simulations, respectively. Two-month-long simulations are performed in this study; a winter month (January 2015) and a summer month (July 2015). The initial conditions for the ocean model are taken from the restart file in the control simulation in Wu et al. (2019). The winter simulation is preceded by a 1-year spin up, while a 17-month spin up precedes the summer simulation. The following analysis is done on 3-hourly output. In general, the impact of the coupling on surface variables can be seen in a very short time (a few hours) whereas the oceanic vertical structure takes longer time (a few days) to adjust. In this study, we mainly focus on the coupling impact on the surface variables.

Data
The impact of OWA interaction processes on the simulations is the focus of this study. We compare against several observational data sources, including remote sensing data (IFREMER CERSAT Global Blended Mean Wind Fields), ERA5 reanalysis data, and in situ measurements. These data sets are briefly introduced in the following.

Remote Sensing Data
IFREMER CERSAT Global Blended Mean Wind Fields are products of remotely sensed surface wind from several data sources, that is, surface wind derived from scatterometers onboard ASCAT-A and ASCAT-B, the SSMIS radiometers, and the WindSat radiometer. The data set covers the global oceans with 0.25 • spatial resolution and 6-hourly temporal resolution.

Reanalyses
ERA5 is the fifth-generation reanalysis by the European Centre for Medium-Range Weather Forecasts, see Hersbach and Dee (2016). ERA5 archives atmospheric parameters with hourly resolution on 31-m horizontal resolution at 137 vertical levels to 0.01 hPa. The ERA5 data were downloaded from this site (https://cds. climate.copernicus.eu). The wave model ECWAM is coupled to the atmospheric model (ECMWF, 2017).

In Situ Measurements
The locations of in situ stations used are shown as dots in Figure 1 and listed in Table 3. The hourly mean wind speed at 9 stations, hourly H s from 12 stations, and hourly SST from 11 stations are compared with our model runs. The data were downloaded from this website (http://marine.copernicus.eu).

Impact of Coupling on the Atmospheric Model
The relationship between the Charnock coefficient and the mean wind speed calculated using equation (19) (the red line) and equation (6) is shown in Figure 2a. The colors represent the normalized distribution (%) of the Charnock coefficient in the same wind bin (1 m/s). The Charnock coefficient calculated from the WW3 model is lower than the one calculated using equation (19) for the condition with U 10 < 4 m/s and U 10 > 15 m/s. Similarly, the drag coefficient from WW3 is lower than the one from the stand-alone WRF model when U 10 > 15 m/s (see Figure 2b). For the other wind conditions, the Charnock coefficient calculated from equation (19) is lower than the one calculated from the WW3 model. Accordingly, the C d is bit higher in the coupled system. One needs to note that the Charnock coefficient calculated from wave spectrum spreads out at the same wind speed. Thus, the roughness length (Charnock coefficient) varies with the wave spectrum even if the mean wind speed is the same.
The impact of the OWA coupling on the mean wind speed (U 10 ) is shown in Figure 3 for January and July 2015. Making the Charnock coefficient sea state dependent in Exp2 reduces the mean U 10 in most areas both in January and July. This is partly caused by the Charnock coefficient and drag coefficient from Exp2 being larger than that from Exp1 for the wind range 4-15 m/s (see Figures 2a and 2b). However, in the northern part of the North Sea, the wave feedback on the atmosphere increases the mean wind speed in January. This is caused by the indirect feedback due to changes in the Charnock coefficient in the coupling domain. The influence of the atmosphere-wave coupling on the mean wind speed is less than 0.2 m/s both in January and July. Adding the atmosphere-ocean coupling processes has a significant impact on the mean wind speed (more than 0.6 m/s in most areas). It decreases the mean wind speed in January in most of the model domain but has the opposite effect along the south coast of Norway and in the Gulf of Finland and the Eastern part of the North Sea. In July, the atmosphere-ocean coupling increases the mean wind simulation in the coastal area ( Figure 3g). This is caused by the higher-resolution SST from the ocean model as well as the influences of current on the wind stress calculation. The wave-current interaction processes (Exp4) increase the mixing in the upper ocean and change the wave simulations of roughness length accordingly. Those processes lead to higher wind speed in January and lower mean wind speed in July (Figures 3d and 3h). Figure 2c shows the

Journal of Advances in Modeling Earth Systems
10.1029/2019MS001761 Figure 3. The impact of ocean-wave-atmosphere coupling processes on the mean 10-m wind, U 10 (m/s), for January (a-d) and July (e-h). The first column represents the mean wind from the control experiment (Exp1). The second to fourth columns represent the differences between experiments and the control experiment, that is, Exp2-Exp1, Exp3-Exp2, and Exp4-Exp3.
normalized distribution of U 10 at different wind speeds (1 m/s bins) for the entire domain. In January, the OWA coupling processes shift the wind distribution slightly down while the opposite occurs in the July run.
In January, the mean U 10 from CERSAT is higher than that from ERA5. The simulation of mean U 10 from the control experiment (Exp1) is larger (smaller) than that from ERA5 (CERSAT) in most of the domain in January (Figures 4a and 4b). In July, Exp1 has a higher mean wind speed compared with ERA5 and CERSAT (Figures 4d and 4e). In the coastal areas, the results from Exp1 overestimate (underestimate) the mean wind speed compared with both ERA5 and CERSAT in January (July). In general, adding the OWA coupling processes (Exp4) reduces the overestimation (underestimation) in the coastal areas in January (July). The low-resolution SST in Exp1 cannot resolve the SST structures in some areas (e.g., along the south coast of Norway and the Gulf of Finland, see Figure 9), which leads to a significant difference in the mean wind speed (Figures 4c and 4f). However, comparing with ERA5, the OWA processes (Exp4) enhance the overestimation of the mean wind speed along the southern coast of Norway in January. In addition to the SST differences, the strong current-wave and current-atmosphere interactions in those areas also contribute to the enhancement due to the strong currents in those areas.
The comparison between the in situ measurements and the model results are shown in Figure 5. Three metrics are used to compare the in situ measurements (the stations are shown in Figure 1): mean absolute error (MAE), root-mean-square deviation (RMSD), and correlation (R). Adding OWA coupling processes (Exp4) improves the simulation in some stations, but worsens the simulations in other stations. In general, adding OWA coupling processes (Exp4) does not have a significant influence on those stations in terms of the three metrics studied in both January and July. The comparison between the model results with in situ measurements can be affected by many factors such as model resolution, the quality of measurements, the station location and the number of stations used in the comparison. Since the magnitude of the physical/dynamical impact of the OWA processes on the simulations is the main aim of this study, rather than an attempt to build a model system ready for forecasting. Coupled systems usually require some tuning of the individual model components  to bring down biases. This has not been done in this study, and it is thus perhaps not so surprising that the coupled system does not significantly improve reduce biases compared with the ocean-only tuned runs.
There is no significant difference between Exp1 and Exp2 in terms of 2-m air temperature (T 2 ) either in January or July (not shown). Compared with Exp3, adding wave-current interaction processes decreases T 2 less than 0.4 • C in July, but does not have a significant influence in January (not shown). The mean T 2 difference between the ERA5 and the simulation Exp1 is shown in Figures 6a and 6c for January and July, respectively. Compared with ERA5, the control experiment (Exp1) underestimates T 2 (by less than 1 • C) in most areas of the Baltic Sea (except in the northern part of the Baltic Sea) in January. In July, the control experiment underestimates T 2 up to 3 • C in the coastal areas. Adding OWA interaction processes enhances the underestimation of the T 2 compared with the results from Exp1 in January except some slight positive impact on the north of the Baltic Sea. But, the OWA interaction reduces significantly the bias in the coastal areas in July and increases the coastal air temperature by more than 1 • C (and by more than 2 • C in the north of the Baltic Sea). The coastal upwelling induced lower SST in the ocean model ( Figure 10) and consequently lower air temperature along the southeastern part of the Gotland and Öland Islands (Figure 6d).

Impact of Coupling on the Wave Model
The mean H s in January is around 2 m in the Baltic Sea and more than 2.5 m in most areas of the North Sea (Figure 7). In contrast, the mean H s in July is much lower (less than 1.4 m in both the Baltic and North Sea) since the wind speed in July is much lower than in January (Figures 2c and 3a-3e). The OWA coupling affects the mean H s , as shown in Figure 7. The general difference pattern between the different experiments is quite similar to the impact on the mean wind field (Figure 3) since waves are mainly determined by the wind. Compared with the atmosphere-wave coupling process and wave-current coupling processes, the atmosphere-ocean coupling processes have a larger impact on the wave height. This is because the higher resolution SST and current feedback on the atmosphere increases the wind, which in turn increases the significant wave height. In contrast to the wind field, the wave-current interaction has a significant influence on the simulation of waves compared to that from the atmosphere-wave coupling (see Figures 7d and 7h). It is clear that the wave-current interaction processes have a significant influence in the high-current areas, such as along the south coast of Norway and in the English channel. In those areas, the wave model is more sensitive to changes to the current field than to changes in the wind field. Since the wave frequency can be shifted by the current due to Doppler effect, the wave height will decrease (increase) when the currents and waves travel in the same (opposite) direction. In addition, the water level changes will also affect the wave field in shallow areas. The comparison between the wave model results and observations of H s is shown in Figure 8. In January, Exp4 (OWA coupling) is worse than Exp1 in terms of MAE and RMSD but it improves the correlation slightly.
In July the results are mostly neutral with a small tendency for Exp4 to improve the MAE, RMSD, and R as well as mean error (not shown). Figure 9 shows the SST at time 2015-07-03T00:00 UTC. Due to the low resolution of the SST in ERA-Interim (79-km resolution), coastal areas are not well resolved. In summer, upwelling events are very common along the coastlines of the Baltic Sea. This leads to very different SSTs compared with the conditions found in the central basin (Lehmann et al., 2012). In Figure 9, we see that NEMO can capture the lower SST caused by upwelling in the coastal areas. These upwelling events obviously affect atmospheric stability and heat flux (Sproson & Sahlee, 2014) and thus they will affect the local atmospheric dynamics in terms of wind speed and air temperature (see Figures 3 and 6).

Impact of Coupling on the Ocean Model
The mean SST in January and July from ERA5 and Exp3 are shown in Figure 10. The mean SST from Exp3 is a little lower (higher) than that from ERA5 in January (July). The wave-current interaction processes increase (decrease) the SST in January (July) due to the enhanced mixing by surface waves. Accordingly, the mixed layer depth is increased by the wave processes in July. In January, the mixed layer depth does not change appreciably from Exp3 to Exp4. The wave-current interaction impact on the SST is quite small in January except along the south coast of Norway (the SST increased by more than 0.2 • C). In contrast, the wave-current interaction has a significant influence on the SST in July, which reduces the SST more than 0.5 • C in the areas away from the coast (see Figure 10). But, in some coastal areas, the wave-current interaction processes increases the SST, such as the Swedish east coast with high upwelling frequency areas. These results agree with the results of Wu et al. (2019) where it was found that the wave-current interaction processes reduce the upwelling intensity in simulations using a wave-forced (uncoupled) ocean model.
The SST comparison between the simulation results and in situ measurements in 11 stations are shown in Figure 11. One can see that the OWA coupling reduces the MAE and RMSD and increases the correlation in some stations, but worsen the simulation in others. The OWA coupling has a greater influence on the summer simulation in terms of SST. In general, the OWA coupling processes do not yield a significant reduction of biases for the model setup used in this study.  6. Discussion

Momentum Balance at the Air-Sea Interface
The ratio between the ocean-side and air-side stress is a function of the sea state, and thus of fetch and duration of the wind. The ratio typically falls between 0.8 and 1.2, but can in extreme cases be both lower and higher under different sea states (i.e., Alari et al., 2016;Wu et al., 2019). Instead of using equation (8) to calculate the ocean-side stress, several studies estimate the ocean-side stress using the simple ratio (e.g., Breivik et al., 2015) (20) on the assumption that the direction of the air-side stress and the ocean-side stress are the same. However, we know that the wave-supported stress, in , and the momentum flux from waves to currents, ds , are determined by the wave spectrum and will in general not coincide exactly with the local wind stress direction, especially in the presence of swell and when the wind is turning. As shown in the study by Semedo et al. (2011), swell waves usually dominate in the extratropical oceans. This is also, albeit to a lesser extent, the case for the North Sea (Semedo et al., 2015). Nonlocal wave systems (swell) and veering winds alter the direction of in and ds . Figure 12 shows one snapshot of the stress direction differences between the airside stress and the other stress at 2015-07-03T00:00 UTC. The distribution of wave age, C p ∕U 10 , is shown in Figure 12a. One can see that there are some high wave age areas along the south part of Gotland and the middle part of the North Sea. In those high wave age areas, the peak direction difference between the mean wind and the peak wave is significant. In some areas, the difference is larger than 90 • (see Figure 12b). The directional difference between a and ds is shown in Figure 12c. In some high wave age areas, the direction difference between a and ds is significant (larger than 80 • in some areas), mainly caused by the contribution of the swell systems. In coastal areas, the directional difference between a and ds is significant. This is mainly caused by the shallow-water bottom interaction processes. Figure 12d shows the directional difference between a and in in the Gulf of Finland. The directional difference is more than 80 • along the north coast. In this case, the wind direction in the Gulf of Finland is parallel with the coastline to the east. The wave direction is more northward than the wind direction. Thus, the influence of the topography and wave reflection by the coastline makes the wave field more complex in the coastal areas, in particular along the north coast. Like the coupled model used by Pianezze et al. (2018), WW3 calculates its own friction velocity from the WRF wind speed through a bulk formulation valid for neutral atmospheric conditions. This causes a slight imbalance between the two models. However, as stated by Pianezze et al. (2018), the difference between the WRF and WW3 friction velocities is generally only a few percent and about 10% under high wind speed. This momentum flux imbalance will be investigated further in a future study.   The physical processes are complex at the air-sea interface under high wind conditions. Sea spray generated by intensive wave breaking can significantly alter the heat and momentum fluxes. It is also one of the possible reasons why the coupled system only has a slightly impact on the simulations in January where there are several low-pressure systems passing through the domain generating high winds. How to implement the fluxes under high wind conditions in the WW3 (wind stress) and WRF model will also be investigated in a future study.

Wave-Current Interaction
In our coupled system the Stokes drift enters the wave-averaged momentum equation through the Coriolis-Stokes forcing (CSF) and the Stokes advection. The Stokes advection is also present in the tracer The surface current, SST, and oceanside stress in the area shown as the red box in Figure 1 at 2015-07-03T00:00 UTC is shown in Figures 13a, 13c, and 13e, respectively (from Exp3). The wave-current interaction processes impact on those parameters are shown in the right column in Figure 13 (the difference between Exp4 and Exp3). At this time step, the surface current is strong (more than 0.25 m/s east of the Stockholm area and north and west of Gotland). Adding wave-current interaction (Exp4) reduces the current in the high current areas, in some areas by more than 30% (Figure 13b). It also reduces the SST by more than 1 • C in some areas. The oceanside stress and the wave-current impact on the stress are shown in Figures 13e and 13f. As can be seen the wave-current interaction has some influence on the stress but less than 10% and only in a quite limited area. This snapshot is representative of the situation throughout the month-long simulations in that the main influence from wave-current interactions is found in areas with high currents and stress. It thus seems likely that the difference in SST and currents are related not only to the stress but also to the Stokes drift (shown in Figure 14a).
Background currents shift the wave frequency through the Doppler effect (Uchiyama et al., 2009). The circulation in the simulation domain is complex with current and waves variously opposing and aligning with each other. In coastal areas, the current-wave interaction can modify the bottom friction (Signell et al., 1990). These processes affect the wave model simulation and indirectly the entire coupled system, as shown in Figure 7. Figure 14 shows the wave-current interaction impact on the waves at 2015-07-03T00:00 UTC. The upper panel shows the surface Stokes drift, significant wave height, peak wave period, and the Charnock coefficient, respectively (from Exp3). The bottom panel shows the wave-current interaction impact on those parameters. It is clear that the wave-current interaction changes the surface Stokes drift by more than 10% in some areas (the wave-current interaction can change the monthly average surface Stokes drift by up to 8% in coastal areas in January and by more than 12% in June). The influence of wave-current interaction on the significant wave height shows a similar pattern as the Stokes drift. Due to the Doppler effect, the peak wave period has changed significantly in some coastal areas. But the spatial pattern is different from the significant wave height and the surface Stokes drift. From equation (6) it is clear that the Charnock coefficient can be altered by any process affecting the wave spectrum. Here the influence pattern of wave-current on the Charnock coefficient is more complex (see Figure 14).
For completeness we note again that we have not investigated the impact of Langmuir turbulence in this study. It is clear from implementations of the Langmuir enhancement factor (McWilliams & Sullivan, 2000) in K-profile parameterizations of upper-ocean turbulence that the effect may be significant, in particular, in altering the mixed-layer thickness of coupled climate models Li et al., 2016Li et al., , 2017. Other effects not investigated here are the nonbreaking wave-induced upper ocean mixing proposed by Qiao et al. (2004), the orbital-motion induced mixing proposed by Babanin (2006) and the radiation stress (Dietrich et al., 2011).

Influence of the Horizontal Resolution of the SST Boundary Conditions
Two more experiments (Exp1h and Exp2h) were devised in order to investigate the impact of the resolution of the SST boundary conditions. In these two experiments, the SST fields from Exp3 are used as boundary

10.1029/2019MS001761
conditions. Otherwise, the experiments are identical to Exp1 and Exp2. It should be noted that the two experiments represent a somewhat unrealistic stand-alone setup since the SST in standalone models are usually taken from a realistic data set (not from a coupled simulation), such as a reanalysis (as we do in Figure 3). Figure 15 shows the impact of the atmosphere-wave coupling (left column) and atmosphere-ocean coupling (right column) on the mean 10-m wind speed. The atmosphere-wave coupling impact on the mean wind speed (Figures 15a and 15c) is of the same order as the difference between Exp2 and Exp1. The spatial pattern is a bit different in January due to the indirect feedback from the SST difference. It is not surprising that the atmosphere-ocean coupling influence is smaller with higher SST boundary conditions, compared with the one using low-resolution ERA-Interim SST (Figures 3c and 3g). In particular, the impact is smaller in coastal areas where the high-resolution SST field captures coastal features (e.g., coastal upwelling) similar to those found in the coupled model. However, the atmosphere-ocean coupling due to diurnal SST patterns and current feedback (Exp3-Exp2h) can still have a significant influence on the results. The difference between Exp2h and Exp1h (Exp3-Exp2h) in terms of T 2 and H s exhibit similar patterns as the U 10 (not shown).
The comparison shows that the coupled and uncoupled with the NEMO model are clearly affected by the coarse resolution of the ERA-Interim boundary conditions. These four experiments (Exp1 to Exp4) demonstrate the impact of going from an uncoupled system with coarse resolution boundary conditions to a fully coupled system. The final two experiments demonstrate that even when a high-resolution (albeit unrealistic) boundary field is used for the stand-alone runs (Exp1h and Exp2h), significant differences persist compared with the fully coupled Exp3 (from which the boundary fields are taken).

Conclusions
In this study, a fully coupled ocean-wave-atmosphere system, UU-CM, is developed. In the coupled system, the momentum flux balance was built considering several different stress terms, that is, airside stress, oceanside stress, wave-supported stress, and stress from waves to currents. The airside stress is estimated from a sea state-dependent Charnock coefficient. The surface wave impact on the upper ocean was added in the coupled system in terms of wave-breaking-induced TKE flux, CSF, Stokes drift advection, and the Stokes drift impact on the mass transport. The current influence on the waves and the airside stress were also added to the system. Four sensitivity experiments were designed to investigate the impact of these OWA interaction processes on the simulations with a domain covering the Baltic and North Seas.
Based on simulations of two months (January and July) in 2015, we find that the atmosphere-wave coupling (sea state-dependent Charnock coefficient) reduces the mean U 10 and H s slightly, which is due to the lower drag coefficient under moderate wind conditions (4 < U 10 < 15 m/s) from the atmosphere-wave coupling model ( Figure 2). The atmosphere-ocean coupling processes as well as the high-resolution SST from NEMO reduce the mean wind speed and H s throughout the domain in January, but increases the mean wind speed and H s in the coastal areas in July. This is mainly because the SST from the coupled system is in general higher (lower) than the one from ERA-Interim in July (January), thus altering the wind and in turn the wave field. The high resolution of the SST in the coupled system captures the influence of local systems such as coastal upwelling, and their influences on the atmosphere. Wave-current coupling processes increase (decrease) the mean U 10 and SST in January (July). In the summer, the wave-induced upper ocean mixing is strong, which leads to a significant influence (decrease) on the SST and accordingly decreasing wind and wave height. However, the ocean-current interaction has a very small impact on the SST in January. The main influence from wave-current interactions on waves may be caused by the direction difference between current and waves and lead to the decrease of H s in January.
Compared with ERA5 and remote sensing data, the coupled system improves the simulation results in terms of mean wind, air surface temperature, and SST. However, the coupled system does not yield a significant reduction in bias compared with the stand-alone control experiment.