Recipes for How to Force Oceanic Model Dynamics

The current feedback to the atmosphere (CFB) contributes to the oceanic circulation by damping eddies. In an ocean‐atmosphere coupled model, CFB can be correctly accounted for by using the wind relative to the oceanic current. However, its implementation in a forced oceanic model is less straightforward as CFB also enhances the 10‐m wind. Wind products based on observations have seen real currents that will not necessarily correspond to model currents, whereas meteorological reanalyses often neglect surface currents or use surface currents that, again, will differ from the surface currents of the forced oceanic simulation. In this study, we use a set of quasi‐global oceanic simulations, coupled or not with the atmosphere, to (i) quantify the error associated with the different existing strategies of forcing an oceanic model, (ii) test different parameterizations of the CFB, and (iii) propose the best strategy to account for CFB in forced oceanic simulation. We show that scatterometer wind or stress are not suitable to properly represent the CFB in forced oceanic simulation. We furthermore demonstrate that a parameterization of CFB based on a wind‐predicted coupling coefficient between the surface current and the stress allows us to reproduce the main characteristics of a coupled simulation. Such a parameterization can be used with any forcing set, including future coupled reanalyses, assuming that the associated oceanic surface currents are known. A further assessment of the thermal feedback of the surface wind in response to oceanic surface temperature gradients shows a weak forcing effect on oceanic currents.

