Surface Currents and Significant Wave Height Gradients: Matching Numerical Models and High-Resolution Altimeter Wave Heights in the Agulhas Current Region

Indeed, all these features should have a large on waves and in particular on wave Still, the quantitative impact of on is not well known due to the complexity of the random and currents that are found in the ocean and the lack observations both and Here, we compare novel satellite altimetry data and state of the art phase-averaged numerical wave models forced both by wind and currents. Currents used are taken from the oceanic model Coastal and Regional Ocean COmmunity, run at high resolution. The influence of current field resolution is investigated by applying Gaussian filters of different width to that same high-resolution current field. We find that a numerical wave model that uses currents with resolutions of km or less and a directional resolution of 7.5° can provide accurate representations of the significant wave height gradients found in the Agulhas current. Using smoother current fields such as derived from altimeters measurements coarse directional spectral resolution or larger directional spread of the wave model underestimates gradients and extreme wave heights. Hence, satellite altimetry provides high-resolution wave height with a gradient magnitude that is highly sensitive to underlying surface current gradients, at resolutions that may not be resolved by today's altimeters measurements.


Introduction
Surface gravity waves generated by wind (hereinafter waves) interact with surface currents at all scales due to a wide range of processes (Phillips, 1977). Except for very short fetch near the coast or for the shortest wave components, the growth of waves in the presence of winds is only significant over large scales, so that the local gradients in the dominant wave properties are generally dominated by current gradients (Phillips, 1984). In the ocean, it appears that refraction, which focuses wave energy in current jets that flow in the wave direction, is probably the dominant source of variations of wave heights at scales 50-200 km with a minimal effect of wind gradients (Ardhuin et al., 2017). For currents speeds much weaker than the waves phase speed, it is the rotational part of the current that is expected to explain the variations in wave direc-Abstract Advances in the understanding and modeling of surface currents have revealed the importance of internal waves, mesoscale and submesoscale features. Indeed, all these features should have a large influence on wind waves and in particular on wave heights. Still, the quantitative impact of currents on waves is not well known due to the complexity of the random wavefields and currents that are found in the ocean and the lack of observations of both currents and waves at scales shorter than 150 km. Here, we compare novel satellite altimetry data and state of the art phase-averaged numerical wave models forced both by wind and currents. Currents used are taken from the oceanic model Coastal and Regional Ocean COmmunity, run at high resolution. The influence of current field resolution is investigated by applying Gaussian filters of different width to that same high-resolution current field. We find that a numerical wave model that uses currents with resolutions of ∼30 km or less and a directional resolution of 7.5° can provide accurate representations of the significant wave height gradients found in the Agulhas current. Using smoother current fields such as derived from altimeters measurements alone, coarse directional spectral resolution or larger directional spread of the wave model generally underestimates gradients and extreme wave heights. Hence, satellite altimetry provides high-resolution wave height with a gradient magnitude that is highly sensitive to underlying surface current gradients, at resolutions that may not be resolved by today's altimeters measurements. This is demonstrated here for relatively steady currents averaged over 3 years.

