The effects of ocean surface waves on global intraseasonal prediction: case studies with a coupled CFSv2.0–WW3 system

This article describes the implementation of a coupling between a global forecast model (CFSv2.0) and a wave model (WW3) and investigates the effects of ocean surface waves on the air–sea interface in the new framework. Several major wave-related processes, including the Langmuir mixing, the Stokes–Coriolis force with entrainment, air–sea fluxes modified by the Stokes drift, and momentum roughness length, are evaluated in two groups of 56 d experiments, one for boreal winter and the other for boreal summer. Comparisons are made against in situ buoys, satellite measurements, and reanalysis data to evaluate the influence of waves on intraseasonal prediction of sea surface temperature (SST), 2 m air temperature (T02), mixed layer depth (MLD), 10 m wind speed (WSP10), and significant wave height (SWH). The wave-coupled experiments show that overestimated SSTs and T02s, as well as underestimated MLDs at mid-to-high latitudes in summer from original CFSv2.0, are significantly improved due to enhanced vertical mixing generated by the Stokes drift. For WSP10s and SWHs, the wave-related processes generally reduce biases in regions where WSP10s and SWHs are overestimated. On the one hand, the decreased SSTs stabilize the marine atmospheric boundary layer and weaken WSP10s and then SWHs. On the other hand, the increased roughness length due to waves reduces the originally overestimated WSP10s and SWHs. In addition, the effects of the Stokes drift and current on air–sea fluxes also rectify WSP10s and SWHs. These cases are helpful for the future development of the two-way CFSv2.0–wave coupled system.


