The impact of horizontal resolution (1/12° to 1/50°; 6 to 1.5 km at midlatitudes) on Gulf Stream separation, penetration, and variability is quantified in a series of identical North Atlantic experiments. The questions the authors seek to address are twofold: 1) Is the realism of the modeled solution increased as resolution is increased? 2) How robust is the modeled mesoscale and submesoscale eddy activity as a function of grid spacing and how representative is it of interior quasigeostrophic (QG) or surface quasigeostrophic (SQG) turbulence? This study shows that (i) the representation of Gulf Stream penetration and associated recirculating gyres shifts from unrealistic to realistic when the resolution is increased to 1/50° and when the nonlinear effects of the submesoscale eddies intensifies the midlatitude jet and increases its penetration eastward, (ii) the penetration into the deep ocean drastically increases with resolution and closely resembles the observations, and (iii) surface power spectra in the 70–250-km mesoscale range are independent of the horizontal resolution and of the latitude and are representative of 2D QG and SQG turbulence.
1 Convergence studies are unusual because parameters are often changed as resolution is increased. Here, we follow in the footsteps of Hurlburt and Hogan (2000), Smith et al. (2000), Oschlies (2002), Bryan et al. (2007), Lévy et al. (2010), Thoppil et al. (2011), Marzocchi et al. (2015), and Biri et al. (2016) by reporting on the impact of horizontal resolution (1/12° to 1/50°) on Gulf Stream separation, penetration, and variability using a series of identical North Atlantic experiments. The questions we seek to address are twofold: 1) Is the realism of the modeled solution increased as resolution is increased? 2) How robust is the modeled mesoscale and submesoscale eddy activity as a function of grid spacing and how representative is it of interior quasigeostrophic (QG) or surface quasigeostrophic (SQG) turbulence?
It is generally recognized that a minimum resolution of 1/10° is required for a proper representation of midlatitudes’ western boundary currents and associated eddies (Paiva et al. 1999; Smith et al. 2000; Maltrud and McClean 2005; Chassignet and Marshall 2008). A grid spacing of 1/10°, however, is not sufficient to resolve the Rossby radius of deformation with two grid points at all latitudes (Hallberg 2013) and therefore does not allow for a proper representation of baroclinic instability and associated eddies throughout the domain. Furthermore, effective resolution is limited by the numerical dissipation range (Soufflet et al. 2016), and models with resolution on the order of 1/10° are now being referred as eddying models, while eddy-resolving models are configurations that truly resolve the first Rossby radius of deformation throughout the domain (CLIVAR 2014). Despite the major improvements observed with 1/10° grid spacing, the solutions remain extremely sensitive to choices in boundary conditions and subgrid-scale parameterizations (Ezer and Mellor 2000; Chassignet and Garraffo 2001; Bryan et al. 2007; Chassignet and Marshall 2008), and there is a continuous need to quantify the value added of increased resolution. Furthermore, horizontal resolution of the O(1) km is required to explicitly resolve submesoscale motions with a horizontal scale of the O(10) km. Submesoscale physics have been shown to play a significant role in the vertical fluxes of mass, buoyancy, and tracers (Thomas et al. 2008; Capet et al. 2008; Fox-Kemper et al. 2008; Klein et al. 2011; Roullet et al. 2012; Capet et al. 2016). However, only a few studies have been able to report on the impact of the submesoscale motions on the large-scale oceanic circulation because of the computing cost associated with numerical simulations with O(1) km grid spacing. Using a hydrodynamic (no active thermodynamics) model, Hurlburt and Hogan (2000) showed a significant improvement in western boundary current pathways with increased resolution (1/16°, 1/32°, and 1/64°, respectively). Lévy et al. (2010) compared the mean characteristics of an idealized, flat bottom, basin-scale seasonally subtropical and subpolar gyres configuration in a suite of numerical experiments varying in horizontal resolution (1°, 1/9°, and 1/54°). They found that the nonlinear effects of the submesoscale eddies that emerge at 1/54° strongly intensify the jet that separates the two gyres, making it more zonal and penetrating farther to the east. The authors state that their results are presumably highly constrained by the idealized geometry of their domain, but we find that many of their results carry over to our series of identical North Atlantic experiments (1/12°, 1/25°, and 1/50°, respectively).
In this paper, we show that (i) the representation of Gulf Stream penetration and associated recirculating gyres shifts from unrealistic to realistic when the resolution is increased to 1/50°, (ii) the penetration of EKE into the deep ocean is drastically different and closely resembles observations, and (iii) surface power spectra in the 70–250-km mesoscale range are independent of the horizontal resolution and of the latitude and are representative of 2D QG and SQG turbulence (k−5 SSH spectra in energetic regions, where k is the wavenumber, not as steep in quiescent regions). The paper is organized as follows: Section 2 describes the model configuration and spinup procedure. The surface mean and turbulent circulation is discussed in section 3 as a function of resolution and in comparison to observations. It also focuses on the impact of ageostrophic motions on the representation of surface eddy kinetic energy (EKE) and the fact that current altimetry measurements likely underestimate surface eddy kinetic energy by as much as 30%. The impact of resolution on the deep-ocean circulation is shown to be quite significant in section 4. Power spectra are used in section 5 to quantify the degree of 2D QG and SQG turbulence present in the numerical simulations. Finally, the results are summarized in section 6.
2. Model configuration and spinup
The Hybrid Coordinate Ocean Model (HYCOM) configuration used in this paper is identical to that of Xu et al. (2010, 2012, 2013, 2014, 2015, 2016) and covers the North Atlantic from 28°S to 80°N (Fig. 1). The vertical coordinate in HYCOM (Bleck 2002) is isopycnal in the stratified open ocean and makes a dynamically smooth and time-dependent transition to terrain following in shallow coastal regions and to fixed pressure levels in the surface mixed layer and/or unstratified seas (Chassignet et al. 2003, 2006). No inflow or outflow is prescribed at the northern and southern boundaries. Within a buffer zone of about 3° from the northern and southern boundaries, the 3D model temperature, salinity, and depth of isopycnal interface are restored to the monthly Generalized Digital Environmental Model (GDEM; Teague et al. 1990; Carnes 2009) climatology with an e-folding time of 5–60 days that increases with distance from the boundary.
The horizontal resolution for the three experiments are 1/12°, 1/25°, and 1/50° (9, 4.5, and 2.25 km at the equator and 6, 3, and 1.5 km in the Gulf Stream region, respectively). The 1/12° model topography is based on the 2′ Naval Research Laboratory (NRL) digital bathymetry database, which combines the global topography based on satellite altimetry of Smith and Sandwell (1997) with several high-resolution regional databases. The 1/25° and 1/50° topographies are linearly interpolated from the 1/12° topography and do not contain additional high-resolution topographic features. In the vertical, the simulation contains 32 hybrid layers with density referenced to 2000 m (σ2): 28.10, 28.90, 29.70, 30.50, 30.95, 31.50, 32.05, 32.60, 33.15, 33.70, 34.25, 34.75, 35.15, 35.50, 35.80, 36.04, 36.20, 36.38, 36.52, 36.62, 36.70, 36.77, 36.83, 36.89, 36.97, 37.02, 37.06, 37.09, 37.11, 37.13, 37.15, and 37.20 kg m−3. The 1/12° configuration resolves the first Rossby radius of deformation up to 60°N, while the 1/50° resolves it everywhere (Hallberg 2013). Furthermore, the 1.5-km grid spacing of the 1/50° configuration resolves up to the fifth internal Rossby radii of deformation at midlatitudes based on the above density distribution of the hybrid layers. As in Chassignet and Garraffo (2001) and Xu et al. (2010), the horizontal viscosity operator is a combination of Laplacian [A2 = max(0.05Δx2 times the deformation tensor A)] and biharmonic (A4 = V4Δx3), where Δx is the grid spacing, A is a constant coefficient, and V4 is the biharmonic diffusive velocity. The viscosity and diffusion parameters are listed in Table 1. The values for the coefficients in the 1/25° decrease as a function of the grid size and are half that of the 1/12°. The values for the coefficients in the 1/50° are kept close to that of the 1/25° in order to isolate the impact of resolving the submesoscale on the solution. The K-profile parameterization of Large et al. (1994) is used for vertical mixing in the surface mixed layer as well as in the ocean interior. The bottom drag is quadratic with a coefficient of 10−3 and a background RMS flow speed of 5 × 10−2 m s−1.
Sea ice–related processes are modeled using an “energy loan” where freezing takes place whenever latent heat is needed to keep the mixed layer temperature from dropping below the freezing level. When the ocean–ice system is being heated, the incoming energy is used to melt the ice before the water temperature is allowed to rise above the freezing level (Semtner 1976).
The three configurations are initialized using potential temperature and salinity from the GDEM climatology and were spun up from rest for 20 years (Fig. 2) using climatological atmospheric forcing from ERA-40 (Uppala et al. 2005), with 3-hourly wind anomalies from the Fleet Numerical Meteorology and Oceanography Center 3-hourly Navy Operational Global Atmospheric Prediction System (NOGAPS) for the year 2003. The year 2003 is considered to be a neutral year over the 1993–present time frame in term of long-term atmospheric pattern the North Atlantic Oscillation (NAO). Heat fluxes are computed using the bulk formulae of Kara et al. (2005). The freshwater flux (evaporation, precipitation, and river runoffs) is treated as a virtual salinity flux, and the sea surface salinity (SSS) is restored to monthly climatology with a piston velocity of 15 m (30 days)−1. The salinity difference (between model and climatology) in SSS restoring is clipped to be 0.5 psu to diminish the damping effect of the restoring on ocean fronts [see Griffies et al. (2009) for a discussion]. There is no tidal forcing.
Figure 2 shows the time evolution of the mean kinetic energy for the three experiments. It takes approximately 5 years for the energy to stabilize, which is the time it takes for the first baroclinic mode Rossby wave to cross the Atlantic basin. The 1/25° is 50% more energetic than the 1/12°, not surprising since the viscosity/dissipation was cut in half when the resolution was increased. The 1/50°, on the other hand, has a mean kinetic energy level similar to the 1/25°. This is expected since the viscosity/dissipation was kept close to that of the 1/25° as the resolution was increased. The slight increase in mean kinetic energy is attributed to the increase in mesoscale and submesoscale activity (discussion in section 3b). Other integrated values such as the Florida Straits transport [30.8, 34.8, and 34.9 Sv (1 Sv ≡ 106 m3 s−1), respectively] and the average overturning streamfunction at 26°N (17.7, 17.8, and 17.5 Sv respectively) are in reasonable agreement with the observed values of ~32 (e.g., Meinen et al. 2010) and ~17 Sv (e.g., McCarthy et al. 2015), respectively.
3. Surface fields
a. Mean circulation
The mean SSH and the mean surface kinetic energy over the last 5 years of the three simulations are compared to the mean AVISO CNES-CLS13 SSH field derived from observations (Rio et al. 2014) over the Gulf Stream region (30°–55°N, 30°–80°W) in Figs. 3 and 4, respectively. Traditionally, a proper representation of the Gulf Stream separation in ocean numerical models has been a challenge [see Chassignet and Marshall (2008) for a review] and still remains an issue for many configurations despite the fact the major improvements are realized when one uses a horizontal resolution on the order of 1/10° (Bryan et al. 2007). In all three simulations, the Gulf Stream separates at the correct location at Cape Hatteras, but its eastward penetration into the interior differs greatly. At 1/12°, the modeled Gulf Stream (Fig. 3e) does not penetrate far into the interior and the recirculating gyre (Fig. 3e), and the highest eddy kinetic energy (Figs. 7c, 8c) is confined west of the New England Seamounts (60°W), results that have already been reported in Chassignet and Garraffo (2001), Haza et al. (2007), and Chassignet and Marshall (2008). The 1/25° simulation does not show a lot of improvement over the 1/12° simulation. It is arguably worse since the Gulf Stream in the 1/25° simulation not only does not extend as a coherent feature past the New England seamounts (Fig. 4d), it exhibits an unrealistically strong recirculating gyre southeast of Cape Hatteras (Fig. 3d) and has excessive surface variability (Figs. 7d, 8d) west of 60°W. It is only when the resolution is increased to 1/50° (Fig. 3c) that the Gulf Stream system appears to settle into a pattern that resembles the observations (Fig. 3a). The Gulf Stream penetration, recirculation gyre, and extension qualitatively compare very well with the latest AVISO CNES–CLS13 mean dynamic topography (MDT; Fig. 3a).
To validate the position of the current main axis as well as the width and intensity of the currents, one needs to be able to perform quantitative model–data comparisons on the scales of interest (i.e., tens of kilometers). There is not of lot of in situ observations that can be directly compared to the model results on those scales, but global climatologies have come a long way by combining altimetry, direct measurements, and surface drifters [see Rio et al. (2014) for a review] and are now able to provide MDT fields with a sharp definition of the western boundary ocean currents and associated fronts. There is, however, still quite a few uncertainties associated with these climatologies as shown by the differences between the 2009 and the 2013 CNES–CLS 1/4° climatologies [Figs. 3a and 3b; see Rio et al. (2011, 2014), respectively]. Improvements in the 2013 CNES-CLS13 MDT from the 2009 CNES-CLS09 MDT arise mostly from using a geoid based on the Gravity and Ocean Circulation Experiment (GOCE) instead of one based on the previous Gravity Recovery and Climate Experiment (GRACE) mission and from removing previously undetected undrogued drifters from the drifting buoy velocities distributed by the Surface Drifter Data Assembly Center (SD-DAC; Grodsky et al. 2011).
Comparison of the model results to these climatologies can be supplemented by comparison to mean dynamic topography constructed along altimetric satellite ground tracks using direct in situ measurements (Carnes et al. 1990; Blaha and Lunde 1992; Chassignet and Marshall 2008). Figure 5a shows the location of bathythermographic data taken during flights under the TOPEX altimeter ground track 932523A in September 1993 (courtesy of the Altimetry Data Fusion Center, Naval Oceanographic Office), where the mean dynamic topography was computed making the assumption that the dynamic topography relative to the geoid can be approximated by the dynamic height relative to a deep pressure surface that is parallel to the geoid (i.e., a level of no motion). Any error in this assumption will lead to different values of mean dynamic topography, but for the spatial scales of interest, as long as the level of no motion is deep enough, choosing a different level just adds or removes a bias. Figure 5b shows the mean SSH along-track 932523A from observations (derived from the bathythermographic data, CNES-CLS09, and CNES-CLS13) and from the three numerical simulations (1/12°, 1/25°, and 1/50°). The latest CNES-CLS13 climatology is closest to the MDT derived from the bathythermographic data, while the CNES-CLS09 profile shows a much stronger southern recirculation, very likely from contamination by the undrogued drifters, which are significantly impacted by wind slippage (Rio et al. 2014). As expected, the 1/50° simulation is closest to the observations with a SSH slope very close to the observations. The slope is weaker in both the 1/12° and 1/25° simulations with the 1/25° exhibiting the biggest departure south of 36°N, where the model is unrealistically energetic.
The 5-yr modeled mean surface velocity can also be compared to the long-term, 20-yr ADCP mean velocity measurements (Fig. 6) made by the Oleander between the U.S. East Coast and Bermuda (Rossby et al. 2014) from 1993 to 2012. As for the mean SSH, the mean position and direction of the 1/50° current vectors is closest to the observations, with the 1/25° being too far south and oriented meridionally and the 1/12° being too zonal. The width of the Gulf Stream is, however, slightly wider in the 1/50° model than in the observations, indicating more fluctuations of the mean axis of the Gulf Stream in the model. The maximum mean speed of the 1/50° modeled Gulf Stream (0.9 m s−1) is also a little less than observed (~1 m s−1; Rossby et al. 2010, 2014).
Figure 7 shows the root-mean-square (RMS) of the SSH for the AVISO observations and the three numerical simulations. The AVISO observations are based on 20 years of altimetry (1993–2012), and the modeled RMS fields are for the last 5 years of the simulations. As mentioned above in section 3a, both the 1/12° and 1/25° eddy variability are confined west of 60°W and the chain of New England seamounts. The 1/25° simulation also shows excessive variability southwest of the Gulf Stream axis, which is reflected in the mean SSH field. The overall distribution of the RMS SSH in the 1/50° is in reasonable agreement with the observations, especially the zonal extent, but it is significantly larger in magnitude, especially around the New England seamounts. We will argue below that the difference is primarily because altimetry only resolves eddies on O(150) km over a 10-day window. The RMS SSH results reported here are consistent with the distribution obtained by Hurlburt and Hogan (2000) using the 6-layer hydrodynamic (no active thermodynamics) NRL Layered Ocean Model (NLOM; Fig. 8). They also found that the eddy variability remained west of the New England seamounts for their 1/16° and 1/32° (equivalent to our 1/12° and 1/25°) and extend south for their 1/32°. It is only at 1/64° (equivalent to our 1/50°) that the variability extends farther eastward. This seems to imply that resolving submesoscale features of the O(10) km is a prerequisite for a correct eastward penetration of the Gulf Stream and the establishment of the recirculating gyres.
Figure 9 displays the observed and modeled surface EKE generated using the geostrophic velocities computed from the 1/50° SSH fields. As for the RMS, the overall distribution of the surface 1/50° EKE is in reasonable agreement with the observations, especially the zonal extent, but it is significantly larger in magnitude. The surface EKE shown in Fig. 9, however, does not take into account ageostrophic motions, and this raises the question as to how significant is the ageostrophic EKE component. Presumably, ageostrophic flows should become more significant as the resolution is increased and submesoscale features arise. Figure 10 displays the difference between the two components. The most striking feature is the asymmetry between the areas north and south of the Gulf Stream main axis, with the ageostrophic contribution being negative (positive) south (north) of the Gulf Stream. This asymmetry can be explained by considering that the largest ageostrophic contribution arises in an area where the flow curvature is significant (meanders, eddies) and where the velocities deviate from geostrophy and satisfy the gradient wind balance, that is, a primary balance between the centripetal acceleration, the Coriolis acceleration, and the horizontal pressure gradient force (Douglass and Richman 2015). Outside the Gulf Stream and high energetic areas, the surface Ekman flow becomes dominant. Ageostrophic motions are approximately 10%–20% of the total velocity in energetic areas of the Gulf Stream and can be as high as 200% in areas dominated by the surface Ekman flow (Fig. 11). A zoom on the region defined by the red rectangle in Fig. 11 shows that the ageostrophic motion in the eddies is always anticyclonic (Fig. 12). This is because the centrifugal force associated with the flow curvature is not taken into account when using the geostrophic balance, which results in the rotational velocities being under (over) estimated in anticyclones (cyclones; Chassignet et al. 1990).
The EKE map displayed in Fig. 9a was derived from along-track measurements optimally interpolated on a regular 1/4° grid by AVISO. By construction, the optimal interpolation filters many of the scales present in nature and is therefore not 100% representative of the observations on space scales less than 150 km (due to measurement noise and errors) and time scales less than 10 days (repeat cycle of the altimeters). To investigate the impact of the subsampling and optimal interpolation on the EKE fields, the 1/50° model outputs were filtered to be more representative of the AVISO gridded outputs. To quantify the impact of the filtering, we compute the SSH wavenumber power spectrum in the 1/50° configuration over the 10° by 20° box shown in Fig. 11. Over this energetic area, the 1/50° modeled power spectrum slope is k−5 in the 70–250-km mesoscale band (red in Fig. 13) and is representative of QG turbulence (see section 5 for a complete description of the power spectra distribution over the domain). Figure 13 illustrates the impact of the filtering on the wavenumber power spectrum. First, the subsampling of the model outputs to the 1/4° grid removes any information below 35 km (green in Fig. 13). Time averaging the outputs over 10 days (dark blue in Fig. 13) or applying a 150-km bandpass filter (turquoise in Fig. 13) bring the power spectrum slope closer to AVISO (black in Fig. 13), but it is only when the two are applied together (purple in Fig. 13) that the modeled wavenumber spectrum resembles most closely that of AVISO. This result is consistent with Biri et al. (2016), who report that the spectrum derived from a 10-day, resampled, 4-km model follows more closely the altimeter spectrum due to aliasing [see also Arbic et al. (2013), who performed a similar exercise but for spectral fluxes instead of spectra]. Applying both the 10-day and 150-km bandpass filters to the model outputs leads to EKE plots (Fig. 14) that resemble more closely the AVISO-derived EKE, suggesting that the AVISO-derived EKE underestimates the observed EKE by approximately 30% in the Gulf Stream region. One should note, however, that even filtered, the modeled surface EKE is higher than observed south of the New England seamounts, suggesting that interactions with the topography may be overemphasized in the model.
4. Interior flows
Most model–data comparisons usually focus on the surface fields because of the scarcity of long time series at depth covering a large spatial area. Scott et al. (2010), expending on the work of Penduff et al. (2006) and Arbic et al. (2009), used a large collection of moored current meter records to assess the ability of three eddying global ocean models [POP, Ocean Circulation and Climate Advanced Model (OCCAM), and HYCOM with 1/10° to 1/12° grid spacing], which resolve the first Rossby radius of deformation throughout most of the domain (i.e., up to 55°–60°N) to simulate the time-averaged total kinetic energy throughout the water column. They found that the models agreed within a factor of 2 above 3500 m and within a factor of 3 below 3500 m. Penduff et al. (2006) suggested that horizontal resolution was probably the most important factor limiting their 1/6° global model in generating realistic eddy kinetic energy, but Scott et al.’s (2010) results were not conclusive in that respect. Thoppil et al. (2011) did, however, find that increasing the model resolution to 1/25° significantly increased the surface and the abyssal EKE and clearly demonstrated that a better representation of upper-ocean EKE is a prerequisite for strong eddy-driven abyssal circulation. In this section, we do show that horizontal resolution has a large impact on the distribution of interior kinetic energy by comparing the three simulations to eddy kinetic energy maps generated from moorings and floats.
In Figs. 15 and 16, vertical sections of modeled zonal velocities and eddy kinetic energy along 55°W are compared to sections based on long-term observations using drifters, floats, and moored current meters (Richardson 1985); 55°W is perhaps one of the most observed sections across the Gulf Stream with measurements first during POLYMODE (Richardson 1985) and then during SYNOP (Bower and Hogg 1996). In this region, most of the deep oceanic variability is generated by the surface currents via vortex stretching. It is very weak at 1/12°, in agreement with the Scott et al. (2010) results. But it does increase significantly as the model resolution is refined, and, at 1/50°, the level and the pattern of the vertical zonal velocities and eddy kinetic energy resemble most closely the observations. This is consistent with Thoppil et al.’s (2011) statement that the surface and abyssal ocean circulation are strongly coupled through the energy cascades that vertically redistribute the energy and vorticity throughout the entire water column.
There are not that many spatial distributions of EKE at depths from observations. The few that exist are based on float measurements (SOFAR or Argo) and vary greatly in coverage and sampling. Figure 17 compares the modeled EKE distribution at 700 m to EKE derived from several years of SOFAR float measurements (Richardson 1993). As for the vertical sections, the 1/50° EKE distribution is closest to the observations, with a 1000 cm2 s−1 maximum west of the New England seamounts and an 800 cm2 s−1 extension past the seamounts. At 1000 m, the Argo float data used by Ollitrault and Colin de Verdière (2014) to derive the EKE do not provide the fine temporal and spatial coverage of the SOFAR floats, and the modeled 1/50° EKE differs significantly (Fig. 18). Filtering the modeled outputs over a 30-day time window and 100 km, as for the Argo data, leads to an EKE distribution that is much closer to the Ollitrault and Colin de Verdière (2014) map (Figs. 18a,c). As for the vertical sections, the magnitude of the 1/12° and 1/25° EKE distribution at 700 and 1000 m are significantly less than in the 1/50° (Figs. 17, 18).
5. Wavenumber spectra distribution
Wavenumber spectra are a useful tool to quantify the energy and variability associated with different scales and regions. In this section, we seek to address the questions of how robust is the modeled mesoscale and submesoscale eddy activity as a function of grid spacing and how representative is it of interior QG or SQG turbulence (Le Traon et al. 2008). The wavenumber spectra distribution is computed over the North Atlantic domain using boxes as in Le Traon et al. (1990), Paiva et al. (1999), Xu and Fu (2012), Richman et al. (2012), Sasaki and Klein (2012), and Dufau et al. (2016), among others. There is some sensitivity to the way spectra are computed, and we follow Sasaki and Klein (2012) by making the box doubly periodic in both the zonal and meridional directions following Lapeyre (2009). Cosine tapered windows are also often used in the literature instead of the doubly periodic method [see, e.g., Richman et al. (2012), which use a 10% cosine taper window], and in order to document the impact of the method, we compare in Fig. 19 the spectra computed using several cosine taper windows (0%, 5%, 10%, 20%, and 40%) and the doubly periodic approach. For a cosine taper window greater than 10%, the results are very close to the doubly periodic (or mirror) approach. The latter method is therefore preferred for its simplicity and reproducibility. The model outputs used to compute the spectra are daily averages, except when discussing the impact of internal waves and inertial motions where hourly snapshots are used.
Figure 20 displays the 1-yr SSH wavenumber spectra over the Gulf Stream energetic region (10° × 20° box defined in Fig. 11) computed from daily averaged outputs of year 20 for the 1/12°, 1/25°, and 1/50°, respectively. Over the mesoscale range of 70 to 250 km, the slope does not vary as a function of the horizontal grid spacing and is k−5.1, k−4.9, and k−5.0, respectively. It is now well documented that there is a strong seasonality associated with enhanced submesoscale activity in the winter mixed layer (Mensa et al. 2013; Sasaki et al. 2014; Callies et al. 2015; Rocha et al. 2016b). There is indeed a change in the slope of the SSH wavenumber spectra between summer and winter at all resolution for scales less than 70 km, but it is more pronounced at 1/12° than at 1/50° (Fig. 20). This does not mean that there is less small-scale instabilities, just that the SSH signature of the high Rossby number submesoscale features is less pronounced in this high EKE region. The seasonality is clearly visible when comparing the relative vorticity distribution between summer and winter in Figs. 21 and 22 and when computing the energy and vorticity spectra as in Figs. 23 and 24. The difference in the number of small-scale coherent features between the 1/12° and the 1/50° is striking, both in the summer and the winter (Figs. 21 and 22). But it is really at 1/50° that one can see an explosion of very small-scale features during the winter months [Fig. 22; see also Fig. 1 of Sasaki et al. (2014)].
Figures 23 and 24 compare the wavenumber spectra for SSH, energy, and relative vorticity between two regions: the highly energetic Gulf Stream region (33°–43°N, 50°–70°W), as shown in Fig. 11, and a low EKE region (20°–30°N, 20°–40°W). In the 70–250-km mesoscale range, the SSH wavenumber spectra slope is k−4.2 in the low EKE region versus k−5 in the high EKE Gulf Stream region. This is consistent with the results put forward by Richman et al. (2012) and Sasaki and Klein (2012), which suggest that highly energetic regions are closer to QG turbulence and low energetic regions are closer to SQG turbulence. Further examination of the SSH, energy, and relative vorticity spectra (Figs. 21 and 22) in the two regions show additional differences. First, in the low EKE region, the SSH wavenumber spectra slope exhibits a stronger seasonal cycle than in the high EKE region, which seems to indicate a stronger SSH signature of the smaller-scale features in winter. Second, the kinetic energy wavenumber spectra in the high EKE region does not differ much when using either the total velocities or the geostrophic velocities. This is consistent with the results of section 3b, which showed that the ageostrophic component forms a small percentage of the total velocity. In the low EKE region, on the other hand, the wavenumber spectra slope in the mesoscale range is reduced by 25% when using geostrophic velocities to compute the spectra. Finally, the biggest difference is in the relative vorticity spectra, which shows a stronger impact of the seasonal cycle in the low EKE region in the 70–250-km mesoscale range.
The spatial distribution of the SSH wavenumber spectra from altimeter observations shows a large spatial latitudinal variability with slopes closer to −5 at midlatitudes and as high as −1 in the tropics (Xu and Fu 2011, 2012; Zhou et al. 2015; Dufau et al. 2016). As in previously published modeling results (Paiva et al. 1999; Richman et al. 2012, Biri et al. 2016), we do not find this latitudinal dependence in the 70–250-km band (Fig. 25): the slope is everywhere between −5 and −4. Several explanations have been put forward to explain the discrepancy including aliasing and noise in the altimetry data (Biri et al. 2016), but it is also conceivable that one may underestimate the impact of high-frequency motions such as internal waves and tides when using daily averages to compute the wavenumber spectra (Richman et al. 2012; Rocha et al. 2016a). To investigate the impact of fast-moving features, the wavenumber spectra of Fig. 25 were recomputed using hourly snapshots (Fig. 26). The biggest difference between daily and hourly spectra is in the tropics, mostly on scales below 70 km where the slopes can be as low as k−1. The impact of fast-moving features on the 70–250-km wavenumber range is much smaller, that is, a decrease in the slope by up to 20% in a couple of 10° × 10° boxes (Fig. 26). Tides have been shown to have a significant impact on the wavenumber spectra on small scales (see Fig. 9 of Rocha et al. 2016a) but does not appear to have an impact on the latitudinal dependence of the wavenumber spectra in numerical models (Richman et al. 2012).
6. Summary and conclusions
In this paper, we quantify the impact of horizontal resolution (1/12° to 1/50°; 6 to 1.5 km at midlatitudes) on a series of identical North Atlantic experiments. First, we find that the representation of the Gulf Stream penetration and associated recirculating gyres shifts from unrealistic to realistic when the resolution is increased to 1/50°. This is consistent with results obtained by Hurlburt and Hogan (2000) using the 6-layer hydrodynamic NLOM model and by Lévy et al. (2010) using the NEMO model in an idealized domain. In all cases, the nonlinear effects of the submesoscale eddies that arise when the resolution reaches 1/50° intensifies the midlatitude jet and increases its penetration eastward. Second, the penetration of the EKE into the deep ocean drastically increases with resolution and closely resembles the observations in the 1/50° configuration. And third, the wavenumber spectra are independent of horizontal resolution and latitudes in the 70–250-km mesoscale range.
Convergence studies such as this one where most parameters are not changed, except for the horizontal resolution, are unusual and the question arises as what one may expect as you continue to increase the horizontal resolution, decrease viscosity, and/or increase the order of the numerical schemes. At some point, for the right amount of internal dissipation and friction, the solution should settle at a level close to the observations. In this paper, we showed that the level of EKE in the 1/50° simulation was comparable to the observations when taking into account the aliasing associated with the altimeter sampling but with one caveat: it was obtained by prescribing absolute winds at the ocean surface. While this is currently the norm for numerical models forced by an atmospheric reanalysis product, this does not allow any ocean feedback to the atmosphere. This feedback takes place via SST [see Small et al. (2008) for a review] and computation of the ocean current/wind shear [see Renault et al. (2016b) for a review]. Earlier studies have shown that the latter can lead to a significant reduction of the surface EKE, and Renault et al. (2016b) demonstrated that, using a coupled model, the current feedback deflects energy from the geostrophic current into the atmosphere and thus dampens eddies. Furthermore, the interaction of ocean eddies and the atmosphere regulates western boundary currents (Ma et al. 2016). In particular, Renault et al. (2016a) showed that the current feedback, through its “eddy killing” effect, can stabilize the Gulf Stream separation, a prerequisite for a proper separation and penetration (Özgökmen et al. 1997). The reduction in surface EKE induced by the surface current feedback can be as high as 30%. This means that future numerical simulations will need to be able to exhibit a higher level of EKE when using relative winds. As stated by Renault et al. (2016b), a bulk-forced oceanic uncoupled simulation should prescribe the surface stress using the relative wind in bulk formulae, which take into account a parameterization of the partial reenergization of the ocean by the atmospheric response. Renault et al. (2016b) propose a surface stress computed with a velocity U that is the wind relative to the current Ua corrected by a current–wind coupling coefficient sw as in U = Ua − (1 − sw)Uo, where Uo is the ocean surface current. The values of the coefficient sw are estimated to be between 0.2 and 0.3 and can vary spatially (R. Abel 2016, personal communication).
In conclusion, it is fair to state that the next threshold for a significant improvement in western boundary currents representation (i.e., Gulf Stream in this paper) is an increase in the horizontal resolution from the eddying 1/10° to submesoscale-enabled 1/50° grid spacing. Not only do the results presented in this paper support some of the results put forward by Hurlburt and Hogan (2000) and Lévy et al. (2010), it also raises the question of the usefulness of intermediate resolutions such as 1/25° or 1/36°, at least for the Gulf Stream region. The computational cost of simulations at 1/50° is extremely large, and while currently available resources do not currently allow for all the sensitivity experiments to numerical choices, stress formulations, vertical resolution, tidal impact, and so on, they will become more common in the future. Finally, we do realize that this paper is somewhat descriptive, but it sets the stage for in-depth analyses of water mass transformations and associated transports that will be presented in a companion paper.
EPC and XX are supported by the Office of Naval Research (Grant N00014-15-1-2594) and the NSF Physical Oceanography Program (Award 1537136). This work is a contribution to SWOT through the NASA Grant NNX16AH79G. The authors thank M. Ollitrault and A. Colin de Verdière for providing the 1000-m float data used in Fig. 18. They also thank B. Arbic, P. Klein, and two anonymous reviewers for their constructive comments. The numerical simulations were performed on supercomputers at the Navy DoD Supercomputing Resource Center, Stennis Space Center, Mississippi, using computer time provided by the U.S. DoD High-Performance Computing Modernization Program.