Multi-nest high-resolution model of submesoscale circulation features in the Gulf of Taranto

Recent oceanographic field measurements and high-resolution numerical modelling studies have revealed intense, transient, submesoscale motions characterised by a horizontal length scale of 100–10,000 m. This submesoscale activity increases in the fall and winter when the mixed layer (ML) depth is at its maximum. In this study, the submesoscale motions associated with a large-scale anticyclonic gyre in the central Gulf of Taranto were examined using realistic submesoscale-permitting simulations. We used realistic flow field initial conditions and multiple nesting techniques to perform realistic simulations, with very-high horizontal resolutions (> 200 m) in areas with submesoscale variability. Multiple downscaling was used to increase resolution in areas where instability was active enough to develop multi-scale interactions and produce 5-km-diameter eddies. To generate a submesoscale eddy, a 200-m resolution was required. The submesoscale eddy was formed through small-scale baroclinic instability in the rim of a large-scale anticyclonic gyre leading to large vertical velocities and rapid restratification of the ML in a time-scale of days. The submesoscale eddy was confirmed by observational data from the area and we can say that for the first time we have a proof that the model reproduces a realistic submesoscale vortex, similar in shape and location to the observed one.


Introduction
Ocean circulation is highly turbulent and occurs over a very wide range of scales, ranging from a few centimetres to thousands of kilometres. As such, it is driven by nonlinear scale interactions that can transfer energy upscale or downscale. Most of the total kinetic energy in the ocean is contained in "mesoscale eddies" (Ferrari and Wunsch 2009), which range from 10 to 300 km in size. Mesoscale flow fields play an important role in the transport and mixing of momentum and tracers across the planet's oceans. Mesoscale ocean eddies are monitored from space using satellite altimeters (Pujol et al. 2012) and are also explicitly resolved in ocean numerical simulations.
New high-resolution oceanographic field measurements and numerical simulations have revealed intense, transient, submesoscale motions characterised by a horizontal length scale of 100-10,000 m. Resolving these submesoscale motions in ocean numerical simulations requires horizontal grid resolutions of O(1 km). Recent very-high resolution models capable of directly resolving submesoscale features have shown significant deviation from eddy-resolving experiments, particularly in terms of the formation of numerous submesoscale eddies, fronts and filamental structures (Mensa et al. 2013;Sasaki et al. 2014;Gula et al. 2016). These experiments suggest that submesoscale physics is an important element in large-scale oceanic circulation. The important contribution of submesoscale processes to the vertical flux of both physical and biogeochemical tracers in the upper ocean has been illustrated by Capet et al. (2008), Thomas et al. (2008), and Lévy et al. (2012). In addition, the stratification of the upper layers by submesoscale processes and the enhancement of connections between the surface and the interior have been demonstrated by Fox-Kemper et al. (2008) and Klein et al. (2008).
In the last two decades, high-resolution numerical models have been used in idealized simulations to analyze submesoscale flow fields (Boccaletti et al. 2007;Lévy et al. 2010;Hamlington et al. 2014). The original idea was to use a 3D high-resolution submesoscale model with highly simplified initial and boundary conditions in order to investigate the impact of different physical processes and their parameterizations on the development of submesoscale structures. Realistic submesoscale-permitting high-resolution simulations have also been carried out in different areas of the ocean (see Table 1). Capet et al. (2008) investigated submesoscale activity in the California Current System using a set of single-nested model simulations with increasing horizontal resolution from 6 km to 750 m. Shcherbina et al. (2013) compared submesoscale statistics from observations with a multiple-nested simulation with up to 500 m resolution in the North Atlantic. Poje et al. (2014) studied submesoscale surface velocity fluctuations in the northern Gulf of Mexico from drifter measurements, including a direct comparison to a realistic model. Haza et al. (2016) analyzed the role of submesoscale motions on Lagrangian transport using a simulation of the Gulf of Mexico circulation with a horizontal resolution of 800m. Jacobs et al. (2016) performed a triple-nested model experiment in the Gulf of Mexico downscaling from 1 km down to 50 m.
In this study, we used realistic flow field initial conditions and a multiple nesting approach to achieve an open ocean horizontal resolution of 200 m, which is suitable for resolving submesoscale processes in this region and we used CTD measurements to confirm our submesoscalepermitting model predictions. This innovative approach will enhance our understanding of submesoscale eddy development and their dynamics. Pinardi et al. (2016) found observational evidence of intense mesoscale and submesoscale variability at the rim of a large scale, semi-permanent gyre in the Gulf of Taranto, Northern Ionian Sea, in the eastern Mediterranean (Fig. 1). The flow field showed a large-scale anticyclonic rim current with intensified jets along the border and an intense mesoscale cyclonic eddy on the eastward side of the gyre. Seven days later, the mesoscale eddy had disappeared and another eddy with a radius of approximately 10 km had formed in the north-western region. In this study, we used our realistic high-resolution model to explore the dynamics of the potentially submesoscale eddy on the eastern corner of the large-scale gyre velocity field over a one week timescale (mettere il periodo).
The paper is organized as follows. Section 2 briefly describes the structured grid component of the SURF platform and the related multi-nesting procedures. Section 3 presents the SURF implementation in the Gulf of Taranto and the main model parameters. The model validation is presented in Section 4. In Section 5, the emergent submesoscale substructures of the anticyclonic gyre are analyzed. The conclusions are presented and discussed in Section 8.  (Madec 2008). It includes a three-dimensional (3D) primitive free-surface ocean equation under hydrostatic and Boussinesq approximations along with turbulence closure schemes and a nonlinear equation of state, which couples the two active tracers (temperature and salinity) to the fluid velocity. The 3D space domain is discretised by an Arakawa-C grid where the model state variables are horizontally and vertically staggered. In the vertical direction, we used stretched z-coordinates distributed along the water column, with appropriate thinning designed to better resolve the surface and intermediate layers. Partial cell parameterisation was used (i.e. the bottom layer thickness varied as a function of position) in order to fit the real bathymetry.
Density was computed according to Jackett and McDougall's nonlinear equation of state (Jackett and Mcdougall 1995). A horizontal biharmonic operator was used for the parameterisation of lateral subgrid-scale mixing for both tracers and momentum. The horizontal eddy diffusivity and viscosity coefficients were parameterised as a function of the parent coarse resolution model. If a 0 is the parent viscosity or diffusivity, the nested model equivalent coefficient was a = a 0 ( x F / x L ) 4 , where x F is the nested grid spacing and x L is the large-scale model grid resolution. The vertical eddy viscosity and diffusivity coefficients were computed following the Pacanowsky and Philander's Richardson number-dependent scheme (Pacanowski and Philander 1981). For cases where unstable stratification was a possibility, a higher value (10 m 2 /s) was used for both the viscosity and diffusivity coefficients.
The Monotonic Upstream Scheme for Conservation Laws (MUSCL) was used for the tracer advection and the Energy and Enstrophy conservative (EEN) scheme was used for the momentum advection (Arakawa and Lamb 1981;Barnier et al. 2006). No-slip conditions on closed lateral boundaries were applied and the bottom friction was parameterised by a quadratic function. To evaluate the surface heat balance, atmospheric fluxes were computed through bulk formulas implemented in MFS (Pettenuzzo et al. 2010).
Two different numerical algorithms were adopted for the open boundary conditions depending on the prognostic simulated variables. For barotropic velocities, the Flather scheme (Oddo and Pinardi 2008) was used, while for baroclinic velocities, active tracers, and sea surface height, the flow relaxation scheme was used Engerdahl (1995). As the parent coarse resolution model only provided the total velocity field, the interpolated total velocity field into the child grid was split into barotropic and baroclinic components. In order to preserve the total transport after interpolation, an integral constraint method was imposed (Pinardi et al. 2003).