Plain Language Summary
Mariners have learned to be wary of severe sea states, especially in strong currents like the Agulhas that flows along the South African coast, where wave heights in the current can be several meters taller than in the surrounding waters. Mariners have also learned to spot currents by watching the water ahead of them. Here, we use satellite measurements of wave heights and a numerical wave model to understand the parameters that control the spatial variation of wave heights across currents. We particularly question the necessary current magnitude and gradient that are required to explain observed wave height gradients. Modeled gradients fade for smooth surface currents like surface currents estimated from satellite measurements of sea level or typical global ocean circulation models. Also, numerical experiments have shown that when incident waves have a narrow range of directions, wave height gradients are sharper. A good wave model should thus resolve both the current features, with a spatial resolution better than 30 km, and the range of wave directions, typically using 48 directions or more. Such a good model can then be used to evaluate the quality of modeled ocean currents by matching the modeled strength of wave height gradients with measurements.
Recent advances in understandings and in ocean modeling of surface ocean dynamic show that the upper ocean is highly energetic not only at the mesoscale, for which the flow is in quasi-geostrophic balance, but also at smaller scales (submesoscales) (McWilliams, 2016). Further, strong ocean currents are associated with sharp and asymmetric velocity fronts, with larger positive vorticity maxima in the Northern Hemisphere (e.g., Gula et al., 2015). Also, the generation of large surface waves has been shown to occur in the presence of strong internal waves (Osborne & Burch, 1980). All these small-scale current features may contain as much surface kinetic energy (KE) as the mesoscales but it is not clear how much they influence the waves. Refraction theory tells us that changes in wave direction for a given wave frequency are the product of the current vorticity magnitude and the scale of the current feature, so that a localized high vorticity may have the same effect as a distributed but lower vorticity. But in practice, ocean waves are random and the different components of their relatively broad spectrum are affected in different ways by the surface vorticity.
The evolution of the wavefield, represented by the wave action spectral densities N(σ, θ), with σ the wave frequency in the frame of reference moving with the local current and θ the wave propagation direction generally follows the wave action equation (Komen et al., 1994;Tolman & Booij, 1998): The contributions of surface currents in Equation 1 come into the advection speeds in longitude   and latitude   , which is the sum of the intrinsic group speed and the surface current, the refraction velocity   , the change of frequency velocity  , and in the right-hand-side source term S because the effective wind velocity that generates waves is the vector difference of wind and surface current velocities (e.g., Ardhuin et al., 2017).
Because the effect of refraction   at position (λ, ϕ) combines with the advection in a new direction θ to produce a change in wave action N at another location (λ′, ϕ′), there is no simple relationship between the current field and wavefield, in other words, surface currents have a nonlocal effect on the distribution of the wave action in the current field. White and Fornberg (1998) have shown theoretically that the spatial distribution of refraction-induced focusing can be predicted for monochromatic waves over a random current with a narrow band spectrum. Still, that does not say much about the spatial distribution of wave heights in this case. The problem is more complex for broad band current spectrum and random waves, for which the significant wave height combines all the spectral components: Guided by these theoretical insights and the solid foundation of the wave action equation (e.g., White, 1999), our understanding of the effects of surface currents on wave height in the real ocean has relied on numerical simulations using Equation 1. These simulations are fairly successful for well-known tidal currents (e.g., Ardhuin et al., 2012), but there are very little data to validate modeled currents and waves in other regions. For example, wave simulations in the Gulf Stream and Drake Passage suggest that the patterns of H s field induced by surface currents are dominated by the refraction (Ardhuin et al., 2017), with a significant impact of small-scale currents. These modeling results could not be validated using standard satellite altimeter data that is dominated by noise for along-track wavelengths shorter than 100 km (Dibarboure et al., 2014). The development of new denoising techniques has revealed a systematic relation between wave height gradients and current vorticity (Quilfen & Chapron, 2019;Quilfen et al., 2018). These filtered data have been compared to preliminary simulations in the Agulhas current using Equation 1 solved by either finite-difference techniques or ray tracing. These comparisons have highlighted the importance of the directional width of the wave spectrum, with stronger H s gradients obtained for narrower incident wave spectra even when only large-scale currents, as derived from gridded altimetry data were used (Quilfen et al., 2018).
These two previous studies by Ardhuin et al. (2017) and Quilfen et al. (2018) have suggested two possible reasons for sharp H s gradient: namely the presence of sharp current gradients or the strong local focalization of waves on a smooth current field. Figure 1 illustrates the first possibility over the Agulhas current, using either large-scale currents of gridded altimetry or a high-resolution modeled current, both described in detail in Section 2.
The present work aims at consolidating these previous analyses and contribute to answering the following questions. What are the parameters controlling the spatial variability of wave heights in a realistic current field? How can these be best reproduced by numerical models? In particular, we focus on the effect of the spatial resolution of the current field and angular discretization of the wave model in relation with the directional spread of wave spectra. Here, we focus on the Agulhas current because of the strong H s signature that is easily captured by satellite altimeters. Further work will be needed for other wave and current regimes.
The numerical model setup and data are presented in Section 2. Results follow in Section 3, with a discussion of the influence of the surface currents resolution in Section 4. Finally, we will conclude this wave-current interactions study in Section 5.