Introduction
Ocean surface gravity waves play an important role in modifying physical processes at the atmosphere-ocean interface, which can influence momentum, heat, and freshwater fluxes across the air-sea interface (Li and Garrett, 1997;Taylor and Yelland, 2001;Moon et al., 2004;Janssen, 2004;Belcher et al., 2012;Moum and Smyth, 2019). For instance, ocean surface waves modify ocean surface roughness to influence the marine atmospheric boundary layer and thus change the momentum, latent heat, and sensible heat transfer (Janssen, 1989(Janssen, , 1991Taylor and Yelland, 2001;Moon et al., 2004;Drennan et al., 2003Drennan et al., , 2005. The breaking waves inject turbulent kinetic energy into the upper ocean, which enhances the mixing process (Terray et al., 1996). Nonbreaking surface waves also affect mixing in the upper ocean by adding a wave-related Reynolds stress (Qiao et al., 2004;Ghantous and Babanin, 2014). The wave-related Stokes drift interacts with the Coriolis force and produces the Stokes-Coriolis force (Hasselmann, 1970). The shear of the Stokes drift is critical for generation of Langmuir circulation, which significantly deepens the mixed layer by strong vertical mixing both at climate scales (Li and Garrett, 1997;Belcher et al., 2012) and at weather scales (Kukulka et al., 2009).
Various wave-related parameterizations have been proposed and used in modeling. The wave-related Charnock parameter (C ch ) defines sea surface roughness and affects wind stress estimates (Pineau-Guillou et al., 2018;Sauvage et al., 2020). There are primarily three methods for defining C ch , considering the wave-induced kinematic stress (Janssen, 1989(Janssen, , 1991, the wave age (Drennan et al., 2003(Drennan et al., , 2005Moon et al., 2004), or the steepness (Taylor and Yelland, 2001). The first two are based on the wind-sea conditions, whereas the latter includes both swells and wind-sea waves. Modifications to these Charnock parameterizations were suggested in recent studies for leveling off roughness under high winds (e.g., Fan et al., 2012;Bidlot et al., 2020;ECMWF, 2020;Li et al., 2021). In the oceanic boundary layer, waves influence upper ocean mixing via wave dissipation and the Stokes drift-induced processes. In Breivik et al. (2015), the wave dissipation-related turbulent kinetic energy flux is found to yield the largest sea surface temperature (SST) differences in the extratropics. The Stokes drift-induced Langmuir turbulence can improve temperature simulation over most of the world oceans, particularly in the Southern Ocean (Belcher et al., 2012;Li et al., 2016). Polonichko (1997), , and Li et al. (2017) indicated that the Langmuir cell intensity strongly depends on the alignment of winds and waves, reaching a maximum when they are aligned. Li et al. (2016) found that the effect of Langmuir cell can be further enhanced by entrainment. In Couvelard et al. (2020), the Stokes drift-related forces can also contribute modestly to deepening of the mixed layer depth (MLD). In the First Institute of Oceanography Earth System Model, Bao et al. (2019) indicated that the nonbreaking wave-induced mixing, the Stokes drift-affected air-sea fluxes, and sea spray are all important for climate estimates.
The wave-related processes at the air-sea interface are complex and important in global coupled systems (e.g., Breivik et al., 2015;Law-Chune and Aouf, 2018;Bao et al., 2019;Couvelard et al., 2020). Most of the coupled models with a wave component at global scale were developed for climate research (e.g., Law-Chune and Aouf, 2018; Bao et al., 2019;Couvelard et al., 2020). Exceptionally, an Integrated Forecasting System (IFS) with fully coupled atmosphere, ocean, and wave components, developed by European Centre for Medium-Range Weather Forecasts (ECMWF) (Janssen, 2004;Bidlot, 2019;Bidlot et al., 2020), has been assembled for global forecasts from medium-range weather scales to seasonal scales (Breivik et al., 2015). The National Centers for Environmental Prediction (NCEP) is establishing its atmosphere-ocean-wave system, in which the Global Forecast System (GFS; the atmosphere module in the Climate Forecast System model version 2.0) is one-way coupled with WAVEWATCH III (WW3).
The effects of wave-related processes are worth further evaluation in different global coupled modeling systems. Since it takes significant time for the wave energy to develop (Janssen, 2004), we investigate the impact of individ- Figure 1. A schematic diagram of the atmosphere-ocean-wave coupled modeling system. The arrows indicate the coupled variables that are passed between the model components. In the diagram, C ch , La t , u s (0), V s , U 10 , and U surf are Charnock parameter (red arrows), turbulent Langmuir number (red arrows), surface Stokes drift velocity (red arrows), Stokes transport (red arrows), 10 m wind (green arrows), and surface current (blue arrows), respectively.
ual wave processes at intraseasonal timescale in a new global atmosphere-ocean-wave system. To achieve this, we couple WW3 to the Climate Forecast System model version 2.0 (CFSv2.0) and then conduct sensitivity experiments in boreal winter and summer for comparison. The effects of upper ocean mixing modified by Langmuir cell, the Stokes-Coriolis force and entrainment, air-sea fluxes modified by surface current and the Stokes drift, and momentum roughness length are evaluated. The CFSv2.0 is a coupled system mainly applied for intraseasonal and seasonal prediction (e.g., Saha et al., 2014). Our work can provide insights into two-way wave coupling of CFSv2.0 and is helpful for the future development of the CFSv2.0-wave coupling system. Two groups of 56 d predictions are conducted for boreal winter and boreal summer. Then, the predictions are compared with observations and reanalysis data. For each group, sensitivity experiments with different wave parameterizations are carried out to evaluate the effects of individual wave-related process.
The rest of the paper is structured as follows: methods and numerical experiments with different parameterizations are described in Sect. 2; the observations and reanalysis data are introduced in Sect. 3, and the results of experiments are evaluated and compared in Sect. 4. Finally, a summary and discussion are given in Sect. 5.  Griffies et al., 2004) as the ocean component. MOM4 is integrated on a nominal 0.5 • horizontal grid with horizontal resolution enhanced to 0.25 • in the tropics and has 40 vertical levels; the vertical spacing is 10 m in the upper 225 m, and it then increases in unequal intervals to the bottom at 4478.5 m. A three-layer sea ice model is included in MOM4 (Wu et al., 2005). GFS uses a spectral triangular truncation of 382 waves (T382) in the horizontal, which is equivalent to a grid resolution of nearly 35 km, and 64 sigma-pressure hybrid layers in the vertical. The time steps of both MOM4 and GFS are 180 s. The ocean and atmosphere components are then coupled at the same rate. In the original two-way coupled system, GFS receives SST from MOM4 and sends fluxes of heat, momentum, and freshwater to MOM4 (black arrows in Fig. 1).
The Chinese Community Coupler version 2.0 (C-Coupler2; Liu et al., 2018) is used to interpolate and pass variables between atmosphere and wave components as well as ocean and wave components. Each component receives inputs and supplies outputs on its own grids. The C-Coupler2 is a common, flexible, and user-friendly coupler, which contains a dynamic 3-D coupling system and enables variables to remain conserved after interpolation.
A schematic diagram of the coupled atmosphere-oceanwave system is shown in Fig. 1. As illustrated, WW3 is twoway coupled with MOM4 and GFS, through the C-Coupler2. WW3 is forced by 10 m wind from GFS (green arrows) and surface current from MOM4 (blue arrows) and generates the wave action density spectrum. Meanwhile, the surface Stokes drift velocity, the Stokes transport, and the turbulent Langmuir number are passed to MOM4 (red arrows; see Sect. 2.3) from WW3, and the surface Stokes drift velocity and the Charnock parameter are passed to GFS (red arrows; see Sect. 2.4 and 2.5). The high-frequency tail assumption for the Stokes drift in WW3 is used with a spectral level decaying as f −5 (frequency). Additionally, the regular ocean surface current velocities from MOM4 are also passed to GFS to calculate the relative wind velocity for the turbulent fluxes together with the surface Stokes drift (blue arrows; see Sect. 2.4).
Both CFSv2.0 and WW3 use warm starts; the initial fields at 00:00 UTC of the first day in each experiment for CFSv2.0 were generated by the real-time operational Climate Data Assimilation System (Kalnay et al., 1996), downloaded from the CFSv2.0 official website (http://nomads.ncep.noaa.gov/ pub/data/nccf/com/cfs/prod, last access: 20 January 2022). To get initial conditions for WW3, a stand-alone WW3 model is set up synchronously (see Sect. 2.2). Since the interactions between waves and sea ice are complicated and beyond the scope of the study, we turn off the coupling between WW3 and CFSv2.0 in areas with sea ice.
In addition, to properly select the coupling frequency between CFSv2.0 and WW3, the root-mean-square errors (RMSEs) of SST, significant wave height (SWH), and 10 m wind speed (WSP10) with different coupling steps for the fully coupled experiment (ALL; details in Sect. 2.6) are calculated and compared (Table S1 of the Supplement). The three components are coupled every time step (180 s) in the 1_STEP_ALL experiment, every 5 steps (900 s) in the 5_STEP_ALL experiment, and every 10 steps (1800 s) in the 10_STEP_ALL experiment. In 10_STEP_WW3, only WW3 is coupled every 10 time steps, whereas GFS and MOM4 remain at the one-time-step (180 s) coupling frequency as in the original settings in CFSv2.0. From Table S1, the 10_STEP_WW3 experiment has a relatively short runtime and small RMSEs. Therefore, the time steps of 10_STEP_WW3 are selected to compromise between computing time and model RMSEs.

Initialization of WAVEWATCH III
In WW3, input of momentum and energy by wind and dissipation for wave-ocean interaction are two important terms (combined as input-dissipation source term) in the energy balance equation (WAVEWATCH III Development Group, 2016), which includes the estimation of the Charnock parameter. Several different packages to calculate the inputdissipation source term (ST) are available in the WW3 version 5. 16, including ST2 (Tolman and Chalikov, 1996), ST3 (Janssen, 2004;Bidlot, 2012), ST4 (Ardhuin et al., 2010), and ST6 (Zieger et al., 2015).
The initial wave fields are generated from 10 d simulation starting from rest in a stand-alone WW3 model. To minimize the biases of initial wave fields, we test simulations with ST2, ST3, ST4, and ST6 schemes and compare the results with Janson-3 observations. Two 10 m wind datasets, the Cross-Calibrated Multi-Platform (CCMP; Atlas et al., 2011)   ment), consistent with findings in Stopa et al. (2016). Thus, the ST4 scheme is chosen to calculate the input and dissipation term and to generate initial wave fields with ERA5 wind forcing for experiments listed in Table 1. The parameters used for the ST4 scheme follow TEST471f from WAVE-WATCH III Development Group (2016), which is the CFSR (CFS Reanalysis) tuned setup and is commonly used at the global scale.

Parameterizations of the Stokes drift-related ocean mixing
The full Stokes drift profile used in MOM4 is obtained by the method by Couvelard et al. (2020), which is based on the work by Breivik et al. (2014Breivik et al. ( , 2016. Breivik et al. (2016) derived the full Stokes drift profile as where u s (0) is the surface Stokes drift velocity, k p = u s (0) 6V s , V s is the Stokes transport, and erfc is the complementary error function. Equation (1) is depth-averaged within each vertical grid interval as where "th" is the thickness of layer k, following Li et al. (2017)

Mixing of Langmuir turbulence
McWilliams and Sullivan (2000) modified the turbulent velocity scale W in K-Profile Parameterization (KPP) for vertical mixing by introducing an enhancement factor ε, to account for both boundary layer depth changes and nonlocal mixing by Langmuir turbulence. Based on their work, Van Roekel et al. (2012) improved the enhancement factor corresponding to alignment and misalignment of winds and waves. Li et al. (2016) evaluated these parameterizations in a coupled global climate model and found that the difference between parameterizations with alignment and with misalignment was not significant, owing to the relatively coarse resolution which cannot accurately represent the refraction by coasts and current features. We use the parameterization from Van Roekel et al. (2012) as well. Because the resolution in our model is relatively coarse too and the angles between winds and waves are less than 30 • in most areas ( Fig. S1i and j in the Supplement), we do not consider misalignment in the study. W (W = ku * /φ, where u * is the surface friction velocity, φ is the dimensionless flux profile, and k = 0.4 is the von Kármán constant) depends on the turbulent Langmuir number; that is, where La t is the turbulent Langmuir number, defined as with u s (0) being the surface Stokes drift velocity. Furthermore, the enhanced W will influence the calculation of boundary layer depth. In KPP the boundary layer depth is determined as the smallest depth at which the bulk Richardson number equals the critical value Ri cr = 0.3; that is, where g is acceleration due to gravity, ρ is density, u is velocity, ρ r is surface density, u r is surface velocity, ρ 0 is the average value of the density, and h is the boundary layer depth. Hence, when W is enhanced, the boundary layer depth h is deepened accordingly.

The Stokes-Coriolis force and associated entrainment
Because the Stokes drift velocity is an increment superimposed on the original current velocity, the Coriolis force Eqs.
(1)- (6), (8) and the Stoke drift together produce an additional so-called Stokes-Coriolis (SC) force (Hasselmann 1970); that is, Here u s (z) is the Stokes drift velocity vector, f is the Coriolis frequency, and z is the vertical unity vector. For consistency, the Stokes drift velocity is also included in advection terms of tracers (e.g., temperature, salinity) and convergence terms (Law-Chune and Aouf, 2018; Couvelard et al., 2020). And the free surface condition for barotropic mode is correspondingly modified to where η is surface elevation, and M curr and M st are the total vertical integral of the regular Eulerian current and the Stokes drift, respectively. To depict the entrainment below the ocean surface boundary layer induced by the Stokes drift, Li et al. (2016) suggested adding the square of the surface Stokes drift velocity (|u s (0)| 2 ) to the denominator of Eq. (7); that is, The boundary layer depth h in KPP from Eq. (10) is then enhanced due to the Stokes drift velocity.

The Stokes drift and sea surface current on air-sea fluxes
At the air-sea boundary layer, the momentum flux (τ ), sensible heat flux (SH), and freshwater flux (E) are calculated as where C d , C h , and C e are surface exchange coefficients for momentum, sensible heat, and freshwater. ρ a is air density. θ and q are potential temperature and humidity differences between air and sea, and V is velocity of air relative to water flow.
In CFSv2.0, V is set to be wind speed (U wind ). However, the effect of ocean surface current should not be ignored. Luo et al. (2005) first indicated that including ocean surface current (U surf ) improves estimates of τ and subsequent ocean response. Renault et al. (2016) further indicated that the improvements of τ by U surf also fed back into the atmosphere. At present, V = U wind − U surf is widely used in coupled ocean-atmosphere models (e.g., Hersbach and Bidlot, 2008;Takatama and Schneider, 2017;Renault et al., 2021). Furthermore, Bao et al. (2019) indicated that as a part of the sea surface water movement with speed magnitude comparable to surface current in mid-to-high latitudes, the surface Stokes drift (u s (0)) should also be included; that is, To account for the effects of the surface currents and of the Stokes drift, Eq. (14) is used in the coupled experiments (Table 1). To complete the coupling, the corresponding modification of the tridiagonal matrix (Lemarié, 2015) has been implemented in CFSv2.0. Note that the direction of Stokes drift is generally consistent with 10 m wind ( Fig. S1i and j in the Supplement), but the directions of surface current and 10 m wind are usually different due to the Coriolis effect ( Fig. S1g  and h). Consequently, the effects of U surf and u s (0) on V depend on the angles between them and U wind .

Parameterizations of momentum roughness
In CFSv2.0, the fluxes of momentum, heat, and freshwater are passed from atmosphere to ocean, and their estimate is critically important. The fluxes are in part determined by surface roughness length, which can be converted to surface exchange coefficients based on Monin-Obukhov similarity theory (Monin and Obukhov, 1954).

The momentum roughness length in GFS
In GFS, the momentum roughness length z 0 has two terms. The first term, z ch , is parameterized by the Charnock relationship (Charnock, 1955), representing wave-resulted sea surface roughness, and the second term, z vis , is the viscous contribution (Beljaars, 1994) for low winds and smooth surface; that is, Here C ch = 0.014 is the constant Charnock parameter, and ν is the air kinematic viscosity. The relation of z 0 in GFS versus 10 m wind speed is shown in Fig. 2 (black line).

The Charnock relationship related to wave state
When ocean surface waves are explicitly considered, the Charnock parameter, C ch , is not a constant (Janssen, 1989(Janssen, , 1991Taylor and Yelland, 2001;Moon et al., 2004;Drennan et al., 2003Drennan et al., , 2005. In the study, we adopt a method developed by Moon et al. (2004), which considered the surface roughness leveling off under extremely high wind speed (Powell et al., 2003;Donelan et al., 2004). Based on observations, Moon et al. (2004) proposed Eq. (16) to estimate the Charnock parameter by the wave age c pi u * (c pi is the peak phase speed of the dominant wind-forced waves) with constant values of a and b changing with 10 m wind speed every 5 m s −1 in the range of 10 to 50 m s −1 .
To obtain continuous values of a and b, we derive a new relationship (Eq. 17 to estimate a and b from 10 m wind speed U 10 by fitting the values in Table 1  Because the observations in Moon et al. (2004) were obtained under tropical cyclones, Eq. (17) is used for U 10 >15 m s −1 , whereas the original Charnock relationship of the WW3 ST4 scheme (Janssen, 1989(Janssen, , 1991 is used for U 10 ≤ 15 m s −1 . The revised parameterization is called ST4-M04. Figure S2 in the Supplement shows the C ch distribution  (17). In general small wind direction variations at low latitudes lead to large wave age and thus low C ch . The situation is opposite at mid-to-high latitudes. The relationships between z 0 and U 10 in GFS, WW3 ST4 scheme (Janssen, 1989(Janssen, , 1991, and ST4-M04 scheme are compared in Fig. 2. The z 0 in GFS increases relatively slowly with increasing wind speed (black). The value of z 0 from ST4 scheme (purple) increases rapidly with wind speed at high winds. In comparison, in ST4-M04 scheme (blue) the rapid increase in z 0 at high wind speed is obviously restrained, although the mean z 0 is slightly higher than that in GFS at wind speed >10 m s −1 due to larger C ch (>0.014 in Fig. S2). Furthermore, since the Charnock number is constant in GFS, the standard deviation (SD) of z 0 at a given wind speed is near zero. Since z 0 is determined only by wind-sea conditions in ST4 and ST4-M04 schemes, the SD at a given wind speed is mainly due to variations in wind fetch and development stage of sea state. The reduced SDs in ST4-M04 scheme, compared to ST4, imply less sensitivity of z 0 to fetch and sea state. Note that ST4-M04 is used in GFS, while z 0 in WW3 is still calculated by the ST4 source term to avoid affecting the balance of adjusted wind input and dissipation.

Set of experiments
A series of numerical experiments is conducted to evaluate the effects of the aforementioned wave-related processes on ocean and atmosphere in two 56 d periods, from 3 January to 28 February 2017 and from 3 August to 28 September 2018 for boreal winter and boreal summer, respectively.
The reference experiment (CTRL) is a one-way coupled experiment, in which CFSv2.0 provides 10 m wind and surface current to WW3, whereas no variable is transferred from WW3 to CFSv2.0. The results of CFSv2.0 in CTRL are consistent with the corresponding CFS Reanalysis data (Saha et al., 2010). For each period, four sensitivity experiments are carried out (Table 1). Based on CTRL, the first is the VR12-AL-SC-EN experiment, in which the Langmuir mixing parameterization is used with the Stokes-Coriolis force and entrainment in MOM4. The second is the Z0-M04 experiment, in which the constant C ch in GFS is replaced by C ch from the WW3 ST4-M04 scheme. The effect of fluxes in GFS generated by V (Eq. 14) is tested in the FLUX experiment. The last experiment is ALL, which includes all three parameterizations.

Data
Due to the availability of in situ and reanalysis data in the simulation periods, only sea surface temperature (SST), ocean subsurface temperature and salinity (T&S), 2 m air temperature (T02), 10 m wind speed (WSP10), and significant wave height (SWH) are used to evaluate the simulation results.

Experimental results
In this section, an evaluation of simulation results is presented. Comparisons are made between model results and observations or reanalysis data. The results in the first 3 d are excluded in the evaluation, since the wave influences bias is approximately 0.32 • C, and the average RMSE is about 1.09 • C from day 4 to day 56 in CTRL (Fig. 3a). The simulated SSTs are generally overestimated, and the large biases (>1.0 • C) are mainly distributed in the Southern Ocean. In Fig. 3b, the global-averaged RMSEs of CTRL (black) increase with time in the first month and then gradually level off. Compared with CTRL, the RMSEs are reduced continuously in VR12-AL-SC-EN and ALL (yellow and red) but not in Z0-M04 and FLUX (purple and blue). To understand the critical process responsible for the bias reduction in ALL, the SST differences are compared across all four experiments (Fig. 3c-f). Clearly, the difference in experiment VR12-AL-SC-EN is similar to that in ALL ( Fig. 3c and f). The spatial correlation coefficient between the SST differences with CTRL of the two experiments ( Fig. 3c and f) is 0.67, significant at 99 % confidence level, and the RMSEs of SST are not different significantly (red and yellow lines in Fig. 3b), indicating the Stokes drift-related parameterizations in VR12-AL-SC-EN mainly contribute to the SST positive bias reduction. This contrasts with Couvelard et al. (2020), where SST overestimations and MLD underestimations are reduced mainly due to the directly modified turbulence kinetic energy scheme. The global mean SST bias in ALL is 0.02 • C with RMSE of 1.03 • C, and in most areas the SST differences compared with CTRL are significant (P ≤ 0.05) (dotted areas in Fig. 3f). Large SST improvements mainly appear in the Southern Ocean, with a regional RMSE decrease from 1.27 to 1.04 • C south of 45 • S ( Fig. 3f and red line in Fig. 5a). The reduction of overestimated SSTs in CTRL (red in Fig. 3a) is because the Stokes drift-related param-eterizations in MOM4 inject turbulent kinetic energy into the ocean, which enhances vertical mixing and subsequently cools the surface waters (Belcher et al., 2012;Li et al., 2016). The modified roughness and relative velocity in Z0-M04 and FLUX also influence upper ocean mixing ( Fig. 3d and e) via changing momentum flux and lead to generally warmer SSTs (purple and blue lines in Figs. 3b and 5a). The effect from Stokes drift-related ocean mixing parameterizations dominates SST changes in ALL.
In boreal summer, the global mean SST bias in CTRL is overestimated approximately 0.29 • C, and the averaged RMSE from day 4 to day 56 is about 1.19 • C. The overestimated SSTs (>1.0 • C) mainly occur in the Northern Hemisphere (Fig. 4a). The global-averaged RMSEs are also generally lower in VR12-AL-SC-EN and ALL than in CTRL (Fig. 4b). The cooling effects in VR12-AL-SC-EN lead to a global mean bias of 0.06 • C, and the large SST improvements mainly occur north of 50 • N ( Fig. 4c and yellow line in Fig. 5b). The changes in SST in Z0-M04 and FLUX (Fig. 4d and e; purple and blue lines in Figs. 4b and 5b) are relatively small. The global mean bias in ALL is 0.04 • C with an RMSE of 1.14 • C (Fig. 4f).
As mentioned before, large improvements of overestimated SST mainly occur at mid-to-high latitudes in local summer. The time series of RMSEs and correlation coefficients of SST between model and observation in the region (0-360 • E, 45-78 • S in boreal winter and 0-360 • E, 50-78 • N in boreal summer) are shown in Fig. 5c-f. The RM-SEs in CTRL (blue in Fig. 5c and d) increase in the first few weeks and then gradually decrease afterward. Compared We also compared T02 from experiments with ERA5 (Fig. 6). Warm biases of T02 appear in both winter and summer in CTRL ( Fig. 6a and b). The changes in T02 in sensitivity experiments (Fig. 6c-j) are generally consistent with the changes in SST in the same experiments (Figs. 3 and 4). The correlation coefficients between the SST and the T02 changes for the ALL experiment in boreal winter and summer (Figs. 3f and 6f and 4f and 6j) are 0.61 and 0.53, respectively, significant at 99 % confidence level. In boreal winter, all wave-coupled experiments except FLUX reduce the T02 mean bias (Fig. 6c-f). VR12-AL-SC-EC has the largest T02 bias reduction compared with CTRL, from 0.55 to 0.17 • C (Fig. 6c). In boreal summer, both VR12-AL-SC-EC and ALL have the largest T02 bias reduction, from 0.29 to 0.08 • C ( Fig. 6g and j). Noticeably, the improvements in RMSEs are not large for all experiments, because the improvements mainly occur in areas with overestimated temperature.

Mixed layer depth (MLD)
To further evaluate the direct effect of the wave-related processes on the upper ocean, we compare the MLD of all experiments with that estimated from Argo profiles in summer. The simulated T&S are interpolated onto the positions of Argo profiles at the nearest time. The MLD is estimated as the depth where the change of potential density reaches the value corresponding to a 0.2 • C decrease in potential temperature with unchanged salinity from surface (de Boyer Montégut et al., 2004;Wang and Xu, 2018).
The time series of MLDs from numerical experiments and Argo south of 45 • S in boreal winter (north of 45 • N in boreal summer) are compared in Fig. 7a (b). The simulated MLDs are generally within the SD of Argo MLDs (shading in Fig. 7). In CTRL, the mean bias (CTRL minus Argo) with SD is −13.15 ± 7.82 m (−6.75 ± 5.29 m) in bo- real winter (summer). The correlation coefficient of MLDs in CTRL with Argo MLDs is 0.55 (0.68) with P ≤ 0.01, and the mean RMSE is 15.30 m (8.55 m) in boreal winter (summer). In ALL, the mean bias (ALL minus Argo) with SD is 7.70 ± 10.42 m (3.30 ± 7.78 m) in boreal winter (summer), and the correlation coefficient of MLDs enhances to 0.63 (0.78). The RMSE south of 45 • S decreases from 15.30 m in CTRL to 12.96 m in ALL. The RMSE north of 45 • N decreases from 6.71 in CTRL to 5.55 m in ALL in the first 6 weeks but the value increases in the last 2 weeks due to overestimation of MLDs. Compared with CTRL (orange in Fig. 7), VR12-AL-SC-EN (yellow) and ALL (dark blue) show significant improvements (P ≤ 0.01) on the underestimated MLDs time series, whereas the MLDs difference between CTRL and Z0-M04 (purple)/FLUX (blue) is nonsignificant. The comparisons of the simulated SWHs in CTRL with the ERA5 also show that the SWHs are overestimated in both winter and summer (Figs. 10a and 11a). In boreal winter, the global mean SWH bias in CTRL is approximately 0.20 m with overestimates (>0.30 m) in the Pacific, the North Atlantic and the Southern Ocean (Fig. 10a), and the average RMSE is about 1.29 m. In boreal summer, the global mean bias in CTRL is approximately 0.17 m with 1.22 m RMSE (Fig. 11a). Similar to WSP10s, the RMSEs of SWHs also increase in the first 2 weeks and then gradually level off (Figs. 10b and 11b). The RMSEs in Z0-M04 and ALL (purple and red) are smaller than in CTRL over most of the time, consistent with changes in WSP10s. The correlation coefficients between changes in WSP10s and changes in SWHs in ALL are 0.77 and 0.73 in boreal winter and summer, re- spectively (Figs. 8f and 10f and 9f and 11f), significant at 99 % confidence level, indicating that the SWHs changes are closely related to changes in wind speeds.
In VR12-AL-SC-EN, the reduction of SST warm biases affects air temperature and stabilizes the marine atmospheric boundary layer (Sweet et al., 1981;O'Neill et al., 2003), and subsequently reduces WSP10s and SWHs with decreased global bias in boreal winter (Figs. 8c and 10c). In Z0-M04, the overestimated WSP10s and SWHs are also reduced (Figs. 8d and 10d) due to the larger z 0 with the ST4-M04 scheme at wind speed >10 m s −1 (Fig. 2). The increase in z 0 enhances wind stress and momentum transferred into the ocean and therefore reduces surface winds (Pineau-Guillou et al., 2018;Sauvage et al., 2020) and consequently reduces SWHs. In FLUX (Figs. 8e and 10e), U surf and u s (0) decrease wind stress and momentum transfer when their directions are consistent with wind directions and vice versa (Hersbach and Bidlot, 2008;Renault et al., 2016). For instance, the angles between wind and current are relatively small (<90 • ) in the northeastern Pacific, reducing the wind stress and thus enhancing WSP10s (Fig. 8e). In contrast, the large angles (>90 • ) between the northwesterlies and the Kuroshio in the northwestern Pacific enhance wind stress, and decrease WSP10s (Fig. 8e). Consequently, improvements occur in areas with misalignment of winds and currents. With all combined effects, the biases of WSP10s and SWHs in ALL in most regions are decreased (Figs. 8f and 10f), with the reduced global RMSEs of 4.17 m s −1 and 1.18 m, respectively. In boreal summer, the improvements of WSP10s and SWHs are relatively small in terms of global averaged RMSEs, because of smaller positive biases in CTRL (Figs. 9a and 11a). In ALL, the global averaged bias of WSP10s (SWHs) is −0.01 m s −1 (0.03 m). The largest reduction primarily appears in the Southern Ocean (Figs. 9f and 11f) to improve the overestimated westerlies and SWHs in CTRL (Figs. 9a and 11a).
Previous studies indicated that ocean surface winds in ERA5 are underestimated in some regions (Belmonte Rivas and Stoffelen, 2019;Kalverla et al., 2020;Sharmar and Markina, 2020). To better demonstrate the effects of waves on WSP10s and SWHs, comparisons of WSP10s and SWHs with the NDBC buoy data are made (Table 2 and Fig. 12). The differences between sensitivity experiments and CTRL are all statistically significant at 95 % confidence level. Buoys are mainly located in the northeastern Pacific, the tropical Pacific, and the northwestern Atlantic oceans (Fig. S3), and buoy identifiers with total numbers, longitudes, and latitudes are listed in Table S3. The method from Hsu et al. (1994) is used to adjust wind speeds from buoy data to the reference height of 10 m.
Compared to the NDBC data, the WSP10s and the SWHs in CTRL are generally overestimated in both winter and summer with positive mean biases (Table 2 and Fig. 12). The reduction of mean biases appears in all experiments except FLUX in boreal winter. The wave-related processes are most effective in areas with positive biases, consistent with previous comparisons with ERA5. In boreal winter, the angles between winds and currents are small. The wind stresses are then reduced in FLUX, and the WSP10s are enhanced, so the positive bias is further enhanced. The improvements in ALL are generally the largest (Table 2), with the WSP10s RMSE of 1.04 m s −1 (1.15 m s −1 ) and the SWHs RMSE of 0.36 m (0.24 m) in boreal winter (summer). As shown in Fig. 12, with the increase in WSP10s and SWHs, the reduction of overestimation in ALL compared with CTRL is more prominent.

Summary and discussion
To investigate the individual role played by wave-related processes on the atmosphere and ocean interface in a coupled global atmosphere-ocean-wave modeling system on an intraseasonal scale, we implement version 5.16 of WW3 into CFSv2.0 over the domain 78 • S-78 • N, using the C-Coupler2. In this coupled system, the WW3 is forced by 10 m wind and surface current generated in CFSv2.0. The Stokes drift-related Langmuir mixing, the Stokes-Coriolis force and entrainment in ocean, air-sea fluxes modified by surface current and the Stokes drift, and momentum roughness length (z 0 ) are considered separately, and the results of sensitivity experiments are compared against in situ buoys, satellite measurements, and ERA5 reanalysis. The effects of waves on intraseasonal prediction are examined in two 56 d cases, one for boreal winter and the other one for boreal summer.  Table 2.
The following key results are found: 1. Overestimated SST, T02, and underestimated MLD in the mid-to-high latitudes in CFSv2.0 are significantly improved, particularly in local summer. Langmuir turbulence, Stokes-Coriolis force and entrainment in VR12-AL-SC-EN lead to enhanced vertical mixing. The enhanced vertical mixing changes the temperature structure in the upper ocean and further affects air temperature. In boreal winter, the regional RMSE of SST (T02) in the Southern Ocean decreases from 1.27 (1.93) in the CTRL experiment to 1.04 (1.67) • C in the ALL experiment. In boreal summer, the effect is weaker because of the relatively smaller ocean areas in the midto-high latitudes of the Northern Hemisphere.
2. In general, all wave-related processes reduce biases for WSP10s and SWHs, particularly in regions where WSP10s and SWHs were overestimated. The decreased SSTs in VR12-AL-SC-EN stabilize the marine atmospheric boundary layer and lead to weakened WSP10s and SWHs. The modified roughness in Z0-M04 generally enhances momentum transfer into the ocean, and so decreases WSP10s and SWHs. The relative windwave-current speed in FLUX also affects wind stress and further influences WSP10s and SWHs. Compared with NDBC buoy observations and ERA5, the ALL experiment shows significant improvements.
In addition to the variables mentioned before, the changes in simulated enthalpy fluxes are also compared, which mainly depend on the WSP10s changes. However, the wave-related effects on enthalpy fluxes are nonsignificant for the 2-month simulation, so the results are not shown. The wave-related parameterizations used in the study mainly improve model biases at mid-to-high latitudes, and SST biases in tropical oceans are only slightly improved (Figs. 3 and 4). Breivik et al. (2015) improved SST as well as subsurface temperature simulations in Nucleus for European Modelling of the Ocean (NEMO) with parameterizations including the wave-related Charnock parameter, modification of water-side stress with wind input and wave dissipation, wave dissipation-related turbulent kinetic energy flux, and the Stokes-Coriolis force. Based on a global NEMO-WW3 coupled framework, Couvelard et al. (2020) modified the Charnock parameter, the Stokes drift-related forces and the Langmuir cell with misalignment of winds and waves, the oceanic surface momentum flux, and the turbulence kinetic energy to reduce SST and MLD biases. In addition, sea spray can enhance air-sea heat fluxes in the tropics (Andreas et al., 2008(Andreas et al., , 2015. We will consider more processes in future studies. Different parameterizations for the same wave-related process also deserve discussion. For ocean surface roughness, the most classic parameterizations are those developed by Janssen (1989Janssen ( , 1991, Taylor and Yelland (2001), and Drennan et al. (2003). The method by Taylor and Yelland (2001) requires the peak wavelength for the total spectrum, whereas that by Drennan et al. (2003) only requires the peak of windsea waves. This difference leads to the fact that the former is more suitable for a mixed sea state, while the latter is more suitable for a young sea state (Drennan et al., 2005). And the effect of Janssen's parameterization (Janssen, 1989(Janssen, , 1991 is similar to that of Drennan et al. (2003), since it is also based on the wind-sea conditions (Shimura et al., 2017).
Our case studies indicate that significant biases in the coupled system remain, probably owing to inaccuracy of the coarse resolution, absence of a coupled wave-ice module, and deficiency of initial fields. In addition, every individual model component could be further improved via new parameter settings or updated parametrization schemes. All these aspects will need to be considered for improving our coupled system.
Author contributions. FX and RS designed the experiments, and RS carried them out. RS developed the code of coupling parameterizations and produced the figures. ZF contributed to the installation and operation of CFSv2.0. LL and HY contributed to the application of C-Coupler2. XL and YZ provided the original code of CFSv2.0. RS prepared the article with contributions from all co-authors. FX and HL contributed to review and editing.
Competing interests. The contact author has declared that neither they nor their co-authors have any competing interests.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Acknowledgements. The authors would like to extend thanks to all developers of the CFSv2.0 model (https://cfs.ncep.noaa.gov/cfsv2/ downloads.html, last access: 20 January 2022). The work was supported by the National Key Research and Development Program of China (no. 2016YFC1401408) and Tsinghua University Initiative Scientific Research Program (2019Z07L01001). We thank Wei Xue and Yao Liu for helping us to test the coupling frequency. We also thank three anonymous reviewers and the handling editor for their constructive comments.
Financial support. This research has been supported by the National Key Research and Development Program of China (grant no. 2016YFC1401408) and the Tsinghua University Initiative Scientific Research Program (grant no. 2019Z07L01001).
Review statement. This paper was edited by Sophie Valcke and reviewed by two anonymous referees.