RESEARCH ARTICLE

10.1029/2018JC014723

**Key Points:**

* Off central Chile subsurface anticyclonic eddies are generated when the Peru-Chile Undercurrent separates from the slope

* After the undercurrent separates, centrifugal instabilities are generated; these instabilities then coalesce to form anticyclonic eddies

* Centrifugal instability is not formed in low-resolution simulations, but EKE is still produced by horizontal Reynolds stresses

**Supporting Information:**

* Supporting Information S1

**Correspondence to:**

M. Contreras,
marcecontreras@udec.cl

**Citation:**

Contreras, M., Pizarro, O., Dewitte, B., Sepulveda, H. H., & Renault, L. (2019). Subsurface mesoscale eddy generation in the ocean off central Chile. *Journal of Geophysical Research: Oceans, 124*, 5700–5722. https://doi.org/10.1029/2018JC014723

Received 29 OCT 2018

Accepted 25 JUL 2019

Accepted article online 31 JUL 2019

Published online 14 AUG 2019

<figure><img src="image_1.png" /></figure>

# Subsurface Mesoscale Eddy Generation in the Ocean off Central Chile

Marcela Contreras¹,² ID, Oscar Pizarro²,³ ID, Boris Dewitte⁴,⁵,⁶,⁷ ID, Hector H. Sepulveda³ ID, and Lionel Renault⁷ ID

¹Graduate Program in Oceanography, Department of Oceanography, University of Concepción, Concepción, Chile,
²Millennium Institute of Oceanography, Concepción, Chile, ³Department of Geophysics, University of Concepcion,
Concepción, Chile, ⁴Centro de Estudios Avanzados en Zonas Áridas (CEAZA), Coquimbo, Chile, ⁵Facultad de Ciencias
del Mar, Universidad Católica del Norte, Coquimbo, Chile, ⁶Millennium Nucleus Ecology and Sustainable Management of
Oceanic Island (ESMOI), Coquimbo, Chile, ⁷LEGOS, Université de Toulouse, IRD, CNRS, CNES, UPS, Toulouse, France

**Abstract** Off the coast of central Chile, subsurface anticyclonic eddies are a salient feature of the oceanic circulation, transporting a significant 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 intensification of anticyclonic vorticity. It reaches magnitudes that are larger than the planetary vorticity while kinetic energy is converted from the PCUC to the eddy flow. 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.

**Plain Language Summary** Subsurface mesoscale eddies are swirling masses of water observed below the surface layer of the ocean (around 100- to 400-m depth). Off central Chile, these eddies have typical diameters of few tens of kilometers. They are formed near the coast, where an intense subsurface poleward flow, namely, the Peru-Chile Undercurrent (PCUC), interacts with the continental slope and the seaward border of the continental shelf. These eddies can travel long distances, toward the open ocean, transporting coastal waters with low dissolved oxygen and high nutrient concentrations and impacting the regional marine ecosystems. We use a high-resolution numerical oceanic model (~0.3 km) to analyze the formation of an eddy near 33.5°S off Chile. We showed that the eddy formation process requires the undercurrent to destabilize and detach from the coast, promoting the generation of submesoscale eddies (diameters ~10 km). This means that in regions of eddy formation, initially, the PCUC drifts offshore transferring momentum to submesoscale eddies. Later on, these eddies begin to coalesce to form an eddy with larger dimensions. We also showed that the spatial resolution of the numerical model can impact the mechanism of transfer of momentum and the vertical structure of eddies.

## 1. 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 (Chelton et al., 2011; 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