Satellite and Modeling Data for Waves in the Agulhas Current
The Agulhas current system is one of the most intense western boundary currents, with velocities exceeding 2.5 m s −1 along the East coast of South Africa, before retroflecting back into the Indian Ocean with large ring eddies shed in the South Atlantic Ocean (Beal et al., 2011;Tedesco et al., 2019). The Agulhas current system is also exposed to very large waves from the Southern Ocean (Young, 1999). are the only available measurement of wave heights. This is particularly the case in strong current regions where moored buoys are more difficult to install. Further, H s measurements along the satellite ground track provide a unique view of the spatial variations of H s , although along one dimension only. Until recently, the analysis of H s variations was limited to wavelengths larger than 100 km, due to the noise associated with the tracking methods used to interpret altimeter waveforms (Ardhuin et al., 2017;Sandwell et al., 2014). The successful application of Empirical Mode Decomposition (EMD) (Huang et al., 1998) to the denoising of H s along-track series now makes it possible to investigate much smaller scales, possibly down to 15 km wavelength or less (Quilfen & Chapron, 2019). Here, we use denoised wave heights from the European Space Agency (ESA) Sea State Climate Change Initiative (SeaState-CCI) version 1 database (Dodet et al., 2020), that uses this denoising technique applied to calibrated Geophysical Data Records from CNES and ESA for the Jason-2, CryoSat-2, and SARAL/AltiKa missions. The analysis of 3 years from 2014 to 2016 in our region of interest gives a total of 4,746 satellite tracks, with one example shown in Figure 1. In the ESA CCI-Sea-States product, Dodet et al. (2020) give a root mean square error (RMSE) on Hs data equal to 0.21, 0.2, and 0.21 m for Jason-2, CryoSat-2, and SARAL/AltiKa missions respectively. The uncertainty of Hs is computed from in situ waves buoys located beyond 200 km offshore. The RMSEs are similar for the three missions both for calibrated and denoised (EMD method) altimeter data ( Figure 10 of Dodet et al. [2020]).

Numerical Wave Model
Our numerical wave model is based on the WAVEWATCH III modeling framework (The WAVEWATCH III® Development Group, 2016) that integrates the action balance Equation 1, discretized on a regular latitude-longitude grid with a resolution of 1/30°. Our baseline configuration uses a spectral discretization into 32 frequencies from 0.037 to 0.7 Hz and 48 directions (Δθ = 7.5°). This model is forced by surface currents, as detailed below, together with operational hourly wind forecasts from the European Centre for Medium-Range Weather Forecasts, at 1/8° resolution. The overall time step used to solve Equation 1 is 390 s, and the solution is obtained with a splitting technique (Tolman, 1992), with a spatial advection step of 130 s, a refraction step of 18 s, and an automatically adjusted source term integration step that can be as short as 10 s. We define the boundaries with three hourly wave spectra from a global model configuration that uses the same wind fields but no current, a spatial resolution of 0.5° and the same spectral discretization as our Agulhas wave model. The wave model grid covers the domain shown in Figure 1, from 40° to 30°S and 16° to 30°W.
The signature of the Agulhas systems is clearly visible in the modeled H s field with a band of larger wave heights. On the example in Figure 1a, one can observe the effect of the main Agulhas current along the coast, including a meander known as a "Natal pulse," located at 29°E, upstream of Port Elisabeth. Large current structures typically have multiple parallel branches caused by the straining of the large-scale field and very sharp boundaries (Figure 1d). In contrast, the H s field computed with the model using surface currents estimated from altimetry measurements (Globcurrent) has blurred patterns (Figure 1b), caused by surface currents with broader features and less intense maxima values (Figure 1e). The large-scale circulation estimated from altimeter data although less energetic is coherent with the Coastal and Regional Ocean COmmunity (CROCO) output snapshot: Agulhas current along the coast, retroflexion, and Agulhas return current. For smaller-scale features, all the 10-100 km structures are missing in the Globcurrent product, including meanders of the Agulhas current along the coast, from 28° to 23°E which play an important role in the current stability (Tedesco et al., 2019). Also, the Agulhas current has a similar transport in both current fields but much sharper gradients and higher maxima, up to 3 m/s in the CROCO model result compared to 2 m/s in Globcurrent.
Altimeter measurements show a narrow H s maximum around 37° in the Agulhas current upstream of the retroflexion (Figure 1c). This narrow peak in H s is closer to the one obtained with the CROCO currents, while the Globcurrent current fields lead to a broad H s maximum.

