Global implications of surface current modulation of the wind-wave field

The influence of ocean surface currents on the global wind-wave field is revisited. State-of-the-art numerical spectral wave model simulations with and without surface currents taken from an eddy resolving global ocean reanalysis were compared. As a global average, simulations forced with currents display significantly better agreement with altimeter derived wave heights. The bias and root mean square error in significant wave heights are mostly reduced when including current forcing, especially in the Southern Ocean. An overall improvement in wave periods and wave direction is also seen when comparing model outputs with the Australian and United States buoy network observations. Including surface ocean current forcing in wave simulations reduces the simulated wave heights in most areas of the world, due to a decreased relative wind given by co-flowing winds and currents. Current-induced refraction generates important changes in wave direction in western boundary current and tropical regions. Furthermore, large and broad changes in friction velocity, atmosphere-to-ocean energy flux, whitecap cover and Stokes drift velocities are observed in equatorial regions. Finally, the importance of the wave model resolution for representing wave–current interactions was tested by comparing results from eddy-permitting (lower resolution) and eddy-resolving (higher resolution) configurations. We conclude that the main patterns of current-induced refraction are well represented in both cases, albeit that the higher resolution simulation represents these in a more detailed manner. Finally, the implications that the observed wave–current interactions have on several ocean processes are discussed.


