Articles | Volume 15, issue 5
Geosci. Model Dev., 15, 2345–2363, 2022
Geosci. Model Dev., 15, 2345–2363, 2022

Model evaluation paper 18 Mar 2022

Model evaluation paper | 18 Mar 2022

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

The effects of ocean surface waves on global intraseasonal prediction: case studies with a coupled CFSv2.0–WW3 system
Ruizi Shi1, Fanghua Xu1, Li Liu1, Zheng Fan1, Hao Yu1, Hong Li1,3, Xiang Li2, and Yunfei Zhang2 Ruizi Shi et al.
  • 1Department of Earth System Science, Ministry of Education Key Laboratory for Earth System Modeling, Institute for Global Change Studies, Tsinghua University, Beijing, 100084, China
  • 2Key Laboratory of Marine Hazards Forecasting, National Marine Environmental Forecasting Center, Ministry of Natural Resources, Beijing, 100081, China
  • 3Department of Atmospheric and Oceanic Sciences, Institute of Atmospheric Sciences, Fudan University, Shanghai, 200433, China

Correspondence: Fanghua Xu (


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.

1 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, 1991; Taylor and Yelland, 2001; Moon et al., 2004; Drennan et al., 2003, 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 (Cch) 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 Cch, considering the wave-induced kinematic stress (Janssen, 1989, 1991), the wave age (Drennan et al., 2003, 2005; Moon 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), Van Roekel et al. (2012), 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 individual 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.

Figure 1A 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, Cch, Lat, us(0), Vs, U10, and Usurf 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.


2 Methods and experiments

2.1 Coupling WAVEWATCH III with CFSv2.0

Version 5.16 of WW3 (WAVEWATCH III Development Group, 2016) developed by the National Oceanic and Atmospheric Administration (NOAA)/NCEP has been incorporated into the CFSv2.0 (Saha et al., 2014) as a new model component. The latitude range of WW3 is 78 S–78 N with a spatial resolution of 1/3; the frequency range is 0.04118–0.4056 Hz, and the total number of frequencies is 25; the number of wave directions is 24 with a resolution of 15; the maximum global time step and the minimum source term time step are both 180 s.

The CFSv2.0 contains two components: the GFS (details are available at, last access: 20 January 2022) as the atmosphere component and the Modular Ocean Model version 4 (MOM4; 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–ocean–wave system is shown in Fig. 1. As illustrated, WW3 is two-way 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 (, 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.

2.2 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 input–dissipation 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) data and the fifth-generation European Centre for Medium-Range Weather Forecasts (ECMWF) Reanalysis (ERA5; Hersbach et al., 2020) data, are used to drive the wave model. Comparing all results, the ST4 scheme with ERA5 wind forcing generates the minimum SWH bias (Table S2 in the Supplement), 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 WAVEWATCH III Development Group (2016), which is the CFSR (CFS Reanalysis) tuned setup and is commonly used at the global scale.

Table 1List of numerical experiments: setups different from control (CTRL) are marked with bold.

Download Print Version | Download XLSX

Figure 2Relationships between momentum roughness length z0 (m) in the coupled system and 10 m wind speed (m s−1); error bars indicate twice the standard deviations for each point.


2.3 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. (2014, 2016). Breivik et al. (2016) derived the full Stokes drift profile as

(1) u s B 16 ( z ) = u s ( 0 ) [ exp 2 k p z - - 2 π k p z erfc ( - 2 k p z ) ] ,

where us(0) is the surface Stokes drift velocity, kp=us06Vs, Vs 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), Wu et al. (2019), and Couvelard et al. (2020).

2.3.1 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 Lat is the turbulent Langmuir number, defined as

(6) L a t = u u s ( 0 ) ,

with us(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 Ricr=0.3; that is,

(7) R i b h = g h [ ρ r - ρ ( h ) ] ρ 0 [ u r - u ( h ) 2 + W 2 ] = R i cr ,

where g is acceleration due to gravity, ρ is density, u is velocity, ρr is surface density, ur 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.

Figure 3The 53 d average SST (C) bias in CTRL (a; CTRL minus OISST), the time series of global-averaged RMSE (b), and the differences between VR12-AL-SC-EN (c) /Z0-M04 (d) /FLUX (e) /ALL (f), and CTRL in January–February 2017 (VR12-AL-SC-EN/Z0-M04/FLUX/ALL minus CTRL). The first 3 d of the simulation is discarded. The dotted areas are statistically significant at 95 % confidence level.

2.3.2 The Stokes–Coriolis force and associated entrainment

Because the Stokes drift velocity is an increment superimposed on the original current velocity, the Coriolis force and the Stoke drift together produce an additional so-called Stokes–Coriolis (SC) force (Hasselmann 1970); that is,

(8) SC Force = u s ( z ) × f z .

Here us(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

(9) η t = - M curr - M st ,

where η is surface elevation, and Mcurr and Mst are the total vertical integral of the regular Eulerian current and the Stokes drift, respectively.

Figure 4The same as Fig. 3 but for August–September 2018.

Figure 5The 53 d averaged latitudinal distribution of SST root-mean-square errors (RMSE), time series of domain-averaged SST RMSE, and correlation coefficient: (a), (b) the latitudinal RMSE in boreal winter/summer compared with OISST; (c), (d) the time series of domain-averaged (0–360 E, 45–78 S/50–78 N) SST RMSE in boreal winter/summer; (e), (f) the time series of domain-averaged (0–360 E, 45–78 S/50–78 N) SST correlation coefficient in boreal winter/summer; differences of RMSE and correlation coefficient time series between VR12-AL-SC-EN/ALL and CTRL are statistically significant at 99 % confidence level, except those in Fig. e​​​​​​​.


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 (|us(0)|2) to the denominator of Eq. (7); that is,

(10) R i b h = g h [ ρ r - ρ ( h ) ] ρ 0 [ u r - u ( h ) 2 + W 2 + | u s ( 0 ) | 2 ] = R i cr .

The boundary layer depth h in KPP from Eq. (10) is then enhanced due to the Stokes drift velocity.

2.4 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 Cd, Ch, and Ce 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 (Uwind). However, the effect of ocean surface current should not be ignored. Luo et al. (2005) first indicated that including ocean surface current (Usurf) improves estimates of τ and subsequent ocean response. Renault et al. (2016) further indicated that the improvements of τ by Usurf also fed back into the atmosphere. At present, ΔV=Uwind-Usurf 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 (us(0)) should also be included; that is,

(14) Δ V = U wind - U surf - u s ( 0 ) .

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 Usurf and us(0) on ΔV depend on the angles between them and Uwind.

2.5 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).

2.5.1 The momentum roughness length in GFS

In GFS, the momentum roughness length z0 has two terms. The first term, zch, is parameterized by the Charnock relationship (Charnock, 1955), representing wave-resulted sea surface roughness, and the second term, zvis, is the viscous contribution (Beljaars, 1994) for low winds and smooth surface; that is,

(15) z 0 = z ch + z vis = C ch u 2 g + 0.11 ν u .

Here Cch=0.014 is the constant Charnock parameter, and ν is the air kinematic viscosity. The relation of z0 in GFS versus 10 m wind speed is shown in Fig. 2 (black line).

Figure 6The 53 d average T02 (C) bias in CTRL (a and b; CTRL minus ERA5), and the differences between VR12-AL-SC-EN (c and g) /Z0-M04 (d and h) /FLUX (e and i) /ALL (f and j) and CTRL (VR12-AL-SC-EN/Z0-M04/FLUX/ALL minus CTRL). The first 3 d of the simulation is discarded. The dotted areas are statistically significant at 95 % confidence level. Panels (a), (c), (d), (e), and (f) are for January–February 2017, and panels (b), (g), (h), (i), and (j) are for August–September 2018.

Figure 7The 53 d time series of domain-averaged (0–360 E, 45–78 S/N) mixed layer depth (MLD, m) in boreal winter/summer: the difference between CTRL and VR12-AL-SC-EN/ALL passes the Student's t test at 99 % confidence level; the time intervals are 6 h; shaded areas indicate twice the standard deviations for Argo.


2.5.2 The Charnock relationship related to wave state

When ocean surface waves are explicitly considered, the Charnock parameter, Cch, is not a constant (Janssen, 1989, 1991; Taylor and Yelland, 2001; Moon et al., 2004; Drennan et al., 2003, 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 cpiu (cpi 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.

(16) C ch = a ( c pi u ) b .

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 U10 by fitting the values in Table 1 of Moon et al., 2004):


Because the observations in Moon et al. (2004) were obtained under tropical cyclones, Eq. (17) is used for U10>15 m s−1, whereas the original Charnock relationship of the WW3 ST4 scheme (Janssen, 1989, 1991) is used for U10≤15 m s−1. The revised parameterization is called ST4-M04. Figure S2 in the Supplement shows the Cch distribution obtained by Eqs. (16)–(17). In general small wind direction variations at low latitudes lead to large wave age and thus low Cch. The situation is opposite at mid-to-high latitudes.

The relationships between z0 and U10 in GFS, WW3 ST4 scheme (Janssen, 1989, 1991), and ST4-M04 scheme are compared in Fig. 2. The z0 in GFS increases relatively slowly with increasing wind speed (black). The value of z0 from ST4 scheme (purple) increases rapidly with wind speed at high winds. In comparison, in ST4-M04 scheme (blue) the rapid increase in z0 at high wind speed is obviously restrained, although the mean z0 is slightly higher than that in GFS at wind speed >10 m s−1 due to larger Cch (>0.014 in Fig. S2). Furthermore, since the Charnock number is constant in GFS, the standard deviation (SD) of z0 at a given wind speed is near zero. Since z0 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 z0 to fetch and sea state. Note that ST4-M04 is used in GFS, while z0 in WW3 is still calculated by the ST4 source term to avoid affecting the balance of adjusted wind input and dissipation.

Figure 8The 53 d average WSP10 (m s−1) bias in CTRL (a; CTRL minus ERA5), the time series of global-averaged RMSE (b), and the differences between VR12-AL-SC-EN (c) /Z0-M04 (d) /FLUX (e) /ALL (f) and CTRL in January–February 2017 (VR12-AL-SC-EN/Z0-M04/FLUX/ALL minus CTRL). The first 3 d of the simulation is discarded. The dotted areas are statistically significant at 95 % confidence level.

Figure 9The same as Fig. 8 but for August–September 2018.

2.6 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 Cch in GFS is replaced by Cch 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.

3 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.

The daily-average satellite Optimum Interpolation SST (OISST) data are obtained from NOAA, with 0.25× 0.25 resolution (Reynolds et al., 2007;, last access: 20 January 2022). The global Argo observational profiles of T&S (Li et al., 2019) are from China Argo Real-time Data Center (, last access: 20 January 2022). The ERA5 datasets of T02, WSP10, and SWH with a spatial resolution of 0.5 are also used (Hersbach et al., 2020;!/dataset/reanalysis-era5-single-levels, last access: 20 January 2022), which assimilated huge amounts of historical data and thus provided reliable hourly estimates. Additionally, the WSP10 and SWH observations from the available National Data Buoy Center (NDBC) buoy data (, last access: 20 January 2022) are used for comparison.

4 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 are weak at the beginning. Compared with observations or ERA5, the general increase in the biases in all experiments is likely a drift from the initial conditions since no data are assimilated.

4.1 Sea surface temperature (SST) and 2 m air temperature (T02)

Figure 3a shows the spatial distribution of 53 d (day 4 to day 56) averaged SST biases in CTRL in boreal winter, defined as SST in CTRL minus OISST. The global mean SST 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.0C) 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).

Figure 10The 53 d average SWH (m) bias in CTRL (a; CTRL minus ERA5), the time series of global-averaged RMSE (b), and the differences between VR12-AL-SC-EN (c) /Z0-M04 (d) /FLUX (e) /ALL (f) and CTRL in January–February 2017 (VR12-AL-SC-EN/Z0-M04/FLUX/ALL minus CTRL). The first 3 d of the simulation is discarded. The dotted areas are statistically significant at 95 % confidence level.

Figure 11The same as Fig. 10 but for August–September 2018.

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 parameterizations 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.0C) 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 RMSEs in CTRL (blue in Fig. 5c and d) increase in the first few weeks and then gradually decrease afterward. Compared with CTRL, RMSEs in VR12-AL-SC-EN (yellow) and ALL (red) are significantly (P≤0.01) reduced by about 0.3 C. The spatial correlation coefficients decrease with time but remain high (>0.90) for all experiments (Fig. 5e and f) with higher values in experiment VR12-AL-SC-EN (yellow).

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.