The mesoscale TFB generates wind and surface stress magnitude, divergence, and curl anomalies (Chelton et al., 2004(Chelton et al., , 2007O'Neill et al., 2010O'Neill et al., , 2012 in response to Sea Surface Temperature (SST) gradients. Small et al. (2008) provides a review of the different processes involved. The mesoscale TFB-induced stress curl anomalies induce Ekman pumping that can have an influence on eddy propagation but not on eddy magnitude (Seo et al., 2016;Seo, 2017). Also, Ma et al. (2016) suggest that the mesoscale TFB, by causing wind velocity and turbulent heat flux anomalies, could regulate the Western Boundary Currents; however, in that study, a large spatial filter (with a spatial cutoff of more than 1000 km) is used.
When estimating the surface stress in a coupled model, the CFB is taken into account by using the wind relative to the oceanic current instead of the absolute wind neglecting surface motion . From an oceanic perspective, coupled simulations are computationally expensive because of the high computational cost of atmospheric models with respect to oceanic ones. Therefore, oceanic simulations are usually forced by an atmospheric product (derived from observations or numerical simulations). However, the recent finding on the role of the CFB in determining the ocean dynamics raises the following questions: • What atmospheric forcing product should be used and how? • When forcing the ocean with absolute wind, hence ignoring the CFB, to what extent do these simulations overestimate the mean currents and the mesoscale activity? • How to accurately take into account the CFB in an uncoupled oceanic simulation? To what extent can the parameterizations proposed in  and (Renault, McWilliams, & Masson, 2017) can mimic a coupled ocean-atmosphere model that includes the CFB?
This study first aims to quantify the biases made on the mean and mesoscale oceanic circulations when forcing an oceanic model with common reanalyses (e.g., CFSR; Saha et al., 2010) or scatterometer products (e.g., QuikSCAT; e.g., Bentamy et al., 2013), with either absolute or relative winds, and then it goes a step further than  and (Renault, McWilliams, & Masson, 2017) by testing the validity of the proposed parameterizations of the wind and stress responses to the CFB. The main goal of this study is to assess to what extent a forced oceanic simulation can mimic the dynamics of a coupled simulation that includes the CFB. Quasi-global ocean-atmosphere coupled simulations (Samson et al., 2017; are first used to mimic the atmospheric forcing from scatterometers and reanalyses by generating synthetic atmospheric forcing. An additional set of uncoupled oceanic simulations is then made forced by this synthetic forcing (Section 2) using or not the parameterizations proposed by  and (Renault, McWilliams, & Masson, 2017).
The following four primary diagnostics are used (Section 2): • magnitude of the vertically integrated current; • coupling coefficient between surface current vorticity and surface stress curl (s ); • eddy wind-work sink of energy induced by the CFB (F e K e ); • mesoscale activity (Eddy Kinetic Energy, EKE).
The oceanic (atmospheric) component uses an Arakawa-C grid, based on a Mercator projection at 1/12 • (1/4 • ) resolution. The geographical domain of this coupled model is a tropical channel extending from 45 • S to 45 • N, with the oceanic grid being a perfect subdivision by three of the atmospheric grid. The ocean vertical grid has 75 levels, with 25 levels above 100 m and a resolution ranging from 1 m at the surface to 200 m at the bottom. The atmospheric grid has 60 eta levels with a top of the atmosphere located at 50 hPa. The WRF default vertical resolution has been multiplied by three below 800 hPa. Thus, the first 33 levels are approximately located below 500 hPa with the first eta base level located at 10 m above the ocean.
As in Samson et al. (2017), the physical package we use in WRF is the longwave Rapid Radiative Transfer Model (RRTM; Mlawer et al. (1997), the "Goddard" Short Wave (SW) radiation scheme (Chou & Suarez, 1999), the "WSM6" microphysics scheme (Hong & Lim, 2006), the Betts-Miller-Janjic (BMJ) convection scheme (Betts & Miller 1986;Janjić, 1994), the Yonsei University (YSU) planetary boundary layer scheme , the unified NOAH Land Surface Model (LSM) with the surface layer scheme from MM5 (Chen & Dudhia, 2001). The planetary boundary layer scheme MYNN2.5 (Nakanishi & Niino, 2006) is also used instead of YSU in an additional specific simulation. Note that in WRF (since Version 3.3.1), the YSU scheme has been modified following Shin et al. (2012) to overcome issues related to a first eta level situated around 10-m. WRF lateral boundary conditions are prescribed from the European Centre for Medium-Range Weather Forecasts ERA-Interim 0.75 • resolution reanalysis (Dee et al., 2011) at 6-hourly intervals.
The ocean physics used in NEMO corresponds to the Upstream-Biased Scheme (UBS; Shchepetkin & McWilliams, 2009) advection for the tracers and the dynamics with no explicit diffusivity and viscosity. The vertical eddy viscosity and diffusivity coefficients are computed from a TKE turbulence closure model (Blanke & Delecluse, 1993). The surface boundary condition on momentum is the surface stress, which is applied as the surface boundary condition on the momentum vertical mixing. The oceanic open boundary conditions are prescribed with an interannual global 1 4 • DRAKKAR simulation (Brodeau et al., 2010). In order to benefit, at a limited cost, from a fully spun-up mesoscale circulation in the initial conditions of the 1 12 • ocean, we first run a 5-year forced ocean simulation initialized from 1 4 • DRAKKAR simulations.

Forced Oceanic Model
To be able to compare the forced simulations to the coupled simulations, it is crucial to have an oceanic model configuration as close as possible to the coupled simulations configuration. Therefore, the forced ocean simulations were performed with NEMO 3.4 using the very same configuration and the exact same initial condition as the coupled simulations. All the forced simulations described hereafter differ only in the atmospheric forcing used or the way they estimate the surface stress. In order to run forced ocean simulations with characteristics as close as possible to the coupled model simulations, we had to modify the bulk formula available in NEMO 3.4. To do so, we used the AeroBulk package (Brodeau et al., 2017, now included in NEMO v4), which offers the possibility to use the bulk formula COARE v3.0 (Fairall et al., 2003), as in our

10.1029/2019MS001715
coupled simulations that use the WRF bulk formulae. We also used hourly forcing fields, which corresponds to the coupling frequency of the coupled simulations. The numerical outputs for the solutions are daily averages.

Surface Wind and Stress
In a forced oceanic model, surface wind stress is either directly prescribed or computed from wind usually taken at 10 m to be compliant with the parameters used in the bulk formulae. Two kinds of prescribed 10-m winds can be defined: • U 10abs : the absolute wind at 10 m, • U 10rel : the 10 m wind relative to the oceanic surface current: Usually, the atmospheric models (such as WRF) provide a 10-m wind diagnostic (U10 and V10 in WRF), which represents in fact U 10rel when ocean surface currents seen by the atmosphere are not set to zero (see also . U 10abs can be simply reconstructed by adding up U o to U 10rel . The surface stress is computed in WRF based on u * , the friction velocity defined in the Monin-Obukhov similarity theory. Because the same vertical scaling, based on the Monin-Obukhov length scale, is used to compute u * and (U10,V10) from the first layer model information, the surface stress can be given as a function of variables given either at 10 m or at the first level of the model. When the surface current is taken into account in WRF, the surface stress is therefore defined by the following equation: Note that to properly take into account the impact of the oceanic surface current in the atmosphere, we must also modify the tridiagonal matrix system solved in the vertical turbulent diffusion scheme (Lemarié, 2015;Renault, Lemarié, & Arsouze, 2019).
If the surface currents are considered as zero in WRF, the same equation applies but with U o = 0, implying U 10rel = U 10abs and so Scatterometer winds are estimated from the pseudo-stress and are also reported as a 10-m equivalent neutral wind relative to the oceanic current, that is, a relative wind to the oceanic current that would exist if the conditions were neutrally stable (Plagge et al., 2012). The 10-m equivalent neutral wind long-term mean is very similar to the 10-m long-term mean (not shown). However, at the mesoscale, because atmospheric stratification can deviate significantly from neutral conditions, the 10-m equivalent neutral wind response to the TFB and to the CFB can be overestimated by 10% to 25% with respect to 10-m wind response (O'Neill et al., 2012;Perlin et al., 2014;Song et al., 2009). Note that an oceanic simulation forced by such a wind product is not affected by this mesoscale overestimation as the simulated eddies are not correlated with those that have been seen by the scatterometers.

Spatial Filtering
As in Seo (2017) and , the mesoscale anomalies are isolated from the large-scale signal by using a spatial filter. The following filter description is derived from  with minor modifications. A field is smoothed using a Gaussian spatial filter with a standard deviation of 4 (12) grid points at 1/4 • (1/12 • ). Gaussian weights of points located at a distance larger than 3 are considered zero. The Gaussian filter is thus applied on a (6 + 1) × (6 + 1) window, which makes 25 × 25 points at 1/4 • , or 73 × 73 points at 1/12 • , and corresponds to 670 km (475 km) at the equator (45 • N) because of the Mercator projection (grid size scales with the cosin of the latitude). One of the properties of the Gaussian filter is that its Fourier transform, and so its spatial frequency response is also a Gaussian function. This allows an analytical definition of the cutoff wave number for a given reduction of the filter response. For example, the power spectrum of the filtered signal will be half of the original signal (−3 dB) at the wave number of about 1/600 km −1 at the equator and 1/420 km −1 at 45 • N. It will be divided by 10 (−10 dB) at the wave number of about 1/330 km −1 at the equator and 1/280 km −1 at 45 • N. Land points are treated as missing data, and the weights of windows including land points are renormalized over the remaining oceanic points. Mesoscale anomalies of a field are then defined as ′ = −[ ], with [ ] the smoothed field. This filter is applied on one coupled simulation (see Section 3) and to the coupling coefficients described hereafter (see Section 2.6).

Mean Oceanic Circulation
The CFB causes a slowdown of the mean oceanic circulation. To characterize this effect in our simulations, the magnitude of the vertically integrated current (m 2 /s) is estimated as where the mean· is defined with respect to long-term temporal averaging (5 years of simulations), u and v are the zonal and meridional currents at each depth, and H the bottom of the ocean. V 5y is the barotropic transport that flows through a section of 1 m long that is perpendicular in any point to the barotropic flow.
To assess the temporal evolution of the mean circulation, the magnitude of the vertically integrated current (m 2 /s) is also estimated using a 3-month running window: where the mean < · > is defined with respect to the 3-month running window. Note that the estimation of the magnitude of the vertically integrated current is nonlinear. As a result, a long-term temporal average of V 3m is systematically stronger than V 5y .

CFB Coupling Coefficients
As in (Renault, McWilliams, & Masson, 2017) and , the coupling coefficient between mesoscale current vorticity and surface stress curl s is defined as the slope of the robust regression (Maronna et al., 2006) between surface stress curl and oceanic current vorticity. Following (Maronna et al., 2006), our M-estimation is based on a biweight function. s is evaluated at each grid point. In order to isolate the stress response due to the mesoscale oceanic surface current, the fields are first temporally averaged using a 29-day running mean to suppress the weather-related variability (Chelton et al., 2007), and the large-scale signal is removed using a high-pass Gaussian spatial filter (as in, e.g., Seo, 2017; see filter description in Section 2d). Note that a, for example, 15-day running mean does not efficiently suppress the weather-related variability.  demonstrate s properly isolates the stress response to the CFB from the TFB. The coupling coefficient s w is estimated as s but using U 10abs instead of the surface stress.

Sink of Oceanic Energy by the Mesoscale CFB
The CFB induces transfers of energy from the mesoscale oceanic currents to the atmosphere. These sinks of oceanic energy are crucial to represent in a oceanic model as they partly control the oceanic mesoscale activity. Temporal filters are commonly used to estimate the eddy work. However, they do not allow proper isolation of the sinks of energy induced by the CFB, because the characteristic lifetime of ocean eddies (usually >90-day) does not allow to filter out wind stress with large spatial scales. In this study, in order to isolate these sinks of energy, we define the mesoscale (eddy) wind work F e K e [m 3 s −3 ] using a spatial filter: where the mean is defined with respect to long-term averaging (5 years of simulations) and ′ is defined as the signal anomalies as estimated using the spatial filter described in Section 2d. Note that this spatial-filter-based definition differs from the usual Reynolds decomposition. For instance, cross-terms do not completely vanish (but remain very weak; see the supporting information). However, it allows better isolation of the oceanic mesoscale signal from large-scale currents that can have a duration of, for example, less than 1 month (e.g., wind-driven current). Nevertheless, similar results can be found using a temporal filter and by estimating difference between a simulation with CFB and a simulation without CFB (Jullien et al., 2020).

Eddy Kinetic Energy
The Eddy Kinetic Energy (EKE, [m 2 /s]) represents the intensity of the mesoscale activity. In this study, it is computed from the geostrophic currents anomalies estimated using the filter described above.

Journal of Advances in Modeling Earth Systems
Again, this spatial-filter-based definition of EKE differs from the usual Reynolds decomposition, and cross-terms such asŪ ′ do not vanish over standing eddies.

Coupled and Forced Simulations
We define the following naming convention to identify the different coupled and forced simulations performed: • The first capital letter defines the type of atmospheric forcing used: C is a coupled simulation; R stands for Reanalysis (as mimicked by a coupled simulation), that is, atmospheric data derived from a coupled simulation; and Sc stands for Scatterometer (as mimicked by a coupled simulation). • The first subscript word indicates whether the coupled model data used for forcing takes into account or not the CFB (i.e., CFB or NOCFB).
The following is the case of forced simulations only: • The second capital letter indicates if the surface stress used to force the the ocean model dynamics is directly prescribed (S) or computed in the ocean model bulk formulae from a prescribed wind speed (W). • In simulations using bulk formulae to compute the wind stress, the following lowerscript, REL or ABS, indicates whether we use a wind relative to currents of the forced ocean simulation or not.
Finally, when relevant, a third capital letter P indicates that a parameterization for the wind response is used in forced simulations, followed by the name of the parameterization in subscript (see Section 3.4 for more details).

Coupled Simulations
As described in Table 1, two reference 5-year coupled simulations are performed over the 1989-1993 period, only differing by their degree of coupling: • In C CFB , WRF gives NEMO hourly averages of freshwater, heat, and momentum fluxes, whereas the oceanic model sends back to WRF the hourly averaged Sea Surface Temperature (SST) and surface currents. The surface stress is estimated using the wind relative to the oceanic current U 10rel . • In C NOCFB , the oceanic model sends back to WRF only the SST. WRF sees zero surface current and the surface stress is a function of U 10abs .
Another simulation has been carried out to briefly assess the mesoscale TFB impact on the mean and mesoscale oceanic circulations: • In C CFB_SMTH , the setup is the same than for C CFB except that the SST sent to WRF is spatially filtered (see description of the filter in Section 2) in order to smooth out the thermal mesoscale structures. C CFB_SMTH thus contains the CFB and the large-scale TFB (as in, e.g., Seo, 2017). This simulation is used to briefly assess the mesoscale TFB impact on the mean and mesoscale oceanic circulations.
Additional coupled simulations have been carried out to assess the sensitivity of the results to the physics taken into account in the model and to the internal variability of the oceanic model (see the supporting information): • The simulation C CFB_MYNN is the same than C CFB except that the MYNN2.5 (Nakanishi & Niino, 2006) Planetary Boundary Layer scheme is used instead of YSU in the atmospheric model. It allows an assessment of the sensitivity of the results to the Planetary Boundary Layer used in the atmospheric model in Supplemental Information. • Two more coupled simulations identical to C CFB have been run on different machines and using different compiler options of optimization (C CFB2 , C CFB3 ). This adds very small perturbations all along the simulation run. This set of three identical coupled simulations defines a small ensemble that is used to provide an estimate of the robustness of the results and of the role of the internal variability in the models (see also Figures S2 and S3). Note that because of the small number of experiments, the internal variability may be underestimated.
All the experiments consider the large-scale thermal feedback.

Mimicking Atmospheric Reanalysis
To mimic the different ways of forcing an oceanic model with a reanalysis, five forced oceanic simulations are carried out (see also Table 1): • R NOCFB S: The model is forced by a surface stress pNOCFB that is directly provided by C NOCFB , that is, a simulation that did not feel the CFB. The surface stress used in this simulation mimics that from a reanalysis such as CFSR or ERA-Interim. • R NOCFB W ABS : The model is forced by U 10abs derived from C NOCFB , mimicking U 10abs from a reanalysis such as CFSR or ERA-Interim. The stress is estimated in the NEMO bulk formulae without taking into account the surface currents of the forced simulation (W ABS ). • R CFB W ABS : The model is forced by U 10abs derived from C CFB , that is, a simulation that felt the CFB. U 10abs mimics the wind that could be provided by future reanalysis (that would take into account the CFB) or by a regional atmospheric simulation that takes into account the CFB. The surface currents of the forced simulation are not taken into account when estimating the stress in the NEMO bulk formulae (W ABS ). • R NOCFB W REL : Same as R NOCFB W ABS , the stress is estimated using U 10abs from C NOCFB , but relative to the surface currents simulated by the oceanic forced simulation. • R CFB W REL : Same as R CFB W ABS , the stress is estimated using U 10abs from C CFB , but relative to the surface currents simulated by the oceanic forced simulation.
An additional simulation forced by the stress from a reanalysis that felt the CFB could be done, but it would correspond to the Sc CFB S simulation described hereafter. Note that all the forced oceanic simulations using NEMO bulk formulae and presented in this subsection use U 10abs .

Mimicking Scatterometer Winds
We saw in Section 2.3 that scatterometer winds can lead to two different errors: an overestimation of the coupling coefficient between mesoscale surface current vorticity and wind curl and a slightly larger long-term averaged wind. However, when forcing an oceanic model with scatterometer wind (and without data assimilation), because the simulated eddies are not correlated with those that have been seen by the scatterometers, the main source of error therefore lies in the fact that the long-term mean of equivalent neutral scatterometer wind is relative to the oceanic motions. In the following experiments, for simplicity and to limit the number of close experiments, U 10rel is used instead of the equivalent neutral wind relative to the oceanic current. The use of the equivalent neutral wind relative to the oceanic motions would lead to similar results.
To mimic the different way of forcing an oceanic model with a scatterometer product, three more simulations are carried out (see also Table 1): • Sc CFB S: The model is forced by a surface stress pCFB that is directly provided by C CFB , that is, a simulation that felt the CFB. The forcing mimics scatterometer surface stress. • Sc CFB W ABS : The model is forced by an U 10rel wind derived from C CFB , mimicking scatterometer wind. This wind is relative to the oceanic current simulated by the C CFB . The surface currents of the forced simulation are not taken into account when computing the wind stress from U 10rel in NEMO bulk formulae (W ABS ). • Sc CFB W REL : Same as Sc CFB W ABS but the surface currents of the forced simulation are taken into account when computing the wind stress from U 10rel in NEMO bulk formulae (W REL ).
Note that although atmospheric products that incorporate scatterometer winds such as CFSR or JRA55-do (Tsujino et al., 2018) should act in a similar way as QuikSCAT, they do not strongly resemble scatterometer data, for example, at ocean fronts, suggesting that the assimilation does not strongly constrain the near-surface JRA55-do (Abel, 2018).  and (Renault, McWilliams, & Masson, 2017) have identified relationships linking the wind and stress response to the oceanic surface currents. To mimic the CFB in a forced oceanic model, they suggest a wind-correction approach or a stress-correction approach. We propose now to test the different approaches in ad hoc parameterizations using atmospheric wind derived from the coupled simulations.  suggest a simple parameterization based on the current-wind coupling coefficient s w (see Section 2.6) estimated from a coupled simulation (C CFB ). s w corresponds to the linear wind response to a given current:

The Wind-Correction Approach
where U o is the surface current and Ua ′ is the wind response to U o . If the wind cannot see the surface current, then s w = 0 (no wind response to the current). If s w = 1, in a forced ocean model, this would be an ideal case, in which there is no loss of energy; that is, the current generates a wind response with a magnitude equal to the current magnitude. When forcing an oceanic model with a product that did not account for CFB, the surface stress should be computed with a wind relative to the current corrected by s w : U 10abs −(1 − s w )U o , which is then used to obtain in the bulk formula (1). Two forced simulations are carried out using this parameterization. The first letter of the simulation name stands for parameterization; the second letter indicates the kind of parameterization used: • R NOCFB W REL Ps w C: a simulation with a constant s w = 0.3, which corresponds to the global mean value of s w as estimated from C CFB . U 10abs is derived from C NOCFB . • R NOCFB W REL Ps w M: a simulation that takes into account monthly and spatial variations of s w . In this simulation, the monthly maps of s w , estimated from C CFB , are provided to the oceanic simulation. U 10abs is derived from C NOCFB .

The Stress-Correction Approach
Alternatively, (Renault, McWilliams, & Masson, 2017) proposed another simple parameterization based on the current-stress coupling coefficient s (see Section 2.6), in a forced oceanic model:

10.1029/2019MS001715
where is the surface stress that includes the CFB effect, 0 is the surface stress that does not include the CFB effect (it can be prescribed or estimated within a bulk formulae), and ′ is the stress response to U o : s can be estimated from a coupled simulation (here C CFB ) or, as demonstrated by (Renault, McWilliams, & Masson, 2017), can be predicted using the magnitude of the wind: with = −2.9 × 10 −3 N s 2 /m 4 , and = 0.008 N s/m 3 . and are derived from the linear regression between |U 10abs | and s and have uncertainties in the range of 20% (see the supporting information). Additionally, for weak wind (i.e., <3 m/s), the regression is not valid anymore and s is taken as a constant of −0.0007 N s/m 3 that corresponds to its value at 3 m/s. Note that (10) can also be estimated from the observations; however, such an estimation suffers from large uncertainties (Renault, McWilliams, & Masson, 2017) because it is based on two different products (AVISO, Ducet et al., 2000, andQuikSCAT, Bentamy et al., 2013), which do not have the same effective spatial resolution. This parameterization (10) represents more physics than using (9) with a time-invariant or climatological s as it allows us to represent the spatial and temporal variations of s and, contrarily to s w , does not suffer the quadratic errors induced by the bulk formula. It is also very flexible as it can be used on top of a prescribed wind stress p or a stress estimated from the absolute wind in a bulk formula. For the main simulations that use a parameterization based on s , the atmospheric stress used to force the model is derived from the C NOCFB simulation, that is, an atmospheric simulation that did not previously feel the CFB. Based on s parameterization, two simulations are carried out (see also Table 1). Again, the first letter of the simulation name stands for parameterization; the second letter indicates the kind of parameterization used: • R NOCFB SPs M: A simulation that is forced by a surface stress pNOCFB and that considers monthly and spatial variations of s derived from C CFB . • R NOCFB SPs P: A simulation that is forced by a surface stress pNOCFB and estimates s using |U 10abs | (10).
Note that estimating the stress in the ocean model using U 10abs instead of a prescribed surface stress would lead to similar results. As detailed in Section 9, a final forced experiment R CFB SPs P has been carried out using the parameterization (10) but with an atmospheric stress forcing pCFB as derived from C CFB . This experiment mimics a reanalysis that felt the CFB (e.g., Coupled-ECMWF reanalyses CERA-SAT and CERA-20C).

Representation of the Surface Stress and Wind Curl in Reanalyses and Scatterometer Products
Scatterometers are fundamentally stress-measuring instruments. They have been extensively used to characterize the air-sea interactions (Chelton et al., 2007;O'Neill et al., 2012; and also to force oceanic models by using their derived wind or stress. Surface stress derived from scatterometer such as QuikSCAT incorporates all the ocean-atmosphere interactions, including the CFB effect on the stress (Chelton et al., 2004;Cornillon & Park, 2001;Kelly et al., 2001;Renault, McWilliams, & Masson, 2017;. As illustrated in Figure 1a and previously shown by Chelton et al. (2001), the surface currents have a significant imprint on the surface stress curl around the Gulf Stream path where a negative current vorticity causes a positive stress curl anomaly. As a result, on the western side of the Gulf Stream, the mean surface stress curl is negative, whereas it is positive on the other side. The wind response to the CFB should have the opposite sign than the stress response . However, consistent with  and as illustrated in Figure 1b, because scatterometer winds are U 10rel , its response to the CFB has the wrong sign, misdiagnosing an increase of the wind over the Gulf Stream instead of a decrease.
Global reanalyses such as CFSR (Saha et al., 2010) or ERA Interim (Dee et al., 2011); or a simulation such as C NOCFB ) are commonly used to force an oceanic model. CFSR couples ocean model SST to the atmosphere, but not ocean currents, whereas ERA Interim is atmosphere only with prescribed SST. Although both reanalyses assimilate data, they do not represent the current imprint on the surface stress nor wind (Belmonte Rivas & Stoffelen, 2019). This is illustrated in Figures 1c and 1d: over the Gulf Stream, because of the lack of CFB, these kind of reanalysis do not represent the imprint of the currents on the wind nor on the stress (but do represent the TFB effect). Some global coupled models (or simulations such as C CFB and recent reanalysis such as CERA-SAT or CERA-20C) take into account the CFB. In such simulations, both surface stress and wind are more realistic as they are marked by the presence of the surface currents: consistent with scatterometer stress, the surface stress curl is marked by dipole of negative and positive stress curl ( Figure 1a). The same behavior can be observed to a lesser extent in wind curl (with opposite signs) with more of a positive sign ( Figure 1e) than a negative sign. Similar results can be found at the mesoscale over, for example, an eddy or a filament.

CFB Effect at the Large Scale
V 5y from C CFB is represented in Figure 2a. The main ocean circulation features are seen in C CFB as, for example, the equatorial currents, the mean gyres, and the Western Boundary Currents. As revealed by the regression between V 5y from C NOCFB and C CFB (Table 2), although C CFB and C NOCFB simulations have a very similar mean circulation pattern, C NOCFB tends to slightly overestimate its strength (see Figure 2b). This is confirmed by Figure 2b that shows the box plots of V 5y as estimated from the coupled simulations. Consistently, from C NOCFB to C CFB , the V 5y distribution is shifted down because the CFB reduces the mean surface stress and slows down the mean currents (Luo et al., 2005;Pacanowski, 1987). As a result, the mean and median V 5y are both overestimated in C NOCFB by ≈ 8% with respect to C CFB (see also Table 3), which is larger than the internal variability (less than 2.5%, estimated from simulations C CFB , C CFB2 , C CFB3 ) Figure 3a represents the monthly evolution of V 3m . Starting from the same initial state, it confirms that the CFB, by reducing the surface stress, induces a slow down of the mean circulation. C CFB and C NOCFB differ after only a couple of months and the difference between V 3m from C CFB , and C NOCFB is relatively stable after ≈10 months: V 3m in C NOCFB is ≈11% larger than C CFB .

CFB Effect at the Mesoscale
The spatial distribution of the long-term average s from C CFB is illustrated in Figure 4a. Figure 5a depicts the box plots of s from C CFB , and C NOCFB .  provide a full assessment of s as simulated by C CFB . s shows temporal and spatial variability, which is driven by the 10-m wind: The larger the wind, the more negative the s (Renault, McWilliams, & Masson, 2017). The larger values of s are therefore mainly situated in the high-latitude regions. On a long-term average, in C CFB , it varies from approximately The box plot of V 5y estimated from C CFB is shown in all the box plots for clarity. The line that divides the box into two parts represents the median of the data, whereas the green dot represents its mean. The length of the box shows the upper and lower quartiles. The extreme lines represent the 95 th percentiles of the distributions. The 5 th percentiles have values so close to 0 (due to large areas very close to zero) that they are located on the x axis and are therefore not visible. Statistical significance of the medians and means shown in the box plots have been tested with a bootstrap method. For all the box plots depicted here and in the following figures, the 95 medians are indistinguishable from the symbol used.
−0.03 to approximately 0 N s/m 3 . Consistent with the literature, in C NOCFB , that is, a simulation that ignores the CFB, s ≈ 0. Figure 6a depicts the spatial distribution of F e K e as estimated from C CFB , and Figure 7a represents the box plots of F e K e as simulated by C CFB and C NOCFB . In agreement with, for example, Xu and Scott (2008) and (Renault, McWilliams, & Masson, 2017), in C CFB F e K e is negative almost everywhere, expressing the large-scale sinks of energy induced by the CFB. In a coupled simulation that neglects the CFB, even if the TFB is taken into account (i.e., C NOCFB ), F e K e ≈ 0. In C CFB F e K e has a global mean and a median of −3.9 × 10 −7 m 3 /s 3 and −1.3 × 10 −7 m 3 /s 3 (Table 3). The largest sinks of energy are confined to the 5 th percentile of the distribution (Figure 7a and Table 4). They reveal locations that are characterized by an intense mesoscale activity and a notably negative s and actually roughly represent the Western Boundary Currents and the Agulhas Return Current. Over these regions F e K e has a mean and a median of −37.7 × 10 −7 m 3 /s 3 and 26.8 × 10 −7 m 3 /s 3 , respectively.
The spatial distribution of the EKE from C CFB is illustrated in Figure 8a; it shows a general agreement with the literature. Figure 9a shows the box plots of the EKE as estimated from C CFB and C NOCFB , and Figure 10a depicts bivariate histogram of the same quantity. The EKE spatial patterns from C CFB and C NOCFB are very similar (see also Table 2), representing, for example, in the same place the eddy-rich regions such as the Western Boundary Currents. However, in agreement with previous studies, C NOCFB systematically overestimates the EKE with respect to C CFB . Consistently, from C NOCFB to C CFB , the mean and median EKE are reduced by ≈27% and 37.5%, respectively (Table 3). These reductions are much larger than the internal variability (less than 2%; see Table 3). The 95 th percentile roughly corresponds to Western Boundary Currents and the Agulhas Return Current, where the mean and median EKE are reduced by ≈15% and ≈20% from C NOCFB to C CFB , which is again larger than the internal variability of the coupled model (less than 2.2%; see Table 4).

Thermal Feedback
The mesoscale TFB causes wind and surface stress magnitude, divergence, and curl anomalies (Chelton et al., 2004(Chelton et al., , 2007O'Neill et al., 2010O'Neill et al., , 2012Desbiolles et al., 2016). Hogg et al. (2009), by using a "crude" parameterization of these surface stress anomalies in an idealized oceanic model, suggest that the mesoscale TFB may weaken the mean circulation by 30% to 40%. However, here, the C CFB_SMTH simulation, which does not include the mesoscale SST feedback, has a mean circulation very similar to that from C CFB (Figure 2a and Tables 2 and 3). Additionally, in Figure 3a, the monthly temporal evolution of V 3m from C CFB and C CFB_SMTH cannot be objectively distinguished. As shown by  and illustrated in Figure 5a, C CFB_SMTH has a s that is very similar to that from C CFB . This is explained by the fact that the surface stress anomalies induced by the mesoscale TFB are barely collocated with the surface currents. As a result, they do not systematically induce conduits of energy between the ocean and the atmosphere. In C CFB_SMTH both F e K e and EKE in C SMTH are very similar to those from C CFB (Figures 7 and 9 and Tables 2, 3, and 4). While spatial differences of F e K e between C CFB and C CFB_SMTH over the whole domain are not distinguishable from differences caused by the internal variability of the model (See Figures S2b and S2d), they reveal sizable difference over Western Boundary Currents and the Agulhas Return Current that are not only due to the internal variability (see the positive values of Figures S2c and S2d). However, they remain much weaker than differences caused by the CFB. This is consistent with the fact that in C NOCFB , that is, a simulation that considers only the TFB, F e K e ≈ 0 (Figure 7a). Note that in C CFB_SMTH , there is also a significant positive F e K e along the coast; however, it is not clear whether it is a filtering artifact or a real signal, but this signal is not present in, e.g., R NOCFB SPs P, which does not consider TFB either. The mesoscale TFB impact on the oceanic dynamics is beyond the scope of this study, but these first diagnostics suggest that in a forced oceanic simulation (with CFB), the lack of mesoscale TFB has a very weak effect on the mean and mesoscale oceanic circulations. Note that when large-scale gradients of SST are also smoothed (as, e.g., in Ma et al., 2016), this may impact the wind, and subsequently, the coupling coefficient s (as it depends on the wind magnitude), and, finally, F e K e .

Forcing an Oceanic Model with an Atmospheric Reanalysis
In this section, the simulations forced by a synthetic reanalysis product (Table 1), namely, R NOCFB S, R NOCFB W ABS , R CFB W ABS , R NOCFB W REL , and R CFB W REL ) are compared to C CFB . The aim of this section is to assess the biases caused by forcing an oceanic simulation with an atmospheric reanalysis in terms of mean circulation, coupling coefficient between the mesoscale surface currents vorticity and stress curl, sinks of energy, and mesoscale activity. Table 2 details the regressions of V 5y and EKE from C CFB and the Reanalysis-like simulations. In all the forced simulations, both V 5y and EKE have a spatial pattern similar to those from C CFB (see also Figure 10 for the EKE). Tables 3 and 4 provide some mean and median values of V 5y , F e K e , and EKE. Figure 2c shows the box plots of V 5y as estimated from R NOCFB S, R NOCFB W ABS , R CFB W ABS , R NOCFB W REL , and R CFB W REL . Not surprisingly, two kind of simulations can be distinguished: the simulations forced with or without CFB. Consistent with the previous section, the oceanic simulations R NOCFB S, R NOCFB W ABS , R CFB W ABS , which ignore the CFB when estimating the surface stress, overestimate the mean and median of V 5y by at least 5.5% (and up to 11.5%; which is larger than the internal variability, see also Table 3). In that sense, they are very similar to C NOCFB (see Figure 2b and Table 3). By contrast, the simulations R NOCFB W REL and R CFB W REL , which take into account the CFB when estimating the surface stress, underestimate V 5y mean Table 3 Mean (median) by 8.2% and 5.3% (7.6% and 3.5%), respectively, again larger than the differences between C CFB and the perturbed coupled simulations (Table 3). This underestimation is explained by the fact that these simulations lack or partially lack a coherent wind response to CFB, which should partly re-energize the mean currents . Figure 3b represents the monthly evolution of V 3m . Consistent with the previous results, after about a year, the first/main temporal adjustment (i.e., a large increase of V 3m ) of the mean circulation is achieved, and the differences between the experiments are robust and clearly visible. It is worth noting that the atmospheric forcing derived from C CFB already contains a mean wind response that is coherent with the C CFB mean surface currents. As a consequence, the mean winds used to force R CFB W ABS (R CFB W REL ) are slightly re-energized in comparison with the forcing used in Figure 3b). The mean oceanic circulation over the 5-year period considered here is not fully deterministic. As a result, the imprint of the C CFB mean surface current on the wind is not fully coherent with those from the forced simulations and none of the forced simulations reproduce the correct value of V 3m . This has an important consequence, as it means that a wind and stress that have an imprint of the surface currents cannot be directly used with a parameterization of the CFB. As discussed in Section 9, they should be first corrected for the influence of the currents on the reanalysis.

Mesoscale
At the mesoscale, again, the simulations forced with and without CFB can be distinguished. On the one hand, the simulations forced without CFB behave in a very similar way as C NOCFB : The resulting stress does not have the imprint of the simulated mesoscale eddies and, thus, does not have any mesoscale correlation with them. As a result, s ≈ 0, F e K e ≈ 0, and global mean and median of the EKE is overestimated by at least 34% and 54%, respectively, compared with C CFB (Figures 4, 5b, 6, 7b, 8, and 9b and Table 3). Consistent with the literature, oceanic simulations forced without CFB overestimate not only the mean circulation but also the mesoscale activity. This overestimation is confirmed by estimating bivariate histogram between, for example, R NOCFB W ABS and C CFB (Figure 10b; similar results are found for the other experiments of the same kind). As C NOCFB , the forced simulations that ignore the CFB largely overestimate the EKE. On the other hand, the simulations forced with CFB, that is, R NOCFB W REL and R CFB W REL , have an opposite behavior. As illustrated in Figure 5b, in this simulation, s is overestimated (i.e., too negative) in term of mean, median, quartiles, and 95 th percentile by roughly 30%, which is consistent with . Note that when the wind is derived from C CFB , it already contains the mesoscale wind response to the mesoscale eddies simulated by C CFB . Even if R CFB W REL is forced by the C CFB winds and uses the same initial conditions, R CFB W REL and C CFB eddies become uncorrelated after a few months. Therefore, the C CFB mesoscale wind anomalies do not properly counteract the surface stress response to the CFB and R CFB W REL and R NOCFB W REL have a very similar s (Figure 5b). The systematic underestimation of the EKE is also confirmed by the bivariate histogram between, for example, R NOCFB W REL and C CFB shown in Figure 10c.
s can be interpreted as a measure of the efficiency of the eddy killing. This is corroborated by Figure 6b that shows the spatial distribution of F e K e as estimated from R CFB W REL (similar results are found using R NOCFB W REL ) and by Figure 7b that depicts the F e K e box plots for these simulations (see also Table 3). When taking into account the CFB, the global mean (median) sinks of energy are overestimated by at least 33% (15%) with respect to C CFB . The most striking features are situated over the Western Boundary Currents (the 5 th percentile of the F e K e and 95 th percentile of the EKE), where the sinks of energy are clearly more intense in R CFB W REL and R NOCFB W REL than C CFB (see also Table 4). As expected, these sinks of energy cause excessive damping of the mesoscale activity everywhere with respect to C CFB , which is much larger than the internal variability of the coupled model (Figures 8a, 8c, and 9b and Tables 3 and 4). This is also consistent with the results of  for the California region.
To conclude this section, when using a reanalysis or a coupled simulation like C NOCFB to force an ocean model, at the large scale, neglecting the CFB leads to an overestimation of the mean oceanic circulation strength. However, when taking into account the CFB by using the relative wind to the oceanic motions in the surface stress estimation, the mean circulation becomes too weak. At the mesoscale, as expected, when ignoring the CFB (by using a prescribed surface stress or an absolute wind), the sinks of energy induced by the CFB are not reproduced, and the EKE is significantly overestimated. By contrast, the use of the wind relative to the oceanic motions causes an overestimation of s and F e K e and, thus, of the damping of the mesoscale activity with respect to C CFB . The prescribed wind does not have a coherent response to the simulated oceanic currents and, thus, does not partially re-energize the mesoscale activity. Note that because the same bulk formulae are used in both WRF and NEMO, the estimated surface stress in WRF is similar to the estimated surface stress in NEMO. As a result, the forced oceanic simulations without CFB are very close to C NOCFB . Forcing an oceanic model with a forcing derived from a reanalysis-like product (with or without CFB) and without a proper parameterization of the CFB is not suitable. However, as discussed in Section 8, such products can be used with a parameterization of the CFB.

Forcing an Oceanic Model with Scatterometer Winds or Stress
In this section, the simulations Sc CFB S, Sc CFB W ABS , and Sc CFB W REL are compared to C CFB . The goal of this section is to quantify the biases induced by forcing an oceanic simulation with scatterometer stress or wind in terms of V 5y , V 3m , s , the sink of energy, and mesoscale activity. Again, as revealed by Table 2, in all the forced simulations, both V 5y and EKE have a spatial pattern similar to those from C CFB . Tables 3 and 4 provide statistics comparison to other simulations.

Large Scale
As for the reanalysis forcing, the simulations forced with absolute wind or stress (i.e., without CFB: Sc CFB S and Sc CFB W ABS ) and with relative wind (i.e., with CFB: Sc CFB W REL ) clearly differ from the others. Figure 2c shows box plots of V 5y as estimated from C CFB and these simulations. When forcing an ocean model with scatterometer stress or absolute wind, the global mean of V 5y is overestimated with respect to C CFB , but only by ≈3%, that is, less than the overestimation observed in C NOCFB but still larger than the internal variability of the coupled model (Table 3). Two main factors can explain such a difference. On one hand, the mean stress and wind used here are derived from C CFB in order to mimic scatterometers; that is, the stress has the imprint of the mean surface currents, and the wind is relative to the mean C CFB surface currents. If the mean surface currents from C CFB were identical to those in Sc CFB S or Sc CFB W ABS , the CFB effect on the large scale would be correctly reproduced when forcing an ocean model with scatterometer data. On the other hand, over the 5-year period considered here, the large-scale currents are only partially deterministic. Therefore, although similar, the mean surface currents in these three scatterometer-like simulations are not identical to that from C CFB . As a result, when forcing an ocean model with a scatterometer-like forcing, the CFB slowdown of the mean circulation is only partly reproduced with respect to that simulated by C CFB . By contrast, when using the relative wind on top of scatterometer-like winds (Sc CFB W REL ), the effect of the CFB on the surface stress is potentially doubly accounted for (not exactly as the mean currents are not totally deterministic). Subsequently, the slowdown effect is overestimated, resulting in an underestimation of V 5y by ≈6.5% with respect to C CFB ; that is, in Sc CFB W REL , V 5y is weaker than that from R CFB W REL (Table 3). Similar results are found when assessing V 3m (Figure 3c). To sum up, from a large-scale circulation perspective, the relative wind cannot be used on top of scatterometer wind as the wind (and stress) already contains the surface stress response responsible of the slow down of the mean circulation.

Mesoscale
As explained in Section 6, at the mesoscale, as for a simulation forced by a reanalysis, neither Sc CFB S nor Sc CFB W ABS can dampen the mesoscale activity as the CFB mesoscale surface stress anomalies present in the prescribed surface stress are not correlated with the eddies simulated by the forced simulations. In both simulations, as illustrated in Figures 5c, 7c, and 9c, and Tables 3 and 4, s ≈ 0, F e K e ≈ 0, and the EKE is overestimated by ≈27% with respect to C CFB , which is again much larger than the internal variability ( Table 3). The bivariate histograms of the EKE between C CFB and Sc CFB S (Sc CFB W ABS ) are x similar to that between C CFB and C NOCFB (not shown). When taking into account the CFB by using the relative wind instead of the absolute wind (Sc CFB W REL ), again because the simulated eddies are not correlated with the prescribed wind, s is overestimated (i.e., too negative), F e K e is too large, and the EKE too weak (in both mean and median, again with differences larger than those cause by the internal variability of the coupled model). Not surprising, the bivariate histogram of the EKE between C CFB and Sc CFB W REL is very similar to that between C CFB and R NOCFB W REL .

Journal of Advances in Modeling Earth Systems
These results on the large scale and on the mesoscale have an important consequence on how to force an oceanic model. Wind or stress forcing derived from scatterometer data may reproduce relatively realistically the mean CFB effect on the large-scale circulation but not the eddy killing effect. It also means that a parameterization of the wind response cannot be used directly on top of a scatterometer forcing as, although it may reproduce the mesoscale effect (see Section 8), it would overestimate the CFB large-scale effect and, thus, the slowdown of the mean circulation. Similar results are found when forcing the oceanic model with wind or the stress.

Parameterizing the CFB in a Forced Oceanic Model
As demonstrated in the previous sections, an oceanic simulation should have a parameterization of the CFB in order to have a realistic representation of the oceanic mean and mesoscale circulations. Because of the absence of a current imprint on the wind or stress, wind reanalyses without CFB are the simplest choice (together with atmospheric uncoupled or coupled models without CFB) when combined with a parameterization of the wind response to the CFB. As described in Section 3, four additional simulations based on reanalysis-like wind (Table 1), namely, R NOCFB W REL Ps w C, R NOCFB W REL Ps w M, R NOCFB SPs M, and R NOCFB SPs P) are carried out to test different kinds of parameterization of the CFB. As for the other forced simulations, in all the parameterized simulations, V 5y and EKE have a similar spatial pattern to those from C CFB (see Table 2 and bivariate histograms in Figure 10defg). Additionally, Tables 3 and 4 summarize the mean and median quantities of V 5y , F e K e , and EKE estimated from these simulations at a quasi-global scale and over the Western Boundary Currents and the Agulhas Return Current.

Large Scale
The box plots of V 5y as estimated from C CFB and the parameterized simulations are illustrated in Figure 2e. Figure 3b represents the monthly evolution of V 3m . All the parameterized simulations correctly reproduce the characteristics of V 5y as simulated by C CFB , although, on average, they have a weak underestimation of mean and median by less than 3%. The temporal evolution of V 3m is also fairly well reproduced in all the simulations. Although more years of simulations should be used to draw a more robust conclusion, a simulation with a parameterization based on s is closer to C CFB than a simulation with a parameterization based on s w (see also Table 3 and 4). Figures 5d, 7d, and 9d represent the box plots of s , F e K e , and EKE from C CFB and from the four parameterized simulations. Tables 3 and 4 provide the mean and median of the F e K e and EKE estimated over the whole domain and over the Western Boundary Currents, respectively. To first order, all the parameterizations improve the realism of the simulation in terms of s , F e K e , and EKE with respect to a classic forced simulation with (or without) CFB (see previous section). The representation of the mean, median, quartiles, and 5 th or 95 th percentile (i.e., the Western Boundary Currents and the Agulhas Return Current) of these variables are clearly improved.

Best Parameterization
The parameterizations based on s lead in general to slightly more realistic results. This is likely because (i) the estimation of s suffers from less uncertainties than the estimation of s w because the stress response to the CFB is relatively larger than the subsequent wind adjustment and, thus, easier to properly isolate  and (ii) because of the quadratic nature of the stress in the bulk formula, the error and uncertainties related to s w are amplified. As a result, the scatterplot between wind curl and current vorticity has a larger spread than that between stress curl and current vorticity (see Figures 8 and 10 of , which may tend to underestimate s w and the partial re-energization of the large-scale circulation. However, at the mesoscale, there is an apparent contradiction. On one hand, s and F e K e are better represented in s -based parameterizations with respect to s w ones. On the other hand, as shown in Figure 9d, in R NOCFB W REL Ps w C, R NOCFB W REL Ps w M, and R NOCFB SPs M, although the sink of energy is either well represented (mean and median) or slightly too large (e.g., in the 5 th percentile), the EKE is too intense both in terms of global means/medians and also over the Western Boundary Currents and the Agulhas Return Current. By contrast, in R NOCFB SPs P, where s is predicted by the wind, the F e K e and the EKE are fairly well represented everywhere, including over the Western Boundary Currents and the Agulhas Return Current. We hypothesize that this counterintuitive result might be explained by nonlinearities in the ocean and during a typical eddy life, as follows: Over an eddy, the use of a constant or statistical s w or s does not allow representation of high-frequency large sinks of energy that happen in the presence of an intense wind burst: The more intense the wind, the larger the sink of energy (Renault, McWilliams, & Masson, 2017). Those large and brief sinks of energy can destabilize an eddy and, thus, diminish its lifetime. Over the lifetime of an eddy, the high-frequency variations of s (represented in C CFB and in R NOCFB SPs P) weakens the transfer of energy from the ocean to the atmosphere while enhancing the dissipation of the eddy energy in the ocean. These high-frequency variations of s and F e K e are not represented in R NOCFB W REL Ps w C, R NOCFB W REL Ps w M, nor R NOCFB SPs M. Indeed, in these simulations, over an eddy, the sink of energy is either constant or it represents the climatological variations of s , slowly diminishing the eddy energy. This effect can be interpreted as an efficiency of the destabilization of an eddy: s in R NOCFB SPs P (varying as a function of U 10abs ) may be more efficient in destabilizing an eddy than in the other parameterized simulations. This conjecture requires further investigation. The better behavior of R NOCFB SPs P in representing the EKE is also confirmed by analyzing the bivariate histograms between the EKE from C CFB and the parameterized experiments (Fig. 10). All the parameterizations improve the represention of the EKE compared to the other forced simulations. However, as revealed by Figure 10g, a parameterization based on a predicted s reduces the spread of the bivariate histogram and allows a better representation of the EKE over its entire energetic range. The remaining scattering in Figure 10g is mainly due to the chaotic nature of EKE as the bivariate histogram between C CFB and R NOCFB SPs P is very similar to that between C CFB and e.g., C CFB2 (see the supporting information).
Using a predicted s (as in R NOCFB SPs P) appears therefore as the best strategy to force an ocean model as it allows to better mimic the behavior of a coupled model, representing more physics and more processes than parameterizations based on simple constant coupling coefficients. Note that these differences among the simulations are even more present over the Western Boundary Currents, that is, over the most eddying region of the world. Figures 4c, 6c, and 8c illustrate the spatial distribution of the mean s , F e K e , and EKE as estimated from R NOCFB SPs P. They can be compared to C CFB and R NOCFB W REL (i.e., a simulation forced with a relative wind but without a parameterization) from which the same quantities are represented on the same Figures. Consistent with the box plot analysis, R NOCFB SPs P is very similar to C CFB and allows us to overcome the overestimation of those quantities when forcing an oceanic model with a relative wind.
In addition, as shown in Section 9, another advantage of R NOCFB SPs P lies in the fact that it can be used on top of a prescribed surface stress or an estimated surface stress. Finally, discrepancies observed between the parameterized simulations and C CFB especially over the Western Boundary Currents can also originate from a slightly different generation of mesoscale activity by baroclinic and barotropic instabilities (Renault, McWilliams, & Penven, 2017). However, this would have a second-order effect on the EKE with respect to the sinks of energy as the ensemble simulations of C CFB give very close results in terms of F e K e and EKE C CFB .

Forcing an Oceanic Model with Future Reanalysis
Global models (reanalysis and climate models) are already or will inevitably take into account the CFB. This is, for example, the case of the latest ECMWF (European Center Medium Weather Forecast) reanalysis.
Paradoxically, this could be problematic from an oceanic modeling forcing perspective. As for a scatterometer, the CFB effects on the stress and wind would have already been included in the wind and surface stress derived from the reanalysis. Forcing an oceanic model with this kind of atmospheric forcing, would result in lower quality of the representation of the deterministic currents (see Section 7). Usually, reanalysis products provide the surface stress and a diagnosed wind at 10 m, which is, if the CFB is accounted for, relative to the surface currents (U 10rel see sections 2). Note that in this case U 10rel already contains the wind response to the CFB. If using U 10rel , the surface stress should first be estimated as usual: rea = a C D |U 10rel |U 10rel . Then, the estimated or provided surface stress ( rea ) should be corrected by removing the influence of the reanalysis surface currents (U orea ) and then by adding up that from the simulated surface currents (U o ), following: If U 10rel is provided, then the surface stress should be estimated following: This procedure should allow the removal of the wind and stress anomalies induced by the oceanic surface currents of the reanalysis and replacing them with those induced by the simulated currents of the oceanic simulation. To test this strategy, a final forced oceanic simulation (R CFB SPs P) is carried out using the surface stress and oceanic surface currents from C CFB and applying equation (11). Figures 2f, 3e, 5e, and 9e and Tables 4, 3, and 2 show the main results at both the large-scale and the mesoscale. In R CFB SPs P the spatial distribution and statistics of both V 5y and V 3m are very similar to the ones simulated by C CFB and R NOCFB SPs P. This indicates that although the stress from the reanalysis already contains the effect of the mean currents, the stress correction (11) allows us to efficiently remove the large-scale CFB imprint of the reanalysis stress and to replace it by one that is coherent with the simulation.
At the mesoscale in R CFB SPs P, the reanalysis eddy imprint on the stress (wind) is also removed by using (11). However, because the eddies of the reanalysis are not in the same place as they are in R CFB SPs P, even without correcting the surface stress, the mesoscale CFB effect should be correctly represented compared to C CFB . This is confirmed by Figure 9f and Tables 3 and 4. In R CFB SPs P s , the sinks of energy and the EKE are also very similar to the ones simulated by C CFB and R NOCFB SPs P. To conclude, in the coming years, such a strategy could be easily applied; however, a sine qua non condition is that the atmospheric fields are provided along with the corresponding oceanic surface currents. A similar strategy could be applied to scatterometer data, but in this case, a coherent surface current measure should be provided (which is not possible yet with the existing satellites).

Summary and Discussion
Several recent studies demonstrate the importance of taking into account the effect of the surface oceanic current on the stress and on the wind (CFB) in order to realistically reproduce the oceanic mean circulation and mesoscale activity. These past studies raise the important question of how to best to force an oceanic model. In this study, using a realistic tropical channel (45 • S to 45 • N) domain for ocean-atmosphere coupled and uncoupled oceanic simulations, we assess to what extent a forced oceanic simulation can reproduce the mean oceanic circulation, the coupling coefficient between the surface oceanic currents and stress, the sinks of energy induced by the CFB, and the mesoscale activity compared to a coupled simulation that includes CFB.
Consistent with previous studies, we first confirm that the CFB slows down the mean circulation and dampens the mesoscale activity. We then assess to what extent an oceanic simulation can be forced by wind or stress as derived from an atmospheric reanalysis or scatterometer products such as QuikSCAT. Figure 11 shows a sketch that sums up the following. Three kinds of forced simulations are distinguished depending on the source of the wind or stress forcing: The wind (stress) can be derived from (1) a reanalysis that did not feel the CFB (labeled a-d, in Figure 11), (2) a reanalysis in which the CFB is accounted for (labeled e-i in Figure 11), and (3) a scatterometer. When using a wind from a reanalysis forcing product that does not include CFB, the forced simulations that use an absolute wind or a prescribed surface stress overestimate the mean circulation strength and the mesoscale activity with respect to a coupled simulation with CFB (a in Figure 11. How to best force an oceanic model? This sketch sums up the main methods. In order to illustrate the possible error made by the different strategies, two main indicators are used: the large-scale circulation and the mesoscale circulation. A blue color indicates an underestimation, the reddish colors gradually highlight an overestimation, and the white color means that the circulation is roughly correctly represented. "Presc." stands for a prescribed coupling coefficient (s w or s ) and "Predic." for a predicted one (s ). The warning sign indicates that coherent surface currents with the atmospheric forcing from the reanalysis-like or from the observations are needed. Such a product is not yet available for the observations (dashed lines). Best way to force an oceanic model is highlighted with bold lines. Figure 11). By contrast, when using the relative wind, they overestimate the slowdown of the mean circulation and the damping of the mesoscale activity induced by the CFB because of their lack of a wind response (d in Figure 11). If the reanalysis previously felt the CFB, because of the already accounted for presence of the large-scale wind response to the CFB, an absolute wind forcing will result in a large overestimation of the mean circulation strength (e in Figure 11) and in an overestimation of the mesoscale activity with respect to a coupled simulation with CFB.
We furthermore test the parameterizations of the CFB based on wind-correction or stress-correction approaches as suggested by  and (Renault, McWilliams, & Masson, 2017) and using an atmospheric forcing derived from a coupled simulation. These approaches are based on the use of a coupling coefficient between wind curl and current vorticity or stress curl and current vorticity. The coefficients can be estimated from a coupled simulation and prescribed to the forced oceanic simulation (c, i, and n in Figure 11) or predicted from the wind magnitude (b, f, and k in Figure 11). When forcing an oceanic simulation with a reanalysis product that does not account for the CFB, while both wind-correction and stress-correction approaches induce a slowdown of the mean circulation and a damping of the mesoscale activity, in contrast to a simulation forced with relative wind and without parameterization, they also allow representation of a partial re-energization of the mean oceanic circulation and of the mesoscale activity caused by the wind response (b and c in Figure 11). The best parameterization is the one based on a predicted s (equation (10); b in Figure 11). Although all the parameterizations obtain relatively similar results in terms of mean and median EKE, the parameterization based on a predicted s has the most realistic results with respect to the fully coupled simulation in eddy-rich regions such as Western Boundary Currents (b in Figure 11). Such a parameterization allows a simulation to represent more physics and to reproduce the daily and the spatial variability of s , avoiding possible compensation of errors when using parameterizations based on prescribed coupling coefficients.
The parameterization based on a predicted s has the advantage of flexibility because it can be used on top of a prescribed or estimated (with absolute wind) surface stress, thus allowing an oceanic model to be driven using climatological forcing. If the CFB effects are absent from both the wind and the stress of an atmospheric product, the user should take this parameterization using the and coefficients found in this study. If the CFB effects are present (as seems likely in future global reanalyses, labeled e-i in Figure 11), the surface currents should be known in a coherent way and filtered out (equations (11) or (12)) as using the parameterization based on the predicted s (labeled f in Figure 11). In such a situation, one could use the and found in this study or estimate and using the surface currents, the surface stress, and the 10-m wind 10.1029/2019MS001715 from the atmospheric product. As shown in the supporting information, because of uncertainties related to the statistical procedure used to extract the mesoscale signal, the coefficients may be increased by ≈20%. This would have the effect of increasing the eddy killing by the coefficient but also the re-energization by the coefficient. However, even by using the coefficients obtained directly from the linear regression, the levels of EKE should be better represented than any other simulations based on the other parameterizations or, obviously, simulations without parameterizations of the CFB (e, g, and h in Figure 11).
Scatterometer products should not be treated as a reanalysis as scatterometer wind and stress already contain the CFB effects but in relation to an unknown current field. As a result, an uncoupled simulation forced by a scatterometer stress (or wind) product (j-n in Figure 11) may result in a realistic mean surface stress and give an accurate estimate of the mean oceanic circulation assuming that mean surface currents are mostly deterministic (with weak eddy-mean flow interactions) and are relatively well reproduced by a model (labels j and l in Figure 11). However, as in any simulation that uses prescribed surface stress without a parameterization of the CFB (or a wind without CFB), the eddies are not correctly damped because the estimated or the prescribed surface stress is not correlated with the simulated eddies. A mesoscale activity too intense may result in a too large an inverse cascade of energy, which might destabilize, for example, Western Boundary Currents (Renault, Marchesiello, et al., 2019). The use of a relative wind without a parameterization will lead, as for the case of a reanalysis, to a poor estimation of both the mean circulation and the mesoscale activity (m in Figure 11). An alternative is to use a scatterometer wind product with a bulk-formula stress and a parameterization of the wind response to avoid overestimation of the eddy damping (k and n in Figure 11). Because the 10-m wind product from a scatterometer represents the wind relative to the oceanic currents, using this wind with a bulk formula and a parameterization based on prescribed coupling coefficients (n in Figure 11) would, therefore, have a realistic eddy killing effect; however, the large-scale effect might be doubly accounted for, weakening the mean circulation. This has been recently shown to be important by Sun et al. (2019) for the North Equatorial Countercurrent. To overcome this issue, one could filter out the CFB imprint (as in f in 11) by using the parameterization based on a predicted s and observed surface currents, for example, from AVISO (k in Figure 11). This should work if and only if the surface currents are known and are observed in a coherent way with the scatterometer wind and stress. However, the AVISO does not have the same effective resolution as QuikSCAT products and suffers from uncertainties. As a result, using the surface current from AVISO to correct, for example, QuikSCAT wind or stress does not allow reconstruction of clean wind or stress product.
Another topic related to this study deals with how to create a data set for forcing ocean models. As done by Large and Yeager (2009) and Tsujino et al. (2018), the time-mean wind speed from NCEP/NCAR and JRA55-do reanalysis, respectively, is adjusted to the time-mean QuikSCAT 10-m neutral relative wind. By doing so, these reanalyses may have the same issue as a scatterometer, that is, a double counting effect of the CFB on the large-scale circulation. To overcome this issue, one could simply add some time-mean ocean surface current to a scatterometer mean wind to estimate time-mean absolute wind using the coupling coefficients. However, again, the time-mean ocean surface current as estimated, for example, from AVISO does not have the same effective resolution as QuikSCAT products and suffers from uncertainties. As a result, adding the time-mean ocean surface current from AVISO to the QuikSCAT mean does not allow reconstruction of a clean time-mean absolute wind. Novels approaches, such as machine learning, could allow to identify oceanic features that imprint the stress by using various data set (such as SST, altimeters) and, thus, to filter out their imprint on scatterometer stress. Future satellite missions (such as SKIM; Ardhuin et al., 2017 andChelton et al., 2018;Rodríguez et al., 2018;WaCM Bourassa et al., 2016) would likely allow measuring 10-m wind, surface stress, and surface currents in a consistent way and, thus, could be another option to construct a data set for forcing ocean models but also to estimate the and coefficients used to predict s ((k) in Figure 11). They would help to better understand and characterize the wind and surface stress responses to both TFB and CFB.
In summary, forced oceanic models should take into account the CFB, which can be done by using the simple parameterizations tested in this study. These parameterizations should also be tested over ideal eddies and regional configurations to assess their effect on, for example, the eddy statistics and Western Boundary Currents dynamics, and we intend to investigate this soon. The wind response to the CFB can also alter the surface heat and freshwater fluxes, although, to our knowledge, this alteration has not been quantified precisely yet. On the one hand, a parameterization based on s w not only allows us to compute the momentum flux response to the CFB but also heat and freshwater fluxes responses in a consistent way. On the other hand, the parameterization based on a predicted s allows a better representation of the mesoscale activity; however, it cannot represent the alteration of the surface heat and freshwater fluxes by the CFB. A possible solution to overcome this limitation is to combine the parameterization based on a predicted s to estimate the surface stress with a parameterization based on s w to take into account the surface heat and freshwater fluxes responses to the CFB. However, this approach does not guarantee that the same winds are used to compute the surface stress, turbulent heat fluxes, and evaporation as was done by Large and Yeager (2009). It therefore comes at the expense of decoupling turbulent heat and freshwater fluxes from the momentum flux. One question that arises from this study is then as follows: Should we favor a more realistic representation of mesoscale activity with an ad hoc approach (i.e., the parameterization based on s ) at the expense of a known relationship among the turbulent fluxes? It may also be noted that even when the wind speed is adjusted in a common manner, turbulent heat fluxes will differ greatly between model simulations, as heat fluxes depend strongly on SST, which is simulation dependent, via the heat-flux feedback (Large & Yeager, 2012). Moreover, in this study, the mesoscale TFB is not parameterized as we show that it does not have a strong impact on the mean and mesoscale circulations.