Subsurface Mesoscale Eddy Generation in the Ocean off Central Chile

Off the coast of central Chile, subsurface anticyclonic eddies are a salient feature of the oceanic circulation, transporting a signi ﬁ cant fraction of coastal water that is rich in nutrients and poor in dissolved oxygen offshore. In this study, the formation mechanism of these eddies is analyzed through a high ‐ resolution (~0.3 km) and low ‐ resolution (~3 km) oceanic model that realistically simulate the regional mean circulation, including the Peru ‐ Chile Undercurrent (PCUC). An analysis of the vorticity and eddy kinetic energy in both simulations indicated that the subsurface eddies can be triggered through a combination of processes that are associated with instabilities of the PCUC. In the high ‐ resolution simulation, we observed that the interaction between the PCUC and topographic slope generates anticyclonic vorticity and potential vorticity close to zero in the bottom boundary layer. The separation of the undercurrent from the slope favors the intensi ﬁ cation of anticyclonic vorticity. It reaches magnitudes that are larger than the planetary vorticity while kinetic energy is converted from the PCUC to the eddy ﬂ ow. These processes set the necessary conditions for the development of centrifugal instabilities, which can form submesoscale structures. The coalescence of submesoscale structures generates a subsurface anticyclonic mesoscale eddy. In the low ‐ resolution simulations (>3 km) centrifugal instabilities are not simulated, and the barotropic conversion of the mean kinetic energy into eddy kinetic energy appears as the main process of eddy formation. We showed that the vertical structure of these eddies is sensitive to the spatial resolution of the model.