4.2 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 boreal 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.

4.3 Wind speed at 10 m (WSP10) and significant wave height (SWH)

Compared with ERA5, the WSP10s in CTRL are generally overestimated in both winter and summer (Figs. 8a and 9a). The global averaged RMSEs of WSP10s in CTRL are 4.25 m s−1 (4.26 m s−1) in boreal winter (summer). The global averaged RMSEs of WSP10s in all experiments increase with time in the first 2 weeks and then gradually level off (Figs. 8b and 9b). The differences of RMSEs between CTRL and other experiments are tiny in the first 10 d, and afterwards the RMSEs in Z0-M04 and ALL (purple and red) become clearly smaller than in CTRL over most of the time.

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, respectively (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 z0 with the ST4-M04 scheme at wind speed >10 m s−1 (Fig. 2). The increase in z0 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), Usurf and us(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.

Table 2The 53 d mean bias with standard deviation (SD) and RMSE for WSP10 and SWH compared with NDBC buoy observation: the bias is calculated as simulation minus NDBC.

Download Print Version | Download XLSX

Figure 12Scatterplots of simulated WSP10/SWH (y axis) vs. buoy WSP10/SWH (x axis): (a) the WSP10 in January–February 2017, (b) the SWH in January–February 2017, (c) the WSP10 in August–September 2018, and (d) the SWH in August–September 2018. The dotted line is y=x. The corresponding mean biases with standard deviations and RMSEs for every experiment are shown in Table 2.


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.

5 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 (z0) 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.

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 mid-to-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 wind–wave-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, 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 (1989, 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 wind-sea 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, 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.

Code and data availability

The code developed for the coupled system can be found under (Shi et al., 2021), including the coupling, preprocessing, run control, and postprocessing scripts. The initial fields for CFSv2.0 are generated by the real-time operational Climate Data Assimilation System, downloaded from the CFSv2.0 official website (, NOAA, 2022b). The daily-average satellite Optimum Interpolation SST (OISST) data are obtained from NOAA (, National Centers for Environmental Information, 2022), and the National Data Buoy Center (NDBC) buoy data are also obtained from NOAA (, NOAA, 2022a). The Argo observational profiles of T&S are available at China Argo Real-time Data Center (, China Argo Real-time Data Center, 2022). The ERA5 reanalysis data are available at the Copernicus Climate Change Service (C3S) Climate Date Store (, Hersbach et al., 2018​​​​​​​).


The supplement related to this article is available online at:

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.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


The authors would like to extend thanks to all developers of the CFSv2.0 model (, 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.


Andreas, E. L., Persson, P. O. G., and Hare, J. E.: A bulk turbulent air-sea flux algorithm for high-wind, spray conditions, J. Phys. Oceanogr., 38, 1581–1596, 2008. 

Andreas, E. L., Mahrt, L., and Vickers, D.: An improved bulk air-sea surface flux algorithm, including spray-mediated transfer, Q. J. Roy. Meteor. Soc., 141, 642–654, 2015. 

Ardhuin, F., Rogers, E., Babanin, A. V., Filipot, J., Magne, R., Roland, A., Der Westhuysen, A. V., Queffeulou, P., Lefevre, J. M., and Aouf, L.: Semiempirical Dissipation Source Functions for Ocean Waves. Part I: Definition, Calibration, and Validation, J. Phys. Oceanogr., 40, 1917–1941,, 2010. 

Atlas, R., Hoffman, R. N., Ardizzone, J., Leidner, S. M., Jusem, J. C., Smith, D. K., and Gombos, D.: A Cross-calibrated, Multiplatform Ocean Surface Wind Velocity Product for Meteorological and Oceanographic Applications, B. Am. Meteorol. Soc., 92, 157–174,, 2011. 

Bao, Y., Song, Z., and Qiao, F.: FIO-ESM version 2.0: Model description and evaluation, J. Geophys. Res.-Oceans, 125, e2019JC016036,, 2019. 

Belcher, S. E., Grant, A. L. M., Hanley, K., Foxkemper, B., Van Roekel, L., Sullivan, P. P., Large, W. G., Brown, A. R., Hines, A., and Calvert, D.: A global perspective on Langmuir turbulence in the ocean surface boundary layer, Geophys. Res. Lett., 39, L18605,, 2012. 

Beljaars, A. C. M.: The parametrization of surface fluxes in large-scale models under free convection. Q. J. Roy. Meteor. Soc., 121, 255–270,, 1994. 

Belmonte Rivas, M. and Stoffelen, A.: Characterizing ERA-Interim and ERA5 surface wind biases using ASCAT, Ocean Sci., 15, 831–852,, 2019. 

Bidlot, J.-R.: Present status of wave forecasting at ECMWF, Workshop on ocean waves, online, 25–27 June 2012, 14784, (last access: 20 January 2022), 2012. 

Bidlot, J.-R.: Model upgrade improves ocean wave forecasts, ECMWF newsletter, online, April 2019, 159, 10, (last access: 20 January 2022), 2019. 

Bidlot, J.-R., Prates, F., Ribas, R., Mueller-Quintino, A., Crepulja, M., and Vitart, F.: Enhancing tropical cyclone wind forecasts, ECMWF newsletter, online, July 2020, 164, 33–37, (last access: 20 January 2022), 2020. 

Breivik, Ø., Janssen, P. A., and Bidlot, J.-R.: Approximate Stokes drift profiles in deep water, J. Phys. Oceanogr., 44, 2433–2445, 2014. 

Breivik, Ø., Mogensen, K., Bidlot, J., Balmaseda, M., and Janssen, P. A. E. M.: Surface wave effects in the NEMO ocean model: Forced and coupled experiments, J. Geophys. Res., 120, 2973–2992,, 2015. 

Breivik, Ø., Bidlot, J.-R., and Janssen, P. A.: A Stokes drift approximation based on the Phillips spectrum, Ocean Model., 100, 49–56, 2016. 

Charnock, H.: Wind stress on a water surface, Q. J. Roy. Meteor. Soc., 81, 639–640,, 1955. 

China Argo Real-time Data Center: Argo Observational Profiles of T&S, China Argo Real-time Data Center [data set],, 2022. 

Couvelard, X., Lemarié, F., Samson, G., Redelsperger, J.-L., Ardhuin, F., Benshila, R., and Madec, G.: Development of a two-way-coupled ocean–wave model: assessment on a global NEMO(v3.6)–WW3(v6.02) coupled configuration, Geosci. Model Dev., 13, 3067–3090,, 2020. 

de Boyer Montégut, C., Madec, G., Fischer, A. S., Lazar, A., and Iudicone, D.: Mixed layer depth over the global ocean: An examination of profile data and a profile-based climatology, J. Geophys. Res., 109, C12003,, 2004. 

Donelan, M. A., Haus, B. K., Reul, N., Plant, W. J., Stiassnie, M., Graber, H. C., Brown, O. B., and Saltzman, E. S.: On the limiting aerodynamic roughness of the ocean in very strong winds, Geophys. Res. Lett., 31, L18306,, 2004. 

Drennan, W. M., Graber, H. C., Hauser, D., and Quentin, C.: On the wave age dependence of wind stress over pure wind seas, J. Geophys. Res., 108, 8062,, 2003. 

Drennan, W. M., Taylor, P. K., and Yelland, M. J.: Parameterizing the sea surface roughness, J. Phys. Oceanogr., 35, 835–848, 2005. 

ECMWF: Official IFS Documentation CY47R1, In chap. PART VII: ECMWF wave model​​​​​​​, ECMWF, Reading, UK,, 2020. 

Fan, Y., Lin, S., Held, I. M., Yu, Z., and Tolman, H. L.: Global Ocean Surface Wave Simulation Using a Coupled Atmosphere-Wave Model, J. Climate, 25, 6233–6252,, 2012. 

Ghantous, M. and Babanin, A. V.: One-dimensional modelling of upper ocean mixing by turbulence due to wave orbital motion, Nonlin. Processes Geophys., 21, 325–338,, 2014. 

Griffies, S. M., Harrison, M. J., Pacanowski, R. C., and Rosati, A.: A technical guide to MOM4, GFDL Ocean Group Tech. Rep, 5, 1–342, (last access: 20 January 2022), 2004. 

Hasselmann, K.: Wave-driven inertial oscillations, Geophys. Astro. Fluid, 1, 463–502,, 1970. 

Hersbach H. and Bidlot, J.-R.: The relevance of ocean surface current in the ECMWF analysis and forecast system. Proceeding from the ECMWF Workshop on Atmosphere-Ocean Interaction, Shinfield Park, Reading, 10–12 November 2008, 61–73, (last access: 20 January 2022), 2008. 

Hersbach, H., Bell, B., Berrisford, P., Biavati, G., Horányi, A., Muñoz Sabater, J., Nicolas, J., Peubey, C., Radu, R., Rozum, I., Schepers, D., Simmons, A., Soci, C., Dee, D., Thépaut, J.-N.: ERA5 hourly data on single levels from 1979 to present, Copernicus Climate Change Service (C3S) Climate Data Store (CDS) [data set],, 2018. 

Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz-Sabater, J., Nicolas, J., Peubey, C., Radu, R., Schepers, D., Simmons, A., Soci, C., Abdalla, S., Abellan, X., Balsamo, G., Bechtold, P., Biavati, G., Bidlot, J., Bonavita, M., De Chiara, G., Dahlgren, P., Dee, D., Diamantakis, M., Dragani, R., Flemming, J., Forbes, R., Fuentes, M., Geer, A., Haimberger, L., Healy, S., Hogan, R. J., Hólm, E., Janisková, M., Keeley, S., Laloyaux, P., Lopez, P., Lupu, C., Radnoti, G., de Rosnay, P., Rozum, I., Vamborg, F., Villaume, S., and Thépaut, J.-N.: The ERA5 global reanalysis, Q. J. Roy. Meteor. Soc., 146, 1999–2049,, 2020. 

Hsu, S. A., Meindl, E. A., and Gilhousen, D. B.: Determining the Power-Law Wind-Profile Exponent under Near-Neutral Stability Conditions at Sea, Applied Meteorology, 33, 757–765,<0757:DTPLWP>2.0.CO;2, 1994. 

Janssen, P. A. E. M.: Wave-induced stress and the drag of air flow over sea waves, J. Phys. Oceanogr., 19, 745–754,<0745:WISATD>2.0.CO;2​​​​​​​, 1989. 

Janssen, P. A. E. M.: The Quasi-linear theory of wind wave generation applied to wave forecasting, J. Phys. Oceanogr., 21, 1631–1642,<1631:QLTOWW>2.0.CO;2, 1991. 

Janssen, P. A. E. M. (Ed.): The interaction of ocean waves and wind, Cambridge University Press, 1st edn., 312 pp., ISBN 0521465400, 2004. 

Kalnay, E., Kanamitsu, M., Kistler, R., Collins, W. D., Deaven, D. G., Gandin, L. S., Iredell, M. D., Saha, S., White, G. H., and Woollen, J.: The NCEP/NCAR 40-Year Reanalysis Project, B. Am. Meteorol. Soc., 77, 437–471,<0437:TNYRP>2.0.CO;2, 1996. 

Kalverla, P. C., Holtslag, A. A. M., Ronda, R. J., and Steeneveld, G.-J.: Quality of wind characteristics in recent wind atlases over the North Sea, Q. J. Roy. Meteor. Soc., 146, 1498–1515,, 2020. 

Kukulka, T., Plueddemann, A. J., Trowbridge, J. H., and Sullivan, P. P.: Significance of Langmuir circulation in upper ocean mixing: comparison of observations and simulations, Geophys. Res. Lett., 36, L10603,, 2009. 

Law-Chune, S. and Aouf, L.: Wave effects in global ocean modeling: parametrizations vs. forcing from a wave model, Ocean Dynam., 68, 1739–1758,, 2018. 

Lemarié, F.: Numerical modification of atmospheric models to include the feedback of oceanic currents on air-sea fluxes in ocean-atmosphere coupled models, INRIA Grenoble-Rhône-Alpes, Laboratoire Jean Kuntzmann, (last access: 20 January 2022​​​​​​​), 2015. 

Li, D., Staneva, J., Bidlot, J.-R., Grayek, S., Zhu, Y., and Yin, B.: Improving Regional Model Skills During Typhoon Events: A Case Study for Super Typhoon Lingling Over the Northwest Pacific Ocean, Frontiers in Marine Science, 8, 613913,, 2021. 

Li, M. and Garrett, C.: Mixed Layer Deepening Due to Langmuir Circulation, J. Phys. Oceanogr., 27, 121–132,<0121:MLDDTL>2.0.CO;2, 1997. 

Li, Q., Webb, A., Foxkemper, B., Craig, A., Danabasoglu, G., Large, W. G., and Vertenstein, M.: Langmuir mixing effects on global climate: WAVEWATCH III in CESM, Ocean Model., 103, 145–160,, 2016. 

Li, Q., Foxkemper, B., Breivik, O., and Webb, A.: Statistical models of global Langmuir mixing, Ocean Model., 113, 95–114,, 2017. 

Li, Z., Liu, Z., and Xing, X.: User Manual for Global Argo Observational data set (V3.0) (1997–2019), China Argo Real-time Data Center, Hangzhou, 37 pp., (last access: 20 January 2022), 2019. 

Liu, L., Zhang, C., Li, R., Wang, B., and Yang, G.: C-Coupler2: a flexible and user-friendly community coupler for model coupling and nesting, Geosci. Model Dev., 11, 3557–3586,, 2018. 

Luo, J. J., Masson, S., Roeckner, E., Madec, G., and Yamagata, T.: Reducing climatology Bias in an ocean-atmosphere CGCM with improved coupling physics, J. Climate, 18, 2344–2360,, 2005. 

McWilliams, J. C. and Sullivan, P. P.: Vertical Mixing by Langmuir Circulations, Spill Sci. Technol. B., 6, 225–237,, 2000. 

Monin, A. S. and Obukhov, A. M.: Basic laws of turbulent mixing in the surface layer of the atmosphere, Contrib. Geophys. Inst. Acad. Sci. USSR, 151, e187, (last access: 20 January 2022), 1954. 

Moon, I., Ginis, I., and Hara, T.: Effect of surface waves on Charnock coefficient under tropical cyclones, Geophys. Res. Lett., 31, L20302,, 2004. 

Moum, J. N. and Smyth, W. D.: Upper Ocean Mixing, in: Encyclopedia of Ocean Sciences, 3rd edn., vol. 1, edited by: Cochran, J. K., Bokuniewicz, J. H., and Yager, L. P., Elsevier, 71–79, ISBN 978-0-12-813081-0, 2019. 

National Centers for Environmental Information: Optimum Interpolation SST, National Centers for Environmental Information [data set],, 2022. 

NOAA: National Data Buoy Center, NOAA [data set],, 2022a. 

NOAA: CFSv2.0 Initial Fields, NOAA [data set],, 2022b. 

O'Neill, L. W., Chelton, D. B., and Esbensen, S. K.: Observations of SST-induced perturbations of the wind stress field over the Southern Ocean on seasonal timescales, J. Climate, 16, 2340–2354,, 2003. 

Pineau-Guillou, L., Ardhuin, F., Bouin, M. N., Redelsperger, J. L., Chapron, B., Bidlot, J. R., and Quilfen, Y.: Strong winds in a coupled wave-atmosphere model during a North Atlantic storm event: evaluation against observations, Q. J. Roy. Meteor. Soc., 144, 317–332, 2018. 

Polonichko, V.: Generation of Langmuir circulation for nonaligned wind stress and the Stokes drift, J. Geophys. Res., 102, 15773–15780,, 1997. 

Powell, M. D., Vickery, P. J., and Reinhold, T. A.: Reduced drag coefficient for high wind speeds in tropical cyclones, Nature, 422, 279–283,, 2003. 

Qiao, F., Yuan, Y., Yang, Y., Zheng, Q., Xia, C., and Ma, J.: Wave-induced mixing in the upper ocean: Distribution and application to a global ocean circulation model, Geophys. Res. Lett., 31, L11303,, 2004. 

Renault, L., Molemaker, M. J., McWilliams, J. C., Shchepetkin, A. F., Lemarie, F., Chelton, D., Illig, S., and Hall, A.: Modulation of Wind Work by Oceanic Current Interaction with the Atmosphere, J. Phys. Oceanogr., 46, 1685–1704,, 2016. 

Renault, L., Arsouze, T., and Ballabrera-Poy, J.: On the Influence of the Current Feedback to the Atmosphere on the Western Mediterranean Sea Dynamics, J. Geophys. Res.-Oceans, 126, e2020JC016664,, 2021. 

Reynolds, R. W., Smith, T. M., Liu, C., Chelton, D. B., Casey, K. S., and Schlax, M. G.: Daily High-Resolution-Blended Analyses for Sea Surface Temperature, J. Climate, 20, 5473–5496,, 2007. 

Saha, S., Moorthi, S., Pan, H., Wu, X., Wang, J., Nadiga, S., Tripp, P., Kistler, R., Woollen, J., and Behringer, D.: The NCEP Climate Forecast System Reanalysis, B. Am. Meteorol. Soc., 91, 1015–1057,, 2010. 

Saha, S., Moorthi, S., Wu, X., Wang, J., Nadiga, S., Tripp, P., Behringer, D., Hou, Y., Chuang, H., and Iredell, M. D.: The NCEP Climate Forecast System Version 2, J. Climate, 27, 2185–2208,, 2014. 

Sauvage, C., Lebeaupin Brossier, C., Bouin, M.-N., and Ducrocq, V.: Characterization of the air–sea exchange mechanisms during a Mediterranean heavy precipitation event using realistic sea state modelling, Atmos. Chem. Phys., 20, 1675–1699,, 2020. 

Sharmar, V. and Markina, M.: Validation of global wind wave hindcasts using ERA5, MERRA2, ERA-Interim and CFSRv2 reanalyzes, IOP C. Ser. Earth Env., 606, 012056​​​​​​​,, 2020. 

Shi, R., Xu, F., Liu, L., Fan, Z., and Yu, H.: The Effects of Ocean Surface Waves on Global Forecast in CFS Modeling System v2.0, Zenodo [code],, 2021. 

Shimura, T., Mori, N., Takemi, T., and Mizuta, R.: Long term impacts of ocean wave-dependent roughness on global climate systems, J. Geophys. Res.-Oceans, 122, 1995–2011,, 2017. 

Stopa, J. E., Ardhuin, F., Babanin, A. V., and Zieger, S.: Comparison and validation of physical wave parameterizations in spectral wave models, Ocean Model., 103, 2–17,, 2016.  

Sweet, W., Fett, R., Kerling, J., and Laviolette, P.: Air-sea interaction effects in the lower troposphere across the north wall of the Gulf Stream, Mon. Weather Rev., 109, 1042–1052,<1042:Asieit>2.0.Co;2, 1981. 

Takatama, K. and Schneider, N.: The Role of Back Pressure in the Atmospheric Response to Surface Stress Induced by the Kuroshio, J. Atmos. Sci., 74, 597–615,, 2017. 

Taylor, P. K. and Yelland, M. J.: The Dependence of Sea Surface Roughness on the Height and Steepness of the Waves, J. Phys. Oceanogr., 31, 572–590,<0572:TDOSSR>2.0.CO;2, 2001. 

Terray, E. A., Donelan, M. A., Agrawal, Y. C., Drennan, W. M., Kahma, K. K., Williams, A. J., Hwang, P. A., and Kitaigorodskii, S. A.: Estimates of Kinetic Energy Dissipation under Breaking Waves, J. Phys. Oceanogr., 26, 792–807,<0792:EOKEDU>2.0.CO;2, 1996. 

Tolman, H. L. and Chalikov, D. V.: Source Terms in a Third-Generation Wind Wave Model, J. Phys. Oceanogr., 26, 2497–2518,<2497:STIATG>2.0.CO;2, 1996. 

Van Roekel, L. P., Foxkemper, B., Sullivan, P. P., Hamlington, P. E., and Haney, S.: The form and orientation of Langmuir cells for misaligned winds and waves, J. Geophys. Res., 117, C05001,, 2012. 

Wang, L. and Xu, F.: Decadal variability and trends of oceanic barrier layers in tropical Pacific, Ocean Dynam., 68, 1155–1168,, 2018. 

WAVEWATCH III Development Group: User manual and system documentation of WAVEWATCH III version 5.16, National Oceanic and Atmospheric Administration, Camp Springs, US, NOAA/NWS/NCEP/MMAB Technical Note 329, 1–361, (last access: 20 January 2022), 2016. 

Wu, L., Staneva, J., Breivik, Ø., Rutgersson, A., Nurser, A. G., Clementi, E., and Madec, G.: Wave effects on coastal upwelling and water level, Ocean Model., 140, 101405,, 2019. 

Wu, X., K. S. Moorthi, K. Okomoto, and H. L. Pan: Sea ice impacts on GFS forecasts at high latitudes, Proceedings of the 85th AMS Annual Meeting, 8th Conference on Polar Meteorology and Oceanography, San Diego, CA, 12 January 2005, 7.4, (last access: 20 January 2022), 2005. 

Zieger, S., Babanin, A. V., Rogers, W. E., and Young, I. R.: Observation-based source terms in the third-generation wave model WAVEWATCH, Ocean Model., 96, 2–25,, 2015. 

Short summary
To better understand the effects of surface waves on global intraseasonal prediction, we incorporated the WW3 model into CFSv2.0. Processes of Langmuir mixing, Stokes–Coriolis force with entrainment, air–sea fluxes modified by Stokes drift, and momentum roughness length were considered. Results from two groups of 56 d experiments show that overestimated sea surface temperature, 2 m air temperature, 10 m wind, wave height, and underestimated mixed layer from the original CFSv2.0 are improved.