Gulf of Taranto model implementation
The model was implemented in the Gulf of Taranto (Fig. 2). A multi-scale conductivity-temperature-depth (CTD) oceanographic campaign was carried out from 1 to 11 October 2014 (Pinardi et al. 2016). The first survey was carried out between 1 and 3 October 2014, during which CTDs were acquired at a large scale (LS1). A second set of surveys on 8 October 2014 and between 8 and 11 October 2014 were conduced at a shelf-coastal scale (SC) and large-scale (LS2), respectively.

Model set-up
A triple nested models experiment was performed (NEST1, NEST2, and NEST3; Table 2). For each nesting, we set the grid spacing ratio to 3, so that the child domain had a grid spacing that was one third of the size of the parent domain. Parent and child models were linked by the initial and lateral boundary conditions. To reduce errors associated with the interpolation procedure, each nested domain was placed within the parent domain such that nested grid cells exactly overlapped the parent cells at coincident cell boundaries. The NEST1 model domain covered an area of On the vertical axis, each of the nested domain levels were the same and consisted of 120 z-levels with a stretching factor of h cr = 30 and a model level with maximum stretching of h th = 100. The locations of the vertical levels were smoothly distributed from 1.4 m to a maximum depth of 2945 m and had level thicknesses that increased with depth from approximately 2.8-90m. The vertical coordinate was defined from the reference coordinate transformation z(k) given by: where the coefficients h sur , h 0 , h 1 , h th and h cr are free parameters to be specified (see Madec et al. 2008 for details). Bathymetry was obtained from the General Bathymetric Chart of the Oceans (GEBCO) datasets by linear interpolation of depth data into the SURF model grid. This dataset contains ocean depths (in metres) at a 30 arc seconds resolution defined on a regular horizontal grid and covering the whole globe.
The initial and lateral boundary conditions for the nested model experiment were extracted from MFS daily mean datasets, which contain temperature, salinity, sea surface height (η) and total velocity (U,V) fields. The MFS model has a horizontal resolution of 1/16 • , contains 72 unevenly distributed layers in the vertical direction (Oddo et al. 2014) and includes a variational assimilation scheme based on the 3D-VAR method of Dobricic and Pinardi (2008). We used MFS data where only the LS1 observations were assimilated by the 3D-VAR scheme.
The atmospheric fields used to force the three consecutive nested models contained wind velocity at 10 m, temperature and humidity at 2 m, as well as total cloud cover and surface pressure from the European Center for Medium-Range Weather Forecasts (ECMWF) operational analyses, which have a 6-h frequency and spatial resolution of 0.125 • . Instantaneous precipitation values were computed from ECMWF operational forecast accumulated precipitations at a frequency of 3 h and a spatial resolution of 0.25 • .