Introduction
Mesoscale eddies play a central role in the transport of heat, momentum, and dissolved substances in the ocean. These eddies affect the primary productivity and biogeochemistry of the ocean and, consequently, are potentially important for the climate of the planet Dong et al., 2014;McGillicuddy et al., 2007). The contribution from mesoscale eddies to the total transport can be comparable to that from large-scale horizontal wind-and thermohaline-driven circulation (Zhang et al., 2014), particularly in regions that are characterized by weak mean circulation, such as the Eastern Boundary Upwelling Systems (EBUSs). In addition to surface mesoscale eddies, which have been rather well documented from satellite data (e.g., Chaigneau et al., 2009;Chaigneau & Pizarro, 2005;Chelton et al., 2011), EBUSs show a distinctive class of subsurface or intrathermocline eddies that are much less well known. These eddies are commonly anticyclonic and have a diameter close to the Rossby radius and a core with a thickness of approximately 500 m (Combes et al., 2015;Hormazabal et al., 2013;Kurian et al., 2011). The anomalies that are generated by this type of eddy have a convex lens shape, with a core that is frequently located in the permanent pycnocline (Collins et al., 2013;Johnson & McTaggart, 2010). These eddies contain relatively homogeneous waters presenting a thermostad-or pycnostad-in its core.
Subsurface anticyclonic eddies have been also observed in other locations in the global ocean. For instance, Richardson et al. (2000)-based on subsurface float measurements-document subsurface eddies in the North Atlantic, which contain water that originated from the Mediterranean Sea (called Meddies). Subsurface mesoscale anticyclonic eddies have also been observed in the Indian Ocean. There, they are sometimes called Reddies because of transport waters from the Red Sea (Shapiro & Meschanov, 1991). Subsurface eddies in the Northwest Pacific off the coast of California are associated with instabilities in the California undercurrent and are called California Undercurrent Eddies or Cuddies (Jerónimo & Gómez-Valdés, 2007;Molemaker et al., 2015). Subsurface eddies in the Southeast Pacific have been described through a combination of data from oceanographic cruises, Argo floats, and gliders (Hormazabal et al., 2013;Johnson & McTaggart, 2010;Thomsen et al., 2016). Generally, these types of eddies in the EBUS could be referred to as Puddies because they originate from the poleward undercurrent (Frenger et al., 2018). Here, we use this name to refer to the subsurface anticyclonic eddies observed off Chile.
Modeling studies in the Pacific EBUS Combes et al., 2015;Kurian et al., 2011) showed that surface eddies are characterized by a maximum vorticity magnitude in the upper 100-m depth and that subsurface eddies have a maximum vorticity magnitude from several tens to a few hundred meters below the surface (~200-m depth). Additionally, surface and subsurface eddies have an asymmetric distribution. Surface eddies of cyclonic vorticity are slightly more frequent than anticyclonic eddies, whereas subsurface eddies are mainly associated with anticyclonic vorticity and would be as frequent as surface anticyclonic eddies Combes et al., 2015).
The eastern South Pacific is characterized by high primary production and by one of the most extensive oxygen-minimum zones of the world ocean (e.g., Fuenzalida et al., 2009). In this region, Puddies are produced by the destabilization of the Peru-Chile Undercurrent (PCUC; Hormazabal et al., 2013;Combes et al., 2015;Thomsen et al., 2016), which flows poleward at the top of the slope and continental shelf and transports Equatorial Subsurface Water (ESSW; Johnson & McTaggart, 2010;Hormazabal et al., 2013). This same water mass is found in the Puddy cores. Off central Chile the ESSW is characterized by relatively warm temperatures, high subsurface salinities (34.2-34.6), high concentrations of nutrients, and very low values of dissolved oxygen (Silva et al., 2009). Puddies can move long distances to the west, transporting ESSW to the subtropical gyre (Frenger et al., 2018). This phenomenon can contribute to the expansion of the oxygen-minimum zones (OMZ) toward well-oxygenated waters, which could have important ecological effects by modulating denitrification and nitrogen fixation (Landolfi et al., 2013). Nearshore, the transport of Puddies can remove subsurface nutrients before they are exposed to the surface by upwelling, reducing both coastal productivity (Gruber et al., 2011;Renault, Deutsch, et al., 2016) and the intensity of the OMZ in this region (Stramma et al., 2012).
Several processes can destabilize the undercurrent and generate subsurface eddies, such as baroclinic or barotropic instabilities, topography effects, or convection (Bosse et al., 2016;D'Asaro, 1988;McWilliams, 1985). For the Puddies that are generated in the EBUS, the formation mechanism must explain the preference toward anticyclonic vorticity. Molemaker et al. (2015) performed numerical simulations off the coast of California and suggested that the mechanism of the Puddy formation is associated with the development of strong anticyclonic vorticity in the bottom boundary layer that exists between the California undercurrent and the slope. When the undercurrent is destabilized, it can separate from the slope, instigating submesoscale instabilities and, thus, submesoscale anticyclonic eddies. These eddies can interact with each other to form a subsurface eddy of a larger scale. Thomsen et al. (2016) described the formation of a Puddy off the coast of central Peru (approximately 12.5°S) based on a comprehensive set of observations, claiming that their results were consistent with the idea from the modeling study by Molemaker et al. (2015).

Journal of Geophysical Research: Oceans
Numerical simulations commonly used to study the mean and mesoscale circulation off central Chile have a too coarse spatial resolution to simulate properly the processes highlighted by Molemaker et al. (2015). The first goal of this study is to describe the mechanisms involved in the Puddy formation, focusing on submesoscale processes and turbulent regimes. To that end, a submesoscale-resolving oceanic simulation is carried out over a domain off central Chile. The second aim of this study is to analyze the effect of the spatial resolution of the model on the Puddies formation by comparing a low-resolution simulation (mesoscale resolving) to the high-resolution simulation.
The paper is organized as follows: The simulation setup is presented in section 2. In section 3, we analyze the origin of the anticyclonic vorticity and potential vorticity that feed the Puddies. A case study of Puddy formation is presented in section 4. In this section, we describe the main processes that are associated with Puddy formation. Furthermore, we analyze the eddy-mean interactions (section 4.2). The analysis of the effect of varying the model resolution and topography on the Puddy formation is analyzed in section 5. Section 6 includes a summary and perspectives for future work.

a MUR Sea Surface Temperature
Satellite-derived sea surface temperature (SST) from the Multi-scale Ultra-high Resolution (MUR) data set were compared with SST from the model. MUR is a daily global data set that is available from 1 June 2002 to the present. This product is a combination of data from the Advanced Very High Resolution Radiometer, Moderate Imaging Spectroradiometer Terra and Aqua, and Advanced Microwave Spectroradiometer-EOS instruments. Ultrahigh resolution (~1 km) is achieved by using the technique called multiresolution variational analysis, which is a statistical interpolation method that is based on wavelet decomposition. Technical details can be found in Chin et al. (1998). These data were provided by the Jet Propulsion Laboratory under support by NASA Making Earth Science Data Records for Use in Research Environments program (https://mur.jpl.nasa.gov/).

b Hydrographic Data
To assess the model's ability to reproduce the observed vertical stratification, we used observations of temperature and salinity obtained with a Conductivity-Temperature-Depth profiler (CTD model SBE-25) during an oceanographic cruise supported by the Fondo de Investigación Pesquera, Chile (FIP2005-01) on board the research vessel R/V Abate Molina between 6 and 17 December 2005. CTD measurements were conducted across five transects between 35.5°and 37.5°S in intervals of 0.5°. Each transect had 10 stations. We used only data from 35.5°S and 36.5°S and from stations between the coast and 75°W. The locations of the stations that were used in this work are shown in Figure 1.

Model Configuration and Forcings
In this study, we used the Regional Ocean Modelling Systems (ROMS; , in its Coastal and Regional Ocean COmmunity (CROCO) version (Debreu et al., 2012). CROCO solves primitive equations by using a terrain-following vertical coordinate with Boussinesq and hydrostatic approximations. Three simulations with successively horizontal grid nesting refinements from a parent grid resolution were performed. The parent simulation used a domain that extended from 10°N to 40°S and from the coast to 90°W with a 1/12°horizontal resolution. This simulation was run for the period 1958-2008 and was analyzed and validated in Dewitte et al. (2012) and Vergara et al. (2016Vergara et al. ( , 2017. This simulation allows for a realistic representation of the PCUC because of its relatively high resolution and usage of a domain that extends up to the equatorial region. To ease the interpretation of the results, the nesting procedure was performed offline from larger to finer scales without feedback from the child grid solution onto the parent grid (Penven et al., 2006). We used the output of the parent simulation (L0) between June 2005 and May 2006 for the initial and boundary conditions of a child simulation with a resolution of~3 km (L1). This year was selected from L0 because it corresponds to a period close to the climatological mean (without an intense El Niño or La Niña event). L1 was run for 15 years. The last year of simulation L1 was used for the initial and boundary conditions for a simulation with 1-km resolution (L2). The period between 15 November 2005 and 15 May 2006 of simulation L2 was used to generate the initial and boundary conditions for a simulation at 0.3-km resolution (L3). The spin-up times of L1, L2, and L3 were assumed to be less than (or equal to) 2 years, 1 year, and 2 months, respectively. The initial and boundary conditions were processed through the ROMS-to-ROMS package by using the methodology in Mason et al. (2010). The Puddy formation was analyzed from the L3 simulation for the period between 14 January and 14 May 2006, which corresponds to a period in the calendar year when the PCUC is more intense (e.g., Combes et al., 2015).
The domains of these three simulations are shown in Figure 1. The topography that was used in L1 and L2 was obtained from SRTM15-plus (Shuttle Radar Topography Mission), which is a global bathymetry data set with 15 arc-sec nominal resolution (~0.5 km; http://topex.ucsd.edu/WWW_html/srtm30_plus.html). To reduce errors that are associated with the pressure gradient, the bottom topography was interpolated onto the model grid and smoothed following Penven et al. (2005). Local smoothing was applied where the steepness of the topography exceeded a factor r ¼ Δh=h ¼ 0:2, where h is the depth and h is the mean depth. The topography of L3 was obtained directly from grid L2.
The model uses terrain-following vertical coordinates (Shchepetkin & McWilliams, 2009), which can select a transition depth between the z level and sigma level (h c ) that does not depend on the minimal depth of the model topography (Lemarié et al., 2012). We chose the transition depth h c = 300, and the parameters that controlled the bottom and surface refinement of the grid were σ b = 2 and σ s = 7, respectively. Simulations L1 and L2 had 50 vertical levels, while L3 had 80, which enabled us to resolve details in the topography and have enough vertical levels in the intermediate ocean layers.
The atmospheric forcing for all the simulations was the same as that for the parent domain. The momentum flux was taken from a downscaled product based on NCEP/NCAR data by using QuikSCAT data for the training period of the statistical model. Over the period that was considered here for the simulation (i.e., June 2005 and May 2006), the momentum flux was almost equivalent to using wind fields from QuikSCAT data (the reader is invited to refer to Goubanova et al., 2011, for details in the statistical model). Heat and freshwater fluxes were derived from the bulk formula (Fairall et al., 2003) by using the air temperature from COADS monthly climatology (da Silva et al., 1994). Short-and long-wave radiation and the relative humidity were taken from COADS. This choice was made initially for the parent simulation because of the large biases in the NCEP atmospheric fluxes near the coast of Peru and Chile. For simplicity and consistency, we used the same approach for the nested domains.
The mechanical interaction between the ocean and the atmosphere (i.e., the current feedback) transfers kinetic energy from the ocean to the atmosphere, inducing surface stress, wind anomalies, and a roughly 30% dampening of mesoscale activity (sometimes called eddy killing; e.g., Renault, Molemaker, et al., 2016). Renault, Molemaker, et al. (2016) suggest a simple parameterization of the wind response to the current feedback that is based on the coupling coefficient between the surface current and the wind (s w ). In the bulk formulae (Fairall et al., 2003) of CROCO, the surface stress can be computed using the relative wind (Ur) to a portion of the oceanic motions following: where Ua is the 10-m wind and Uo is the surface current. Although such a parameterization presents some limitations (e.g., in this study, s w does not present any spatial or temporal variations and is taken as 0.3; Renault et al., 2019), it has the advantage of producing more realistic levels of eddy kinetic energy (EKE) compared to altimetry.
CROCO allows the use of different vertical mixing parameterizations to indirectly include small-scale and microscale processes that significantly affect the movement of larger scales. As in the parent simulation, the vertical mixing of the child simulations was parameterized by using the KPP scheme (K-profile parameterization; Large et al., 1994).
The child simulations generally reproduced the main characteristics of the SST off central Chile. The SST from L3 was compared to MUR data averaged between 14 January and 14 May 2006 ( Figure 2). Along the coast, L3 presented a cold bias by~2°with respect to satellite observations ( Figure 2c). L1 and L2 showed similar spatial patterns to L3 (not shown). While the magnitude of this cool bias was uncertain because of limitations of the satellite data (see Dufois et al., 2012), such a bias is a common feature of regional simulations in EBUS when using satellite scatterometer winds (e.g., Mason et al., 2011;Penven et al., 2005). Scatterometers have a blind zone in coastal regions (within 25 km of the coastline), so winds that are used to force models are extrapolated toward the coast. Therefore, scatterometer winds do not adequately represent the shoreward decreases of alongshore wind speeds near the coast (Astudillo et al., 2017;Capet et al., 2004) because of orography, surface roughness, and air-sea interactions . The poor representation of the wind drop-off yields an overestimation of Ekman transport and an underestimation of Ekman pumping but also tends to enhance vertical mixing and vertical advection (Renault et al., 2012), which also contributes to the cold bias. The wind drop-off can also modify the structure of coastal currents. Consistent with the Sverdrup dynamics, a weak wind drop-off results in a weak wind stress curl that deepens the undercurrent and intensifies the surface current (Capet et al., 2004;Song et al., 2011). Consequently, the current shear is intensified and the isopycnal tilt is increased in the presence of a weak drop-off. These conditions favor the development of baroclinic instabilities that generate eddies. The shape of coastal wind profiles that are used to force models could thus have a relevant role in the mechanism of Puddy formation through modifying the speed and depth of the PCUC core, which should be considered when interpreting the results.
Despite the limitation in the wind forcing mentioned above, our model simulated a relatively realistic mean vertical stratification. Figures 3 and 4 show vertical sections of the temperature and salinity, respectively, from simulation L2 and in situ data at 35.5°S and 36.5°S between 7 and 12 December 2005. In situ data were obtained for December 2005 at the previously mentioned latitudes ( Figure 1). The simulation realistically simulated the vertical structure of the temperature at the different latitudes of interest ( Figure 3). The largest bias was observed above 150 m, where the simulation overestimated the mean temperature by~2°C (Figures 3c and 3f) compared to the observations. The model also tended to underestimate the salinity in the upper 200 m ( Figure 4). However, the maximum salinity at 250 m was reasonably reproduced by the simulation (Figure 4). Importantly, this maximum salinity was associated with the core of the ESSW waters. The difference in the salinity between the observations and L2 showed a similar spatial pattern to the difference in the temperature. Additionally, the vertical structure of the density (sigma-t) from the simulation coincided with that from the observations. As expected, the largest difference in the density between the  in situ data and L2 was located in the upper 150 m. The bias in the stratification could be explained by either spurious diapycnal mixing from advection error (Marchesiello et al., 2009), mesoscale processes that were not reproduced by the model, or the biases in the atmospheric forcing. Overall, the agreement between the model and observations was reasonable considering that the model did not assimilate data.

Anticyclonic Vorticity Along the Slope
Modeling studies in the eastern South Pacific showed that most subsurface mesoscale eddies are anticyclonic Combes et al., 2015;Frenger et al., 2018). According to Molemaker et al. (2015), the origin of the anticyclonic vorticity that is associated with these types of eddies originates from the bottom boundary layer of the poleward undercurrent that flows over the slope. Along the Chilean coast, this mechanism can be related to the PCUC, which is rapidly reduced in the sloping boundary layer, creating both horizontal (∂V/∂x) and vertical (∂V/∂z) velocity shear components. In the above expressions, V is the southward (alongshore) flow and x and z are the cross shore and vertical coordinates, which are positive east and upward, respectively. Thus, the vertical component of the vorticity in the boundary layer (∂V/∂x) is anticyclonic.
To analyze the generation of anticyclonic vertical vorticity when the PCUC interacts with the slope, we calculated the mean vertical vorticity distribution for a 4-month period (14 January to 14 May 2006) in the higher-resolution simulation (L3). Figure 5 shows the mean meridional current and mean vertical Journal of Geophysical Research: Oceans vorticity in different cross-shore sections and the mean flow at 150-m depth. In the upper slope and continental shelf, where the PCUC core is close to the slope (Figures 5a, 5c, and 5d), a narrow band of anticyclonic vorticity is observed near the bottom (Figures 5e, 5g, and 5h). In the 33°30′S section, where the PCUC is detached from the coast (Figure 5b), a vorticity intensification in the slope boundary layer was not observed (Figure 5f). The values of the vorticity in Figure 5 were normalized by f (the Coriolis parameter, which is negative in the Southern Hemisphere). The sections centered at 33°S and 34°S presented more intense poleward flow and consistently larger anticyclonic vorticity values in the boundary layer. Important alongshore variability existed in the poleward flow ( Figure 5m) and the associated boundary layer vorticity, which may have been related to changes in the topography and corresponding changes in the form drag (Gula et al., 2015;Molemaker et al., 2015).
Observations in the core of Puddies that formed in the eastern South Pacific indicated that the Ertel potential vorticity (PV) was close to zero (Thomsen et al., 2016).
is the velocity vector with components in the zonal, meridional, and vertical axes; and b = −gρ/ρ o , g, ρ, and ρ o are the buoyancy, gravity, density, and a reference density, respectively. If PV is of opposite sign to f (i.e., fPV < 0), then the fluid is unstable (Hoskins, 1974). For an easy interpretation, the fPV quantity is used. A possible source of fPV close to zero for the Puddies is the bottom boundary layer. Previous studies showed that the bottom boundary layer may act as a source or sink for the potential vorticity over the slope depending on the direction of the current that interacts with the slope (Williams & Roussenov, 2003). Off the coast of central Chile, the undercurrent is oriented in the direction of the Kelvin wave propagation, so downslope Ekman flow advects lighter water under denser water, driving diabatic mixing and reducing fPV magnitude (Benthuysen & Thomas, 2012). This relationship was verified in the simulation (Figures 5i-5l): When the undercurrent interacted with the slope, a strip of fPV with values close to zero was observed. Importantly, this source of fPV close to zero was not the only source because other diabatic processes could have been involved in the Puddy formation, for example, mixing that was associated with submesoscale instabilities, as shown in section 4.2.

Formation of Subsurface Eddies: A Case Study
To identify the main processes that were involved in the generation of subsurface mesoscale eddies off central Chile, we analyzed the generation of a well-developed eddy that was observed in the highest-resolution (L3) simulation during April 2006. Figure 6 shows the evolution of the vertical component of relative vorticity normalized by f (i.e., ζ z /f), emphasizing the region where the Puddy formed (see the circles in the thick dotted line). Prior to the formation of the Puddy (30 March), the vorticity was smaller to the north of 33°10′S. Between 33°10′S and 33°35′S, a strip of anticyclonic vorticity was observed along the shelf break and upper slope, reaching magnitudes larger than f (Figure 6a; yellow contours in this figure indicate anticyclonic relative vorticity equal to|f|). This strip of anticyclonic vorticity weakened on 5 April (Figure 6b). To the north of 33°20′S, the magnitude of the anticyclonic vorticity increased. This intensification was also observed on 9 April and was associated with the formation of submesoscale coherent structures (Figure 6c). These structures interacted to produce a larger eddy (Figure 6d), which finally separated from the coast and moved slightly southward (24 April, Figure 6e). To the north of the studied eddy (33°24′S; 72°30′W), an intense submesoscale vortex also formed at the coast but did not merge with our main eddy. Later, the main eddy continued moving to the southwest, transporting the vorticity that was generated at the coast. On 14 May, the Puddy was observed at 34°50′S and 74°W (Figure 7). The perturbation of the isopycnals in the eddy core had a typical convex shape, with positive temperature and salinity anomalies (corresponding to ESSW) and fPV close to zero (Figures 7b-7f). These properties characterize Puddies off central Chile (Combes et al., 2015;Hormazabal et al., 2013). Another Puddy was present near 33°42′S and 74°24′W (Figure 7). The formation mechanism of this last Puddy cannot be described because it was not formed in this simulation but rather inherited from the northern boundary condition. This eddy showed similar characteristics to our case study eddy although its diameter is somewhat larger and its core deeper (Figures 7g-7k).
According to our hypothesis, the vorticity flux from the boundary layer that fed the eddy was associated with the separation of the PCUC from the slope. We identify and describe this separation process in the next section.

Undercurrent Separation
The main source of vorticity for Puddies is the bottom boundary layer, but feeding these eddies requires that the boundary layer separates from the slope, transporting the vorticity away from the coast. Both numerical and observational studies showed that the poleward undercurrent in the EBUS can separate from the coast, commonly near abrupt changes in topography Thomsen et al., 2016). Different factors can contribute to coastal current separation: the so-called β effect, the stretching of the vortex tubes, and changes in the bottom topography (Marshall & Tansley, 2001). For Puddies off California, the California Undercurrent separates from the slope because of form drag and bottom-pressure torque from changes in the bottom topography . This phenomenon causes the undercurrent to slow down near the slope and drift off the coast. The separation of the PCUC and its role in the formation of Puddies off the coast of central Chile is analyzed below. However, the causes of this separation are not discussed in this study, but we suggest that the separation is likely produced by similar processes to those for the coast of California as the separation is associated with bottom pressure torque and positive bottom pressure anomaly (which were calculated following Molemaker et al. (2015) and Song & Wright (1998); see Supporting Information, Text S1 and Figure S1).
We focus on the region where a Puddy formed in April 2006 (near 33°25′S). Before this Puddy formation (30 March), the undercurrent had a magnitude of~25 cm/s following the continental shelf break and the upper slope, which are represented by the 200-and 500-m isobaths, respectively (Figure 8a). When the Puddy began to form (5 April) in a region with changes in the topography (33°25′S, indicated by the green line in Figure 8b), the core of the undercurrent moved offshore and substantially weakened, and a positive (northward) current was observed at the coast. This slight coastal current became later a component of the Puddy structure ( Figure 8c). On 14 April, this eddy began to move to the southwest (Figure 8d). After the Puddy formed, the poleward undercurrent returned to the coast in the northern half of the domain, increasing its magnitude near the slope and generating another intense strip of anticyclonic vorticity (24 April; Figure 8e). This current separation and Puddy formation process could also be observed in the vertical

10.1029/2018JC014723
Journal of Geophysical Research: Oceans sections centered at 33°25′S (Figure 9; the section is marked with a horizontal green line in Figures 6 and 8).
As observed in Figure 8, prior to the formation of the Puddy, the poleward undercurrent was more intense and was located close to the slope, with a core near 200-m depth. Below the core of the undercurrent, the isopycnals abruptly bent downward (Figure 9a). When the Puddy began to form, the intensity of the undercurrent decreased and the current displaced offshore. Simultaneously, a positive current developed along the upper slope and continental shelf and the isopycnals flattened (Figure 9b). Subsequently, the currents showed a clear Puddy structure with a core centered at 150-m depth (Figure 9d).
During the Puddy formation, the fPV presented a similar spatial pattern to the vertical component of the relative vorticity between 100-and 500-m depth. Regions with anticyclonic vorticity (Figures 9f-9j) were closely related to regions with fPV close to zero or even slightly negative magnitudes (Figures 9k-9o). Before the Puddy formation, the PCUC interacted with the slope and the slope had anticyclonic vorticity greater than f ( Figure 9f) and a negative fPV (Figure 9k). The region that was dominated by negative fPV extended beyond the boundary layer. When the fPV was negative and far from the boundary, the flow was unstable and inertial instability (also called symmetric instability by Hoskins, 1974) could be generated (e.g., Hoskins, 1974;Thomas et al., 2013). Inertial instability (or negative fPV) was observed in the latitudinal section (33°25′S) in Figure 9 and along the coast of the study region, where the anticyclonic vorticity was close to or greater than |f| and where the undercurrent tended to detach from the slope (Figures 6 and 8).
When the undercurrent separated from the slope, an anticyclonic vorticity core was observed off the slope at approximately 150-m depth (Figures 9h and 9i). The Puddy could also be clearly recognized from the distribution of fPV magnitudes close to zero. This characteristic is typical of Puddy formation, which allows the eddy core to contain relatively homogeneous water. The origin of fPV close to zero could be directly associated with friction in the bottom boundary layer (section 3.1) or the effects of submesoscale instabilities.

Eddy Kinetic Energy Budget
To identify the sources of eddy kinetic energy that is involved in the generation of Puddies, we analyze the terms that are associated with the transfer of kinetic and potential energy between the mean and eddy flows. These terms include the horizontal (HRS) and vertical (VRS) Reynolds stresses and the vertical buoyancy flux (VBF). By using Reynolds decomposition, the velocity and buoyancy can be expressed as u; is the temporal mean of the variable over the period of the Puddy formation (i.e., 30 March to 30 April 2006 and (′) is the fluctuation with respect to the mean. Then, the Reynolds stresses and vertical buoyancy fluxes can be defined as follows (e.g., Gula et al., 2016;Harrison & Robinson, 1978;Kang & Curchitser, 2015): The overbars in (2)  The integrated terms between 100-and 500-m depths are shown in Figure 10. This depth range corresponds to the region of interaction between the PCUC and the slope and is not directly influenced by processes in the surface mixed layer (the integrated terms below 500-m depth did not show a significant difference with Figure 10). The most relevant terms were the HRS and VBF (Figures 10a and c). However, the VBF averaged over the region of the Puddy formation (shown in Figure 10) had smaller magnitudes than the HRS (5 · 10 −2 and −1.36 · 10 cm 3 /s 3 , respectively), suggesting that the HRS was the main source of energy for the EKE. The VRS was an order of magnitude smaller compared to the other two terms (see Figure 10b, where the values were multiplied by 10). The same order of magnitude for the HRS and VRS was observed off the coast of California . In the surface layer, the main transfer of energy was from EPE to EKE

10.1029/2018JC014723
Journal of Geophysical Research: Oceans (VBF > 0), instead, in the subsurface, VBF was roughly an order of magnitude smaller than in the surface layer.
The HRS was characterized by alongshore variations (Figure 10a) that could be explained by the effect of the topography on the energy conversion terms. The shear increased when the flow pressed against the slope because of a combination of bathymetric curvature and inertia effects (Gula et al., 2015). To the north of the topography change (33°30′S), where the eddy that was analyzed in the previous sections was generated, the HRS presented positive values, that is, the PCUC transferred kinetic energy to turbulent flows. Regions where HRS > 0 could be associated with submesoscale instabilities (Gula et al., 2016).
We redefined the Reynolds decomposition of the velocity and buoyancy to analyze the temporal evolution of the energy conversion terms toward the EKE during the Puddy formation. Now, the operator ( − ) consisted of a 10-day running average centered at time t and (′) denotes deviations from this average. We chose this period based on the temporal power spectral density (Zhang et al., 2016) that were computed from the horizontal velocity at 150-m depth between 32°S and 35°S and from the 150-m isobath to 40 km offshore. In the power spectral density, a change in slope occurred at around 10 days (Supporting Information, Figure S2). Thus, the averages were associated with the mean and low-frequency changes (temporal scale longer than 10 days), while the perturbations contained a shorter temporal scale were related to the submesoscale (on the order of f −1 ; e.g., Thomas et al., 2008). By using these new definitions for ( − ) and (′), we estimated the evolution of the terms of the energy transfer (2), with the horizontal and vertical Reynolds stresses and the vertical buoyancy flux called HRS*, VRS*, and VBF*, respectively. We only show the results for HRS* because VRS* and VBF* presented much smaller magnitudes than those of HRS*.
During the Puddy formation, regions where HRS* was positive over the 200-and 500-m isobaths existed ( Figure 11). However, the magnitudes of HRS* in these regions presented temporal and spatial variations. These variations were mainly associated with the interaction of the PCUC with the slope. In regions where the undercurrent interacted with the slope but was close to separate, HRS* was positive (e.g., at 33°25′S on 30 March; Figures 8a and 11a). We also observed HRS* > 0 at 33°35′S on 14 April (Figure 11d), but this result was produced by the interaction between the topography and the northward current of the Puddy. During the Puddy formation in the region where the eddy was located, a dipole of HRS* > 0 and HRS* < 0 was observed, which was associated with the shear of the Puddy currents.
In the regions upstream of the PCUC separation, HRS* was positive and matched with regions of large anticyclonic vorticity, whose magnitudes could exceed the planetary vorticity (Figure 6), and negative fPV ( Figure 11). These regions are indicated in Figure 11 by rectangles in red lines. If fPV ≤ 0 in a region far from the bottom boundary layer, inertial instabilities can occur (Hoskins, 1974). Inertial instabilities that extract the energy of the mean current through horizontal shear (i.e., HRS* > 0) are called centrifugal instabilities (e.g., Thomas et al., 2013). The generation of these instabilities by the interaction between currents and

10.1029/2018JC014723
Journal of Geophysical Research: Oceans topography has been studied numerically in the California Undercurrent Jiao & Dewar, 2015;Molemaker et al., 2015) and Gulf Stream (Gula et al., 2016). This type of instability is associated with energy dissipation, producing small-scale turbulence and diapycnal and isopycnal mixing  and modifying the fPV in turn. Diapycnal mixing and subsequent geostrophic adjustment are the main mechanisms for submesoscale coherent structure generation (McWilliams, 1985). D'Asaro (1988) showed evidence of the formation of these structures in association with the interaction between the current and topography. As these submesoscale structures merge together, the centrifugal instability is stabilized, creating a larger eddy . The eddy has close to zero fPV at its core because centrifugal instability increases the anticyclonic relative vorticity and it promotes vertical mixing that weakens the stratification Jiao & Dewar, 2015).
The regions where VBF* varies do not necessarily correspond to where the Puddy forms (Supporting Information, Figure S3). Negative VBF* values can be expected where centrifugal instabilities are present; nevertheless, we did not observe a significant relationship between the centrifugal instabilities and VBF* < 0. Note that positive VBF* values are also present along the slope. They can be associated with submesoscale baroclinic instability. This instability develops in weakly stratified layers over sloping topography and can be generated by an Ekman adjustment of the undercurrent on a slope (Wenegrat et al., 2018). In summary, the transfer of kinetic energy from the PCUC to the turbulent flux (HRS*) was the main source of energy for the Puddies.

Meridional Velocity, Relative Vorticity, and Potential Vorticity
To investigate the sensitivity of the above results to the bottom topography, we ran the L1, L2, and L3 simulations by using the same configuration in section 2 but deriving the model grids from the bottom topography of the parent L0 grid (i.e., 1/12°resolution) instead of the SRTM15-plus bathymetry. We call these new

10.1029/2018JC014723
Journal of Geophysical Research: Oceans simulations L1s, L2s, and L3s. To compare the intensity of the PCUC in all the simulations inside the region of Puddy formation (between 32°33′S and 33°51′S), we estimated the distribution of the meridional speed at 150-m depth over the outer continental shelf and upper slope (i.e., in the region where the bottom depth was deeper than 150 m and shallower than 500 m). All the simulations showed a unimodal distribution with mean values between −6.6 and −9.8 m/s, consistent with a predominant poleward current (Figure 12a). The calculations were performed for the period between 14 January and 14 May 2006, corresponding to the analyzed period in simulation L3. L1s and L3s distributions showed a larger amount of poleward (negative) speeds than L3 distribution, but the average speed in L3 (−7.7 ± 0.15 cm/s) was between L1s (−9.8 ± 0.13 cm/s) and L3s (−6.6 ± 0.12 cm/s) mean speeds (the errors were calculated by a bootstrap procedure; Efron & Tibshirani, 1994). According to Figure 12a, the changes in the model resolution had a greater effect on the probability density function (PDF) of the poleward current than changes in the details of the bottom topography. Nevertheless, the changes in the model resolution (i.e., comparing L1 s and L3) seemed to have a slightly larger effect on the mean intensity of the undercurrent, with reduced amplitude in the model with finer resolution and more realistic bottom topography.
In contrast to the above results on the poleward undercurrent, the vertical component of the relative vorticity was more sensitive to changes in model (and topography) resolution (Figure 12b), as also was the fPV. In the different simulations, the distribution of the vorticity (normalized by f) was skewed to anticyclonic values, with a skewness of −0.93, −1.05, and −0.13 for L1 s, L3s, and L3, respectively. Both L3 and L3s had relative vorticity values greater than |f| (i.e., values of ζ z /f < -1), while L1s only had values smaller than |f|. Notably, the ζ z /f distribution had larger tails for the L3 simulation, consistent with an increasing relative vorticity when both the topography and model resolution increased. A better representation of the slope may, therefore, contribute to the generation of vorticity, as discussed in section 3. The fPV was predominantly positive in all the simulations. In fact, L1 and L1s showed only positive fPV values. In contrast, L3s and L3 showed some negative fPV values, particularly L3 (Figure 12c). Regions with negative fPV may have been associated with submesoscale centrifugal instabilities. Because L1s (and L1) did not show relative anticyclonic vorticity larger than |f|, this simulation could not replicate centrifugal instabilities, which would be related, in turn, to negative fPV. In summary, higher model resolution and a more realistic slope contribute to the generation of large anticyclonic relative vorticity and reduced fPV, which contribute to the formation of centrifugal instabilities.

Puddy Formation
In L1s, during the Puddy formation, consistent with the high-resolution simulation, anticyclonic vorticity over the slope are generated and the undercurrent separates from the slope. However, striking differences are observed during the formation. In L1s, the interaction between the PCUC (−0.15 m/s) and the slope

10.1029/2018JC014723
Journal of Geophysical Research: Oceans generated anticyclonic vorticity (as observed on the vertical section at 33°S, not shown), but the latter had a lower intensity (~0.25f) compared to what was observed in L3 (Figures 5a, 5e, and 5i); In addition, the fPV magnitudes reduced less in L1s than in L3 (not shown). During the Puddy formation in L1s (between 18 March and 13 April 2006), the PCUC separated from the coast, transporting the anticyclonic vorticity offshore and forming a vorticity filament (Figure 13a). No background strain could have stabilized this filament, so this structure became unstable and rolled up (Figure 13 and c) into a mesoscale vortex (Figure 13d; Hide & Titman, 1967). Because the vorticity amplitude in the low-resolution simulation did not reach values that allowed centrifugal instabilities to develop, the undercurrent separation process in that case is not associated with submesoscale turbulence.
As a further analysis of the effect of resolution on the mechanism of Puddy formation, we estimated the terms of energy transference (2) in L1s, where the temporal mean for the Reynolds decomposition corresponded to the period of the Puddy formation (i.e., 15 March to 15 April 2006) in L1s. In the region where the Puddy formation began (33°10′S and from the coast to 72°W), the transfer of energy was mainly driven by barotropic instabilities (HRS; Figure 14a). VBF remains a second order driver as its average over the region shown in Figure 14 is an order of magnitude smaller than HRS (5.08 · 10 −2 and 4.75 · 10 −1 cm 3 /s 3 , respectively). The VBF dipole revealed in Figure 14c is the signature of the Puddy, which was not observed in L3 (Figure 10c). Possibly, the signature of the Puddy in L3 is hidden by other mechanisms, which are relevant for VBF but not necessarily associated with the Puddy formation, for example, internal waves (Barkan et al., 2017;Bühler & McIntyre, 2005;McWilliams, 2016;Vanneste, 2013). This topic requires further analysis and is beyond the scope of the present paper.

Structure of the Puddies
The vertical structure of the Puddies had weaker anticyclonic vorticity and higher fPV in the low-resolution simulations (Figures 13h and 13i) than the high-resolution simulations (Figures 6e, 6f, 6j, and 6k). Thomsen et al. (2016) studied a set of observations off the coast of Peru and showed that the core of a Puddy had an anticyclonic vorticity of approximately 0.75f and potential vorticity that was close to zero. These values are consistent with those from L3 during Puddy formation and values from the study by Molemaker et al. (2015). Additionally, the temperature and salinity in the Puddy core were not affected by the resolution (Figures 6b, 6c, 6g, and 6h and 13e and 13f).
In both L3 and L1s, the total deformation rate (S) and Okubo-Weiss parameter (W) of the Puddies are defined as where S s = ∂v/∂x + ∂u/∂y is the shearing deformation rate and S n = ∂u/∂x-∂v/∂y is the stretching deformation rate. These terms are relevant in the eddy kinematics, as they allow to determine the existence of a vortex, which occur where rotation dominates over the deformation rate (i.e., W < 0; Okubo, 1970;Weiss, 1991). Puddies simulated in L1s and L3 present the same order of magnitudes of S and W (O(10 −4 ) and O(10 −9 ), respectively). However, the core of the Puddy from L3 is characterized by larger magnitude of W than L1s (Supporting Information, Figure S7).
To describe the main kinematics properties of the Puddies, we estimated the diameter, thickness, rotation period, and swirl velocities of the Puddies from L1s ( Figure 13) and L3 (E1 in Figure 7). The boundaries of the Puddy were identified using the Okubo-Weiss parameter (using a threshold of W = 10 −12 /s 2 ). Puddies from L1s and L3 did not show a significant difference in its kinematics properties. Vertical velocity has a relevant role within Puddies as they may transport tracers vertically (Klein & Lapeyre, 2009). In Puddies, it has been observed that vertical velocities are generated by a combination of horizontal deformation and advection of vertical relative vorticity by ageostrophic vertical shear (Barceló-Llull et al., 2017). When the formation of the Puddy ends (see Figures 6e and 13d), the vertical velocities as simulated by L3 are 3 times larger than those in L1s (Supporting Information, Figure S8). This means that ageostrophy processes are more intense in L3. Barceló-Llull et al. (2017) showed that in the main pycnocline, the vertical velocities are characterized by a dipolar structure. This structure is also found in L1 s (Supporting Information,Figures S8b and S8d) and L3 (Supporting Information,Figures S8a and S8c) although the vertical velocities of the L3 simulation are somewhat noisy at 200-m depth.

Summary and Discussion
In the previous literature, off central Chile, Puddies have been described using in situ observation from oceanographic cruise (Hormazabal et al., 2013) and Argo floats (Johnson & McTaggart, 2010), and numerical models Combes et al., 2015). These studies suggested that the Puddy formation is associated with the PCUC, but the specific formation mechanisms were not analyzed. Off California, Molemaker et al. (2015) described for first time the formation mechanism of the poleward undercurrent eddies. They show that the formations result from a combination of processes associated with the destabilization of the undercurrent, which gives rise to submesoscale instabilities. Numerical models commonly used to study the mean and mesoscale circulation off central Chile have had a too coarse spatial resolution to represent these processes. In this study, a Puddy formation was analyzed by carrying out a set of oceanic simulation that resolves submesoscale processes. The sensitivity of the mechanism associated with the formation of the Puddy is assessed by comparing a submesoscale with a mesoscale resolving simulation.
In the highest-resolution simulation, we showed that the processes that were associated with Puddy formation were similar to those for other regions Thomsen et al., 2016). We observed that the interaction between the PCUC and the slope produce anticyclonic vorticity in the bottom boundary layer and fPV close to zero mainly because of the generation of anticyclonic vorticity by lateral shear over the slope (note that the reduction of stratification in the bottom layer also contributed to reduce fPV). Consistently with previous studies, the Puddy formation began when the poleward undercurrent was separated from the slope, transporting the anticyclonic vorticity and fPV close to zero that was generated in the bottom boundary layer to the open ocean. When the undercurrent was separated, the Puddy core began to form. The eddy increased its size through the coalescence of submesoscale structures that formed when the PCUC became unstable. During this process, the relative vorticity and potential vorticity in the bottom boundary layer were advected offshore.
The undercurrent separation was associated with submesoscale processes and the development of instabilities. To analyze these processes, an energy budget was performed to estimate the sources and sinks of the rate of change of EKE. In contrast to what is usually observed in the surface layer in EBUS systems (Leth & Shaffer, 2001;Marchesiello et al., 2003), the EKE in the subsurface was produced by horizontal Reynolds stress that extracted kinetic energy from the PCUC and not by vertical buoyancy; such a balance is sensitive to the horizontal resolution.
Over the analyzed period of Puddy formation in the regions where the PCUC separated from the slope, the anticyclonic vorticity was greater than the Coriolis parameter and the submesoscale kinetic energy was extracted from the PCUC. These conditions are necessary for the development of centrifugal instabilities, which are associated with a loss of balance. During the Puddy formation, centrifugal instabilities induce a transfer of kinetic energy from the PCUC to the Puddy and a diapycnal mixing, which is also responsible of the formation of submesoscale and mesoscale structures. The mixing generated by these instabilities is also associated with dissipation, which causes a direct energy cascade. When the submesoscale structures begin to merge, the dynamic becomes governed by geostrophic turbulence. These results are verified through a wave number spectrum analysis (Supporting Information, Figure S4), which reveals the evolution of the turbulent processes associated with the transfer of kinetic energy (Ferrari & Wunch, 2009). During the development of centrifugal instabilities the spectrum slope is equal to −2 suggesting ageostrophic processes at work in the energy cascade (Boyd, 1992;Callies & Ferrari, 2013;Capet et al., 2008;Klein et al., 2008). By contrast, when the submesoscale structures merged, the observed spectrum slope is close to −3, which is consistent with interior quasi-geostrophic turbulence (Charney, 1971;Callies & Ferrari, 2013;Klein et al., 2008).

Journal of Geophysical Research: Oceans
The mechanism of Puddy formation is furthermore analyzed in a low-resolution simulation that does not resolve submesoscale processes. In this simulation, the relative vorticity in the slope was less intense and fPV was larger than in the high-resolution simulations. This result implies that the necessary conditions for the development of centrifugal instabilities were not met in the low-resolution simulations, such as L1s (described in section 5), resulting in a more stable current system than in the high-resolution simulations. Conversely to the high-resolution simulations, where centrifugal instabilities developed, Puddy formation was triggered by barotropic instabilities in the low-resolution simulations. A wave number spectral analysis during the Puddy formation from L1s reveals that the magnitudes of the slopes are larger than 3, which means that submesoscale processes are not reproduced by L1s (Supporting Information, Figure S5).
The vertical structures of Puddies are affected by the spatial resolution of the simulations. We expect that the Puddies from simulations with a lower resolution than 3 km present a greater difference with the Puddies from simulations analyzed in this study. From the parent simulation (L0; simulation used to forces the child simulations) we analyzed the vertical structure of a Puddy at 33°36′S and 73°12′W (Supporting Information, Figure S6). The Puddy core shows lower anticyclonic vorticity (ζ z~− 0.2 f) and higher magnitudes of potential vorticity (fPV~0.5 · 10 −13 ) in L0 than in L1s and L3 Puddy simulations. We also calculated the deformation rate and the Okubo-Weiss parameter of the Puddy in L0 simulation, which were roughly five times smaller than in the L1s and L3 Puddy simulations (Supporting Information, Figure S7). This means that Puddies from L0 tend to be less coherent than Puddies from the highest-resolution simulation. These results could suggest that Puddies from L0 have less chance to survive when interacting with other oceanic flows than Puddies in the high-resolution simulation. Furthermore, previous studies have shown that the lifetime of the eddies is shorter when they have smaller magnitude of relative vorticity compared to the neighboring flows, allowing the water masses to leave the core (Early et al., 2011;McWilliams, 1985). However, this topic requires further analysis because eddy survival also depends on other factors such as its spatial characteristics, strain, among others (McWilliams, 1985).
Here, only one case of Puddy formation was analyzed in each simulation. We cannot confirm that these cases correspond to the formation mechanism dominant in the region. Long-term high-resolution simulations (at least 10 years and <1 km) are needed to confirm the representativeness of the formation mechanisms described, as well as assess their mean properties and seasonal dependence and estimate the role of internal variability (dispersion). These simulations can also contribute to describe the impact of the Puddy formation on exchange of water mass between the coast and open ocean, which is not only produced by the eddy-induced across-shelf velocities Thomsen et al., 2016) but also by the mesoscale stirring that produce small-scale structures (Thomsen et al., 2016).
We now discuss the implications of our results for understanding some aspects of the circulation off central Chile. First, the PCUC variability is tightly linked to remote equatorial forcing through intraseasonal (between 20 and 80 days) coastal-trapped waves (CTWs; Shaffer et al., 1997;Pizarro et al., 2002;Combes et al., 2015). Our results question the extent to which CTWs could have a role in the temporal variability of Puddy formation; the downwelling CTWs would intensify the undercurrent for periods that are similar to the duration of Puddy formation. Nearshore negative wind stress curls can also force the PCUC (Chaigneau et al., 2013), which can intensify the PCUC during the austral summer-autumn (Combes et al., 2015). We expect that Puddy formation is more probable during this period. This issue could be addressed from long-term simulations with our model setup.
The PCUC has been observed to intensify during El Niño events, possibly because of the downwelling CTWs that are generated by equatorial Kelvin waves (Colas et al., 2008;Dewitte et al., 2012). Combes et al. (2015) showed that the number of Puddies 250 days after the El Niño peak (when the sea level was higher) decreased along the coasts of Peru and Chile. This decrease could be associated with a weakening of the PCUC that was observed seven months after the El Niño peak (Pizarro et al., 2001). This weakening would be produced by a geostrophic current toward the equator that was generated by the offshore propagation of the sea level because of the radiation of extratropical Rossby waves during El Niño (Ramos et al., 2008). This issue regarding the equatorial control of Puddy formation is also relevant for understanding the OMZ because Puddies can transport coastal biogeochemical properties to the open ocean, particularly dissolved oxygen. For instance, the advection of nitrogen-deficient waters (which represents the accumulated nitrogen loss of the past; Gruber & Sarmiento, 1997) can occur from the bottom boundary layer of the slope into the open ocean during Puddy formation (Thomsen et al., 2016). Thomsen et al. (2016) suggested that Puddies might be crucial in resupplying NO − 3 to the productive continental margin at depth, where nitrogen loss and organic matter exportation are thought to be highest (Kalvelage et al., 2013). Puddies are thus a conduit by which the OMZ can fluctuate at a variety of timescales, motivating the investigation of the mechanisms by which Puddies can be influenced by climate models such as the ENSO. This topic is a direction for future work. Technology Council (CONICYT). We are thankful for the supported by FONDECYT Grant 1181872 and the Millennium Science Initiative (IC 120019). We thank Jonathan Gula for help with section and the analysis of bottom pressure anomalies (Supporting Information). Cruise data used in Figures 3 and 4 are freely available under request to the Chilean National Center for Hydrographic and Oceanographic Data (Centro Nacional de Datos Hidrográficos y Oceanográficos de Chile; http://www. shoa.cl/n_cendhoc/). CROCO model is available at http://www.croco-ocean. org. All input data set and configuration of our CROCO simulations are described in section . Powered@NLHPC: This research was partially supported by the supercomputing infrastructure of the NLHPC (ECM-02).