Currents Fields Used for Forcing the Wave Model
Given the large influence of surface current details, we have designed a series of simulations with currents at different resolutions. These current fields are based on surface current estimates from the CROCO model (Debreu et al., 2012) without data assimilation nor tidal forcing with a resolution of 1/36° both in latitude and longitude. The CROCO model domain is larger than the WW3 model domain that is shown in Figure 1 and covers 15.1°-33.7°E and 40.4°-27.2°S. This CROCO model configuration is expected to produce surface currents that are statistically consistent with the real ocean and has been used for several process studies (Tedesco et al., 2019). However, for any particular time and location, the variable current structure is not expected to reproduce the stochastic behavior of the ocean as no data assimilation is used within the model domain. The CROCO model has been forced at the surface by the ERA-interim reanalysis and boundaries have been forced by a global reanalysis GLORYS. We have also used low-pass-filtered CROCO currents as an input forcing for the wave model. These are obtained by applying an isotropic two-dimensional Gaussian filter on both zonal and meridional components of the current velocity vector. This filter is defined by its standard deviation σ c (Figure 2a). We emphasize that the alternative approach of rerunning CROCO at different resolutions may produce very different results and would require some tuning of each model configuration that is beyond the scope of the present work.
The filtered current fields effective resolution is the result of the convolution of the Gaussian filter and the original current field. Theoretically, the spectrum of the filtered current is the product of the original current spectrum and the spectrum of the Gaussian filter. In practice, it means that the current spectrum rolls of sharply for wavelengths shorter than L c = 4σ c or an effective resolution of 2σ c . An example of omnidirectional KE (KE = <U> 2 + <V> 2 , <.> 2 denotes the variable's variance) spectra are given in Figure 2b. These spectra are from an azimuthal integration of a two-dimensional spectral analysis applied on surface current fields at a specific time. Details of the spectral method are available in Section 2.2 of Ardhuin et al. (2017). These 1-D spectra represent the distribution of variance across spatial scales. The clear drop of variance for filtered surface currents at high wavenumbers shows that the Gaussian filtering process has removed small spatial scales in the Agulhas region.
Seven surface current fields have thereby been created, with effective resolutions ranging from 10 to 100 km.  The filtering of the current field results in the removal of small-scale structures, including small mesoscale eddies and filaments, as well as the smoothing of the large-scale structures. Alternatively, we also used a surface current forcing taken from the Globcurrent product (Rio et al., 2014). This Globcurrent product has a spatial resolution of 1/4° both in latitude and longitude and is temporally resolved at 1 day. It provides the geostrophic component of the total surface currents estimated from the sea surface height (SSH) measured by altimeters, and a mean dynamic topography that combines other data sources (Rio et al., 2014). A similar spectral analysis described above has been applied on Globcurrent product and revealed its effective resolution 150 km. The 60 km resolution-filtered CROCO current has scales similar to those in the Globcurrent field, with a lower surface currents intensity for filtered surface current (due to filtering process). We note that the surface relative vorticity ζ of the filtered current ( Figure 3) is similar to the ones presented in Figure  17c of Chelton et al. (2019) in the Coastal California current for similar resolution (few kilometers, 20 and 80 km). Figures 3a-3d illustrate how the wave height patterns follow the surface vorticity patterns as already shown in Figure 13 of Quilfen et al. (2018). Figure 3 (left) shows a H s maximum where the normalized vorticity is positive in the main stream of the Agulhas (southwestward) and also show that the H s gradient is sharp for WW3 results forced with high-resolution currents and become blurred for poorly resolved surface current. We have run our wave model during 3 years, from 2014 to 2016, with the appropriate surface currents (fully resolved from CROCO model, filtered and estimated by altimetry), wind and boundary conditions forcings.

Spatial Variability of H s in Realistic Surface Currents Field
Wave-current interactions have been simulated in the Agulhas current from 2014 to 2016. Filtered altimetry data have been studied for the same time frame and all model outputs have been interpolated in time and space on those altimeters tracks. One example of model-satellite comparison is displayed in Figure 1c. Except for the topographically trapped flow patterns, the high-resolution CROCO model is not expected to have current features in the same place as the real features, but it may still have realistic eddy sizes and meander shapes. We will thus compare the statistical properties of modeled and measured H s .
In particular, we consider the statistical properties of the along-track H s gradient defined as with dr the along-track distance between successive 1 Hz measurements (dr is typically 7 km), and ΔH s the difference between successive H s measurements taken 1 s apart. Statistics of ∇H s have been interpolated on a regular grid with a resolution of 1/8° by 1/8° in longitude and latitude. The mean values are shown in Figure 4, ranging from 0 to 3 cm per km.
A few high values of the H s gradient right at the coast are clearly visible for the simulation without current.
These high values can be explained by partial sheltering caused by headlands, all the large gradients appear in regions of strong current gradients, and specifically in the main Agulhas current, from 29°E 33°S to 17.5°E 39°S. The values of the mean ∇H s measured in the main Agulhas branch are in the range of 1.5-3 cm/km (Figure 4i) which is remarkably high, and corresponds to the maximum values shown in Figure 1. These persistent maximum gradients are located exactly where the model has the strongest current, and where the largest H s gradients are also predicted in Figure 4a. This is the well-known region of strong focalization of waves caused by wave refraction over the current (Gutshabash & Lavrenov, 1986;Kudryavtsev et al., 2017;Quilfen & Chapron, 2019). Indeed when propagating against a current that is uniform in the flow direction, waves of a given period and direction can be trapped: when coming from the center of the current toward its edge they turn back toward the center at the location where the current reaches a certain value (Kenyon, 1971). The waves behavior is similar to the propagation of light waves along an optical fiber where light waves are trapped and propagate within a range of specific refraction's index values that depends on their initial incidence angle. Quilfen and Chapron (2019) have demonstrated with ray tracing and assuming the wave action is conserved along the ray, that where waves are trapped, strong ∇H s are measured. Figure 4 shows that the maximum ∇H s signal is upstream 26°E, where the main Agulhas current is known to be stable. Downstream of 26°E, the current is bimodal with occasional disturbances known as Natal pulses (Lutjeharms & Roberts, 1998;Paldor & Lutjeharms, 2009).
Around 22°E, the Agulhas current comes off the Agulhas Bank and the current direction veers to the south, which probably explains the lower values of ∇H s as the current direction is less favorable for trapping the dominant south-westerly waves, resulting in this lower gradient of wave heights. Beyond that point, ∇H s increases again but it is more spread out in the north-south direction.
MARECHAL AND ARDHUIN 10.1029/2020JC016564 8 of 16 Nowhere does the much coarser and weaker current in the Globcurrent product produces H s gradients larger than 2 cm/km ( Figure 4g). Yet, the Globcurrent product leads to modeled gradients in the retroflexion region, around 38° S, 25 °E, that are similar to those given by the CROCO model, both weaker than observed. ∇H s in the main Agulhas current are similar for CROCO filtered at 60 km and Globcurrent, as shown in Figure 3 through the H s field. As the effective current resolution is degraded from 10 to 60 km, the mean H s gradient progressively vanishes with a particularly clear drop from 60 km (Figure 4e) to 70 km ( Figure 4f). The magnitude of the gradients can be quantified by different percentiles, as shown in Figure 5. For the 95th percentile and above, we find that 60% of the H s gradient is obtained for effective current resolutions of 30 km or less. The uncertainty of the ∇H s has been estimated from the known uncertainty of the altimeters Hs measurements (see Section 2.1). Assuming that the error over the ∇H s is uncorrelated (only contaminated by a random noise), we have perturbed the Hs measurements independently with a white noise characterized by an amplitude equal to the RMSE of each satellite missions. The computed standard deviation (not shown here) is very large and disagrees with our numerical model estimations. Hence, the RMSE (estimated only from in situ wave buoy) of the denoised data from the ESA CCI-Sea-States product seems to be overestimated.

Spectral Analysis
In order to obtain a more quantitative analysis, we perform the same spectral analysis on the model and satellite data. We use overlapping windows following Welch (1967), with the Fourier transform computed after detrending and applying a Hanning window. Results are presented in Figure 6. In order to help with the interpretation, the surface current velocity (  2 2 U V ) was also analyzed along the same tracks. One spectrum is computed for each track. All spectra have been averaged to obtain one averaged spectrum for each numerical simulation for each surface currents forcing field.  The H s spectra (Figure 6a) show that between resolutions of 200 and 30 km, and even down to the smaller resolved scale, the resolution of the surface currents drive the H s variability. For wavelengths between 50 and 100 km, simulations forced by the Globcurrent surface currents show a H s variability higher than simulations forced with surface currents filtered at 60, 70, and 100 km, whereas surface currents from Globcurrent have an effective along-track resolution around 150 km. This along-track resolution is consistent with the 150-250 km resolution of SSH gridded altimeter data in the Agulhas region (Ballarotta et al., 2019). Using a wave model forced by different surface current fields, Figure 6b reveals what was already reported by Ardhuin et al. (2017), that is, the lower the surface currents KE (<U> 2 + <V> 2 , <.> 2 denotes the variable's variance), the lower the H s spectrum. Surface KE spectrum computed from surface current taken from Globcurrent fields shows a level of variability for wavelengths in the range 50-200 km that is similar to the 40-km filtered current.
For all simulations, the shape of the spectrum of the modeled H s is very similar to the KE spectrum, and slightly steeper, around k −3.4 for H s compared to k −3.0 for the KE spectrum (exponents have been computed through a linear regression) for scales smaller than 100 km. The same behavior was found for realistic simulations in Gulf Stream and Drake Passage (Ardhuin et al., 2017). As the spectral level in the current forcing is reduced, the H s spectrum is reduced in the same proportion until it reaches a background level. For a wavelength of 100 km, this background level is around 0.08 m 2 /cycle/km, which is very close to the variability associated with the wind field in the analysis by Ardhuin et al. (2017). This parallel behavior of the H s and KE spectra may be due to the dominant balance between propagation and refraction terms in the action balance Equation 1.

Surface Current Resolution and Gradients of H s
In the ocean, surface currents are energetic at mesoscales and submesoscales, with features such as fronts, eddies, and filaments. Waves interact with those features, and refraction explains the spatial redistribution of the wave action density that results in a change of H s . In the Agulhas system, numerical wave simulations forced with highly resolved surface currents, rich in mesoscale structures, show that the small features and sharp gradients are important for simulating realistic ∇H s , statistically consistent with filtered altimeter MARECHAL AND ARDHUIN 10.1029/2020JC016564 10 of 16 Figure 6. Left panel (a), averaged significant wave height spectra from model and altimetry data. Right panel (b), averaged surface kinetic energy spectra. All spectra have been obtained by averaging all along-track spectra (4,746 tracks) from altimeters measurements (black solid line) and interpolated simulated data (in colors). The associated surface currents resolution are given in the legend. λ is the wavelength. data ( Figure 5). We find that an effective resolution of 30 km, which resolves features with wavelengths larger than 60 km is necessary to reproduce most of the wave height gradients, which can be quantified by its median value or higher percentiles shown in Figure 5. Given that the high-resolution CROCO model that provides our forcing current does not assimilate observations, its features other than the largest scales of the Agulhas current are not expected to be in the right places at the right time, it is difficult to define a wave-gradient based metric that could be used to further validate the CROCO model for different regions or scales. Quilfen et al. (2018) argued that using a finite-difference numerical wave model to solve the action balance Equation 1 generally underestimate the ∇H s , showing marked differences between finite differences and ray-tracing solutions. Here, we find that it is the choice of a large-scale current from Globcurrent that explains the relatively weak modeled H s gradient.

Directional Resolution in Wave Models
In the limit of a large number of directions and a fine spatial resolution, the solution to the wave action equation obtained here with third-order finite-difference refraction and advection schemes (Leonard, 1991;Tolman, 2002) should be identical to the one obtained with backward ray tracing (Ardhuin & Herbers, 2005;Booij et al., 1999;Longuet-Higgins, 1957;O'Reilly & Guza, 1993). In practice, the number of discrete model directions is limited by the cost in memory storage and computation time, and most wave model implementations use 24-36 directions. Given the importance of refraction in the presence of current gradients (Ardhuin et al., 2012;Holthuijsen & Tolman, 1991), we used 48 directions in the analysis presented above. We examine here the importance of the directional resolution and how the numerical solution is smoothed by the use of a small number of directions. We have thus repeated our simulations (same forcing files and same boundary conditions) different directional resolutions (Δθ), using 24 (Δθ = 15°), and 180 (Δθ = 2°) directions instead of 48 (Δθ = 7.5°). The refraction time step Δt r has been changed in proportion to keep a constant ratio Δt r /Δθ. We have further checked than reducing the other time steps had minimal effects on the solution. The spectral analysis described in Section 3.2 has been repeating for those new simulations and presented in Figure 7a. Because the Δθ = 2° simulation is extremely costly, the wave model has been run for 4 months only, from the January 1 to the April 30, 2015. The altimeters track have been extracted for the same time frame and the model outputs have been interpolated on those tracks.
Spectral analysis shows that the model setup with a finer directional resolution (N θ = 48 instead of 24) has a larger variability of H s at all scales, with an increase of the Power Spectral density (PSD) by about a factor of 2, similar to what was found for Drake Passage by Ardhuin et al. (2017). In addition, for scales smaller than 100 km, H s variability is stronger for simulations forced with higher resolution currents. Further refining the directional resolution to 180 directions gives a further increase in H s variability. When the narrow directional discretization is combined with high-resolution currents, the modeled H s spectrum is within 30% of the satellite measurements for all scales shorter than 100 km.
A typical example of spatial variability along a transect is shown in Figures 7b and 7c, with a much sharper peak of H s in the model runs using 180 or 48 directions.

Influence of Incident Waves Directional Spreading (σ θ )
We generally expect that a fine directional resolution is most important when the directional wave spectrum is very narrow. In these conditions, wave energy can be focused in a small area, as predicted by the analysis of monochromatic waves with rays traced with parallel directions outside the current region (White & Fornberg, 1998). In contrast, broad wave spectra have focal points in different locations for the different spectral components, which effectively smears the regions of maximum H s .
In order to quantify that effect in realistic conditions, we have rerun the model with modified boundary conditions. Instead of taking the directional wave spectra E(f, θ) straight from a global hindcast, we now make these spectra broader or narrower in directions, without changing the spreading along the frequency nor the mean direction at each frequency. The details of the method are given in Appendix A. The conservation of the total variance and mean direction between all original spectra and new spectra has been verified. At each frequency, the original directional spreading has been changed by ±30%. Examples of the resulting H s fields are displayed in Figures 8a-8c. Figure 8 illustrates how a decrease of σ θ induces an increase in the number of small H s structures and an amplification of structures already existing and vice versa. This is better quantified along a track that is close to the upwave (western) boundary. The left peak at 39.5°S in Figure 8d has a variation of H s from 3.45 m with a broader spectrum to 3.85 m with a narrower spectrum. This 25% change in wave energy is a typical order of magnitude. Besides the peak, some fluctuations of H s between 37° and 39°S are much reduced for the broader spectra.
Following the method used previously, we now look at the averaged H s spectra for each 1-year long simulation, with different boundary conditions. The result shows higher variability, by about 50%, at all scales for incident waves with lower values of the directional spread σ θ . The shape of the H s spectra is very similar for all simulations with a steeper slope for wavelengths shorter than 125 km.
Our simulations have confirmed that over a real current system like the Agulhas, the spatial variability is sensitive to the spectral width of the wavefield, and to the numerical resolution used in models with narrower spectra and finer resolution producing stronger gradients. Unfortunately, the directional spread is one of the worst modeled parameters (Stopa et al., 2016). More directional data, such as provided by the SWIM instrument on the China France Ocean Satellite (Hauser et al., 2017), may help design better model parameterizations and can be used for data assimilation with important impact in strong current regions. The performance of data assimilation for directional SWIM data has already been proved in the Southern Ocean, particularly for wave age and H s (Aouf et al., 2021).

Conclusion
Surface currents modify the wavefield in a complex way that is not just local (Ardhuin et al., 2017;Kudryavtsev et al., 2017;White & Fornberg, 1998), creating a spatial pattern of wave properties that can be important for applications and that may reveal properties of the ocean currents that are otherwise difficult to obtain. MARECHAL AND ARDHUIN 10.1029/2020JC016564 12 of 16 Figure 7. a) Averaged significant wave height spectra for altimeters measurements (in black) and for modeled data (colors). Blue spectra are for modeled wave height forced with surface current from Globcurrent (Glob.) and red spectra for high resolved (HR) CROCO forcing. (b) Instantaneous simulated significant wave height field highly resolved in directions (180 dirs). (c) An example of modeled wave heights interpolated along an altimeter track for different directional resolution, the location of the track is in black line on panel (b). λ is the wavelength. CROCO, Coastal and Regional Ocean COmmunity.
Large mesoscale current systems such as the Agulhas current are places where particularly strong H s gradients are found (Lavrenov, 1998;Quilfen & Chapron, 2019). Combining state of the art of wave modeling and novel-filtered altimetry data, we have investigated the factors that lead to these large gradients, and under which conditions they can be reproduced by numerical models. The present work shows that model forced with realistic and high resolved surface currents, statistically consistent with the real upper ocean dynamics and sufficiently discretized in direction, is able to capture sharp significant wave height gradient measured by satellite altimeters. These sharp gradients are much reduced in the results of wave models that are forced by surface currents derived from a combination of mean dynamic topography (Rio et al., 2014) and sea level anomalies derived from these same altimeters that measure the wave heights. This low resolution of satellite-derived currents (Ballarotta et al., 2019;Chelton et al., 2019) is related to the sparse tracks of existing and planned nadir altimeters, but it is also due to the along-track noise level in the processing used today for altimeter data.
Besides the structures of the forcing current, the numerical implementations of wave models will typically miss part of the true gradients of the wavefield due to numerical diffusion. Here, we find that high spectral resolutions, using 48 or more directions systematically produces finer details, in a way that is statistically consistent with altimeter data. Thanks to wave simulations, we have shown that this effect is most pronounced for waves coming from the Southern Ocean with a small directional spreading. Indeed, simulated waves in the Agulhas region show sharper wave height spatial structures when waves spectra forcing at the boundaries are narrow, for identical surface currents in the model parameterization.
Reproducing realistic wave height gradients is important not only for marine safety but also for studying upper ocean processes driven by wave breaking. It is also a necessity to capture sea states biases in ocean remote sensing of wide range of variables, from sea level (Minster et al., 1991) to sea surface salinity (Reul & Chapron, 2003) or surface currents Marié et al., 2020).
We found that the gradients of significant wave heights can be quantified in satellite altimeter data in a way that is useful to make a statement on the quality of the ocean currents, in the context of numerical wave MARECHAL AND ARDHUIN 10.1029/2020JC016564 13 of 16 Figure 8. Up figures, two-dimensional significant wave height field snapshot (November 4 00:00 UTC) for (a) the unchanged directional spreading (σ θ ) boundary spectra, (b) the extended σ θ boundary spectra (+30%), and (c) the reduced σ θ boundary spectra (−30%). The solid line is the footprint of one altimeter track for the same date, the significant waves height simulated are displayed in (d) for unmodified (black line), extended (red line) and reduced (purple line) σ θ . (e) The averaged simulated H s spectra over 1 year for the three simulations. modeling. We can imagine that many future developments will further constrain the currents by using (1) more information about the wavefield than just the wave height, (2) measurements over a broader area than the narrow pencil beam of nadir altimeters, and (3) different analyses and techniques. For the first type of future developments, we can mention the use of directional measurements provided by the China France Ocean Satellite (Hauser et al., 2017), launched in October 2018 and the understanding of directional spectral evolution in currents provided by Villas Bôas and Young (2020). For the second aspect, we are expecting a wealth of data, including wave measurements, from the soon-to-be-launched Surface Water Ocean Topography mission (Morrow et al., 2019). As for the third aspect, it can involve the use of different metrics. For example, Villas  showed that the magnitude of the wave height gradient was also related to the slope of the current KE spectrum, which is an interesting quantity for diagnosing the upper ocean dynamics (Le Traon et al., 2008).
Our modification of the boundary conditions is done by a modification of D(f, θ), without changing E(f).
There can be an infinite number of ways to modify D(f, θ). Here, first compute the directional moments a 1 (f), b 1 (f), a 2 (f), and b 2 (f) are computed from D(f, θ) following O' Reilly et al. (1996). These are the discrete Fourier coefficients of the directional distribution D(f, θ).
From the modified parameters, a new directional distribution D′(f, θ) is estimated using the Maximized Entropy Method (Lygre & Krogstad, 1986).