Introduction
Since the pioneering studies of Longuet-Higgins and Stewart (1960, 1961, 1964, it has been understood that there are complex non-linear interactions and energy transfers between ocean wind-waves and the mean flow. These are two-way interactions, i.e., wind-waves can influence surface currents and in turn currents can modify wave properties. In this study, we focus on the influence of ocean surface currents on wind-waves across the globe. Indeed, there are various ways in which currents can affect waves: waves can change their properties, such as their steepness or frequency, due to exchanges of energy with the mean flow via the radiation stress; besides, non-uniform currents can refract wave trains in a manner analogous to the way in which changes in depth induce wave refraction. Further, the transfer of energy from the atmosphere to the ocean can be modified in the presence of currents (the effective wind that acts on waves propagating on a moving medium might be different than the real wind). While these effects have been investigated generally at local scales, such as the effects of tidal circulation on coastal wave fields (e.g., Jones, 2000;Ardhuin et al., 2012;Rapizo et al., 2015;Lewis et al., 2019), and the understanding of more theoretical aspects of this interaction is in constant progress of San Clemente Island, Munk et al. (1963) measured the arriving swell waves and, through backtracking techniques and assuming swell waves follow great circle paths, inferred the position of the storm that had generated them. They concluded that intense storms in the Southern Ocean produced waves that travelled across the Pacific Ocean in 5 to 15 days before reaching the coasts of California. However, on occasions, the method of Munk et al. (1963) yielded errors of up to 1000 km in the position of the originating storms, placing them in the Antarctic Continent. Kenyon (1971) studied the wave refraction produced by the shear in ocean currents and, using an idealized model of the Antarctic Circumpolar Current (ACC), found it to be plausible that the ACC could deflect the direction of swell waves and thus explain the errors in Munk's measurements. Gallet and Young (2014) took up this idea by using modern satellite-derived ocean currents products to demonstrate that mesoscale vorticity in the ocean can significantly deflect the direction of waves, and thereby providing a compelling explanation for Munk's errors.
Many studies, using different methodologies and across a range of spatial scales, have shown the indubitable sensitivity of waves to ocean currents. Using SAR data from the SIR-B mission, Irvine and Tilley (1988) gave initial accounts on how waves can be refracted, and even trapped, by the Agulhas Current and how this could translate into extreme wave generation. Ardhuin et al. (2017), carrying out wave simulations with and without current forcing, showed that wave height variability in the open ocean is modulated by currents at scales of 10-100 km, and that the spectrum of spatial variations of currents is proportional to that of wave heights. Wandres et al. (2017), using a coupled ocean circulation and wave model, showed that the Leeuwin Current, especially its eddies and meanders, has a significant impact on the wave climate of Western Australia, being able to produce a change of ±25% in significant wave height and ±20 • in mean direction. More recently, Rapizo et al. (2018) showed that a significant improvement in the performance of global wave simulations can be achieved by including surface current forcing from a global reanalysis. They showed that the main effect of including currents is the reduction of wave heights due to a diminished relative wind in the presence of co-flowing current fields (this is especially evident in the Southern Ocean where the westerly winds blow in the same direction as the ACC). However, there are also current-induced refraction and energy convergence processes that occur in localized areas, such as the Agulhas Current or the Equatorial Current and Counter-Current systems. When comparing model results with observations taken from a moored buoy in the Southern Ocean, Rapizo et al. (2018) observed that the current-forced wave simulation performed better than the simulation without currents in terms of wave heights and period. In terms of wave direction, they found a small improvement for peak direction and a slight worsening for mean direction. They attributed this worsening to the coarse resolution of their simulation and the current forcing product (0.5 • ), presumably incapable of effectively capturing the ocean circulation features that modify wave direction through current-induced refraction. Likewise, Quilfen et al. (2018) found that their wave model simulations underestimated the magnitude of the wave height variability, which they also attributed to their model resolution (0.25 • ) being too coarse, as well as to errors in the ocean current data. Overall, these recent studies suggest that increased wave model resolution is an important factor to consider in order to properly represent wave-current interactions.
This paper revisits and challenges the current view on the importance of model resolution for wave-current interaction research. It builds on previous studies extending and improving some of their aspects. Specifically, we analyse a three-year high resolution (0.1 • ) wave model simulation across the Southern Ocean, with and without current forcing, to assess the importance of using an eddy-resolving resolution to properly represent wave-current interaction. We investigate the role of wave-current interaction on other processes such as energy and momentum fluxes between the atmosphere and the ocean, Stokes drift velocities (Stokes, 1847, of critical importance for search and rescue activities and tracking of particles in the ocean), whitecap coverage, among others. The paper is organized as follows: Section 2 describes the configuration of the wave and current models and the altimeter and buoy datasets used to validate the wave simulations; Section 3.1 presents a comprehensive validation of wave parameters from the simulations with and without currents against altimeter and moored buoys observations; Section 3.2 presents an analysis of the effects that ocean currents have on different wave parameters; and Section 3.3 discusses the importance of wave model resolution by comparing results from two sets of wave simulations, with high and low resolutions. Finally, Section 4 discusses key aspects of the main results, together with the main conclusions extracted from this study.

Wave model configuration
A multi-grid implementation of WAVEWATCH III v5.16 was utilized (Tolman, 2009), consisting of two grids: a global grid, from 78 • S to 78 • N with 0.4 • spatial resolution (i.e., at the margins of being eddypermitting); and a regional sub-grid spanning the Southern Ocean around the globe from 65 • S-28 • S, with a 0.1 • spatial resolution (eddyresolving). The atmospheric forcing for the wave model was the 10 m (above the mean sea level) wind data from the Climate Forecast System Reanalysis (CFSRv2, Saha et al., 2010), which has a spatial resolution of approximately 0.2 • in the horizontal and a temporal resolution of 1 h. The model considers the temporally varying sea ice concentration (also taken from CFSR), following the approach described by Tolman (2003). Default threshold concentrations were considered (< 0.25: ice has no effect on wave propagation; > 0.75: ice is treated as land; > 0.25 and < 0.75: dampening of the wave energy depending on the ice concentration value). The obstruction masks used to account for small islands not resolved by the model grid (Tolman, 2003) were constructed using the Global Self-consistent Hierarchical High-resolution Shoreline (GSHHS) dataset. Water level effects were not incorporated in the wave simulations in this study. Some of the most relevant features of the model's configuration are: the Ardhuin et al. (2010) source term parametrizations (ST4), specifically configured for utilization with the CFSR winds; the Discrete Interaction Approximation (DIA, Hasselmann et al., 1985) for nonlinear wave-wave interactions; the third-order Ultimate Quickest propagation scheme (Leonard, 1979(Leonard, , 1991, including the garden sprinkler effect correction (Tolman, 2002); JONSWAP bottom friction; and Battjes and Janssen (1978) shallow water depth breaking. The model is spectrally discretized in 29 frequencies, varying from 0.035 to 0.5 Hz with an increment factor of 1.1, and directions taken every 15 • . With this configuration, we run two sets of simulations: one including ocean surface current data taken from the Bluelink Reanalysis (BRAN, Oke et al., 2013) as an additional forcing for the wave model, and another simulation without any surface current forcing. BRAN is an eddy-resolving ocean reanalysis (more specifically, BRAN is 0.1 • x 0.1 • in the horizontal) that uses data assimilation to produce a realistic quantitative description of the three-dimensional ocean circulation for the last three decades. We have used the surface current product from BRAN (at 2.5 m depth). By comparing the differences in these two sets of simulations (with and without current forcing) we assess the sensitivity of wind-wave modelling to the inclusion of surface currents in both eddy-permitting and eddy-resolving model grid configurations. We have run three years of wave simulations from January 2014 until December 2016 (December 2013 was run but considered as a spin-up month). As mentioned earlier, one of the main effects of the inclusion of surface current forcing in wave simulations is a modification of the relative wind. WAVEWATCH III incorporates this process in a rather simple fashion, balancing wind and current vectors. The relative wind vector (from which wind speed and direction are determined) in the presence of currents is computed as: where ⃖⃖⃖⃖ ⃗ represents the absolute wind vector (with units of m/s), ⃖⃖ ⃗ is the ocean current vector (in m/s) and is a dimensionless coefficient. WAVEWATCH III considers = 1, which implies the largest effect of currents (e.g., if the measured 10 m absolute wind speed and direction are the same as the surface currents, the relative wind input to the model is null).

Altimeter wave dataset
Altimetric satellites orbit the Earth measuring significant wave heights with an almost global coverage. Altimeter observations represent an outstanding data source to assess a global wave model's performance (Rascle and Ardhuin, 2013;Durrant et al., 2014;Roland and Ardhuin, 2014;Hemer et al., 2017). Here, we have used altimeter derived significant wave heights to evaluate the performance of the current-forced and non-current-forced wave simulations. The altimeter dataset was calibrated, quality-controlled and assessed by Ribal and Young (2019). It comprises all the available altimeter missions from 1985 to the present. We have selected the along-track significant wave height data of those missions that operated during the period 2014-2016: JASON-2, CRYOSAT-2, HY-2 A, SARAL, JASON-3 and SENTINEL-3 A. We have mapped these data onto the global 0.4 • resolution grid of our model. That is, for each bin of our model grid, we selected the satellite data that fell into that bin and averaged it over the course of an hour. In this way, we compared the significant wave height at each grid point of the model with all the available satellite observations that corresponded to that grid point. With this approach, the comparisons were localized, and excessive averaging avoided. However, the availability of altimeter data is inhomogeneous, i.e., it is higher close to the satellite tracks but there are areas in which the amount of data points is much lower.

Wave buoy data
Wave buoy measurements are usually regarded as the most accurate wave information available, and they are often used to validate altimeter wave height products (Zieger et al., 2009;Ribal and Young, 2019). For this study, we have selected a representative set of buoys from Australia and the United States (US), publicly available through the Australian Ocean Data Network (AODN: https://portal.aodn.org. au/) and US National Data Buoy Center (NDBC: https://www.ndbc. noaa.gov/) portals, respectively. All the buoys have significant wave height information. However, some of them compute it from a waveto-wave analysis and others from spectral analysis. Most of the buoys also measure wave direction (peak direction for some of them, and peak direction for sea and swell separately for others). In every case, the WAVEWATCH III parameters were computed via spectral analysis, and the necessary processing/filtering (such as removing bad data or averaging buoy observations with a higher sampling frequency than the model output) was applied to match the observations with the model's results. The wave parameters selected for comparison (significant wave height, ; peak period, ; and peak direction, ) were computed from modelled WAVEWATCH III directional wave spectra as follows: with ( , ) being the wave spectral density as a function of frequency and direction (2D spectrum), the peak frequency computed from the 1D spectrum (frequency at which the highest energy level is attained), = ∫ 2 0 cos ( ) ( , ) and = ∫ 2 0 sin ( ) ( , ) . It is important to note that WAVEWATCH III outputs the spectrum in terms of intrinsic (and not absolute) frequency. Therefore, for a correct comparison against buoy measurements, a Doppler shift correction should be applied to the output of the current-forced simulation to compute the absolute frequency. The conversion from intrinsic to absolute frequency was done using the dispersion relationship of wind-waves, With and the absolute and intrinsic (angular) frequency, respectively, ⃖ ⃗ the wavenumber vector and ⃖⃖ ⃗ the surface current vector. The inner product in Eq. (5) was approximated using the mean wavelength and mean wave direction.

Comparison against observations
In the first instance, the significant wave heights ( ) simulated by the model were evaluated against those derived from the satellite altimeter observations. The following statistical parameters were used to evaluate the model runs as indicators of the degree of similarity between observed ( ) and modelled ( ) values: the Pearson linear correlation coefficient R, the mean bias ( 1 ∑ ( − )), and the root mean square error RMSE . The normalized bias was computed as 100 * . For wave direction, the linear correlation coefficient formula was taken from Jammalamadaka and Sengupta (2001), and the bias and RMSE were computed keeping in mind that wave direction is an angular variable. Treating the direction values as points in a circumference of unit radius, the difference between two given direction values is defined not as the arithmetic subtraction between the two (as we would do for wave height or period), but as the smaller arclength along the circumference (this way, the difference between 350 • and 10 • is +20 • , not +340 • ; conversely, the difference between 10 • and 350 • is −20 • ). Likewise, the mean direction was computed treating the direction values as unit vectors and calculating the direction of their resultant vector (i.e., dividing the direction values dataset in x and y components, computing the average of each component separately and forming the resulting vector with them); thus, the mean direction between 10 • and 350 • would be 0 • , and not 180 • . When comparing with buoy measurements, the bias in wave direction is defined as the average of all the differences in direction values (model minus observations). Therefore, a positive bias in direction shows that the model represents waves in a more clockwise direction relative to north than the observations. Fig. 1 shows the Pearson's linear correlation coefficient (left panels), mean bias (middle panels) and RMSE (right panels) of for the simulation without currents (a, d, g), for the current-forced simulation (b, e, h), and the difference between those statistics (c, f, i). The correlation coefficient between altimeter derived and modelled significant wave heights is > 0.9 at high latitudes of both hemispheres and it is lower (∼0.7) at low latitudes and in coastal areas. The bias and RMSE distributions observed in Fig. 1 are a common feature of global wave models. The mean bias is mostly positive, meaning that the model overestimates the observed wave heights, except in the western equatorial Pacific where the model simulates lower wave heights than the observed ones. The highest absolute errors are found at high latitudes, especially in the Southern Ocean. These errors are related to several factors, such as a positive bias in the reanalysis winds used to force the model, the fact that the wave model does not incorporate information about drifting icebergs which block the wave energy, among others. There are also high errors close to the ice edge, presumably related to the parameterization of wave propagation over areas partially covered with ice being too simple. By visual comparison of the performance of both simulations (with and without currents), it can be seen that incorporating current forcing in the wave model significantly increases the correlation at low latitudes, and reduces the mean bias and RMSE, especially in the Southern Ocean. The bottom panels show the difference in correlation, bias and RMSE, computed as the value of the current forced simulation minus that of the simulation without currents. In this way, a positive value in Fig. 1.c means that the correlation coefficient in the current-forced simulation is higher than in the simulation without currents. Indeed, there are substantial improvements in correlation (i.e., higher values in the current forced simulation) at low latitudes, especially at the eastern side of ocean basins, in the south-eastern African coasts and, to a lesser extent, in the Antarctic Circumpolar Current and East Australian Current regions. Regarding the bias, because it is a signed quantity, when comparing bias values from different simulations, we decided to take the difference in their absolute values (e.g., a bias of −0.2 m is regarded more accurate than a bias of +0.3 m). Wherever the difference in bias or RMSE is negative (blue in Fig. 1.f and/or Fig. 1.i), it means that the absolute bias or RMSE of the current-forced simulation is lower than that of the simulation in the absence of currents. In general, most areas of the world show a reduction of the mean bias or the RMSE, and hence an improvement of the wave model performance, when including current forcing. However, there are also some areas where the mean absolute bias or the RMSE of the current-forced simulation is higher than that of the simulation without currents. These areas are mainly the western Pacific and the waters surrounding the Antarctic shelf. As a global average, including current forcing in the model reduces the bias from 0.21 m to 0.17 m, reduces the RMSE from 0.29 m to 0.27 m, and increases the correlation between modelled and observed wave heights from 0.906 to 0.911.
On the other hand, similar results are obtained comparing the observed and simulated extreme wave heights: Figure S1 of the Supporting Information section (SI) shows the observed 90th percentile (top panel) and the differences in simulated 90th percentile from the wave simulations with and without currents (bottom panel). As expected, the extreme wave heights are reduced in most areas of the world, particularly in the Southern Ocean, due to co-flowing winds and currents lowering the energy transfer from the atmosphere. However, including current forcing increases the 90th percentile in the Agulhas Current, East Australian Current, Equatorial Counter Currents and Gulf Stream regions, likely due to exchanges of energy between the waves and the main flow in conditions of opposing waves and currents. The same comparison against altimeter data was carried out for the 90th percentile ( Figure S2). The correlation values of 90th percentile wave heights are significantly lower than for the total . Further, although there is an overall improvement in the correlation values of the 90th percentile wave heights, some areas reveal a decrease in correlation for the current-forced simulation ( Figure S2.c). In terms of bias and RMSE, there is a substantial improvement in the current-forced simulation, especially at high latitudes, as well as in the Agulhas Current, East Australian Current and Equatorial Counter Current regions. However, the mean absolute bias and RMSE of the 90th percentile increase (i.e. become worse) with the inclusion of current forcing in the Tropics (except in the Equatorial Counter Current areas), particularly in the Western Pacific basin. Fig. 2 shows the performance of both wave model simulations in terms of R, bias and RMSE throughout the year and separated by regions. The bias and RMSE are expressed here as percentages of the mean observed wave heights. The correlation values for the currentforced simulation are consistently higher than those of the simulation without currents for each month of the year. Likewise, the mean bias and RMSE are consistently lower throughout the seasons. The wave simulation in the absence of currents performs reasonably well globally, with total R-values > 0.85, a mean bias < 10% of the average observed wave heights, and RMSE < 20%. The current-forced simulation always performs better, with correlation values approximately 0.005-0.01 larger, mean bias values ∼1.5% smaller and RMSE ∼1% smaller. The performance of the current-forced simulation is also improved in each ocean basin, with the improvement being slightly greater for the Southern Hemisphere basins. While these differences are small, it should be considered that these statistics are computed over long periods of time and over large areas, while the effects of wavecurrent interaction can be more localized, as observed in Fig. 1 (e.g., if only a portion of the Southern Ocean or of the Agulhas Current is selected, the differences will be larger). It is also worth noticing that the improvement in estimates is seen at each ocean basin throughout the year ( Figure S3).
Next, a comparison against buoy data is presented, with the aim of evaluating the performances of both wave simulations (with and without currents) in representing the observed wave heights, periods and directions close to the coast. First, we show the evaluation of the WAVEWATCH III runs against the Australian buoys. A total of 15 Australian buoys was selected, three in Western Australia (WA), one in South Australia (SA), one in Tasmania (TAS), eight in New South Wales (NSW), two in Queensland (QLD), and finally one in deep waters of the  (middle) and RMSE (bottom) of the simulations without currents (red line) and with currents (blue line), for each month of the year (left panels) and separated by regions (right panels). The bias and RMSE are expressed as percentages of the mean observed wave heights by the altimeters. All the correlation values for the current-forced simulation are statistically higher (at the 95% confidence interval) than those of the simulation without currents. . The green dots show the differences in correlation, computed as the value in the simulation with currents minus that of the simulation without currents. The red squares and blue triangles show the differences in mean bias and RMSE respectively, computed as the value in the current-forced simulation minus that of the simulation without currents. Units are in metres for bias and RMSE differences. (c) Same as (b) but for peak direction. Units are in degrees for bias and RMSE differences. (d) Same as (b) but for peak period. Units are in seconds for bias and RMSE differences. Southern Ocean. Table S1 (see SI) shows the correlation, bias and RMSE values for , , and between the Australian buoy measurements and WAVEWATCH III output (for both simulations, with and without currents). Good agreement was found in terms of , with R∼0.9, mean biases ranging between 2 cm and 30 cm, and RMSE values generally between 30 cm and 40 cm. The has correlation values ∼0.5, a mean bias ∼1.5s for the WA, SA and TAS buoys but lower than 0.8s for the NSW and QLD buoys, and RMSE∼2.4s. In terms of , we found R∼0.6, the mean bias ranges from −16 to 26 • across locations, and RMSE∼40 • . A similar performance of the model is found when comparing against NDBC buoys (Table S2).
Here, only the differences in the performances of both simulations are shown to highlight the potential improvement achieved by considering current forcing in the wave simulations. Panel (a) of Fig. 3 shows the Australian buoys that were selected for this comparison, whereas panels (b), (c) and (d) show the differences in correlation, . The green dots show the differences in correlation, computed as the value in the simulation with currents minus that of the simulation without currents. The red squares and blue triangles show the differences in mean bias and RMSE respectively, computed as the value in the current-forced simulation minus that of the simulation without currents. Units are in metres for bias and RMSE differences. (c) Same as (b) but for peak direction. Units are in degrees for bias and RMSE differences. (d) Same as (b) but for peak period. Units are in seconds for bias and RMSE differences. absolute value of the bias and RMSE in , and respectively. As in Fig. 1, a positive difference in correlation, or a negative difference in absolute bias or RMSE, represents an improvement in the wave model performance when including currents. As before, the difference in the absolute value of bias is presented in Fig. 3. In terms of , there is a statistically significant (at the 95% confidence level) improvement in correlation in all the NSW buoys (except Tweed Heads) and in Brisbane. The RMSE improves for all locations except for Cape du Couedic, Tweed Heads and Gold Coast. The absolute bias (average of the differences) improves only for the WA buoys, Cape Sorell, Crowdy Head, Brisbane and the Southern Ocean Flux Station (SOFS) buoy. Nevertheless, the average wave height for the whole period is better represented by the current-forced simulation in most locations (not shown). In terms of peak period, there is a general statistically significant (at the 95% confidence level) improvement in the WA, NSW and QLD buoys (except Tweed Heads), and a broad reduction of the bias and RMSE across all sites. Finally, in terms of peak direction, there is a significant improvement for Rottnest Island, Albany, all the NSW buoys (except Eden, Byron Bay and Tweed Heads) and Brisbane, whereas in Eden and Gold Coast there is a statistically significant worsening. The RMSE of the current forced simulation is substantially lower for the NSW buoys. The mean bias improves for 8 of the 14 stations and worsens for Rottnest Island, Cape Naturaliste, Crowdy Head, Coffs Harbour, Tweed Heads and Gold Coast. Fig. 4 presents a comparison analogous to that in Fig. 3 but for the NDBC buoys. In terms of , there is an overall increase in the correlation values when adding currents (although not statistically significant at the 95% confidence level), and a substantial decrease in the absolute value of the bias (except for buoys 46078 and 46025) and RMSE. In terms of peak period, there is a statistically significant improvement in the correlation values (the difference in correlation is not statistically significant for buoys 46073, 46012, 42055, 42002 and 42057), and a substantial decrease in the absolute bias and RMSE (except for buoys 46073 and 42055). The differences in the performance statistics of estimates are more variable: while there is an overall increase in the correlation values (although statistically significant only for 6 buoys), some locations show a decrease in correlation (not significant); the absolute bias in peak direction decreases at 8 buoy locations and increases at 10 locations; whereas the RMSE shows a general decrease in the current-forced simulation.
Although differences in bulk statistics between simulations appear small, the inclusion of currents as a forcing of the wave model produces major changes during specific wave events (e.g., like the one shown in Fig. 5). During the last week of March 2016 in Sydney, the measured significant wave heights were consistently > 1.5 m, and > 3.5 m in the last day of the month. The current-forced wave simulation agrees much better with the observations and has substantial differences with the simulation without currents of up to 1.2 m in wave heights, 5s in period and 80 • in direction during that event. During that day, the meandering southward flow of the East Australian Current generated two very intense anticyclonic eddies centred around (152 • E; 36 • S) and around (152 • E; 39 • S), which had a significant impact on , the mean periods and direction estimates in the Tasman Sea, with the maximum differences being near Sydney. The directional wave spectrum at the Sydney buoy location for 2016-03-27T12:00Z of the simulation without currents shows a northward propagating wave mode with a frequency of around 0.07 Hz (period of ∼14 s) and a westward propagating mode with a frequency of around 0.09 Hz (period of ∼11 s). For the simulation with currents, the energy in the northward propagating mode is greatly increased, while the westward propagating mode decreases its intensity. The one-dimensional spectra reveal that including current forcing shifts the peak in the spectrum towards lower frequencies, due to the increased energy in the northward propagating mode. The spatial patterns of differences in wave height, period and direction show a positive difference in areas where the currents and the northward propagating wave mode travel in opposite directions, and negative differences when they travel in the same direction.
While undoubtedly the ocean circulation at the coast (e.g., tidal currents, not accounted for in the Bluelink reanalysis) will have a significant impact on the measured wave properties nearer the coast, the results presented in Fig. 5 allow us to conclude that the larger scale ocean circulation (as seen by the Bluelink reanalysis and presumably other ocean circulation reanalysis products) can also have a significant influence on the wave climate at the coast. Therefore, current forcing should be an important consideration for coastal wave modelling studies in East Australia.

Changes in wave properties due to currents
In this section, a comparison between the output from simulations with and without currents is presented. Differences in several wave parameters were analysed to further understand the effects of ocean surface currents on the global wind-wave climate. Results here should be interpreted in light of those presented also in Section 3.1, that is, any conclusion drawn in this section is valid in those areas where the inclusion of currents as a forcing in the wave model translates into an improvement of the model's performance (these include most areas of the world with the exception of the western Pacific and waters surrounding the Antarctic continent and the Arctic). For brevity and convenience, we use the following acronyms for the major ocean currents: Antarctic Circumpolar Current=ACC, Agulhas Cur-rent=AC, Madagascar Current=MC, East African Coast Current=EACC, East Australian Current=EAC, Indonesian Throughflow=ITF, Equatorial Currents=EC (i.e. the North Equatorial Current and South Equatorial Current), Equatorial Counter-Current=ECC, Brazil Current=BC, North Brazil Current=NBC, Caribbean Current=CC, Gulf Stream=GS and Kuroshio Current=KC. Fig. 6 outlines the regions associated with these ocean currents together with the spatial structure in their average speeds from 2014-2016, computed with data from the BRAN reanalysis. In addition, Fig. 6 presents the time-averaged relative vorticity for that period for the major ocean current systems described in the paper. Fig. 7 shows the time-averaged percental changes in significant wave height ( ), mean period ( 01 ) and changes in mean direction ( ) for the three-year period in which the simulations were run (2014)(2015)(2016). In each case, the difference in wave parameters was computed as the value of the variable in the current-forced simulation minus that of the simulation without currents. The inclusion of current forcing in the model leads to an overall decrease in the significant wave height in most areas of the world. This is mostly due to a diminished relative wind in areas with co-flowing winds and currents. The highest absolute changes in are observed in the Southern Ocean (not shown), although they represent less than 8% of the average wave heights in that region ( Fig. 7.a). Here, the ACC flowing in the same general direction as the predominant westerlies makes the relative wind lower than the real wind (Rapizo et al., 2018). In the AC region, east of the African coast, there is a dipole of increasing wave heights close to the coast and decreasing heights further offshore. A similar feature is observed in the CC region, and it is also present in the MC and EAC, although with less intensity. This pattern of increase and decrease in could be related to both exchanges of energy between the waves and the currents and/or to convergence and divergence of wave energy due to current-induced refraction (however we do not identify the relative contribution of each processes to the observed differences in wave parameters in this study). In the Agulhas retroflection (a prominent turnabout of the southern end of the Agulhas Current system, where the current direction reverses and flows back into the Indian Ocean), a pattern of alternating positive and negative differences in wave height arises, which is a consequence of the intense mesoscale eddy activity that characterizes this area. Eddies can refract the incoming wave trains in very complex manners (since the change in wave direction is determined by the sign of the current's shear, it is different on each side of the eddy), hence redistributing the wave energy and generating energy convergence and divergence areas, which in turn translate into higher and lower wave heights, respectively (e.g., Mathiesen, 1987). The reduced relative wind in the Southern Ocean generates waves of lower amplitude. These waves will propagate away from their generation area, and since they have lower heights in the Southern Ocean, they will also present reduced wave heights in the Southern Hemisphere ocean basins. That is why, on average, the wave heights are reduced in most parts of the ocean basins when incorporating currents into the model. Important differences are observed in the equatorial region, with an average decrease in wave heights in the North and South EC of around 8% of the mean and a slight increase in wave heights in the ECC region. In the Northern Hemisphere, little or no change is observed in the KC and the GS. However, there is an observed decrease in in the North Pacific and North Atlantic Oceans centred around 50 • N, also related to a lower relative wind. The most conspicuous percental changes in wave height are observed in the NBC, CC, ITF and at the eastern African coasts. The NBC and CC Fig. 7. (a) Percental average difference in significant wave height, computed as the value in the current-forced simulation minus that of the simulation without currents, normalized by the mean wave height in the simulation without currents . Units are in %. (b) Same as (a) but for mean period. Units are in %(c) Same as (a) but for mean direction. Units are in degrees. Most grid points present a statistically significant difference (at the 95% confidence level, based on a t-student test), except those locations where the percental difference is approximately between ±1% for and 02 . For wave direction, it is considered that the differences are insignificant when the mean direction in the current-forced simulation lies between the 95% confidence interval of the mean direction of the simulation without currents. Significant differences are attained in those locations with a mean difference of more than ±1 • . significantly reduce by up to 15% of the mean values, and in the NBC retroflection area there is a slight increase in . The circulation in the Gulf of Mexico increases the wave heights at the eastern side of the Gulf. In addition, the Florida Current increases by more than 15% of their mean values. The ITF is an ocean current system located in an area with extremely complicated geography. It is composed of a complex suite of currents that bifurcate and converge, but collectively describe a net east to west transport from the Pacific into the Indian Ocean. As such, there are alternating increases and decreases of wave heights due to the inclusion of currents. The most important changes occur north of Sulawesi Island, where is reduced by more than 15% of the mean in this area. The southernmost part of the ITF affects the wave climate off northwest Australia, increasing the wave heights close to the coast and decreasing them further offshore. The South Indian Ocean EC flowing northward of Madagascar significantly reduces by more than 15% of the mean values in this region. Close to the coast, the flow is split into the southward flowing AC (which increases the wave heights close to the coast) and the EACC (which decreases the wave heights). The changes in wave height due to the inclusion of current forcing in the wave simulation can explain the improvements in the simulated observed in Fig. 1.f. For example, the mean bias in for the simulation without currents in the western Pacific is mainly negative (Fig. 1.d). Therefore, if the inclusion of currents increases the wave heights in this area, this will translate into an improvement of estimates. Indeed, Fig. 7.a shows that the wave heights are increased in a narrow band in the western equatorial Pacific, where a reduction of the mean bias and RMSE is observed (Fig. 1.f and i). Similarly, in the NBC region close to the coast, the mean bias in for the simulation without currents is negative (Fig. 1.d). Including currents reduces the wave height in this area (Fig. 7.a), therefore exacerbating the negative bias, which is deemed as a worsening of estimates in this area (Fig. 1.f). Nevertheless, the correspondence between Fig. 1.f and Fig. 7.a is not absolute (possibly due to altimeter under-sampling in the equatorial region).
To better understand the conditions that lead to an increase or decrease in wave height (and mean period), Fig. 8 shows the timeaveraged changes in and mean period ( 01 ), as a function of the current speed and the wave propagation relative to the currents' direction. A value of 0 • in the x-axis of Fig. 8 indicates that waves propagate in the same direction as the currents, and a value of −180 • or 180 • indicates that waves and currents are propagating in opposite directions. Fig. 8.a and b show the differences in significant wave height for the December-January-February and June-July-August months, respectively. When wave propagation is in the same direction as the currents, there is an overall decrease in wave height. Conversely, when waves propagate against the currents, there is an overall increase in wave heights. Something analogous happens with the mean wave period: there is a general decrease (increase) in wave period when waves propagate with (against) the currents. These results are well in agreement with those presented by Barnes and Rautenbach (2020) for the South African region. However, it remains unclear what are the physical mechanisms that lead to these differences: while a reduced relative wind (for winds and currents flowing in the same direction) should decrease the wave period, the change in intrinsic wavenumber magnitude (the ''concertina'' effect) should increase the wave period for waves going in the same direction as the currents. In addition, currentinduced refraction can significantly alter the wave climate of regions otherwise unaffected. Fig. 7.b shows the percental differences in 01 between the simulations with and without currents. There is an average decrease of mean wave periods in the Southern Ocean, which could also be related to a less intense relative wind (e.g., a JONSWAP spectra will become narrower in frequency and the peak will move to lower frequencies with increasing wind speeds Hasselmann et al., 1973). Further, these differences are propagated throughout the ocean basins, producing lower-period waves in the South Indian, Pacific and Atlantic Oceans. Similar to the differences, there is a dipole of increasing wave periods close to the coast and decreasing periods off-shore in the AC, MC, CC and EAC regions. Fig. 8.c and d shows that, on average, waves tend to decrease (increase) their periods when they propagate in the same direction as (opposite to) the currents. This would be particularly true in the AC, MC and EAC regions which receive southerly waves from the Southern Ocean. Furthermore, this is supported by Figure S4, which shows the average differences in mean wave period off-shore of eastern Australia for the June-July-August months, together with the one-dimensional frequency spectra of both wave simulations for two distinct locations: one where the waves propagate in the same direction as the current (where there is a decrease in the wave energy in the spectrum and the peak shifts towards higher frequencies), and another one where waves oppose the flow (where there is an increase of energy in the spectrum with a shift towards lower frequencies). These differences in wave spectral density levels could be produced by exchanges of energy with the mean flow (mediated by the radiation stress, Longuet-Higgins and Stewart, 1964), or by current-induced refraction, affecting areas that, in the absence of currents, would present a different wave climate. The substantial increase in wave period for current speeds of 1-1.75 m/s for waves propagating approximately between −50 • and −120 • of the currents' direction observed in Fig. 8 can be traced back to the interaction of Southern Ocean swell waves with the EC (see Figure  S6). In the equatorial regions, there is a time-average increase in the mean wave period. Here, the wave climate is highly complex, with multiple wave modes presenting different frequencies and directions, and being affected differently by the ocean circulation ( Figures S5 and  S6). For example, the inclusion of currents increases the wave period of the Southern Ocean swell at location (120 • W; 0 • ) mostly in the June-July-August months, whereas it reduces the period of the south-easterly equatorial waves during this time. On the other hand, the Northern Hemisphere generated waves increase their periods due to currents throughout the year (Figures S5 and S6). In this case, the swell waves propagate opposite to the currents' direction, and the south-easterly waves travel in the same direction as the local currents. In the North Indian Ocean, there is also a decrease in wave periods that can be traced back to the lower period waves coming from the Southern Ocean. However, at the western side of the Arabian Sea and the Bay of Bengal, there is a significant increase in wave periods of around 7%. These changes seem to reach the coasts of Pakistan, India and Myanmar. In addition, a reduction of more than 10% of the mean period is observed to the north of Sulawesi Island. Off northwest Australia, the wave period is increased by around 8%.
Importantly, there are significant changes in mean wave direction due to the inclusion of currents in the model. Here, a positive difference in direction means that the waves of the current-forced simulation propagate in a more clockwise direction than in the simulation without currents (for example, the wave direction in the simulation without currents could be 90 • (to the E) and the direction in the case with currents 135 • (to the SSE)). Likewise, a negative difference implies an anti-clockwise rotation of waves with the inclusion of currents. Kenyon (1971) stated that waves which propagate against a variable current will be deflected in the direction of increasing current speed. Conversely, waves propagating in the same direction as the currents will be deflected in the direction of decreasing current speed. In the Southern Ocean, there are positive differences in wave direction, which means that currents refract the eastward propagating waves on average a few degrees to the south. There are very strong directional differences close to the southeast African and east Australian coasts, once more evidencing the dipole of positive differences (clockwise rotation) to the right of the current's flow direction and negative differences (anticlockwise rotation) to the left of the current. This pattern is also evident in the BC offshore of Uruguay and south Brazil, in the KC and, to a lesser extent, in the GS and NBC. Given the exposure of the south African coasts to Southern Ocean generated swell, the AC can induce a clockwise rotation of these waves if they are close to the coast, and an anti-clockwise rotation if they are further offshore. This could translate into a swell wave focusing in the central axis of the AC, as well as in the EAC. In deep waters at around 30-40 • latitude, on average there are insignificant changes in mean direction, hence current-induced refraction is less important in these areas. However, waves in the Tropics are significantly refracted by the EC and ECC. In the eastern equatorial Pacific, there are time-averaged changes in mean wave direction of around +10 • . As shown in Figures S5 and S6, the equatorial Pacific is an area with a very complex and multi-modal wave climate, and therefore examining changes in the direction of individual wave modes will be more accurate. Nevertheless, we can conclude that significant wave refraction occurs in this area. In the western Pacific, complex refraction patterns are observed, with clockwise and anticlockwise rotations induced by the currents. However, this is an area where the performance of the wave model worsens when including currents. Therefore, the changes in wave direction might not be a real feature but may be a consequence of a poorer representation of the ocean circulation in this area by the Bluelink Reanalysis. At high latitudes of the Northern Hemisphere, although small, the main changes in mean direction are negative, meaning that the currents refract waves to a more northward direction. The patterns of differences in mean wave direction approximately follow the patterns of relative vorticity in western boundary currents and the ACC. Negative (positive) vorticity values match with clockwise (anticlockwise) differences in wave direction. In the AC or EAC regions, the main currents flow poleward and the prevalent waves propagate in the opposite direction (equatorward). According to Kenyon (1971), waves will be deflected in the direction of increasing current speed, hence producing clockwise rotations in the western side of the AC or EAC, and anticlockwise rotations in their eastern sides (Fig. 7.c). Fig. 9 shows the time-averaged percental differences in wind friction velocity (a), wave energy flux (CgE, b), Stokes drift speed (c), and the changes in Stokes drift direction (d). The friction velocity is a parameter used to characterize the stress transmitted from the wind to the ocean surface, and it can be modified under different sea surface roughness states (e.g., waves). The parameterization of wind friction velocity in WAVEWATCH III depends on the source term physics choice. For ST4, it is an adaptation of Janssen (1991) that includes a ''sheltering term'' ( ), designed to reduce the drag coefficient at high winds. Including current forcing in the wave model reduces the friction velocity in most areas of the world, particularly the EC, ACC, EACC, NBC and CC, and increases it mainly in the ECC, Florida Current and NBC retroflection area (Fig. 9.a). The signs of the changes in in these  areas approximately match those of the friction velocity. For example, in the EC regions, a diminished relative wind due to co-flowing winds and currents reduces the surface stress, and consequently the friction velocity values are lower. As a result, the energy and momentum fluxes from the atmosphere to the ocean decreases ( Figure S7.b), and with it the generated wave heights (Fig. 7.a). The opposite occurs in the ECC regions, where the currents and the wind flow in opposite directions and therefore the relative wind computed by the wave model (and hence the wind friction velocity and the transfer of energy from the atmosphere to the ocean) is greater. In addition, the whitecap coverage is also reduced ( Figure S7.c), due to a decreased wind stress acting on this area (whitecaps arise from waves breaking in deep waters, entraining air bubbles and forming patches of foam in the surface of the ocean, Leckler et al., 2013;Scanlon et al., 2016). In short, the wind friction velocity, atmosphere-to-ocean energy and momentum flux, and whitecap cover changes due to the inclusion of currents, follow the same pattern. Their most salient features are a decrease in the North and South EC, an increase in the ECC and Florida Current, and a decrease in the central North Indian Ocean and in the EACC region.
The changes in friction velocity range between ±10% ( Fig. 9.a) in the Tropics, whereas the changes in whitecap cover and energy flux from the atmosphere vary between ±50% of their mean values ( Figure S7). We find that the wave energy flux (CgE, Fig. 9.b) decreases in most areas of the world, especially in the AC, EACC, NBC and CC regions. However, there is a slight increase in the ECC region and close to the southeast African and Australian coasts. This increase can be linked to the increases in the spectral energy levels (and therefore in significant wave height) and wave period observed in these areas (e.g., Fig. 7, S4, S5 and S6). The wave energy flux is an important variable for assessing the potential force of the waves on coastal processes and on coastal or offshore infrastructure. We find that ignoring current forcing in wave simulations can lead to inaccurate estimates of the wave energy flux in many coastal regions around the world (for example, > +25% mean difference in the coast of Mozambique, > +10% in the east coast of Australia, Madagascar and India, < −20% difference in the coasts of Somalia). The decrease in CgE in the southern coast of Australia and the increase in the eastern coast due to the inclusion of current forcing could help to improve the bias in CgE observed by Hemer et al. (2017).
Finally, differences in Stokes drift speed and direction are shown in panels 9.c and 9.d, respectively. As mentioned earlier, the Stokes drift is an important consideration for search and rescue activities and tracking of particles in the ocean, such as plastics, plankton, or drifting debris (Van Den Bremer and Breivik, 2018;Dobler et al., 2019). Recent theoretical developments have been made regarding the solution for Stokes drift velocities in the presence of ocean currents (Henry, 2019).
Here, we analyse the differences in the Stokes drift representation by our simulations (with and without current forcing). Interestingly, the time-averaged differences in Stokes drift speed present an analogous pattern to that of the friction velocity: decreasing values in the ACC (a maximum absolute average change of 0.25 cm/s), EACC, central region of the North Indian Ocean, EC, NBC and CC, and an increase in the ECC and Florida Current. The changes in Stokes drift speed range between ±30% in low latitudes. In terms of Stokes drift direction, there are significant differences in the AC, MC, ITF, EAC, NBC, KC, GS and especially in the equatorial region. As seen in the mean wave direction, in this case there is also a dipole of positive differences in Stokes drift direction (clockwise rotation) to the right of the mean flow direction, and negative differences (anti-clockwise rotation) to the left of the current's direction in the major western boundary currents (AC, MC, EACC, EAC, KC BC, BC, NBC, CC ad GS). The differences in Stokes drift direction can also be linked to the differences in the wind friction velocity direction ( Figure S7.a). However, these are Eulerian differences, in the sense that they are computed as the time-average at each grid point; a Lagrangian difference (following the particle as it moves) might have a significantly different impact on the tracking of particles in the ocean. Rapizo et al. (2018) showed that a significant improvement in estimates is attained by considering current forcing in wave simulations. However, when comparing with buoy observations (SOFS), they concluded that despite the improvements in estimates of wave height and period, only a slight improvement in peak wave direction estimates and a slight worsening in mean direction were found for the currentforced simulation. They attributed this to the current forcing product and model resolution being too coarse (0.5 • ), presumably incapable of representing current-induced refraction properly: the change in wave direction is proportional to the magnitude of the current shear in a direction perpendicular to the waves propagation (Dysthe, 2001), which would be better represented in an eddy-resolving simulation (rather than eddy-permitting). Rapizo et al. (2018), using two distinct surface current products (CFSR currents with a resolution of 0.5 • , and HYCOM currents with a resolution of 1/12 • ), studied the spatial patterns of differences in between simulations with and without currents. They found that the differences in were significantly greater in the simulation with HYCOM currents, pointing to the importance of the resolution and the accuracy of the current forcing dataset in correctly representing wave-current interactions. Here, we use the Bluelink Reanalysis surface current dataset (with a resolution of 0.1 • ) as forcing for our simulations. To test the importance of the wave model resolution on representing wave-current interaction processes, apart from the simulations described in Section 2.1 (i.e., a global grid of 0.4 • resolution with another grid of 0.1 • resolution across the Southern Ocean), another set of simulations (with and without currents) was undertaken using only a global grid with a 0.4 • spatial resolution but otherwise keeping the same configuration described in Section 2.1. The significant wave height output from the model in both currentforced simulations (with and without the higher resolution grid in the Southern Ocean) was compared against altimeter observations as in Section 3.1. In each simulation, we found a significant improvement in the estimations with the inclusion of currents. However, when we explored the result sensitivity to increased resolution in the Southern Ocean, any improvements due to this higher resolution were much less obvious. Fig. 10 shows the difference in performance metrics (linear correlation coefficient, mean bias and RMSE) between the simulations with and without the higher resolution grid (in each case, surface currents were included as a forcing of the model).

Importance of wave model resolution
Analogously to Fig. 1, a positive value in the correlation differences ( Fig. 10.a) means that the correlation is increased when using the higher-resolution grid, and a negative value in Fig. 10.b or c represents a reduction in the mean bias or RMSE and hence an improvement in the model performance. As expected, there are no changes in the model performance in the Northern Hemisphere. In the Southern Ocean, the average of the mean bias (absolute value) across all grid points is reduced from 29 cm to 27.5 cm, and the RMSE from 35.5 cm to 34.9 cm. Many regions in deep waters of the Southern Ocean show a reduction of the mean bias by up to ∼10 cm and of RMSE by up to ∼5 cm. However, in various coastal areas and downstream of them, the (current-forced) simulation with the higher resolution grid performs worse than the simulation with the global grid only. This issue could be related to the fact that the higher resolution grid reaches shallower areas, where presumably our global model can incur errors and misrepresent bathymetric features and transformations of shoaling waves. On the other hand, errors in the current data from the Bluelink Reanalysis in coastal areas could also contribute to this effect. Nonetheless, deep-water areas show an overall slight improvement in estimates. Rapizo et al. (2018) previously found the inclusion of current forcing improved the estimates of peak wave direction only slightly and worsened the mean wave direction estimations at the SOFS buoy location, and they related this issue to the accuracy and resolution of the current forcing. The rationale is that a wave model with a spatial resolution that is too coarse will incorporate current data in an eddypermitting manner (it would only represent eddies when the Rossby radius of deformation is much larger than the grid spacing), whereas an eddy-resolving model would be expected to be more accurate. While the accuracy of the current forcing is paramount, the wave model resolution is also an important consideration. This would be especially true for representing current-induced refraction and changes in wave direction.
Here, we tested the differences in the performance of wave models forced with currents, but with model grids of differing resolution (0.4 • and 0.1 • ). First, we selected the same storm event in Sydney shown in Fig. 5 and analysed how both simulations (with and without the high-resolution Southern Ocean grid) represented the differences in wave height and direction when including currents, in an area surrounding the south part of Australia (Fig. 11). This was the most intense wave-current interaction event at the Sydney buoy location for 2014-2016. Fig. 11 shows that significant changes in , higher than ±1 m, occurred in the Tasman Sea and the Southern Ocean that day, as well as complex changes (although of lower amplitude) off western Australia. Off eastern Australia, the current-induced refraction was particularly intense that day, producing changes in wave direction of more than ±20 • . The ACC and the Leeuwin Current also produced changes of wave direction, although of lower magnitude. While undoubtedly these features are resolved in a more detailed manner in the higher-resolution simulation, the lower resolution simulation does not fail to represent the general current-induced refraction patterns, nor the changes in wave height. Besides, the intensity of the changes is very similar in both cases. Fig. 11 shows the differences in the representation of current-induced refraction between the higher and lower resolution simulations, in terms of (in panel (c)) and (in panel (f)): the differences in wave direction are generally smaller than 4 • , and differences in are well below ±15 cm; for current-induced changes in wave height of more than ±1 m, this represents a maximum difference of ∼10% (for the most intense wave-current interaction event at the Sydney location). A comparison of these current-forced simulations (with and without the higher resolution grid in the Southern Ocean) with the Australian buoy network reveals that while some improvements are attained with the high-resolution simulation (especially in terms of wave period statistics, and also correlation in wave direction off eastern Australia), there are also areas where there is a worsening of wave estimates ( Figure S8). In short, depending on the application, running a (current-forced) wave simulation with a 0.4 • grid provides a reasonable representation of wave-current interaction processes. Some areas will result in improvements in certain wave parameters when using a high-resolution grid.

Discussion and conclusions
In this study, we have shown that a significant improvement in the performance of wave simulations can be achieved by including ocean surface current data from a global reanalysis as an additional forcing to the wave model. The correlation values for the current-forced simulation are statistically larger (at the 95% confidence level) than those of the simulation without currents for every month of the year and for every ocean basin ( Fig. 2 and S3). We found there is an important reduction of the mean bias and RMSE in most areas of the world, and particularly in the Southern Ocean: here, global wave models tend to overestimate the observed wave heights, in part due to biases in the wind data used as input, drifting icebergs that block the wave energy but not accounted for in the model, and not considering surface current forcing. Since the westerly winds blow in the same direction as the ACC, the relative wind that acts on the surface of the ocean is lower than the real wind. WAVEWATCH III incorporates this effect by balancing current and wind vectors, assuming a proportionality coefficient of 1 ( in Eq. (1)). In reality, the wind, waves and currents will interact and modify the boundary layer profile and hence the relative wind will be different to the one given by Eq. (1). In this implementation, the choice of = 1 most likely produces an overestimation of the relative wind in most conditions. A fully coupled ocean-wave-atmosphere model would be necessary to represent these interactions more accurately (see, for example, Renault et al., 2016). Nevertheless, this simplistic approach partially captures this process and helps to reduce the bias in wave heights. Other areas with improvement in the wave model performance are the AC and EACC regions and the south Indonesian islands. A large area eastward of the Drake Passage also exhibits a substantial improvement in wave heights. In the Agulhas retroflection area, there are alternating bands of improvement and worsening (in terms of mean bias and RMSE), which suggest errors in the spatiotemporal location of the mean current, meanders and eddies in the Bluelink reanalysis data used as input. The model performance in the western Pacific region is also worse when incorporating currents. This is a very challenging area for ocean modelling because of its complex geography, with multiple atolls and small islands covering this region. A global simulation with a relatively coarse grid resolution is unable to capture the complex geographic and bathymetric features of this area. Nevertheless, as a global average, the current-forced simulation performs better than the simulation without currents in terms of correlation, bias and RMSE throughout the year. Moreover, the improvement is also seen if separating the data into different ocean basins. This improvement can yield major benefits for a broad range of practicalities that involve wave modelling, such as wave forecasts, industry-based applications, coastal studies, future wave projections, among others. In addition, many wave modelling studies have historically disregarded the influence of ocean currents (e.g., coastal studies: Harley et al., 2009;swell climate analysis: Semedo et al., 2011; assessments of projected changes of wave climate: Morim et al., 2019), although that trend has been reversing in recent years (Ardhuin et al., 2017;Rapizo et al., 2018;Hegermiller et al., 2019).
A comparison was carried out between the simulation outputs and buoy observations from the Australian and NDBC buoy networks. Although most buoys are in intermediate and even shallow waters, and the spatial resolution of the global grid used in this study (0.4 • globally and 0.1 • in the Southern Ocean) is too coarse to properly capture the wave processes and transformations that occur at the coast, this comparison is helpful to evaluate how the large-scale wave-current interactions influence the coastal wave climate. In addition, the model suitably represents the observed wave parameters, with average correlation values of 0.91, 0.57 and 0.57 for , and respectively, for the Australian buoy network (Table S1 of the Supporting Information section (SI)). For the NDBC buoy network, the average correlation values are 0.94, 0.58 and 0.68 for , and respectively (Table  S2 of the SI). The main conclusion we draw from this comparison is that introducing current forcing into wave simulations produces a substantial improvement in the model estimates of wave height, period and direction, albeit some statistics are made worse in a few locations. It has been previously shown that tidal currents (not included in the Bluelink reanalysis) can modulate the coastal wave climate. Here, we see that the large-scale circulation of the East Australian Current also has a significant impact on the wave properties measured at the east coast of Australia. However, as noted by Oke et al. (2013) and Chiswell and Rickard (2014), the ocean surface currents from the Bluelink reanalysis present significant differences with the observations in the EAC and AAC. These errors in the current forcing dataset can be a factor that contributes to the worsening in the representation of wave parameters observed for some buoys: e.g., Figure 17 of Chiswell and Rickard (2014) shows large discrepancies in the ocean surface velocities in south-east Australia close to the Eden buoy location, where we observe a statistically significant worsening in the wave direction estimates for the simulation with currents.
On the other hand, while the average improvement is small and not statistically significant for some locations, currents can have a significant impact on waves on short intervals of time. Fig. 5 shows the , and measured by the Sydney buoy during the last days of March 2016. At this location, between 2014 and 2016, the wave heights of the simulation without currents were changed by < 5% of their value during 50% of the time due to the inclusion of current forcing. For 80% of the time, the changes in were < 10%. However, during the event shown in Fig. 5 in the last days of 2016, there were significant changes in wave height of up to 1.2 m (2.8 m in the current-forced simulation, 75% higher than 1.6 m in the simulation without currents), 5s in peak period (a 50% change from ∼11s to ∼16s) and 80 • in direction. Further, the current-forced simulation output values are substantially closer to the observed wave conditions than the simulation without currents. This has direct implications for coastal management studies. Not considering currents led to a significant underestimation of wave heights during this event, which can impact other wave derived quantities, such as wave energy fluxes or estimations of wave setup. In the long term, wave-current interactions could affect estimates of wave contributions to coastal sea-levels, maximum and return period wave heights, design parameters of coastal structures, among many others. Importantly, wave direction is substantially misrepresented, and this variable plays a critical role in many important coastal processes, such as alongshore sediment transport (Komar, 1971). Moreover, Harley et al. (2017) investigated the June 2016 Collaroy Beach storm (the most damaging storm in the last 40 years in that area) and observed that the anomalous wave direction (instead of wave heights) was primarily responsible for making this storm so destructive. Here we show that the meandering flow of the EAC can significantly influence the wave climate in the Tasman Sea and modulate the wave direction at the coast.
Considering that an overall improvement in wave simulation is achieved by introducing current forcing, the differences in various wave parameters were analysed to further understand wave-current interaction processes and their implications. Fig. 7 shows the percental changes in wave heights, mean periods and directions. Wave heights are reduced in most parts of the world, except in the AC and EAC close to the coast, the ECC and in some regions of the North Indian Ocean, and the northern section of the CC (there is also a slight increase in wave heights close to the Antarctic Continent, but since there is not an improvement of the model performance by including currents in this area, this feature is disregarded). The well-known overestimation of observed wave heights by the model is reduced by accounting for the effect of ocean currents on waves. Waves also have shorter periods in the current-forced simulation in most ocean basins, and larger periods in the Tropics, at the east and northwest coasts of Australia and southeast coast of Africa. A correct representation of wave periods is fundamental to improve the propagation of swell waves by numerical wave models. Because the wave propagation speed in deep-waters depends on the wave period, the longer wave periods of the simulation without currents on most areas of the world (Fig. 7.b) would translate into faster-propagating waves, and therefore wave models would tend to predict early swell arrival times (Jiang et al., 2016). Fig. 8 shows that, in general, waves tend to increase their wave heights and periods when they propagate against the main flow, and vice versa. This agrees with results shown in Fig. 7: In the Agulhas Current or East Australian Current regions, where wave heights and periods are increased (decreased) where the predominant southerly waves propagate against (with) the current. The changes in wave properties in the equatorial region are more challenging to scrutinize, given the complex and multi-modal wave climate of the region. However, the largest increases in wave period seem to arise from the Northern and Southern Hemisphere swell waves ( Figures S5 and S6). There are also large percental changes in wave period in the Arabian Sea and Bay of Bengal. Wave-current interactions in this area deserve a more detailed analysis, given the complexity of the main wave climate as well as the strong seasonality in the (reversing) ocean circulation in the area, both strongly influenced by the Asian Monsoon (e.g., Shankar et al., 2002). Although the consideration of current forcing leads to an overall improvement in wave period estimates, the physical mechanisms that lead to the observed differences in wave period in many areas of the world are not yet fully understood. Importantly, there are considerable average differences in wave direction, particularly in the Tropics, the ITF and the AC and EAC regions. The positive/negative differences at the right/left sides of the current's main flow direction in the AC, MC, EAC, BC, NBC and GS represent clockwise/anticlockwise changes in wave direction, in agreement with results from Kenyon (1971). This suggests that waves propagating opposite to the currents' main direction may focus and/or be trapped in the central axis of the currents. This is most likely to occur in the AC, MC and EAC which are exposed to Southern Ocean generated swell waves. The changes in wave direction in the AC, MC and EAC seem, at this scale, to reach the coast. Further, Figs. 3 and 5 show that an overall improvement in wave direction is achieved with the current-forced simulation along eastern Australia. Given that wave direction has historically been challenging to simulate in wave models, these results may be beneficial to improving estimates of wave direction in future studies, as well as providing improved understanding and representation of coastal processes, for which wave direction plays a critical role (e.g., along-shore sediment transport). Furthermore, this will undoubtedly affect global wind-wave climate projections studies, which so far have been carried out using surface wind fields taken from the Coupled Model Intercomparison Project (CMIP), with no ocean surface current forcing (Hemer et al., 2013;Morim et al., 2019). The changes in wave heights, periods and direction due to the inclusion of currents shown in Fig. 7 are approximately of the same magnitude (and in some cases greater) than the projected changes in these wave parameters calculated by Hemer et al. (2013) and Morim et al. (2019). Wave-current interactions could be particularly relevant for wave projections near and along the east coast of Australia, given the significant wave refraction induced by the EAC (Figs. 5 and 7.c) and the on-going and projected increase in EAC transport and poleward extension as climate changes into the future (Oliver and Holbrook, 2014;Oliver et al., 2015).
The changes in wave conditions by the inclusion of currents also modify the ocean surface roughness, which has consequences for the coupled climate system. The whitecap cover is reduced by 0.02-0.04 in the ACC, EACC, NBC, CC and EC, and it is increased by 0.02 in the ECC region. The global average whitecap coverage is estimated to be around 3.6% (Blanchard, 1963). In the equatorial region, considering current forcing changes the whitecap cover by more than 25%. The friction velocity and atmosphere-to-ocean energy flux show similar patterns, although the changes in friction velocity range between ±10% and the energy flux into the ocean changes by ±50% (it should be noted that the atmosphere-to-ocean momentum flux changes are equivalent to those of the energy flux). The generation of whitecaps in deep water plays an important role in the enhancement of heat, momentum, gas and particle transfer between the ocean and the atmosphere, as well as being high sunlight-reflectance areas that contribute to increase the planetary albedo and are therefore fundamental in global radiation budget calculations (Cavaleri et al., 2012). The results shown in Figs. 7, 9 and S7 point to the critical importance of wave-current interaction for climate modelling in Tropical areas. The inclusion of currents leads to reduced ocean surface roughness and whitecap coverage in the EC region, and reduces the energy and momentum transfers from the atmosphere. The contrary happens in the ECC area. Since waves can also influence surface currents and winds, better representation of these processes could be achieved by implementing a fully-coupled atmosphere-wave-ocean model. In addition, a comparison of the relative wind computed by the wave model with measured winds (either in situ or from satellites) is necessary to better assess the validity of these results. Furthermore, significant changes in Stokes drift speed (∼±30%) and direction (∼±10 • ) in western boundary currents and particularly in the equatorial region are attained by incorporating current forcing. The changes in Stokes drift direction follow a very similar pattern to the changes in friction velocity direction ( Figure S7.a). Currents are a fundamental consideration in particle tracking studies (e.g., Onink et al., 2019), which often use independent ocean circulation and Stokes drift data (as opposed to fully coupled models). Given the inherent turbulence of the upper ocean, particle tracking studies must use a large number of particles to provide statistical consistency to their results. Hence, if the change in Stokes drift speed or direction can advect a particle to a slightly different position, this could have significant impacts on the final trajectory of the particle.
Finally, the sensitivity of the wave-current results to wave model resolution was explored in Section 3.3. Current-forced wave simulations with and without the higher resolution (0.1 • ) grid in the Southern Ocean were compared. Fig. 10 shows that by incorporating an eddyresolving model configuration (instead of an eddy-permitting one), slightly better estimates of significant wave height in deep waters were achieved. This is likely a product of a more accurate representation of the surface wind field from the CFSR reanalysis (of 0.2 • spatial resolution). Based on Fig. 11 of Rapizo et al. (2018), we hypothesized that the patterns of current-induced refraction would be radically different in both simulations (high and low resolution). Instead, we found ( Fig. 11) that changes in the significant wave height and mean wave direction follow the same patterns in both simulations -albeit that the higher resolution simulation represents these in a more detailed manner. The changes in wave period due to the inclusion of current-forcing are also equivalent in both simulations (not shown). The most extreme wavecurrent interaction event at the Sydney buoy location for 2014-2016 was chosen for Fig. 11. However, other time periods were also selected and analysed, all yielding similar results. Figure S8 compares these two simulations against the buoy observations in Australia, and shows that the improvements achieved when incorporating the high-resolution grid are mostly in terms of peak wave period. Estimates of wave height and direction are improved for some locations but worsened for others. While the introduction of current forcing greatly improved the performance of global wave simulations, using a high-resolution grid did not necessarily translate into significant further improvements. Rather, we found that current-induced changes in wave parameters were relatively well represented using the coarser (0.4 • ) model grid, in comparison with the eddy-resolving 0.1 • -grid configuration in the Southern Ocean. This is appreciated by examining the representation of wave direction against observations from the Australian buoys (within the 0.1 • grid domain) and NDBC buoys (within the 0.4 • global grid): Figures 3 and 4 show that there is an improvement in wave direction with the inclusion of currents for most of the buoys. While here we have only compared wave simulations with resolutions of 0.4 • and 0.1 • , ocean surface currents of smaller spatial scales can also have a substantial impact on waves (e.g., Ardhuin et al., 2017). In addition, an assessment of the sensitivity of wave modelling to the inclusion of surface current forcing from different ocean reanalysis products could yield a greater understanding of global and regional wave-current interaction processes.
In conclusion, surface currents impact wind-waves in various ways, and WAVEWATCH III certainly captures these processes in a one-way coupling implementation, broadly improving simulation performance. Undoubtedly, there will be some errors in the current forcing data taken from the BRAN reanalysis, and a fully coupled modelling approach will be more beneficial to capture the complex interactions between the ocean, atmosphere and waves and provide more realistic outputs. Nevertheless, we have shown that global deep-water and coastal wave parameter estimates can be significantly improved with the methodology applied in this study, and that major benefits can be achieved based on this approach for the wave modelling community.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.