Submesoscale Impacts on Mesoscale Agulhas Dynamics

Mesoscale dynamics of the Agulhas Current system determine the exchange between the Indian and Atlantic oceans, thereby influencing the global overturning circulation. Using a series of ocean model experiments compared to observations, we show that the representation of mesoscale eddies in the Agulhas ring path improves with increasing resolution of submesoscale flows. Simulated submesoscale dynamics are validated with time‐mean horizontal‐wavenumber spectra from satellite sea surface temperature measurements and mesoscale dynamics with spectra from sea surface height. While the Agulhas ring path in a nonsubmesoscale‐resolving (1/20)° configuration is associated with too less power spectral densities on all scales and too steep spectral slopes, the representation of the mesoscale dynamics improves when the diffusion and the dissipation of the model are reduced and some small‐scale features are resolved. Realistic power spectral densities over all scales are achieved when additionally the horizontal resolution is increased to (1/60)° and a larger portion of the submesoscale spectrum is resolved. Results of an eddy detection algorithm applied to the model outputs as well as to a gridded sea surface height satellite product show that in particular strong cyclones are much better represented when submesoscale flows are resolved by the model. The validation of the submesoscale dynamics with sea surface temperature spectra provides guidance for the choice of advection schemes and explicit diffusion and dissipation as well as for further subgrid‐scale parameterizations. For the Agulhas ring path, the use of upstream biased advection schemes without explicit diffusion and dissipation is found to be associated with realistically simulated submesoscales.


Introduction
In the greater Agulhas region, the subtropical supergyre (Speich et al., 2007) strongly exchanges water masses between the Atlantic and the Indian Ocean. This has important implications for the global thermohaline circulation and in particular the Atlantic meridional overturning circulation (Beal et al., 2011;Speich et al., 2007;Weijer et al., 2002). The turbulent upper ocean inflow of Indian Ocean waters into the Atlantic, the "Agulhas leakage," strongly controls the South Atlantic heat and salt budgets and thus its stratification (De Ruijter et al., 1999;Weijer et al., 2002). A part of this water reaches the South Atlantic western boundary and is transported into the equatorial current system where it is an important driver of decadal sea surface temperature (SST) variability (Lübbecke et al., 2015;Rühs et al., 2018). In respect to the far-field influence, the Atlantic multidecadal oscillation is found to correlate with the strength of the Agulhas leakage with a lag of 15 years . This is consistent with the transit times from the Agulhas region into the North Atlantic derived by Rühs et al. (2013).
The Agulhas Current (AC) is the western boundary current of the subtropical gyre in the South Indian Ocean. It flows along the continental slope, separates from the Agulhas Bank south of Africa, retroflects to the east, and forms the Agulhas Return Current (ARC; Lutjeharms & Ansorge, 2001). With a transport rate of 70 Sv (1 Sv = 10 6 m 3 /s) and a core velocity of more than 2 m/s measured at the Agulhas Current Time-series (ACT) section (Beal et al., 2011), the AC is the strongest western boundary current of the Southern Hemisphere. At the retroflection, large rings are shed that are referred to as Agulhas rings. With observed diameters of more than 500 km (Arhan et al., 2003) and available potential energy of up to 71 PJ (Goni et al., 1997), Agulhas rings are among the largest and most energetic ocean eddies. After their shedding, they decay quickly to less than half of their initial size, mixing the bulk of their water mass properties into the surrounding Cape Basin De Ruijter et al., 1999). The Agulhas rings propagate northwestward into the South Atlantic partly due to the effect and partly due to the advection by the Benguela drift (De Ruijter et al., 1999). They transport mass, momentum, heat, and salt from the Indian into the Atlantic Ocean and The region south of 34 • S is a key region for the AC system where small-scale dynamics can have a large impact on the Agulhas leakage and thus on the global ocean. At the northern boundary of the AC, shear-edge instabilities lead to the formation of meanders, shear-edge eddies, filaments, plumes, and submesoscale vorticies (Krug et al., 2017;Lutjeharms et al., 1989;Lutjeharms, Boebel, et al., 2003;Lutjeharms, Penven, et al., 2003). The warm and salty plumes advected from the northern flank of the AC move either onto the Agulhas Bank or drift northwestward into the South Atlantic Ocean as "Agulhas filaments" (Lutjeharms & Cooper, 1996). The filaments transport only 0.1 Sv (Lutjeharms, 2007) but provide about 13% of the total inter-ocean salt flux (Lutjeharms & Cooper, 1996) and thus are a minor but important contribution to the leakage of salt from the Indian Ocean. Meanders and shear-edge eddies grow in size to scales of 50-100 km traveling downstream (Lutjeharms et al., 1989). Shear edge and Natal Pulse related cyclones are observed to be either shed southwestward into the open ocean or to feed their vorticity into a lee eddy on the western side of the Agulhas Bank Lutjeharms, Boebel, et al., 2003;Penven et al., 2001). Lutjeharms, Penven, et al. (2003) further found indication in a model that the shear-edge eddies and other passing cyclonic features are trapped between the eastern side of the Agulhas Bank and the AC, providing a prereservoir of cyclonic vorticity that irregularly leaks into the lee eddy. The lee eddy is observed to occasionally detach from the shelf (Penven et al., 2001) and to travel southwestward . In conclusion, the formation of Agulhas cyclones (shedded shear-edge and Natal Pulse-related cyclones as well as detached lee eddies) and filaments is strongly linked to small-scale activity. Agulhas cyclones interact with anticyclonic Agulhas rings  and sometimes even trigger their shedding . This indicates also a dependence of anticyclonic Agulhas rings and their shedding on small-scale flows. Further cyclones and anticyclones are found to form under strong interaction with smaller-scale flows in the frontal zone of the southern Benguela upwelling system and travel westward where they interact with the Agulhas rings (Lutjeharms & Stockton, 1987;Rubio et al., 2009).
Recent high-resolution model studies found that submesoscale dynamics play an essential role in the mixing of Agulhas ring waters into the Cape Basin and thus in setting the properties of the water masses that are exported into the Atlantic (Capuano, Speich, Carton, & Blanke, 2018;Sinha et al., 2019). On the basis of a high-resolution regional model ((1/108) • ) for a very limited domain, Capuano, Speich, Carton, & Blanke (2018) found that submesoscale dynamics in the Cape Basin drive a strong exchange of both upper and intermediate waters of Atlantic and Indian Ocean origin. They conclude that the simulation of a realistic Indo-Atlantic exchange requires the resolution of a large portion of the submesoscale spectrum. Supporting In this study, we investigate the impact of submesoscale flows on the dynamics of mesoscale eddies in the Agulhas region via model comparison. We use a nested modeling approach to achieve a submesoscale-permitting horizontal resolution: the implementation of two-way interacting grid refinements ("nests") for the greater Agulhas region into global ocean general circulation models with coarser resolution ("hosts") Holton et al., 2017;Schwarzkopf et al., 2019). Several of such configurations have been developed within the INALT family (Schwarzkopf et al., 2019) with horizontal grid refinements of down to (1/10) • and (1/20) • (i.a., INALT20r), as well as a double-nested configuration with (1/60) • horizontal resolution (INALT60). In this study, we compare a nonsubmesoscale-permitting INALT20r simulation to parallel experiments where we gradually increase the vertical resolution, decrease the diffusion and the dissipation of the model, and increase the horizontal resolution to resolve more and more small-scale and submesoscale flows. The model development is lead by a systematic model validation down to the smallest resolved scales and results in a further developed submesoscale-permitting INALT60 that is shown to be in very good agreement with the observations. Commonly, ocean general circulation models are validated by comparing the simulated sea surface height (SSH) anomaly to those of gridded products derived from along-track satellite altimetry as well as by comparing simulated and observed vertical sections through the large-scale currents, such as the ACT vertical section through the Agulhas (Under) Current (Beal et al., 2015). This common validation covers only the simulated large-scale circulation and integrated eddy activity. For a scale-dependent validation of the simulated mesoscale dynamics, we additionally compare the time-mean horizontal-wavenumber spectra of SSH from the model simulations with those from JASON-2 along-track measurements as well as of the SSH spectral slope from the simulations with predictions by quasigeostrophic theory.
Submesoscale flows are associated with Rossby numbers of order 1 (Thomas et al., 2008). The Rossby number is a measure of the influence of Earth's rotation on the flow defined as Ro = | ∕f|. Here = v∕ x− u∕ y is the vertical component of the relative vorticity (u and v are the flow components in zonal [x] and meridional [y] direction) and = 2Ω sin is the planetary vorticity with the Earth's rotation rate Ω = 7.2921×10 −5 rad/s and the latitude . The Rossby number as well as the normalized relative vorticity ∕f can thus be used to quantify and compare the occurrence of submesoscale currents in the model experiments. The validation against observations is, however, a challenging task. In situ observations with a focus on submesoscale dynamics have only been carried out in limited regions of the Agulhas region. Submesoscale eddies with diameters of 10-20 km were observed at the inshore front of the AC with RAFOS floats  and with seagliders (Krug et al., 2017). They were attributed to shear instabilities of the AC jet (Krug et al., 2017). The signal of present-day satellite altimetry exceeds the associated noise at scales larger than about 50 km in the greater Agulhas region (Dufau et al., 2016). Dufau et al. (2016) further point out that in winter time, when the submesoscale dynamics are most active (Callies et al., 2015), the rough sea state shifts this signal-equals-noise scale to even larger scales. Present-day altimetry is thus not able to resolve submesoscale dynamics. Very high resolution (order of tens of meters) satellite observations of reflectance and derivable products such as chlorophyll concentration are hard to interpret in terms of the underlying ocean dynamics and are not comparable to models without biogeochemical components. One of the few remaining possibilities today for an area-covering validation of submesoscale mixed-layer dynamics with satellite products is the analysis of its imprint into the SST field. For this study, we apply spatial spectral analysis to Moderate Resolution Imaging Spectroradiometer (MODIS) data with a horizontal resolution of about 1 km and compare the time-mean spectra to those from the models. The spectral approach applied to both SST and SSH provides a scale-dependent model validation and comparison from the largest simulated scales down to the smallest submesoscales resolvable with a (1/60) • model. Moreover, the validation with SST spectra can justify or guide the choice of explicit diffusion and dissipation, advection schemes as well as the development of new subgrid-scale parameterizations.
The paper is structured in the following ways: In section 2, the model configurations and experiments are described in detail and evaluated with respect to the simulation of submesoscale flows with Rossby numbers of order 1. Subsequently, the simulated ocean dynamics are compared and validated against observations: the submesoscale dynamics with horizontal-wavenumber SST spectra from swath data in section 3, the large-scale circulation and the integrated eddy activity with a gridded SSH product and a vertical section