Multi-nesting time concatenation and model spin-up time
Limited area ocean models require an initial spin-up time in order to produce dynamically adjusted fields after initialization from the interpolation of coarser ocean model fields (Simoncelli et al. 2011). We set two goals to define the spin-up time and design the multi-nesting time concatenation. On the one hand, we wanted to provide a realistic initial flow field for the model simulation. However, we also wanted to provide a forecast at different spatial resolutions for the period 9-11 October 2014 (i.e. corresponding to the LS2-CTD data collection days).
To attain both of these goals, the first nesting simulation started on 5 October 2014 at 00:00 and ran until 12 October 2014 at 24:00 (Fig. 3). The initial interpolated flow fields for this simulation were provided by the MFS model with LS1-CTD data assimilated. The second and third nesting simulations started at 00:00 on 6 and 7 October 2014, respectively, and ran until 12 October 2014 at 24:00. They were initialised from the interpolation of the parent coarse-grid model fields (the NEST1 and NEST2 models, respectively) after a spin-up time of one day. This multinesting time concatenation produced dynamically adjusted fields to the higher-resolution nested grid models within the time period of LS2 CTD data collection (i.e. 9-11 October 2014).

Horizontal smoothing with a Shapiro filter
After the first nesting, we found that the model generated gridpoint-scale noise. This consisted of oscillations Fig. 3 Multi-nesting time concatenations adopted by the structured and unstructured grid relocatable ocean platform for Forecasting (SURF) model in the MREA experiment. LS1 and LS2 denote the periods of the first and second large-scale surveys. CS denotes the data collection period for the coastal survey of the same magnitude as the background flow which appeared sporadically in many areas of the model domain and grew to large amplitudes in less than one day (Fig. 4). Numerical diffusion can be used to suppress these unphysical gridpoint-scale oscillations in models and is achieved by applying a spatial filter to the fields (i.e. a diffusive term that is applied separately to the variable each N time steps). Previous computational ocean dynamics studies have implemented the Shapiro filter to control the gridpointscale numerical oscillations (e.g. Shops and Loughe 1995;Klinger et al. 2006). This filter was introduced in the 1970s by Shapiro (1970) and Shapiro (1975). It is a high order linear filter that efficiently removes gridpoint-scale noise without affecting the physical structures of a field. The Shapiro filter of the 2N accuracy order applied to a variable based on the expression: wherew i is the filtered value of variable w at point x i , I is the identity operator and δ 2N is the even composition This filter is a discrete symmetric operator with a (2N + 1) point stencil. It acts as a low-pass filter that preserves the low frequency content (i.e. largest wavelengths) and totally dissipates the high frequency content (i.e. shortest wavelengths) from the the original field. In our triple nested models, a 4th order Shapiro filter (N = 2) was applied at all grid points and to the temperature, salinity, sea surface height and 3D velocity fields. The filter was applied separately to each layer. We adopted the simplest technique where fields were extrapolated to land before applying the one-dimensional (1D) filter in the zonal direction followed by another filter in the meridional direction. The filter involved was implemented ten times, once per day, which successfully eliminated the unwanted 2 x waves and significantly reduced the amplitudes of other poorly-resolved short waves, especially the 3 x and 4 x waves that tended to accumulate energy during model integration.