<footnote>©2019. The Authors.
This is an open access article under the
terms of the Creative Commons
Attribution License, which permits use,
distribution and reproduction in any
medium, provided the original work is
properly cited.</footnote>
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 (Colas et al., 2012; 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 (Colas et al., 2012; 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).
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.

## 2. Methodology and Data Description

### 2.1. Data

#### 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.

### 2.2. Model Configuration and Forcings

In this study, we used the Regional Ocean Modelling Systems (ROMS; Shchepetkin & McWilliams, 2003, 2005), 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. (2016, 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
<figure><img src="image_2.png" /></figure>

**Figure 1.** Topography (m) of the domain of the simulation L1 (~3-km resolution) from the SRTM15_PLUS data set. The areas that are delineated by thick yellow lines show the domain of simulations L2 (~1-km resolution) and L3 (~0.3-km resolution). The 200- and 1,000-m isobaths are delineated by black lines. The red dots show the stations of an oceangoing cruise carried on in December 2005 used to evaluate model temperature and salinity. The bottom left corner of the figure shows the L0 domain (parent simulation with 1/12° resolution), with the areas of the child simulations delineated by blue lines.

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 = \Delta h / \bar{h} = 0.2$, where $h$ is the depth and $\bar{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 $\sigma_b = 2$ and $\sigma_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
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Mean sea surface temperature (°C) between 14 January and 14 May 2006 from Multi-scale Ultra-high Resolution (MUR) data (a) and L3 (0.3 km, b). (c) Difference between MUR and L3.</figcaption></figure>

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:

$$Ur = Ua - (1 - s_w)Uo, \qquad (1)$$

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 (Renault, Hall, & McWilliams, 2016). 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
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Vertical sections of the temperature (shaded, °C) and density (sigma-t, black contours) during December 2005 at 35°30′S (upper) and 36°30′S (bottom): from (left) in situ data, (center) a model (L2 simulation), and (right) the difference between the observations and model.</figcaption></figure>

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

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Vertical sections of the salinity (shaded, psu) and density (sigma-t, black contours) during December 2005 at 35°30′S (upper) and 36°30′S (bottom): from (left) in situ data, (center) a model (L2 simulation), and (right) difference between the observations and model.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Vertical sections of the meridional velocity (a–d), vertical component of the relative vorticity (normalized by the Coriolis parameter; e–h), and the product of the potential vorticity and Coriolis parameter (i–l) at (1) 33°S, (2) 33°30'S, (3) 34°S, and (4) 34°36'S. The scale of the x axis in each section is different. The gray contours show the meridional velocity at 0, -5, and -10 cm/s (a–d), the relative vorticity (normalized by the Coriolis parameter) at -0.33 (e–h), and the product of the potential vorticity and Coriolis parameter at 0.8 s⁻⁴ (i–l). (m) Meridional velocity at 150-m depth. The black lines show the location of the vertical sections. The arrows show the velocity at 150-m depth. These variables were averaged between 14 January and 14 May 2006.</figcaption></figure>

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.

### 3. Anticyclonic Vorticity Along the Slope

Modeling studies in the eastern South Pacific showed that most subsurface mesoscale eddies are anticyclonic (Colas et al., 2012; 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 ($\partial V/\partial x$) and vertical ($\partial V/\partial z$) velocity shear components. In the above expressions, *V* is the southward (along-shore) 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 ($\partial V/\partial 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
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). Here, $PV = (\zeta + f z) \times \nabla b$, where $\zeta = \nabla \times u$ is the relative vorticity; z is a unit vertical vector; $u = (u, v, w)$ is the velocity vector with components in the zonal, meridional, and vertical axes; and $b = -g\rho/\rho_o$, g, $\rho$, and $\rho_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.

## 4. 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., $\zeta^2/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.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Temporal evolution of the vertical component of the relative vorticity (normalized by the Coriolis parameter) at 150-m depth, during Puddy formation. The black contours indicate the 200- and 500-m isobaths. Figure 6e shows a different region than that in a–d. The region in a–d is delineated by the rectangle with thin black lines in Figure 6e. The horizontal green lines indicate the latitude of the vertical sections in Figure 9. The circles in the dotted black lines in a–d indicate the zone of Puddy formation, while the circles in the dotted black line in e indicate the location of a well-developed Puddy. The contours with yellow lines correspond to the vertical component of the relative vorticity (normalized by the Coriolis parameter) and equal to −1.</figcaption></figure>

**4.1. 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 (Molemaker et al., 2015; Thomsen et al., 2016). Different factors can contribute to coastal current separation: the so-called $\beta$ 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 (Molemaker et al., 2015). 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
<figure><img src="image_8.png" /><figcaption>**Figure 7.** (a) Vertical component of the relative vorticity (normalized by the Coriolis parameter) at 150-m depth. The area that is delineated by black lines corresponds to the region in Figure 6e. The black thicker lines denote the locations of the vertical sections of Puddies E1 (b–f) and E2 (g–k). The vertical sections show the salinity (b, g), temperature (c, h), meridional velocity (d, i), vertical component of the relative vorticity (normalized by the Coriolis parameter; e, j), and the product of the potential vorticity and the Coriolis parameter (f, k). The black contours denote the potential density, with intervals of 0.1 kg/m³. Puddy E1 corresponds to the formation that was analyzed. The black line with dots shows the track since the start of E1 formation to 14 May 2006.</figcaption></figure>

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 *f*PV 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 *f*PV 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 *f*PV (Figure 9k). The region that was dominated by negative *f*PV extended beyond the boundary layer. When the *f*PV 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 *f*PV) was observed in the latitudinal
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Temporal evolution of the meridional velocity at 150-m depth, during the analyzed Puddy formation. The arrows show the velocity at 150-m depth. Black contours show the 200- and 500-m isobaths. The green lines denote the locations of the vertical sections in Figure 9. The black circles indicate the zone of Puddy formation.</figcaption></figure>

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 *f*PV magnitudes close to zero. This characteristic is typical of Puddy formation, which allows the eddy core to contain relatively homogeneous water. The origin of *f*PV close to zero could be directly associated with friction in the bottom boundary layer (section 3.1) or the effects of submesoscale instabilities.

## 4.2. 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, v, w) = (u, v, w) + (u', v', w') and b = b + b', where (·) 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):

$$ \text{HRS} = -\overline{u'}^2 \frac{\partial \overline{u}}{\partial x} - \overline{u'} \overline{v'} \frac{\partial \overline{u}}{\partial y} - \overline{v'}^2 \frac{\partial \overline{v}}{\partial y} - \overline{u'} \overline{v'} \frac{\partial \overline{v}}{\partial x}, \quad (2a) $$

$$ \text{VRS} = -\overline{u'w'} \frac{\partial \overline{u}}{\partial z} - \overline{v'w'} \frac{\partial \overline{v}}{\partial z}, \quad (2b) $$
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Vertical sections at 33.4°S (green lines in Figures 6 and 8) that show the temporal evolution of the meridional velocity (a–e), vertical component of the relative vorticity (normalized by the Coriolis parameter) (f–j), and product of the potential vorticity and the Coriolis parameter (k–o) during the Puddy formation. The black contours denote the potential density, with intervals of 0.1 kg/m³.</figcaption></figure>

$$VBF = \overline{w'b'}. \qquad (2c)$$

The overbars in (2) also indicate the average over the period between 30 March and 30 April 2006. If the HRS, VRS, and VBF are positive, kinetic energy is transferred from the mean current (MKE)—or from eddy available potential energy (EPE) in the case of VBF—to EKE. Conversely, if the magnitudes of the HRS and VRS (VBF) are negative, then the transfer is from EKE to MKE (EPE).

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⁻² and −1.36 · 10 cm³/s³, 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 (Dewar et al., 2015). In the surface layer, the main transfer of energy was from EPE to EKE
<figure><img src="image_11.png" /><figcaption>**Figure 10.** (a) HRS, (b) VRS multiplied by 10, and (c) VBF integrated between 100- and 500-m depths. The terms were calculated considering the decomposition of the variables as the sum of the temporal mean over the period between 30 March and 30 April 2006 and the fluctuation with respect to this mean. The black contours show the 200- and 500-m isobaths.</figcaption></figure>

(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 anticlonic vorticity, whose magnitudes could exceed the planetary vorticity (Figure 6), and negative *f*PV (Figure 11). These regions are indicated in Figure 11 by rectangles in red lines. If *f*PV ≤ 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
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Temporal evolution of HRS\* integrated between 100- and 500-m depths during the Puddy formation. The thin black contours show the 200- and 500-m isobaths. The black circles show the zone of Puddy formation. The meridional velocity and potential vorticity (multiplied by the Coriolis parameter) at 150-m depth are shown in thick black (15 cm/s) and thin green (0) contours, respectively. The rectangles with red lines indicate the regions where centrifugal instabilities developed.</figcaption></figure>

topography has been studied numerically in the California Undercurrent (Dewar et al., 2015; 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 (Dewar et al., 2015) and modifying the *f*PV 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 (Molemaker et al., 2015). The eddy has close to zero *f*PV at its core because centrifugal instability increases the anticyclonic relative vorticity and it promotes vertical mixing that weakens the stratification (Dewar et al., 2015; 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.

**5. Sensitivity to Model Resolution and Topography**

**5.1. 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
<figure><img src="image_13.png" /><figcaption>**Figure 12.** PDFs of the (a) meridional velocity, (b) vertical component of the relative vorticity (normalized by the Coriolis parameter), and (c) potential vorticity multiplied by the Coriolis parameter at 150-m depth. These PDFs were computed between 32°33′S and 33°51′S and from the 150-m to 500-m isobaths for the period between 14 January and 14 May 2006. The green, orange, and blue lines show the PDFs from L1s (~3-km resolution and low-resolution bottom topography), L3s (~0.3-km resolution and low-resolution bottom topography), and L3 (~0.3-km resolution and high-resolution bottom topography), respectively. The standard errors of the mean are indicated in parentheses on top of each panel and are estimated by using a bootstrap procedure (Efron & Tibshirani, 1994).</figcaption></figure>

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 *f*PV. 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 *ζ*²/*f* < −1), while L1s only had values smaller than |*f*|. Notably, the *ζ*²/*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 *f*PV was predominantly positive in all the simulations. In fact, L1 and L1s showed only positive *f*PV values. In contrast, L3s and L3 showed some negative *f*PV values, particularly L3 (Figure 12c). Regions with negative *f*PV 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 *f*PV. In summary, higher model resolution and a more realistic slope contribute to the generation of large anticyclonic relative vorticity and reduced *f*PV, which contribute to the formation of centrifugal instabilities.

## 5.2. 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
<figure><img src="image_14.png" /><figcaption>**Figure 13.** (a–d) Temporal evolution of the vertical component of the relative vorticity (normalized by the Coriolis parameter) at 150-m depth during the Puddy formation from simulation L1s (~3-km resolution and low-resolution bottom topography). The arrows indicate the velocity at 150-m depth. The green circles indicate the zone of Puddy formation. The black line in (d) denotes the location of the vertical sections of the (e) salinity, (f) temperature, (g) meridional velocity, (h) product of the potential vorticity and the Coriolis parameter, and (i) vertical component of the relative vorticity (normalized by the Coriolis parameter). The contours in (e) to (i) indicate the potential density. The contour interval is every 0.1 kg/m³.</figcaption></figure>

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 *f*PV 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⁻² and 4.75 · 10⁻¹ cm³/s³, 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.
<figure><img src="image_15.png" /><figcaption>**Figure 14.** (a) HRS, (b) VRS multiplied by 10, and (c) VBF integrated between 100- and 500-m depth. The terms were calculated from L1s (~3-km resolution and low-resolution bottom topography) between 15 March and 15 April 2006.</figcaption></figure>

### 5.3. 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

$$S = \sqrt{S_s^2 + S_n^2}, \qquad (3a)$$

$$W = S_s^2 + S_n^2 - \zeta^2, \qquad (3b)$$

where $S_s = \partial v/\partial x + \partial u/\partial y$ is the shearing deformation rate and $S_n = \partial u/\partial x - \partial v/\partial 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. Puddies from L1s and L3 present similar diameters (36.74 and 40.45 km, respectively). The shallower and deeper boundaries of Puddies from L1s (L3) were -60.75 (-99.35) m, and -274.94 (-344.87) m, respectively. The swirl velocities ($V_s$) were calculated by averaging the magnitudes of velocities ($[u^2 + v^2]^{1/2}$) along the $10^{-12} s^{-2} W$ contour. The swirl velocities of Puddies from L1s and L3 were 23.72 and 21.49 m/s. The rotation periods (estimated as $T = 2\pi R/V_s$ where R is the radio) of the Puddy from L1s and L3 were 5.89 and 6.97 days. We did not estimate the properties of the Puddies that involve their trajectories (e.g., transport and life cycle) because the period and domain of L3 simulation are not enough to analyze the Puddy displacement.

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.

**6. Summary and Discussion**

In the previous literature, off central Chile, Puddies have been described using in situ observation from ocea-
nographic cruise (Hormazabal et al., 2013) and Argo floats (Johnson & McTaggart, 2010), and numerical
models (Colas et al., 2012; Combes et al., 2015). These studies suggested that the Puddy formation is asso-
ciated 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 destabilili-
zation 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 for-
mation 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 forma-
tion were similar to those for other regions (Molemaker et al., 2015; 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 instabil-
ities. 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).
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 ($\zeta^{z} \sim -0.2 f$) and higher magnitudes of potential vorticity ($fpV \sim 0.5 \cdot 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 (Colas et al., 2012; 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⁻₃ 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 time-scales, 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.

## Acknowledgments

Marcela Contreras was supported by a master scholarship in Chile from National Chilean Research and 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).

## References

Astudillo, O., Dewitte, B., Mallet, M., Frappart, F., Rutllant, J. A., Ramos, M., & Illig, S. (2017). Surface winds off Peru-Chile: Observing closer to the coast from radar altimetry. *Remote Sensing of Environment*, **191**, 179-196. https://doi.org/10.1016/j.rse.2017.01.010

Barceló-Llull, B., Pallàs-Sanz, E., Sangrà, P., Martinez-Marrero, A., Estrada-Allis, S. N., & Aristegui, J. (2017). Ageostrophic secondary circulation in a subtropical intrathermocline eddy. *Journal of Physical Oceanography*, **47**(5), 1107-1123. https://doi.org/10.1175/JPO-D-16-0235.1

Barkan, R., Winters, K. B., & McWilliams, J. C. (2017). Stimulated imbalance and the enhancement of eddy kinetic energy dissipation by internal waves. *Journal of Physical Oceanography*, **47**(1), 181-198. https://doi.org/10.1175/JPO-D-16-0117.1

Benthuysen, J., & Thomas, L. N. (2012). Friction and diapycnal mixing at a slope: Boundary control of potential vorticity. *Journal of Physical Oceanography*, **42**(9), 1509-1523. https://doi.org/10.1175/JPO-D-11-0130.1

Bosse, A., Testor, P., Houpert, L., Damien, P., Prieur, L., Hayes, D., & Karstensen, J. (2016). Scales and dynamics of Submesoscale Coherent Vortices formed by deep convection in the northwestern Mediterranean Sea. *Journal of Geophysical Research: Oceans*, **121**, 7716-7742. https://doi.org/10.1002/2016JC012144

Boyd, J. P. (1992). The energy spectrum of fronts: Time evolution of shocks in Burgers' equation. *Journal of the Atmospheric Sciences*, **49**(2), 128-139. https://doi.org/10.1175/1520-0469(1992)049<0128:TESOFT>2.0.CO;2

Bühler, O., & McIntyre, M. E. (2005). Wave capture and wave-vortex duality. *Journal of Fluid Mechanics*, **534**, 67-95. https://doi.org/10.1017/S0022112005004374

Callies, J., & Ferrari, R. (2013). Interpreting energy and tracer spectra of upper-ocean turbulence in the submesoscale range (1–200 km). *Journal of Physical Oceanography*, **43**(11), 2456-2474. https://doi.org/10.1175/JPO-D-13-063.1

Capet, X., McWilliams, J. C., Molemaker, M. J., & Shchepetkin, A. F. (2008). Mesoscale to submesoscale transition in the California Current System. Part I: Flow structure, eddy flux, and observational tests. *Journal of Physical Oceanography*, **38**(1), 29-43. https://doi.org/10.1175/2007JPO3671.1

Capet, X. J., Marchesiello, P., & McWilliams, J. C. (2004). Upwelling response to coastal wind profiles. *Geophysical Research Letters*, **31**, L13311. https://doi.org/10.1029/2004GL020123

Chaigneau, A., Dominguez, N., Eldin, G., Vasquez, L., Flores, R., Grados, C., & Echevin, V. (2013). Near-coastal circulation in the Northern Humboldt Current System from shipboard ADCP data. *Journal of Geophysical Research: Oceans*, **118**, 5251-5266. https://doi.org/10.1002/jgrc.20328

Chaigneau, A., Eldin, G., & Dewitte, B. (2009). Eddy activity in the four major upwelling systems from satellite altimetry (1992-2007). *Progress in Oceanography*, **83**(1-4), 117-123. https://doi.org/10.1016/j.pocean.2009.07.012

Chaigneau, A., & Pizarro, O. (2005). Eddy characteristics in the eastern South Pacific. *Journal of Geophysical Research*, **110**, C06005. https://doi.org/10.1029/2004JC002815

Charney, J. G. (1971). Geostrophic turbulence. *Journal of the Atmospheric Sciences*, **28**(6), 1087-1095. https://doi.org/10.1175/1520-0469(1971)028<1087:GT>2.0.CO;2

Chelton, D. B., Schlax, M. G., & Samelson, R. M. (2011). Global observations of nonlinear mesoscale eddies. *Progress in Oceanography*, **91**(2), 167-216. https://doi.org/10.1016/j.pocean.2011.01.002

Chin, T. M., Milliff, R. F., & Large, W. G. (1998). Basin-scale, high-wavenumber sea surface wind fields from a multiresolution analysis of scatterometer data. *Journal of Atmospheric and Oceanic Technology*, **15**(3), 741-763. https://doi.org/10.1175/1520-0426(1998)015<0741:BSHWSS>2.0.CO;2

Colas, F., Capet, X., McWilliams, J. C., & Shchepetkin, A. (2008). 1997-1998 El Niño off Peru: A numerical study. *Progress in Oceanography*, **79**(2-4), 138-155. https://doi.org/10.1016/j.pocean.2008.10.015

Colas, F., McWilliams, J. C., Capet, X., & Kurian, J. (2012). Heat balance and eddies in the Peru-Chile current system. *Climate Dynamics*, **39**(1-2), 509-529. https://doi.org/10.1007/s00382-011-1170-6

Collins, C. A., Margolina, T., Rago, T. A., & Ivanov, L. (2013). Looping RAFOS floats in the California Current System. *Deep Sea Research Part II: Topical Studies in Oceanography*, **85**, 42-61. https://doi.org/10.1016/j.dsr2.2012.07.027

Combes, V., Hormazabal, S., & Di Lorenzo, E. (2015). Interannual variability of the subsurface eddy field in the Southeast Pacific. *Journal of Geophysical Research: Oceans*, **120**, 4907-4924. https://doi.org/10.1002/2014JC010265

da Silva, A. M., Young, C. C., & Levitus, S. (1994). *Atlas of Surface Marine Data 1994, Vol. 1: Algorithms and Procedures*, NOAA Atlas NESDIS, **6** (p. 83), 20910-3282. Washington, DC: U.S. Government Printing Office.

D'Asaro, E. A. (1988). Generation of submesoscale vortices: A new mechanism. *Journal of Geophysical Research*, **93**(C6), 6685-6693. https://doi.org/10.1029/JC093iC06p06685

Debreu, L., Marchesiello, P., Penven, P., & Cambon, G. (2012). Two-way nesting in split-explicit ocean models: Algorithms, implementation and validation. *Ocean Modelling*, **49-50**, 1-21. https://doi.org/10.1016/j.ocemod.2012.03.003

Dewar, W. K., McWilliams, J. C., & Molemaker, M. J. (2015). Centrifugal instability and mixing in the California Undercurrent. *Journal of Physical Oceanography*, **45**(5), 1224-1241. https://doi.org/10.1175/JPO-D-13-0269.1

Dewitte, B., Vazquez-Cuervo, J., Goubanova, K., Illig, S., Takahashi, K., Cambon, G., & Ortlieb, L. (2012). Change in El Nino flavours over 1958-2008: Implications for the long-term trend of the upwelling off Peru. *Deep Sea Research Part II: Topical Studies in Oceanography*, **77-80**, 143-156. https://doi.org/10.1016/j.dsr2.2012.04.011

Dong, C., McWilliams, J. C., Liu, Y., & Chen, D. (2014). Global heat and salt transports by eddy movement. *Nature Communications*, **5**(1), 3294. https://doi.org/10.1038/ncomms4294

Dufois, F., Penven, P., Whittle, C. P., & Veitch, J. (2012). On the warm nearshore bias in Pathfinder monthly SST products over Eastern Boundary Upwelling Systems. *Ocean Modelling*, **47**, 113-118. https://doi.org/10.1016/j.ocemod.2012.01.007
Early, J. J., Samelson, R. M., & Chelton, D. B. (2011). The evolution and propagation of quasigeostrophic ocean eddies. *Journal of Physical Oceanography*, 41(8), 1535–1555. https://doi.org/10.1175/2011JPO4601.1

Efron, B., & Tibshirani, R. J. (1994). *An introduction to the bootstrap*. New York: Chapman & Hall/CRC Press.

Fairall, C. W., Bradley, E. F., Hare, J. E., Grachev, A. A., & Edson, J. B. (2003). Bulk parameterization of air–sea fluxes: Updates and verification for the COARE algorithm. *Journal of Climate*, 16(4), 571–591. https://doi.org/10.1175/1520-0442(2003)016<0571: BPOASF>2.0.CO;2

Ferrari, R., & Wunsch, C. (2009). Ocean circulation kinetic energy: Reservoirs, sources, and sinks. *Annual Review of Fluid Mechanics*, 41(1), 253–282. https://doi.org/10.1146/annurev.fluid.40.111406.102139

Frenger, I., Bianchi, D., Stührenberg, C., Oschlies, A., Dunne, J., Deutsch, C., & Schütte, F. (2018). Biogeochemical role of subsurface coherent eddies in the ocean: Tracer cannonballs, hypoxic storms, and microbial stewpots? *Global Biogeochemical Cycles*, 32, 226–249. https://doi.org/10.1002/2017GB005743

Fuenzalida, R., Schneider, W., Garcés-Vargas, J., Bravo, L., & Lange, C. (2009). Vertical and horizontal extension of the oxygen minimum zone in the eastern South Pacific Ocean. *Deep Sea Research Part II: Topical Studies in Oceanography*, 56(16), 992–1003. https://doi.org/10.1016/j.dsr2.2008.11.001

Goubanova, K., Echevin, V., Dewitte, B., Codron, F., Takahashi, K., Terray, P., & Vrac, M. (2011). Statistical downscaling of sea-surface wind over the Peru–Chile upwelling region: Diagnosing the impact of climate change from the IPSL-CM4 model. *Climate Dynamics*, 36(7-8), 1365–1378. https://doi.org/10.1007/s00382-010-0824-0

Gruber, N., Lachkar, Z., Frenzel, H., Marchesiello, P., Münnich, M., McWilliams, J. C., & Plattner, G. K. (2011). Eddy-induced reduction of biological production in eastern boundary upwelling systems. *Nature Geoscience*, 4(11), 787–792. https://doi.org/10.1038/ngeo1273

Gruber, N., & Sarmiento, J. L. (1997). Global patterns of marine nitrogen fixation and denitrification. *Global Biogeochemical Cycles*, 11(2), 235–266. https://doi.org/10.1029/97GB00077

Gula, J., Molemaker, M. J., & McWilliams, J. C. (2015). Gulf Stream dynamics along the southeastern US seaboard. *Journal of Physical Oceanography*, 45(3), 690–715. https://doi.org/10.1175/JPO-D-14-0154.1

Gula, J., Molemaker, M. J., & McWilliams, J. C. (2016). Topographic generation of submesoscale centrifugal instability and energy dissipation. *Nature Communications*, 7(1), 12,811. https://doi.org/10.1038/ncomms12811

Harrison, D. E., & Robinson, A. R. (1978). Energy analysis of open regions of turbulent flows—Mean eddy energetics of a numerical ocean circulation experiment. *Dynamics of Atmospheres and Oceans*, 2(2), 185–211. https://doi.org/10.1016/0377-0265(78)90009-X

Hide, R., & Titman, C. W. (1967). Detached shear layers in a rotating fluid. *Journal of Fluid Mechanics*, 29(1), 39–60. https://doi.org/10.1017/S002211206700062X

Hormazabal, S., Combes, V., Morales, C. E., Correa-Ramirez, M. A., Di Lorenzo, E., & Nuñez, S. (2013). Intrathermocline eddies in the coastal transition zone off central Chile (31–41 S). *Journal of Geophysical Research: Oceans*, 118, 4811–4821. https://doi.org/10.1002/jgrc.20337

Hoskins, B. J. (1974). The role of potential vorticity in symmetric stability and instability. *Quarterly Journal of the Royal Meteorological Society*, 100(425), 480–482. https://doi.org/10.1002/qj.49710042520

Jerónimo, G., & Gómez-Valdés, J. (2007). A subsurface warm-eddy off northern Baja California in July 2004. *Geophysical Research Letters*, 34, L06610. https://doi.org/10.1029/2006GL028851

Jiao, Y., & Dewar, W. K. (2015). The energetics of centrifugal instability. *Journal of Physical Oceanography*, 45(6), 1554–1573. https://doi.org/10.1175/JPO-D-14-0064.1

Johnson, G. C., & McTaggart, K. E. (2010). Equatorial Pacific 13 C water eddies in the eastern subtropical South Pacific Ocean. *Journal of Physical Oceanography*, 40(1), 226–236. https://doi.org/10.1175/2009JPO4287.1

Kalvelage, T., Lavik, G., Lam, P., Contreras, S., Arteaga, L., Löscher, C. R., & Kuypers, M. M. (2013). Nitrogen cycling driven by organic matter export in the South Pacific oxygen minimum zone. *Nature Geoscience*, 6(3), 228–234. https://doi.org/10.1038/ngeo1739

Kang, D., & Curchitser, E. N. (2015). Energetics of eddy–mean flow interactions in the Gulf Stream region. *Journal of Physical Oceanography*, 45(4), 1103–1120. https://doi.org/10.1175/JPO-D-14-0200.1

Klein, P., Hua, B. L., Lapeyre, G., Capet, X., Le Gentil, S., & Sasaki, H. (2008). Upper ocean turbulence from high-resolution 3D simulations. *Journal of Physical Oceanography*, 38(8), 1748–1763. https://doi.org/10.1175/2007JPO3773.1

Klein, P., & Lapeyre, G. (2009). The oceanic vertical pump induced by mesoscale and submesoscale turbulence. *Annual Review of Marine Science*, 1(1), 351–375. https://doi.org/10.1146/annurev.marine.010908.163704

Kurian, J., Colas, F., Capet, X., McWilliams, J. C., & Chelton, D. B. (2011). Eddy properties in the California current system. *Journal of Geophysical Research*, 116, C08027. https://doi.org/10.1029/2010JC006895

Landolfi, A., Dietze, H., Koeve, W., & Oschlies, A. (2013). Overlooked runaway feedback in the marine nitrogen cycle: The vicious cycle. *Biogeosciences*, 10(3), 1351–1363. https://doi.org/10.5194/bg-10-1351-2013

Large, W. G., McWilliams, J. C., & Doney, S. C. (1994). Oceanic vertical mixing: A review and a model with a nonlocal boundary layer parameterization. *Reviews of Geophysics*, 32(4), 363–403. https://doi.org/10.1029/94RG01872

Lemarié, F., Kurian, J., Shchepetkin, A. F., Molemaker, M. J., Colas, F., & McWilliams, J. C. (2012). Are there inescapable issues prohibiting the use of terrain-following coordinates in climate models? *Ocean Modelling*, 42, 57–79. https://doi.org/10.1016/j.ocemod.2011.11.007

Leth, O., & Shaffer, G. (2001). A numerical study of the seasonal variability in the circulation off central Chile. *Journal of Geophysical Research*, 106(C10), 22,229–22,248. https://doi.org/10.1029/2000JC000627

Marchesiello, P., Debreu, L., & Couvelard, X. (2009). Spurious diapycnal mixing in terrain-following coordinate models: The problem and a solution. *Ocean Modelling*, 26(3-4), 156–169. https://doi.org/10.1016/j.ocemod.2008.09.004

Marchesiello, P., McWilliams, J. C., & Shchepetkin, A. (2003). Equilibrium structure and dynamics of the California Current System. *Journal of Physical Oceanography*, 33(4), 753–783. https://doi.org/10.1175/1520-0485(2003)33<753:ESADOT>2.0.CO;2

Marshall, D. P., & Tansley, C. E. (2001). An implicit formula for boundary current separation. *Journal of Physical Oceanography*, 31(6), 1633–1638. https://doi.org/10.1175/1520-0485(2001)031<1633:AIFFBC>2.0.CO;2

Mason, E., Colas, F., Molemaker, J., Shchepetkin, A. F., Troupin, C., McWilliams, J. C., & Sangrà, P. (2011). Seasonal variability of the Canary Current: A numerical study. *Journal of Geophysical Research*, 116, C06001. https://doi.org/10.1029/2010JC006665

Mason, E., Molemaker, J., Shchepetkin, A. F., Colas, F., McWilliams, J. C., & Sangrà, P. (2010). Procedures for offline grid nesting in regional ocean models. *Ocean Modelling*, 35(1-2), 1–15. https://doi.org/10.1016/j.ocemod.2010.05.007

McGillicuddy, D. J., Anderson, L. A., Bates, N. R., Bibby, T., Buesseler, K. O., Carlson, C. A., & Hansell, D. A. (2007). Eddy/wind interactions stimulate extraordinary mid-ocean plankton blooms. *Science*, 316(5827), 1021–1026. https://doi.org/10.1126/science.1136256
McWilliams, J. C. (1985). Submesoscale, coherent vortices in the ocean. *Reviews of Geophysics*, **23**(2), 165-182. https://doi.org/10.1029/RG023i002p00165

McWilliams, J. C. (2016). Submesoscale currents in the ocean. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences*, **472**(2189). https://doi.org/10.1098/rspa.2016.0117

Molemaker, M. J., McWilliams, J. C., & Dewar, W. K. (2015). Submesoscale instability and generation of mesoscale anticyclones near a separation of the California Undercurrent. *Journal of Physical Oceanography*, **45**(3), 613–629. https://doi.org/10.1175/JPO-D-13-0225.1

Okubo, A. (1970). Horizontal dispersion of floatable particles in the vicinity of velocity singularities such as convergences. *Deep-Sea Research*, **17**(3), 445–454. https://doi.org/10.1016/0011-7471(70)90059-8

Penven, P., Debreu, L., Marchesiello, P., & McWilliams, J. C. (2006). Evaluation and application of the ROMS 1-way embedding procedure to the central California upwelling system. *Ocean Modelling*, **12**(1-2), 157–187. https://doi.org/10.1016/j.ocemod.2005.05.002

Penven, P., Echevin, V., Pasapera, J., Colas, F., & Tam, J. (2005). Average circulation, seasonal cycle, and mesoscale dynamics of the Peru Current System: A modeling approach. *Journal of Geophysical Research*, **110**, C10021. https://doi.org/10.1029/2005JC002945

Pizarro, O., Clarke, A. J., & Van Gorder, S. (2001). El Niño sea level and currents along the South American coast: Comparison of observations with theory. *Journal of Physical Oceanography*, **31**(7), 1891–1903. https://doi.org/10.1175/1520-0485(2001)031<1891:ENOSLA>2.0.CO;2

Pizarro, O., Shaffer, G., Dewitte, B., & Ramos, M. (2002). Dynamics of seasonal and interannual variability of the Peru-Chile Undercurrent. *Geophysical Research Letters*, **29**(12), 1581. https://doi.org/10.1029/2002GL014790

Ramos, M., Dewitte, B., Pizarro, O., & Garric, G. (2008). Vertical propagation of extratropical Rossby waves during the 1997–1998 El Niño off the west coast of South America in a medium resolution OGCM simulation. *Journal of Geophysical Research*, **113**, C08041. https://doi.org/10.1029/2007JC004681

Renault, L., Deutsch, C., McWilliams, J. C., Frenzel, H., Liang, J. H., & Colas, F. (2016). Partial decoupling of primary productivity from upwelling in the California Current system. *Nature Geoscience*, **9**(7), 505–508. https://doi.org/10.1038/ngeo2722

Renault, L., Dewitte, B., Marchesiello, P., Illig, S., Echevin, V., Cambon, G., & Ayers, J. K. (2012). Upwelling response to atmospheric coastal jets off central Chile: A modeling study of the October 2000 event. *Journal of Geophysical Research*, **117**, C02030. https://doi.org/10.1029/2011JC007446

Renault, L., Hall, A., & McWilliams, J. C. (2016). Orographic shaping of US West Coast wind profiles during the upwelling season. *Climate Dynamics*, **46**(1-2), 273–289. https://doi.org/10.1007/s00382-015-2583-4

Renault, L., Masson, S., Oerder, V., Jullien, S., & Colas, F. (2019). Disentangling the mesoscale ocean-atmosphere interactions. *Journal of Geophysical Research: Oceans*, **124**, 2164–2178. https://doi.org/10.1029/2018JC014628

Renault, L., Molemaker, M. J., McWilliams, J. C., Shchepetkin, A. F., Lemarié, F., Chelton, D., & Hall, A. (2016). Modulation of wind work by oceanic current interaction with the atmosphere. *Journal of Physical Oceanography*, **46**(6), 1685–1704. https://doi.org/10.1175/JPO-D-15-0232.1

Richardson, P. L., Bower, A. S., & Zenk, W. (2000). A census of Meddies tracked by floats. *Progress in Oceanography*, **45**(2), 209–250. https://doi.org/10.1016/S0079-6611(99)00053-1

Shaffer, G., Pizarro, O., Djurfeldt, L., Salinas, S., & Rutlant, J. (1997). Circulation and low-frequency variability near the Chilean coast: Remotely forced fluctuations during the 1991–92 El Nino. *Journal of Physical Oceanography*, **27**(2), 217–235. https://doi.org/10.1175/1520-0485(1997)027<0217:CALFVN>2.0.CO;2

Shapiro, G. I., & Meschanov, S. L. (1991). Distribution and spreading of Red Sea Water and salt lens formation in the northwest Indian Ocean. *Deep Sea Research Part A. Oceanographic Research Papers*, **38**(1), 21–34. https://doi.org/10.1016/0198-0149(91)90052-H

Shchepetkin, A. F., & McWilliams, J. C. (2003). A method for computing horizontal pressure-gradient force in an oceanic model with a nonaligned vertical coordinate. *Journal of Geophysical Research*, **108**(C3), 3090. https://doi.org/10.1029/2001JC001047

Shchepetkin, A. F., & McWilliams, J. C. (2005). The regional oceanic modeling system (ROMS): A split-explicit, free-surface, topography-following-coordinate oceanic model. *Ocean Modelling*, **9**(4), 347–404. https://doi.org/10.1016/j.ocemod.2004.08.002

Shchepetkin, A. F., & McWilliams, J. C. (2009). Correction and commentary for “Ocean forecasting in terrain-following coordinates: Formulation and skill assessment of the regional ocean modeling system” by Haidvogel et al., J. Comp. Phys., 227, pp. 3595–3624. *Journal of Computational Physics*, **228**(24), 8985–9000. https://doi.org/10.1016/j.jcp.2009.09.002

Silva, N., Rojas, N., & Fedele, A. (2009). Water masses in the Humboldt Current System: Properties, distribution, and the nitrate deficit as a chemical water mass tracer for Equatorial Subsurface Water off Chile. *Deep Sea Research Part II: Topical Studies in Oceanography*, **56**(16), 1004–1020. https://doi.org/10.1016/j.dsr2.2008.12.013

Song, H., Miller, A. J., Cornuelle, B. D., & Di Lorenzo, E. (2011). Changes in upwelling and its water sources in the California Current System driven by different wind forcing. *Dynamics of Atmospheres and Oceans*, **52**(1-2), 170–191. https://doi.org/10.1016/j.dynatmoce.2011.03.001

Song, Y. T., & Wright, D. G. (1998). A general pressure gradient formulation for ocean models. Part II: Energy, momentum, and bottom torque consistency. *Monthly Weather Review*, **126**(12), 3231–3247. https://doi.org/10.1175/1520-0493(1998)126<3231:AGPGFF>2.0.CO;2

Stramma, L., Prince, E. D., Schmidtko, S., Luo, J., Hoolihan, J. P., Visbeck, M., & Körtzinger, A. (2012). Expansion of oxygen minimum zones may reduce available habitat for tropical pelagic fishes. *Nature Climate Change*, **2**(1), 33–37. https://doi.org/10.1038/nclimate1304

Thomas, L. N., Tandon, A., & Mahadevan, A. (2008). Submesoscale processes and dynamics. *Ocean modeling in an Eddying Regime*, **17**, 17–38. https://doi.org/10.1029/177GM04

Thomas, L. N., Taylor, J. R., Ferrari, R., & Joyce, T. M. (2013). Symmetric instability in the Gulf Stream. *Deep Sea Research Part II: Topical Studies in Oceanography*, **91**, 96–110. https://doi.org/10.1016/j.dsr2.2013.02.025

Thomsen, S., Kanzow, T., Krahmann, G., Greatbatch, R. J., Dengler, M., & Lavik, G. (2016). The formation of a subsurface anticyclonic eddy in the Peru-Chile Undercurrent and its impact on the near-coastal salinity, oxygen, and nutrient distributions. *Journal of Geophysical Research: Oceans*, **121**, 476–501. https://doi.org/10.1002/2015JC010878

Vanneste, J. (2013). Balance and spontaneous wave generation in geophysical flows. *Annual Review of Fluid Mechanics*, **45**(1), 147–172. https://doi.org/10.1146/annurev-fluid-011212-140730

Vergara, O., Dewitte, B., Montes, I., Garçon, V., Ramos, M., Paulmier, A., & Pizarro, O. (2016). Seasonal variability of the oxygen minimum zone off Peru in a high-resolution regional coupled model. *Biogeosciences*, **13**(15), 4389–4410. https://doi.org/10.5194/fg-13-4389-2016

Vergara, O., Dewitte, B., Ramos, M., & Pizarro, O. (2017). Vertical energy flux at ENSO time scales in the subthermocline of the Southeastern Pacific. *Journal of Geophysical Research: Oceans*, **122**, 6011–6038. https://doi.org/10.1002/2016JC012614

Weiss, J. (1991). The dynamics of entrophy transfer in two-dimensional hydrodynamics. *Physica D: Nonlinear Phenomena*, **48**(2-3), 273–294. https://doi.org/10.1016/0167-2789(91)90088-Q
Wenegrat, J. O., Callies, J., & Thomas, L. N. (2018). Submesoscale baroclinic instability in the bottom boundary layer. *Journal of Physical Oceanography*, 48(11), 2571–2592. https://doi.org/10.1175/JPO-D-17-0264.1

Williams, R. G., & Roussenov, V. (2003). The role of sloping sidewalls in forming potential vorticity contrasts in the ocean interior. *Journal of Physical Oceanography*, 33(4), 633–648. https://doi.org/10.1175/1520-0485(2003)33<633:TROSSI>2.0.CO;2

Zhang, Z., Tian, J., Qiu, B., Zhao, W., Chang, P., Wu, D., & Wan, X. (2016). Observed 3D structure, generation, and dissipation of oceanic mesoscale eddies in the South China Sea. *Scientific Reports*, 6(1), 24,349. https://doi.org/10.1038/srep24349

Zhang, Z., Wang, W., & Qiu, B. (2014). Oceanic mass transport by mesoscale eddies. *Science*, 345(6194), 322–324. https://doi.org/10.1126/science.1252418

Reference From the Supporting Information

Kundu, P. K., & Cohen, L. M. (1990). *Fluid mechanics*, (p. 638). Calif: Academic Press.