Model Configurations and the Simulation of Submesoscale Dynamics
The hydrostatic and Boussinesq-approximated primitive equations are solved with the Nucleus for European Modelling of the Ocean Version 3.6 (Madec & the NEMO team, 2014). The Nucleus for European Modelling of the Ocean model configurations used for this study are part of the INALT family (Schwarzkopf et al., 2019). In this section, we give an overview of the used model configurations. For more details, we refer to Schwarzkopf et al. (2019). The INALT configurations are global ocean general circulation models with horizontal grid refinements for the greater Agulhas region and adjacent ocean basins. The global ocean is simulated building on the configuration ORCA025 (Barnier et al., 2006) and consists of an ocean general circulation model coupled to the viscous-plastic sea ice model Louvain-la-Nueve Ice Model Version 2 (Fichefet & Maqueda, 1997). For the discretization, a time step of 20 min and an Arakawa C-grid (Arakawa & Lamb, 1977) are used. Horizontally, the host grid has an eddy-permitting resolution of (1/4) • distributed on a tripolar grid with poles located at the South Pole and over Siberia and Canada to avoid singularities at the geographical North Pole. Vertically, 46 and 120 z levels were used with a partial cell at the bottom (25-m minimum vertical extent) to improve the influence of bathymetry on the dynamics of the ocean (Barnier et al., 2006). The vertical distribution of the z levels is shown in Figure 4, and the data set used for the ORCA025 bathymetry is described in Barnier et al. (2006). ORCA025 uses free-slip lateral momentum boundary condition for the component parallel to the boundary and no-normal flow condition.
For tracer advection, the total variance dissipation (TVD) scheme (Zalesak, 1979) has been used and for momentum advection, the vector invariant (VI) form with energy-and enstrophy-conserving scheme (Arakawa & Hsu, 1990) for the vorticity term. As the original energy-and enstrophy-conserving scheme has been shown to be associated with symmetric instabilities of the computational kind (Ducousso et al., 2017), a Hollingsworth-corrected version was used (Hollingsworth et al., 1983). The numerical diffusion and dissipation associated with both advection schemes are assisted by explicit horizontal diffusion and dissipation terms in the primitive equations. For tracers, an isoneutral Laplacian scheme with a nominal horizontal eddy diffusivity of 300 m 2 /s was used. For momentum, a horizontal bilaplacian scheme with a nominal horizontal eddy viscosity of −1.5 × 10 11 m 4 /s was applied. Vertical diffusion and dissipation are represented with a 1.5-level turbulent kinetic energy scheme (Blanke & Delecluse, 1993). The downward flux of horizontal momentum at the bottom is parameterized as C D u h,btm √ u 2 btm + v 2 btm + with the bottom drag coefficient C D = 0.001, the horizontal velocity vector in the bottom grid cell u h,btm = (u btm , v btm ), and = 0.0025 m 2 /s accounting for unresolved currents due to tides, internal wave breaking, and others. A diffusive scheme with a horizontal mixing coefficient of 1,000 m 2 /s is used for the parameterization of the effect of the bottom boundary layer dynamics on tracers.
For the configuration INALT20r, a horizontal grid refinement of five and a temporal refinement of 3 with respect to the host grid is applied, resulting in a horizontal resolution of (1/20) • (≈4.5 km) and a nest time step of 400 s for the greater Agulhas region (6-50 • S, 20 • W to 70 • E, see Figure 1). The nesting is done with the two-way nesting scheme Adaptive Grid Refinement in FORTRAN (Debreu et al., 2008) enabling an active interaction between both grids and a two-way exchange of signals through the boundary of the nest. In contrast to the global host grid, a no-slip boundary condition is used for the horizontal momentum within the nest. The no-slip condition is the more physical one as the velocities vanish at the resting seafloor. The free-slip condition was nevertheless used for the coarse-resolution host grid as else the simulated boundary current transports are too weak. For the high-resolution nest of a comparable (1/20) • configuration, Schwarzkopf et al. (2019) found only minor changes of the AC and AUC transports when the sidewall boundary condition is switched from free to no slip. Thus, for this study, the no-slip lateral boundary is used for the nests. Eddy diffusivity is linearly scaled down from the host grid values to 60 m 2 /s, and viscosity is quadratically scaled down to −6 ×10 9 m 4 /s. The bathymetry for the nest is generated via interpolation from the ETOPO1 bathymetry (www.earthmodels.org/data-and-tool/topography/etopo). The initial conditions, for the experiments of this study, are taken from a 30-year spin-up of INALT20r that was initialized with climatological temperature and salinity fields from Steele et al. (2001). It was spun up from 1980 to 2009 under atmospheric forcing and bulk formulas developed for the Coordinated Ocean-Ice Reference  Stewart et al. (2017) to derive the requirements for the vertical grid spacing as shown in Figure 4. All land backgrounds that are shown in this paper are taken from Blue Marble Next Generation (Stöckli et al., 2005).
Experiments 2 Large & Yeager, 2009). The forcing variables are prescribed with a 2 • spatial resolution. Atmospheric temperature, wind speed, and humidity are prescribed 6-hourly, radiation daily, and precipitation monthly. A vertical axis with 46 z levels (L46) and a resolution of 6 m near the surface to 250 m in the deep ocean has been used for this integration ( Figure 4). All of the following four experiments start on 1 January 2010 from the spin-up fields and are integrated under JRA55-do(v1.3) forcing (Tsujino et al., 2018) from 2010 to 2017. JRA55-do provides surface forcing in a temporal and horizontal resolution of 3 hr and 0.5 • . Since a large part of the spectrum of submesoscale flows evolve at time scales shorter than a day (McWilliams, 2016), its representation in the simulation is expected to improve with a better resolved daily cycle in the forcing. The fast evolution of submesoscale dynamics is also accompanied by the need for high-frequency model output, which is written as 4-hr means for the surface grid cells. For the horizontal-wavenumber spectra, an additional sea surface snapshot is written every 5 day at noon (UTC) into the model output.