Validation of model predictions with CTD data
In order to test and quantify the improvement obtained by the higher resolution child model compared with the coarse resolution parent model, we evaluated the root mean square error (RMSE) between the quantities simulated by each nested model ψ m and the observed quantities ψ o , defined by: where N is the total number of CTD data points confined within the nested model domain and ψ stands for either temperature or salinity. For each nesting, we interpolate both child and father results over the CTD data depths and then computed the RMSE between data and simulations at 5-m intervals along each CTD cast. The RMSE resulting from the comparison between SURF results at different depths and CTD data is shown in Fig. 5  In order to demonstrate the similarity between the measured currents (Fig. 1) and the model results, we computed geostrophic velocities at 10 m depth for the NEST3 model. The geostrophic velocities are obtained from the dynamic height computed using the model temperature and salinity fields at 00:00 on 11 October 2014 (Fig. 6) and 100 m reference level. The geostrophic flow field shows a submesoscale cyclonic vortex with a diameter of ∼ 4 km in the northwestern border of the anticyclonic gyre. With respect to Pinardi et al. (2016), the center of the small-scale vortex in the NEST3 is shifted northward and eastward by approximately 5 km. Its diameter might seem slightly smaller than in Fig. 1 but the experimental submesoscale has been mapped only by one CTD station and its observational size depends on the objective analysis correlation length scale. We argue that this comparison demonstrates that the model is capable to reproduce the observed submesoscale to a degree of accuracy that has not been demonstrated before.

Submesoscale features associated with the anticyclonic gyre circulation
In this section, we analyze the emergent submesoscale structures of the central anticyclonic gyre of the Gulf of Taranto. The submesoscale dynamics are examined through a comparison of the various horizontal resolution model predictions from the coarser resolution 'eddy-resolving' MFS ( x ∼ 6 km) model in which submesoscale motions are not resolved, to the NEST1 ( x ∼ 2 km) and NEST2 ( x ∼ 0.7 km) models. The later are able to resolve the full range of mesoscale dynamics and part of the submesoscale regime, up to the finer resolution 'sub-mesoscale resolving' NEST3 ( x ∼ 0.2 km) model which is able to resolve almost the full range of spatial scales of the submesoscale dynamics. Root mean square error (RMSE) between the structured and unstructured grid relocatable ocean platform for forecasting (SURF)-child solutions and conductivity-temperature-depth (CTD) data (red dots) and between SURF-father results and CTD data (blue dots) for temperature ( • C; left panels) and salinity (psu; right panels) as a function of depth. The top panels show the large-scale Mediterranean Forecasting System vs. NEST1. Central panels show NEST1 vs. NEST2. Bottom panels show NEST2 vs. NEST3

Submesoscale structures in the anticyclone gyre border
To illustrate the effects of increased resolution, we considered horizontal sections at a 10 m depth of instantaneous temperature and current (Fig. 7), relative vorticity (Fig. 8), and vertical velocity together with potential density fields ( Fig. 9) at 00:00 on 11 October 2014 for the three embedded nests. Our goal was to detect the small, O(5) km submesoscale eddy in Fig. 1 in some of our nested model grids.
We found that at 5, 4, and 3 days after initial conditions for the NEST1, NEST2, and NEST3, respectively, the dominant mesoscale patterns could be recognised across all solutions. However, additional features emerged at smaller scales when permitted by increases in the grid resolution. In the MFS model (Figs. 7-9, top-left panels), the flow field showed a large-scale anticyclonic rim current with intensified jets along the border. The NEST1 model (Figs. 7-9, top-right panels) showed a density distribution with sharper gradients then the parent MFS model. In addition, smallerscale temperature patterns were visible along the border of the anticyclonic gyre, which were not present in the parent model. At this horizontal resolution, the model was not able to resolve most of the submesoscale regime, and submesoscale activity was absent. As the resolution was increased to 700 m (NEST2 - Figs. 7, 8, and 9, central-right panels) and 200 m (NEST3 - Figs. 7, 8, and 9, bottom-right panels), coherent eddies emerged leading to a denser and well-defined vortex population, covering a wide range of scales. A submesoscale cyclonic vortex with a diameter of 4 km was found in the northwest region of the central anticyclonic gyre located at 40 • 17 N, 17 • 4 E for both NEST2 and NEST3, which was much weaker in the first, and well defined in the latter. This submesoscale vortex was confirmed by observational data collected in the study area (Fig. 1).
The similar but different structure of the vorticity, density, and velocity fields in the area of the submesoscale in NEST2 and NEST3 is interesting. Starting with NEST2, the vorticity field (Fig. 8) became generally filamentous but particularly in the area of the unstable rim current. The latter   (Fig. 7) generating an intensified convergence downstream of the split, where the submesoscale forms. The main difference between the flow field in NEST2 and NEST3 was the intensity of the rim split jet and the secondary jet that only formed in NEST3. With respect to vorticity, in NEST2, the vorticity was threadlike, while in NEST3 the vorticity threads thickened, at the location of the submesoscales and in convergent areas. Overall, the complexity and scales of the vorticity in NEST3 were unique with respect to the other resolutions, showing that the 200 m grid spacing gives a dynamically different regime. Regarding the density and vertical velocity field, NEST2 and NEST3 were again structurally different from the other nestings and MFS, as shown in Fig. 9. Sharper fronts and larger vertical velocity areas developed along the unstable rim current, and in the case of NEST3, the isopycnal surfaces with a cold density anomaly contour detachment folded at the location of the submesoscale. In NEST3, the vertical velocity amplitude in the submesoscale frontal regions reached values of 90 m day −1 .