The Nonsubmesoscale-Resolving (1/20) • Experiment (INALT20r.L46.HighDiff)
The first experiment is performed with the described configuration with 46 vertical levels and is called INALT20r.L46.HighDiff to highlight the horizontal and vertical resolution as well as the diffusion and dissipation setup. This reference experiment is strongly eddying but permits almost no submesoscale flows. shows that even in the Southern Hemisphere winter, when the submesoscales should be most active (Callies et al., 2015), the surface flow is mainly associated with low Rossby numbers and small-scale features are almost absent from the simulation (Figure 3a). There are no submesoscale flows simulated in the Agulhas  . The vertical grid spacing as a function of depth for the 46 vertical levels (L46) and the 120 vertical levels (L120). Crosses show the requirements on the vertical grid spacing to resolve the first, second, and third baroclinic mode as derived by Stewart et al. (2017) from CTD measurements at the locations shown in Figure 1c. The right panel shows a zoom into the near-surface box marked in the left panel.
ring path. Only some smaller-scale currents are found in the retroflection and occasionally east of the AC. The time series of Ro > 1 shows only a very weak seasonal cycle ( Figure 2b).

Increasing the Vertical Resolution (INALT20r.L120.HighDiff)
Increasing horizontal resolution provides the potential for resolving smaller-scale currents. To achieve this, it is necessary that the vertical resolution resolves at least the vertical structure of the horizontal flows (Stewart et al., 2017). Under the assumption that baroclinic modal basis functions represent the vertical structure of the horizontal flows, requirements for the vertical resolution to resolve each mode can be derived from CTD measurements by demanding that the model grid at the location of the measurement has to be associated in the vertical with at least three cells between each pair of modal function zero crossings as well as between the surface (bottom) and the nearest zero crossing (Stewart et al., 2017). The respective requirements for each mode derived by Stewart et al. (2017) from World Ocean Circulation Experiment observations are shown in Figure 4 for the measurement within the core Agulhas region (Figure 1c). If the vertical grid spacing at a specific depth is smaller than the requirements, the respective mode is resolved by the vertical grid. Only the first baroclinic mode is resolved with L46 ( Figure 4). With reference to the requirements derived by Stewart et al. (2017), we construct a new vertical grid spacing with 120 vertical levels (L120) that resolves the third baroclinic mode almost everywhere.
Near the surface and in the deep ocean, the vertical grid spacing is chosen independently from the baroclinic mode approach. A further requirement for the vertical grid spacing is that it shall permit the processes in the mixed-layer even in summer when the mixed-layer depth is occasionally less than 20 m in the Agulhas region (de Boyer Montégut et al., 2004). Thus, we placed 10 vertical levels in the upper 20 m with a minimum grid spacing of 1 m at the surface. Capuano, Speich, Carton, & Laxenaire (2018) argue that a uniform vertical grid spacing of about 50 m-in particular at intermediate depths-is needed for the realistic simulation of the exchange of Indian and Atlantic Ocean intermediate waters. As a compromise with respect to computational costs and the near-surface focus of this study, we gradually increase the vertical grid spacing to a maximum of 100 m in the deep ocean. The L120 experiment is named INALT20r.L120.HighDiff. Ro > 1 is found to occur as rarely as with L46 (Figure 2a and 2b) and small-scale features show up again only occasionally in the respective winter time snapshot of near-surface Ro ( Figure 3b). As also the results of almost all other analyses, we present in this paper, are similar for L46 and L120, the results for INALT20r.L46.HighDiff are not discussed explicitly but are shown in the figures throughout the paper. The similar results do not imply that the increase in vertical resolution is not important for the simulation of small-scale flows. We hypothesize that a too high model diffusion and dissipation inhibit the simulation of small-scale flows and thus masks the effect of an increase in vertical resolution.