Evolution of submesoscale structures
The submesoscale structure that formed along the northwestern region of the central anticyclonic gyre is illustrated in Fig. 10. This figure shows the sequence of snapshots of the horizontal distribution of temperature (left panels), relative vorticity (central panels), and vertical velocity (right panels) at a 10 m depth every 6 h from October 10 2014 00:00 (top panels) to October 10 2014 18:00 (bottom panels). A meander visible initially at 40 • 15 N, 16 • 55 E on the northwestern border of the central gyre was formed by splitting the rim current into two jets that re-join after a few tens of kilometers, forming a large meander. In the downstream region of the split, a convergence zone had been created, indicated by larger areas of downwelling vertical velocities and a region of high density intrusion (26.3 entrapped density versus 26.1 sigma along the rim current). After 12 h, the positive vorticity thread evolved into patches that define the cyclonic sub-mesoscale. The vorticity signature of the submesoscale was non-symmetric with a larger positive vorticity area at the interface with the main rim current. After 18 h, a 26.3 density patch had detached from the previous high density intrusion and larger upwelling velocities had developed near the cyclonic submesoscale eddy center. To sum up, the birth of a submesoscale seems to be the small-scale equivalent of a meandering mesoscale jet and eddy formation which occurs at an accelerated rate in several hours and over space scales of a few tens of kilometres. The dimensions of the detached density patches and of the vorticity threads were in the order of a few kilometres, and vertical velocities were larger than in the mesoscale equivalent cases (Robinson et al. 1988;Ikeda et al. 1989;Stern 1989).

Surface mixed-layer
The submesoscale activity was associated with an intense vertical velocity in the upper ocean which increased with the increasing resolution (Fig. 11). The longitude crosssection was taken along latitudes 40 • 16 48 N, see section indicated in Fig. 10, where the small-scale meandering perturbations evolve into submesoscale vortices, as disussed in Section 5.2. The NEST1 model showed vertical velocities of an order of a few meters per day, typical of eddyresolving models. As the resolution was increased to 700 m (NEST2) and 200 m (NEST3), intense vertical velocities developed in the upper ocean up to ∼ 100 m day −1 , which are typically one order of magnitude larger than those associated with the mesoscale. The vertical velocity distributions also showed alternate upwelling and downwelling regions alternately, and high vertical velocities values were found in the proximity of the vortices. The intense vertical velocities in the mixed layer may communicate with the mesoscale upand down-welling associated with the geostrophic frontal meander scale. Thus, submesoscale processes are instrumental in transferring properties and tracers, vertically, between the surface ocean and the interior. They enhance, for example, the nutrient supply and the exchange of dissolved gases with the atmosphere (Brannigan 2016).
Another impact of submesoscale activity was the restratification of the mixed layer (ML). The contour lines in Fig. 10 Submesoscale structures in NEST3 temperature ( • C) and current data (left panels), relative vorticity (s −1 ; central panels), and vertical velocity (m/day) together with contour lines of potential density (contour intervals are 0.1kg m −3 ; right panels) at z = 10 m. The time interval between consecutive rows is 6 hours starting from 00:00 on 10 October 2014 (top panels) and ending at 18:00 on 10 October 2014 (bottom panels). 'S' denotes the position of submesoscale features. Black arrows, whose length is proportional to velocity, denote the ocean current direction Fig. 11 represent the evolution of the potential density for various horizontal resolutions. At 00:00 on 10 October 2014 (top panels), the ML was weakly stratified with numerous lateral density gradients. During the course of the simulation, we observed a reduction in the strength of the lateral density gradient, highlighting that submesoscale turbulence contributes to a rapid restratification by the slumping of the horizontal density gradient in the ML. The 2-km resolution model failed to contain the influence of the submesoscale dynamics on the restratification of the ML. However, as the horizontal resolution increased from 700 to 200 m, an increase in the restratification rate was clearly observed, reconfirming the importance of submesoscale dynamics in the restratification phase.

Discussion: submesoscale eddy generation
It has been suggested that mesoscale-driven surface frontogenesis energizes submesoscale flows (Lapeyre and Klein 2006;Roullet et al. 2012). The theory of frontogenesis at the ocean surface is well understood and the essential physics can be understood with quasi-geostrophic (QG) dynamics (Stone 1966). A mesoscale strain field sharpens lateral buoyancy gradients at the surface more effectively than in the interior of the ocean. Ageostrophic circulation develops in response to the increasing lateral buoyancy gradient, as described by the omega equation (e.g. Hoskins et al. 1978). In the interior, this circulation weakens the lateral buoyancy gradient and consequently light water downwells on the dense side and dense water upwells on the light side of the gradient. However, at the surface, vertical velocity disappears and ageostrophic circulation cannot counteract the increase in the lateral buoyancy gradient. Therefore, the mesoscale strain field is left unopposed, leading to the formation of strong submesoscale surface fronts.
Another mechanism that may energize submesoscale flows in the ML are the small-scale baroclinic or frontal instabilities, sometimes called 'mixed-layer instabilities'. These instabilities allow perturbations to extract the potential energy stored in the lateral buoyancy gradients generated by mesoscale stirring or by spatial variations in atmospheric forcing. In idealised ML models, where ML modes grow on a prescribed front (Boccaletti et al. 2007;Fox-Kemper et al. 2008), at a few hundred metres of depth, ML disturbances grow on horizontal scales to an order of 1-10 km and over timescales of 1 day. In the real ocean, ML instabilities occur in the presence of energetic mesoscale eddy fields. ML modes can therefore grow on mesoscale buoyancy gradients and can be shared by mesoscale strain fields. Eddy edge waves can propagate both along the surface and along the sharp increase in stratification at the base of the ML.
These two mechanisms produce distinct submesoscale flow characteristics and vertical fluxes (Callies et al. 2016). The main differences compared with mesoscale-driven surface frontogenesis are that ML instabilities energize the entire depth of the ML and produce larger vertical velocities. The net relative vorticity small scale increases could be in line with the development of MLIs, which release submesoscale eddy kinetic energy (EKE) extracted from available potential energy (APE) through slumping of the isopycnals, which contributes to rapid restratification.
The submesoscale eddy generation mechanism involves both of the above mentioned processes during part of the submesoscale formation period. The formation of the submesoscale in this paper was generated by baroclinic instability of the rim current of a large scale anticyclonic gyre (also indicated as frontogenesis). This instability leads to regions of strain which sharpened lateral buoyancy gradients. As a consequence, sharp density fronts (Fig. 9) and vorticity threads or filaments (Fig. 8) are produced. In this paper, we document for the first time the single submesoscale eddy generation mechanism because we have observational data at the local scale. It is shown that the rim current becomes unstable and splits into two jets. A sharp front develops where the submesoscale eddy is generated (Fig. 10), probably due to the secondary instability of the split jet (local frontogenesis). A density patch occlusion occurred which forms the cyclonic submesoscale eddy with large downwelling/upwelling velocities and thickening of the vorticity threads.