Reducing the Diffusion and the Dissipation of the Model (INALT20r.L120.LowDiff)
The diffusion and the dissipation of the model are further constraints for the simulation of the small-scale flows that are potentially resolved by the given horizontal and vertical grid spacing. The model diffusion (dissipation) is the combined effect of numerical and explicit diffusion (dissipation) in the primitive equations. The diffusion (dissipation) removes tracer variance (kinetic energy) mainly from the smaller scales. The numerical diffusion (dissipation) depends on the choice of the advection schemes for tracer (momentum).
For the combination TVD/VI, the numerical contribution has to be assisted by additional explicit diffusion and dissipation to avoid the excitation of spurious numerical modes. In the configurations above, the explicit horizontal diffusion (dissipation) is realized as an Laplacian (Bilaplacian) operator, applied to the model solution multiplied by a constant diffusion (dissipation) coefficient. The coefficients have to be large enough to prevent numerical instabilities everywhere in the model domain. Short sensitivity experiments with INALT60 showed that the small scales are still strongly damped with TVD/VI when the coefficients are tuned to a minimum value protecting the simulation from numerical instability (not shown). By contrast, a sensitivity test with a flux-form third-order upstream biased scheme (UBS; Farrow & Stevens, 1995;Madec & the NEMO team, 2014;Webb et al., 1998) for both tracer and momentum without explicit diffusion and without explicit dissipation showed no exertation of numerical modes and a conservation of the small-scale energy down to scales of about 10 km. UBS schemes do not require additional explicit diffusion and dissipation to keep the simulation numerically stable. We perform an experiment where we switch in the nest from TVD/VI with explicit diffusion and dissipation to UBS without explicit diffusion and dissipation. On the host grid we keep the TVD/VI setup as described above. We refer to the experiment as INALT20r.L120.LowDiff. The model time steps had to be decreased to 15 min in the host grid and to 150 s in the nest to run the simulation.
In the winter time snapshot of near-surface Ro, small eddies, sharp fronts, and other smaller-scale currents occur more often compared to the previous configurations-especially in regions where mesoscale activity is strongest (Figure 3c). Small-scale activity is now also found in the Agulhas ring path, mainly between the rings, as well as south of the ARC. Ro > 1 still occurs on average in only 0.27% of the domain area (Figures 2a  and 2b) and shows up only occasionally in the winter time snapshot (Figure 3c), indicating that most of these smaller-scale flows are close to geostrophic balance. The time series of the area coverage of Ro > 1 shows, although weak, the typical seasonal cycle of submesoscale activity (Figure 2b).

Increasing the Horizontal Resolution (INALT60)
Within INALT20r, a secondary nest with a horizontal resolution of (1/60) • (≈1.5 km) is implemented covering the core Agulhas region (0-40 • E, 25 • -45 • S, see Figure 1), together constituing the configuration INALT60. The INALT60 experiment presented by Schwarzkopf et al. (2019) was a first sensitivity study with respect to the horizontal resolution simulating the 1980s. Here, we build on their work and use a further developed INALT60 with L120, UBS advection schemes, JRA55-do forcing, and an integration period from 2010 to 2017. The submesoscale-permitting INALT60 simulation is the main experiment for this study and performed as the fourth parallel simulation. Again, two-way nesting is provided by Adaptive Grid Refinement in FORTRAN and the bathymetry is interpolated from ETOPO1 which has also a resolution of (1/60) • . From INALT60, we write both 4-hr means and 5-day snapshots for the upper 250 m of the secondary nest.
In INALT60, small-scale motions are simulated all over the greater Agulhas region (Figure 2f). Strongest small-scale activity is again found, where the mesoscale eddy activity is strongest. The distribution of ∕f shows the typical skewness toward cyclonic motions at higher Ro (Rudnick, 2001). Ro > 1 occur much more often than in the (1/ For all experiments, the 3-D averaged kinetic energy density (here not shown) within the region of the secondary nest reaches a new stable level after a 2-year secondary spin-up. Thus, we mainly analyze the model period from 2012 to 2017 in this study.

Validating Submesoscale Dynamics with SST Spectra
In this section, we validate the simulated submesoscale dynamics by comparing its time-mean horizontal-wavenumber SST spectrum and its characteristics to those derived from MODIS    twice a day with a horizontal resolution of about 1 km. The MODIS data are classified to be either "best," "good," "questionable," "bad," or "not processed." The latter two categories mainly represent the presence of clouds and are excluded from the analysis. A swath with extraordinary good quality taken on 30 October 2017 of the northeastern Agulhas ring path is partly shown in Figure 5a and compared to an INALT60 snapshot (Figure 5b). Submesoscale features such as fronts, filaments, and small-scale eddies can be seen in both the model solution and the observations. Submesoscale cold-core eddies are found in both data sets, for example, around 15.0 • E, 35.8 • S (marked by the arrows in Figures 5a and 5b). The respective absolute value of the SST gradient | h SST| (Figure 5c and 5d) is a measure of the horizontal-wavenumber SST power spectral density (PSD). It is of comparable strength in the simulation and the observations. The MODIS swath, however, shows features of even smaller scales than in INALT60. A week earlier on 22 October 2017, a swath with very large cloud-free areas was taken. The SST gradients are shown in Figure 5e and compared to INALT60 (Figure 5f). The model SST gradients are too strong in the retroflection region as well as along  Figure 6a). We expect that the model is associated with too large PSD in these regions.
To investigate the spatial distribution of the spectral characteristics, the Agulhas region is subdivided into subregions of 4 • × 4 • (squares in Figure 6a and Figure 5). For both the model and the observations, the data of each subregion are treated separately. First, the noise of the MODIS data is reduced by applying a two-dimensional running mean with a box length of five data points to each swath within the subregion. Subsequently, for each subregion and each swath (model snapshot), the geographical coordinates are transformed into Cartesian coordinates using the European Petroleum Survey Group Geodesy numbers 4326 (World Geodetic System 1984) and 3395 (Cartesian). After that, the data are linearly interpolated onto a regular grid with 1-km grid spacing and scanned for cloud-free lines longer than 200 km in both quasi-zonal  Figure 6. The horizontal scale is identified as half the wavelength . and quasi-meridional directions. The data are not used for the spectral analysis in the respective subregion, if there are less than 100 of such lines in one swath (model snapshot) within the subregion (and thus if more than ≈90% of the subregion is cloud covered). In particular south of 41 • S only about 400 useful swaths are found for the period 2012-2017 for each subregion (Figure 6a). If there are more than 100 cloud-free lines longer than 200 km in the swath within the subregion, in both horizontal directions a one-dimensional horizontal-wavenumber spectrum is computed for each of these lines. For the computation of the each spectrum, mean and trend are subtracted and the data are multiplied by a Hanning window and an associated amplitude correcting factor of 2. To average spectra based on data lines with different length L, the data line is zero padded to a total length of 3,000 km before, which is not exceeded within the analysis. Afterward, the one-dimensional discrete Fourier transform is computed using the python numpy fast Fourier transform algorithm. The PSD is then computed as the squared absolute value of the first half of the Fourier transform normalized by the original length L. Values for wavelengths > L∕2 are removed as they are strongly reduced by the Hanning window. For the model experiments, the spectral analysis is applied to snapshots every 5 days at noon (UTC). Similar to the observations, the spectral analysis is applied to quasi-zonal and quasi-meridional data lines on the regular Cartesian grid, the model output is interpolated on. The spatial-mean spectra within each subregion for both the models and the observations are averaged over the period 2012-2017. To avoid a potential temporal bias in the observations due to a seasonal cycle in cloud cover, monthly averages are calculated before the total time average. We refer to half the horizontal wavelength as the horizontal scale, as in the Fourier analysis the amplitude of quasi-elliptic features, such as eddies, imprint to a wavelength that is about twice as large as the horizontal extent of the features.
Maps of the mean spectral slope for the 20-to 50-km scale band show a comparable pattern for INALT60 and the observations (Figure 6, left): steeper slopes of about k −2.5 are found for the AC, the retroflection, the upwelling region west of Africa, the ARC, and the area east of the AC and shallower slopes close to k −2.1 in the Agulhas ring path. The slopes in the Agulhas ring path are found to be slightly shallower in the model than in the observations, while in the region of the southwest Indian Ocean subgyre slightly steeper slopes are found in the model. The largest differences are found in the southernmost part of the domain south of 41 • S, where the cloud cover is very large (Figure 7a). Maps of the mean PSD in the same scale band show again a similar pattern for INALT60 and the observations (Figure 6, right): large PSD in the upwelling region, the AC and ARC as well as the retroflection and south of the Agulhas ring path. However, the PSD is much larger in INALT60 in the region of the retroflection as well as along the subtropical front, while it is associated with the same amplitude in other regions such as the Agulhas ring path or the subgyre region ( Figure 6, bottom right panel). This is consistent with the comparison of the SST gradient snapshots (Figures 5c-5f). The Agulhas ring path is a region where the spectral characteristics of INALT60 are found to agree very well with the observations. Consistently, the time-mean spectrum of INALT60, derived for a 9 • × 8 • region in the Agulhas ring path with the same procedure as described above for 4 • × 4 • subregions, is almost identical to the one from MODIS (Figure 7a). Both spectra follow an almost straight line down to scales of 15 km with a slope of k −2 , where k is the wavenumber. At 8 km, the spectrum of INALT60 drops strongly due to the dominant effect of the diffusion and the dissipation of the model. Below these scale, the spectrum from MODIS is associated with more PSD than INALT60 consistent with that smaller features are . Orange contours show cross section "northeastward" velocities with a contour interval of 1 cm/s, black and violet contours mark "southwestward" velocities with a contour interval of 5 and 10 cm/s. The 5-cm/s isotach from the observations in the lower right panel is also shown in the plots for the models with a cyan line. covered by the observations. In INALT20r.L120.HighDiff almost no small-scale flows are simulated in the Agulhas ring path. Consistently, its spectrum is below 100 km associated with a too steep slope and thus too less PSD. When the diffusion and dissipation of the model are reduced (INALT20r.L120.LowDiff), the spectrum is found to be in good agreement with the observations down to 25 km. This scale can further be identified as the scale, where the SST spectra of the (1/20) • simulations drop due to the dominant effect of the model diffusion and dissipation. INALT20r.L120.LowDiff 70 ± 25 6.6 ± 5.5 7 1± 29 6.6 ± 5.4 INALT60.L120.LowDiff 73 ± 26 6.5 ± 5.2 6 7± 27 6.5 ± 5.8 Observations (Beal et al., 2015) 78 ± 31 7.0 ± 5.1 -- The region of the retroflection is found to be a region of steeper spectral slopes (≈ k −2.5 ) and larger PSD.

Differences in Large-Scale Dynamics and Integrated Eddy Activity
For all four model experiments, the pattern of the simulated large-scale near-surface AC system compares well with observations. Simulated SSH anomalies, averaged from daily data for the period 2012-2017, show a very good agreement with the satellite-derived daily L4 product provided by the Copernicus Marine and Environment Monitoring Service CMEMS (Figure 8, left). Location and strength of the SSH gradient associated with the AC, its separation and the ARC compare well with the observations. Moreover, the meanders of the time-mean ARC are more or less at the observed locations and the Agulhas ring path, identifiable as a region of enhanced SSH anomaly, has the observed extent. While in INALT20r.L120.LowDiff and INALT60, the amplitude of the SSH anomalies in the Agulhas ring path is closer to the observations, the time-mean retroflection extends in both experiments too far south.
The temporal standard deviation of the SSH is a measure of the integrated geostrophic eddy activity. Its overall patterns are generally in good agreement with the observations, as expected from the good agreement of the time-mean SSH (Figure 8, middle). Deviations from the observations are found in the retroflection where the (1/20) • experiments show too much eddy activity, while the amplitude in INALT60 compares well. In the Agulhas ring path, the SSH variability of INALT60 is strongest and closest to the observations. Also very regional features, such as the observed small region of enhanced eddy activity south of the Agulhas Bank around 35 • S, 24 • E, that is driven by instabilities of the inshore front of the AC (Lutjeharms, Penven, et al., 2003), is captured very well in size and amplitude in INALT60. North of the ARC, all experiments slightly underestimate the observed SSH variability. The excess southward extent of the retroflection, in particular in INALT60, is consistently associated with more eddy activity compared to the observations. Due to the short integration period it is unclear whether this and other regional differences are systematic or only due to the strong nonlinearities in the AC system. Further, the along-track satellite measurements are interpolated both in space and time for the daily gridded product, which is thus only an approximation for the true daily means. Both arguments highlight the limitations of the validation of the eddy field with the standard deviation of gridded SSH and the need for a statistically more robust approach applied directly to along-track data (see section 5).
The comparison of the time-mean cross-section velocity of the AC at the ACT section with observations by Beal et al. (2015) for the period from May 2010 to February 2013 shows for all experiments a very good  Figure 10. For the greater Agulhas region (a), the spectra were first averaged in time and for each subregion (rectangulars in Figure 10) before a spatial average over all subregions. The spectra for the model experiments have been computed from snapshots each fifth day and for the observations from the measured tracks each tenth day. For the observations, the original spectrum is shown, as well as the one where the mean noise of the 15.0-to 37.5-km scale band is removed from the whole spectrum. Thin straight lines show k −5 and k −4 slopes, where k is the wavenumber. Dark gray lines mark the 79-to 150-km scale band for the computation of the spectral slope for Figure 10. The horizontal scale is identified as half the wavelength . For the Agulhas ring path (b), all spectra for the measurement tracks that range over more than one subregion within the black box in Figure 10 are averaged. The y axis is the same as in (a). A zoom with linearly scaled axes is shown in the subwindow for the marked box.
agreement with respect to the strength and the vertical extent of the AC (Figure 8, right). The 30-cm/s isotach extends down to about 1,000 m, and the 5-cm/s isotach down to about 1,800 m. The horizontal extent of the AC is not comparable for such a short time period, as it strongly depends on the nonlinear recirculation from the ARC into the AC . The time-mean AC transports in the simulations have values between 70 and 84 Sv in the respective time period and compare well with the observed value of 78 Sv (Table 1). The AC transport is here defined as the "southwestward" barotropic cross-section transport integrated vertically and from the coast to the distance where the time-mean barotropic cross-section transport changes its sign (Beal et al., 2015). Due to the nonlinearities of the Agulhas system, the transports of the model experiments and the observations are not correlated in time (not shown). INALT20r.L120.LowDiff and INALT60 are associated with temporal standard deviations of the AC transport of about 25 Sv, which is higher than values of about 20 Sv for the high-diffusion/dissipation experiments but smaller than 31 Sv derived from the observations ( Table 1). The AC transport in INALT60 is closest to the observations.
For the representation of the AUC below the AC, a gradual improvement is found from INALT20r.L120.HighDiff to INALT60. Largest improvements are found when the diffusion and the dissipation are reduced while only a minor change is found when the horizontal resolution is increased. INALT20r.L120.HighDiff is associated with a weak time-mean AUC of less than 4 cm/s that is only found in the first 150 km from the coast (Figure 8c). With reduced diffusion and dissipation, near the bottom several AUC cores are found. Their extent from the coast of 200 km is in agreement with the observations. The simulated multicore structure does not have to be unrealistic, owing to the interpolation of the observations. The main core at the slope is associated with time-mean velocities of more than 5 cm/s, which are below the observed 9 cm/s. The reduction of the diffusion and the dissipation is nevertheless associated with about a doubling of the AUC transport averaged over the observed period as well as its standard deviation from 3.6 ± 2.6 to 6.6 ± 5.5 Sv which is close to the observed 7.0 ± 5.1 Sv ( Table 1). The AUC transport is here defined as the "northeastward" cross-section transport below 2,000-m depth integrated vertically and from the coast to a distance of 200 km. In INALT60, both the horizontal extent and the maximum time-mean velocities of more than 8 cm/s (Figure 8i), as well as an AUC transport of 6.5 ± 5.2 Sv (Table 1) are very close to the observations.  Figure 9b) and for the eddy detection ( Figure 11).