Comparison with previous studies
Several studies have been carried out in the last few years to analyze submesoscale appearance using realistic very highresolution simulations as listed in Table 1. These studies have been performed in several areas of the ocean and they differ from the present study partly in terms of their research questions and scope but they are a useful benchmark to compare with our results. Capet et al. (2008) investigated the submesoscales in the California Current System. They used realistic highresolution simulations with the Regional Oceanic Modeling System (ROMS). They used an off-line, one-way nesting technique to downscale from the large-scale ROMS configuration (USW12) with 12-km horizontal grid spacing to different resolution up to 0.75 km horizontal resolution. They demonstrated that a regime transition occurs at the 1.5 to 0.75 km resolution, looking at the change in the spectrum slopes with varying resolution. The key process that generated the submesoscales was found to be frontogenesis which sharpens surface density fronts down to a horizontal scale of few kilometers making them baroclinically unstable. The fronts exhibit small-scale meanders indicative of submesoscale frontal instability process occurring in the surface layer. Shcherbina et al. (2013) presented vorticity, divergence, and strain rate statistics of submesoscale turbulence in the North Atlantic Mode Water region using current velocity data obtained with synchronous ADCP sampling. The observations were statistically compared with numerical model predictions in this region using multiple-nested realistic ROMS simulations, from about 6 km resolution down to 0.5 km. They confirmed that submesoscales are more often cyclonic as it is the case of our vortex. Furthermore they also indicate the atmospheric forcing as responsible for the initiation of submesoscale activity which might be also our case. Pinardi et al. (2016) document that the week before the submesoscale was found in the observational data, a large wind and precipitation event occurred, probably enhancing frontogenesis of the Gulf of Taranto temperature gradients. In addition, here we show that the submesoscale eddy is connected instability of the Gulf of Taranto gyre border which might be what Shcherbina et al. [2013] mention as an important source of submesoscale turbulence in terms of nonlinear mesoscale eddies interactions. Poje et al. (2014) studied the structure of submesoscale surface velocity fluctuations in the northern Gulf of Mexico using high-frequency position data provided by the nearsimultaneous release of 300 surface drifters. Their data showed that the submesoscale fluctuations were setting the local dispersion properties for drifter clusters at 100 m scales, not captured by the large scale eddy-resolving data assimilating models. Our results show the localization of the submesoscale turbulence along the gyre border at 200 m resolution and this might indicate a general agreement with the drifter data experiment. Haza et al. (2016) and Jacobs et al. (2016) analyzed the role of submesoscale motions on Lagrangian transport in the Gulf of Mexico circulation. They used realistic highresolution simulations with the HYbrid Coordinate Ocean Model (HYCOM). They showed that the main processes affecting the drifter clustering include mesoscale frontogenesis which is also our case. Furthermore, they found that drifters released inside a mesoscale vortex could leak out of it by the submesoscale turbulence that might be also true for our case as shown by the vertical velocity inside the submesoscale that is indicative of horizontal divergence in the surface velocity field.
More investigations are needed to understand the role of submesoscale eddies and their formation. In our paper, we show for the first time that a realistic single submesoscale vortex is generated via a frontogenetic process at the border of a gyre/mesoscale eddy vortex as probably documented by the previous literature but non shown at the local, single submesoscale eddy process.

Summary and conclusions
In this study, the submesoscale motions associated with a large-scale anticyclonic gyre in the central Gulf of Taranto were examined using submesoscale-permitting simulations with realistic flow field initial conditions and multiple nesting techniques. The hierarchy of nested models highlighted one-way consecutive horizontal grid nesting from a parent grid with a resolution of x ∼ 6000 m, and with child grids at x ∼ 2000 m, x ∼ 700 m, and x ∼ 200 m.
The flow field showed a large-scale anticyclonic rim current with intensified jets and additional structural complexity (e.g. elongated filaments and small-scale cyclonic vortices) emerging in higher resolution nests. To generate submesoscale eddies, a 200-m resolution was required. A submesoscale cyclonic vortex with a diameter of 4 km was found in the northwest region of the central anticyclonic gyre located at 40 • 17 N, 17 • 4 . This eddy was confirmed by observational data collected in the study area. Increasing the resolution highlighted an increase in the vertical velocity field of up to 100 m/day near the surface, confirming that submesoscale turbulence is associated with intense vertical movements in the upper ocean. During the course of the simulation, in higher resolution nests, a reduction in the strength of the lateral density gradient was observed, highlighting that submesoscale turbulence contributes to a rapid restratification by slumping the horizontal density gradient in the ML.
Our submesoscale eddy generation mechanism was probably due to the instability of the rim current of a large-scale anticyclonic gyre where both the density front sharpening and vorticity threads or filaments are produced. The rim current split as in a typical mixed barotropic/baroclinic instability process. The development of MLIs produced larger vertical velocities and rapid restratification by slumping the horizontal density gradient.
We believe that this work is original because we show the growth of submesoscale structures around an anticyclonic gyre that dominates the circulation in the Gulf of Taranto. Submesoscale activities in this region have never been investigated before with realistic very high-resolution simulations. In addition, the stratification of this region produces a small Rossby deformation radius of 10-12 km which is different from that of the Atlantic (40-50 km). We directly compared our high resolution geostrophic current with the observational evidence and the comparison results in a major validation of our model output. Thus, we can say that for the first time we have a proof that the model reproduces a realistic submesoscale vortex, similar in shape and location to the observed one.