SSH Spectra
In this section, we validate the time-mean horizontal-wavenumber spectra of SSH from the model experiments against those from Native GDR-D Jason-2 L3 along-track SSH anomalies from 2012-2015. The Ssalto/Duacs altimeter products were produced and distributed by CMEMS (http://www.marine. copernicus.eu). The SSH anomalies are available every 10 days with a resolution of about 5.1 km and are referenced to the mean SSH product CLS01 (Hernandez & Schaeffer, 2001). We excluded the data from 2016 and later from the analysis as the satellite orbit has been changed in October 2016. The location of the observed tracks is shown in Figure 10. To investigate the regional distribution of the spectral slope, we divide the Agulhas region into 4 × 13 subregions oriented at the crossings of the measurement tracks ( Figure 10). The geographical coordinates of the data from observations and models are transformed into Cartesian coordinates, and the data are interpolated onto a regular 1-km spacing. One spectrum is computed for every measured track within each subregion. In the rare case of data gaps, for each segment longer 200 km a spectrum is computed and afterward all segment spectra are averaged. The time-mean spectrum for each subregion is derived by averaging the time-mean spectra of both diagonals. Averaging over all subregions gives one time-mean SSH anomaly spectrum for the Agulhas region (Figure 9a).
The spectrum from JASON-2 for the Agulhas region decreases toward smaller scales exponentially with a slope close to k −5 , where k is the wavenumber (Figure 9a). At scales of 90 km, the spectrum turns to almost white noise due to the instrument noise. Xu and Fu (2012) suggested to identify the instrument noise as the mean PSD of the 25-to 35-km wavelength band (12.5-to 17.5-km scale band) and to remove that from the whole spectrum. However, at the small scale end below 13 km, a drop in the spectrum is found-probably due to the interpolation technique-that would lead to a too low estimate of the mean noise. Thus, we remove the mean 15.0-to 37.5-km noise from the entire spectrum, where it is associated with almost white noise. The corrected spectrum drops with ≈ k −5 down to scales of 75 km and changes to a k −4 slope below.
For the models, we interpolate the 5-day snapshots of SSH anomaly of the period 2012-2015 onto the mean track of the satellite data and proceed similarly as with the measured data. The SSH anomalies are referenced to the 2012-2015 mean SSH. The mean track is identified as the linear regression through all measured data points on the respective diagonal in each subregion.
The SSH spectra of the high-diffusion/dissipation experiments show too less PSD on all scales and decrease with an almost constant slope down to scales of about 30 km (Figure 9). At smaller scales, they drop due to the reduction of tracer variance (kinetic energy) by the model diffusion (dissipation).
INALT20r.L120.LowDiff is associated with more PSD in the small-scale range consistent with more simulated small-scale features as demonstrated in sections 2 and 3. The spectrum from INALT60 shows very similar slopes compared to the observations above scales of 75 km (k −5 ) as well as below (k −4 ). The PSD is at all scales much higher than in INALT20r.L120.LowDiff and very close to the observations. The remaining small lack of PSD might be related to the absence of tidal induced motions in INALT60.
The map of the mean spectral slope of the 79-to 150-km scale band reveals for the observations that regions of large mesoscale activity such as the Agulhas ring path, the retroflection, and AC/ARC and adjacent regions are associated with steeper slopes close to k −5 . Open ocean regions such as southwest of the Agulhas ring path and north and south of the ARC east of 36 • E show shallower slopes close to k −4 (Figure 10e). Away from the Agulhas ring path, all model configurations are overall in good agreement with the observations with respect to the spectral slope. For the open ocean k −4 slopes are found and for the retroflection, the ARC and adjacent regions k −5 slopes ( Figure 10). In the Agulhas ring path, too steep slopes of around k −6 are found for the high-diffusion/dissipation configurations (Figure 10a and 10b). Reducing the diffusion and the dissipation leads to shallower spectral slopes in the Agulhas ring path that are closer to the observations but for some subregions still steeper than k −6 . In INALT60, the spectral slopes in the Agulhas ring path are found to be close to k −5 and compare very well with those from the observations.
To identify the impact of the submesoscale dynamics on the large mesoscales in the Agulhas ring path and to validate their simulation, we extend the length of the data lines and perform the computation of the spectra for a large region in the Agulhas ring path (26.88-41.17 • S, 1.42-12.74 • E, Figure 10). Within this region, a spectrum is computed for each measured track that ranges over more than one of the small subregions identified before. The spectral analysis is performed analogous to the method described above. Subsequently, all spectra are averaged to one time-mean SSH spectrum for the Agulhas ring path and presented in Figure 9b. For scales smaller than 200 km, the results are similar to those of the whole Agulhas region. The shallowing of the spectral slope, when the diffusion and the dissipation of the model are reduced, is more pronounced in the Agulhas ring path consistent with the maps of the spectral slope. For scales larger than 200 km, the PSD for the high-diffusion/dissipation experiments is found to be almost only half of the observed. This indicates a too weak mesoscale eddy field. Reducing the diffusion and the dissipation of the model leads to a small improvement in PSD. If also the horizontal resolution is increased, the PSD shows very good agreement with the observations.

Agulhas Eddy Detection
In the Agulhas ring path, an increase of the mesoscale SSH spectral density, as shown in the previous section, can mainly be attributed to either an increase in the number of eddies or an increase in the strength of the eddies, or both. To specify this, the SSH-based automated eddy detection algorithm presented by Chelton et al. (2011) is applied to the daily CMEMS gridded L4 SSH product as well as to the daily model outputs interpolated onto the CMEMS grid for the same region in the Agulhas ring path (black box in Figure 10 and the period 2012-2017). The data sets are temporally detrended and the domain average is removed for each time step before the eddy detection. Anticyclones (cyclones) are detected as connected sets of minimum 8 to maximum 1,000 pixels with SSH anomalies above (below) a given threshold, a local maximum (minimum) with a maximum amplitude larger than 1 cm and a maximum distance of any pair of pixels of 700 km. The connected pixel set is thereby divided into interior and perimeter pixel, where the latter are associated with boundaries to the surrounding of the pixel set. The maximum amplitude A max of an eddy is then computed as the SSH difference of the value of the local interior extrema and the mean value of the perimeter pixels. For each time step, the SSH fields are scanned for connected pixel sets that fulfill the above-mentioned conditions for a threshold of −1 m (1 m) to 1 m (−1 m) with an increment of 0.01 m for anticyclones (cyclones). If a connected region is identified to be an eddy, its interior pixels are removed from the further detection. The scale of an eddy is identified to be the diameter of a circle with the same area as the detected pixel set. A single eddy is detected at each daily mean model output time step as long as it remains in the investigated domain. In the subsequent analysis we do not investigate single eddies but analyze the eddy characteristics in a statistically average sense. We present histograms of the average number of eddies in dependence on their amplitude, scale, and direction of rotation for the Agulhas ring path region and the period 2012-2017. For a specific amplitude or scale range, the number of detected eddies with the respective properties is therefore summed up for all time steps and then divided by the total number of time steps.
The average eddy number does not increase substantially throughout the model experiments and is in good agreement with the observations. On average, about nine anticyclones as well as nine cyclones are detected in the domain for all experiments as well as for the observations. Only a small increase in average cyclone number from 8.4 to 9.2 is found when the diffusion and the dissipation of the model are reduced and a very small increase in anticyclone number from 8.7 to 9.1 when the horizontal resolution is increased. Both values are very close to those of the observations (9.3 and 9.2). Moreover, the distribution of the scale of all detected eddies is found to be very similar for all models and observations (Figure 11b and 11f).
Mesoscale eddies and in particular cyclones are found to strengthen when the diffusion and the dissipation are reduced and especially when the horizontal resolution is increased. Weak anticyclones (cyclones) with amplitudes between 1 and 7 cm occur on average more often in the INALT20r configurations (in sum ≈ 4.1 (≈5.5)) than in INALT60 (3.7 (3.8)) and the observations (3.3 (4.1)) ( Figure 11a and 11d). The amount of medium amplitude anticyclones (cyclones) with amplitudes between 7 and 20 cm is too small in the high-diffusion/dissipation experiments (≈ 2.3 (≈ 2.6)), larger in INALT20r.L120.LowDiff (2.8 (3.6)), and similar in INALT60 (3.3 (4.5)) and the observations (3.7 (4.2)). The increase in medium amplitude eddy number is restricted to scales below about 300 km (not shown). Two strong amplitude anticylones with amplitudes larger than 20 cm are detected on average in all simulations and in the observations. The average distribution of their scale is very similar (Figure 11c). About one strong cyclone is detected on average in INATL60 and the observations with a similar scale distribution (Figure 11g). In contrast to the anticyclones, on almost all scales too less strong cyclones are detected on average in INALT20r.L120.LowDiff (0.5) and in particular in the high-diffusion/dissipation experiments (≈0.3). Consistently, the maximum amplitude of

Conclusion and Discussion
In this study, we compare a series of ocean model experiments to observations to show that the representation of mesoscale eddies in the Agulhas ring path improves with an increasing resolution of submesoscale flows. Commonly, ocean models are validated by comparing the time-mean and standard deviation of simulated SSH with gridded satellite altimeter products as well as vertical sections through important currents with respective measurements. Here, we extend the common validation of the mesoscale dynamics by comparing time-mean horizontal-wavenumber spectra of simulated SSH with those from JASON-2 along-track altimetry as well as the results of an SSH-based eddy detection algorithm applied to the model outputs as well as to a gridded SSH satellite product. The simulated small-scale mixed-layer dynamics are validated with SST horizontal-wavenumber spectra from the model experiments and MODIS swath measurements.
We start with the nonsubmesoscale-resolving reference configuration INALT20r.L46.HighDiff with (1/20) • horizontal resolution, 46 vertical levels, and relatively high model diffusion and model dissipation. The validation shows that, in the Agulhas ring path, the simulation is associated with too weak PSDs on all scales, too steep spectral slopes, and too weak mesoscale eddies. In particular, strong cyclones are not represented, while strong anticyclones are represented well. We hypothesize that the latter is due to a good representation of the Agulhas ring shedding in the simulation. Three parallel simulations are performed where we gradually increase the vertical resolution to 120 vertical levels (INALT20r.L120.HighDiff), decrease the diffusion and the dissipation of the model (INALT20r.L120.LowDiff), and increase the horizontal resolution to (1/60) • for the core Agulhas region (INALT60). An increase of only the vertical resolution does not lead to a substantial improvement in small-scale, mesoscale, and large-scale dynamics. An additional reduction of the diffusion and the dissipation of the model leads to an improvement of the AC variability as well as AUC variability and strength to a good agreement with the observations by Beal et al. (2015). In the Agulhas ring path, the reduced diffusion and dissipation lead to the simulation of some small-scale flows, larger PSD on all scales, shallower spectral slopes, and a better representation of cyclone strength and number. A small increase in cyclone number is found that might be attributable to the improved representation of the AC/AUC variability. Most of the mentioned results nevertheless show still a large gap to the observations. This gap is closed with the implementation of the secondary (1/60) • nest. In INALT60, much more and much stronger submesoscale flows occur all over the domain and the mesoscale dynamics are well represented with respect to SSH variability, spectral slope, PSD, and eddy number and strength. It has to be clarified why substantial changes in SSH spectral slope are found only for the Agulhas ring path, while there is a good agreement with the observation throughout the model experiments away from the Agulhas ring path, consistent with similar temporal standard deviations of SSH. A limitation of the eddy detection method is the cutoff at the boundaries of the investigated region. The scale of detected eddies that continue outside the domain is underestimated. However, this affects models and observations equally so that it does not bias the comparison.
Comparing the time-mean SST horizontal-wavenumber spectra from the simulations and from MODIS satellite measurements is one of a few possible methods available today for a systematic validation of the mixed-layer dynamics and further for a justification of and guidance for the choice of advection schemes and explicit diffusion and dissipation. To our knowledge, this is the first study applying this approach. The 20-to 50-km mean spectral slope shows a similar spatial pattern for INALT60 and the observations with shallower slopes of about k −2.1 in the Agulhas ring path and steeper slope of about k −2.5 in the regions of the AC, the retroflection, the ARC, and the subgyre. Further, the PSD agrees very well for the Agulhas ring path and the subgyre region east of the AC. In the Agulhas ring path, the spectrum for INALT60 (INALT20r.L120.LowDiff) is in very good agreement with the observations down to 8 km (25 km). This justifies the use of upstream biased advection schemes without explicit diffusion and dissipation for this region. Below 8 km (25 km), the spectra drop due to the effect of the model diffusion and dissipation. The respective wavelengths are close to 10 times the grid spacing in accordance with the effective resolution identified by Soufflet et al. (2016). In regions of strong large-scale temperature gradients as for the region of the retroflection and along the subtropical front, INALT60 is associated with too large PSD for all comparable scales. The associated lack of SST damping across all scales might be attributable to nonrepresented atmosphere-ocean interaction in ocean-only models that are most important in regions of strong SST gradients. Comparing an example wintertime snapshot with large cloud-free areas with the respective INALT60 snapshot shows too strong temperature gradients in the model. This gives some evidence that the offset is not an artifact of the validation method.
A major limitation of the validation via infrared derived SST is the nonavailability of satellite data for cloudy regions. Largest offsets are found in regions with large cloud cover where the lowest amount of data is available. Systematic changes of the absolute SST for cloudy and cloud-free regions do not affect the spectral analysis. However, a dependence of submesoscale turbulence on clouds as well as a bias in the orientation of the cloud-free lines to either zonal or meridional direction can lead to biases in the validation. A further bias can arise from uncertainties in the estimation of the bulk SST from the skin SST during the postprocessing of the MODIS products (Donlon et al., 2002;Kilpatrick et al., 2015;Schluessel et al., 1990) used here. A part of the uncertainties might be clarified by a spectral comparison to 4-km × 8-km resolution satellite-based radar-derived surface velocities (Rouault et al., 2010) or new high-resolution satellite altimeter products from the SWOT mission that will start in 2020 and are predicted to cover scales down to 10 km in the Agulhas region (Dufau et al., 2016).
In contrast to SSH, the model SST spectra show similar PSD at the larger scales (>100 km). We attribute this to the fact that the SST variability is to some extent prescribed at the larger scales by the atmospheric forcing which is associated with a coarse horizontal resolution of 0.5 • (≈50 km). The results are, however, not contrary to the ones for SSH as SST amplitudes are no measure of the geostrophic flow. The Agulhas region is associated with very strong surface heat fluxes that are in the retroflection directed toward the atmosphere throughout the year (Walker & Mey, 1988;Yu & Weller, 2007). These heat fluxes lead to a fast reduction in SST gradients, for example of recently shed Agulhas rings (Olson et al., 1992) compared to their surrounding.
Simulated mesoscale dynamics are validated and compared with time-mean SSH anomaly horizontal-wavenumber spectra against JASON-2 along-track measurements. The spectral SSH slopes can further be compared to theoretical predictions. The balanced component of the flow can be decomposed into components due to interior potential vorticity anomalies and surface density anomalies (Lapeyre & Klein, 2006). The two extrema, which one of both components is neglectable, predict spectral slopes of k −5 for interior-quasigeostrophic (IQG) dynamics (Charney, 1971) and k −11/3 for surface-quasigeostrophic (SQG) dynamics (Blumen, 1978;Held et al., 1995), where k is the wavenumber. In SQG theory, advection by the ageostrophic flow is neglected. If included, it leads to a theoretical k −4 slope (Boyd, 1992;Xu & Fu, 2011). Semigeostrophic theory predicts a k −8/3 slope for the kinetic energy spectrum (Andrews & Hoskins, 1978) that is close to the IQG prediction of k −3 . For the greater Agulhas region, the SSH spectral slope has been derived from satellite altimetry to be close to the SQG predicted value of k −3.7 for the 100-to 300-km wavelength band (scales of 50-150 km) by Le Traon et al. (2008) and for the 70-to 250-km wavelength band (scales of 35-125 km) by Xu and Fu (2011). In the core Agulhas region (the region with (1/60) • resolution in INALT60), Xu and Fu (2011) found steeper slopes of k −4 consistent with SQG theory including ageostrophic advection. These estimates, however, did not remove the measurement noise from the spectra before the computation of the spectral slope. Xu and Fu (2012) showed that the spectral slope estimates for the Agulhas region are associated with values of k −4 − k −4.5 , when the mean white noise of the 25-to 35-km wavelength band is subtracted from the whole spectra before that. Regions of strong mesoscale activity are found to be associated with steeper spectral slopes than quieter regions (Chassignet & Xu, 2017;Dufau et al., 2016;Richman et al., 2012;Sasaki & Klein, 2012;Xu & Fu, 2012). In this study, we find slopes from JASON-2 along-track observations that are very close to SQG prediction including ageostrophic advection in regions of low mesoscale activity (southwest of the Agulhas ring path as well as north and south of the ARC east of 36 • E) and close to the IQG and semigeostrophic prediction in more energetic regions (everywhere else in the Agulhas region) for scales that are larger than about 75 km. A scale of 75 km corresponds to a wavelength of 150 km, which is close to the wavelength 2 R d of the Rossby radius of deformation R d that is about 25 km in the Agulhas region (Chelton et al., 1998). For smaller scales, just above the noise floor, shallower slopes are found that follow the SQG including ageostrophic advection prediction. Unbalanced motions might also contribute to the shallowing of the slope (Callies & Ferrari, 2013;Qiu et al., 2018). However, a separation of the SSH spectrum into a balanced and an unbalanced contribution, as applied by Qiu et al. (2018) to a global (1/48) • horizontal resolution model, showed that balanced motions in the core Agulhas region dominate the SSH spectrum down to scales of 25 km or even smaller. Further research is needed to identify, which theoretical model is appropriate to describe the underlying dynamics. The good agreement of the spectral slope from the observations and INALT60 gives some confidence that the shallowing of the spectral slope at scales of about 75 km we found for the observations is not an artifact of the instrument noise. Our results thus suggest that the spectral slope changes within the scale band that Xu and Fu (2012) used for the computation of the slope. This can explain why they estimated shallower slopes than we find here. The steep slopes of k −5 are only identified when they are averaged over a scale band that is located above the wavelength of the Rossby radius.
Idealized gyre simulations with integration periods of 100 years by Lévy et al. (2010) showed that the simulated gyre circulation strengthened when the horizontal resolution is increased from (1/9) • to (1/54) • . Assuming that their results are adaptable to the Indo-Atlantic supergyre, they indicate that there might be a submesoscale impact on the strength of the AC, and thus on the Agulhas ring shedding, that is not captured in our simulations. A much longer simulation, resolving the submesoscale dynamics for the full Indo-Atlantic supergyre, would be needed to address this.
In conclusion, the increasing resolution of submesoscale flows is found to be accompanied by a better representation of the mesoscale eddies in the simulated Agulhas ring path. The eddies are found to strengthen, the more submesoscale flows are resolved by the model. The model solution converges to the observations in our (1/60) • ocean simulation INALT60. Our results indicate that submesoscale dynamics drive a strong inverse kinetic energy cascade in the Agulhas ring path. This will be investigated in a follow-on study.