Oxygen Variability During ENSO in the Tropical South Eastern Pacific

The Oxygen Minimum Zone (OMZ) of the Tropical South Eastern Pacific (TSEP) is one of the most intensely deoxygenated water masses of the global ocean. It is strongly affected at interannual time scales by El Nino (EN) and La Nina (LN) due to its proximity to the equatorial Pacific. In this work, the physical and biogeochemical processes associated with the subsurface oxygen variability during EN and LN in the period 1958 – 2008 were studied using a regional coupled physical-biogeochemical model and in situ observations. The passage of intense remotely forced coastal trapped waves caused a strong deepening (shoaling) of the OMZ upper limit during EN (LN). A close correlation between the OMZ upper limit and thermocline depths was found close to the coast, highlighting the role of physical processes. The subsurface waters over the shelf and slope off central Peru had different origins depending on ENSO conditions. Offshore of the upwelling region (near 88°W), negative and positive oxygen subsurface anomalies were caused by Equatorial zonal circulation changes during LN and EN, respectively. The altered properties were then transported to the shelf and slope (above 200 m) by the Peru-Chile undercurrent. The source of nearshore oxygenated waters was located at 3°S – 4°S during neutral periods, further north (1°S – 1°N) during EN and further south (4°S – 5°S) during LN. The offshore deeper (< 200 – 300 m) OMZ was ventilated by waters originating from ~ 8°S during EN and LN. Enhanced mesoscale variability during EN also impacted OMZ ventilation through horizontal and vertical eddy fluxes. The vertical eddy flux decreased due to the reduced vertical gradient of oxygen in the surface layer, whereas horizontal eddy fluxes injected more oxygen into the OMZ through its meridional boundaries. In subsurface layers, remineralization of organic matter, the main biogeochemical sink of oxygen, was higher during EN than during LN due to oxygenation of the surface layer and in spite of a reduced primary production. Sensitivity experiments highlighted the larger impact of equatorial remote forcing with respect to local wind forcing during EN and LN.

The Oxygen Minimum Zone (OMZ) of the Tropical South Eastern Pacific (TSEP) is one of the most intensely deoxygenated water masses of the global ocean. It is strongly affected at interannual time scales by El Niño (EN) and La Niña (LN) due to its proximity to the equatorial Pacific. In this work, the physical and biogeochemical processes associated with the subsurface oxygen variability during EN and LN in the period 1958-2008 were studied using a regional coupled physical-biogeochemical model and in situ observations. The passage of intense remotely forced coastal trapped waves caused a strong deepening (shoaling) of the OMZ upper limit during EN (LN). A close correlation between the OMZ upper limit and thermocline depths was found close to the coast, highlighting the role of physical processes. The subsurface waters over the shelf and slope off central Peru had different origins depending on ENSO conditions. Offshore of the upwelling region (near 88 • W), negative and positive oxygen subsurface anomalies were caused by Equatorial zonal circulation changes during LN and EN, respectively. The altered properties were then transported to the shelf and slope (above 200 m) by the Peru-Chile undercurrent. The source of nearshore oxygenated waters was located at 3 • S−4 • S during neutral periods, further north (1 • S−1 • N) during EN and further south (4 • S−5 • S) during LN. The offshore deeper (<200-300 m) OMZ was ventilated by waters originating from ∼8 • S during EN and LN. Enhanced mesoscale variability during EN also impacted OMZ ventilation through horizontal and vertical eddy fluxes. The vertical eddy flux decreased due to the reduced vertical gradient of oxygen in the surface layer, whereas horizontal eddy fluxes injected more oxygen into the OMZ through its meridional boundaries. In subsurface layers, remineralization of organic matter, the main biogeochemical sink of oxygen, was higher during EN than during LN due to oxygenation of the surface layer. Sensitivity experiments highlighted the larger impact of equatorial remote forcing with respect to local wind forcing during EN and LN.

INTRODUCTION
Oxygen Minimum Zones (OMZ) are large water bodies of the open ocean with low concentrations of dissolved oxygen (DO < 22 µmol kg −1 ) at subsurface depths (∼50-900 m; Levin, 2003;Karstensen et al., 2008). Low DO has a strong impact on the nitrogen cycle: a significant portion (35%) of the nitrogen loss in the global ocean occurs in these regions (Devol et al., 2006) due to denitrification and anammox (Babbin et al., 2014). The oxycline, which separates the well-mixed surface waters with the lowoxygenated subsurface waters, also impacts marine life by limiting the habitat of several species in the water column (Diaz and Rosenberg, 2008;Gutiérrez et al., 2008;Bertrand et al., 2011). These features, along with the recently observed increase of their vertical extent (Stramma et al., , 2010Ito and Deutsch, 2013;Schmidtko et al., 2017;Breitburg et al., 2018) have made OMZ studies an oceanographic hot topic in the recent decades.
Complex physical and biogeochemical processes are involved in the TSEP OMZ formation and maintenance. Vertical mixing between the shallow, oxygenated, surface mixed layer and subsurface OMZ waters is reduced due to the relatively weak alongshore winds. A sharp permanent pycnocline also mitigates vertical mixing and prevents local ventilation of subsurface waters (Fiedler and Talley, 2006). Furthermore, it is located in the poorly ventilated "shadow zone" of the Subtropical South Eastern Pacific, out of reach of the ventilated water masses pathways which subduct in the thermocline (Luyten et al., 1983). Indeed, the residence time of OMZ waters can be very long: water masses may recirculate within the OMZ during more than a year (e.g., ∼15 months at 8 • S, Czeschel et al., 2011). On the other hand, as it is relatively close to the equatorial region, the northwestern boundary of the OMZ is ventilated by the equatorial current system and is sensitive to its variability. In particular, equatorial subsurface undercurrents (SSCCs) or Tsuchiya jets (Tsuchiya, 1975) transport relatively DO-rich waters eastward (Montes et al., 2014). The nearshore OMZ is then ventilated by the poleward Peru-Chile Undercurrent (PCUC), which connects the offshore equatorial region to the Peruvian coastal upwelling system (PCUS) (Zuta and Guillén, 1970;Codispoti et al., 1989;Chaigneau et al., 2013). The PCUC is recognized as the main source of the upwelling waters (Huyer et al., 1987). It is mainly fueled by the primary (2 • S−6 • S, pSSCC) and secondary (6 • S−10 • S, sSSCC) subsurface undercurrents (Montes et al., 2010), and to a much lesser extent by the Equatorial Undercurrent (2 • N−2 • S, EUC; Wyrtki, 1967). In addition, the TSEP OMZ is embedded in a region of intense mesoscale activity (Chaigneau et al., 2008) due to the presence of unstable boundary currents Colas et al., 2012). These eddies generate horizontal (e.g., Vergara et al., 2016) and vertical eddy fluxes which impact the OMZ structure.
The very high PCUS productivity also impacts the OMZ. The large amount of organic matter (OM) eventually sinks into the OMZ, and oxygen is consumed by microbial respiration during OM remineralization (e.g., Paulmier et al., 2006;Cavan et al., 2017). Overall, the temporal variability of the OMZ extension is assumed to be mainly driven by physical processes (e.g., Xu et al., 2015), but its maintenance is due to the oxygen respiration (Paulmier et al., 2006). The OMZ is variable over a wide range of temporal and spatial scales. At seasonal time scales, its variability is related to the local wind forcing, remineralization and mesoscale circulation (eddies and filaments; Thomsen et al., 2016;Vergara et al., 2016). At interannual time scales, the OMZ is mainly affected by the warm (El Niño, EN) and cold (La Niña, LN) phases of the El Niño Southern Oscillation (ENSO). During EN, intraseasonal Equatorial Kelvin Waves (IEKW) are generated by wind anomalies in the Equatorial Pacific region (e.g., Kessler et al., 1995). They propagate eastward and reach the western coasts of America. There, they trigger coastal trapped waves (CTWs) which propagate poleward and impact the vertical structure of physical and biogeochemical variables alongshore off Peru (Gutiérrez et al., 2008;Echevin et al., 2014;Graco et al., 2017) and Chile (Ulloa et al., 2001). The water column thus becomes oxygenated during EN and a deeper oxycline (>100 m) is observed (Gutiérrez et al., 2008(Gutiérrez et al., , 2016, while during LN (e.g., 1999LN (e.g., -2000, a shallower oxycline is described (Graco et al., 2007(Graco et al., , 2017. These DO changes are strongly linked with the physical variability, as both the thermocline and oxycline are displaced vertically during the passage of downwelling (warm) and upwelling (cold) CTWs Graco et al., 2017).
Besides, the DO content of the PCUC could also be modified during ENSO. As it is fed mainly by the pSSCC and the EUC during LN and by the sSSCC during EN (Montes et al., 2011), variability in DO concentrations in the equatorial region could impact differently the ventilation of the OMZ. Few modeling works have described the OMZ variability at interannual time scales. Using a regional model, Mogollón and Calil (2017) simulated the impact of the strong 1997-1998EN and 1999-2000 LN events on the OMZ. They found a deepening (∼150 m) of the OMZ during 1997-1998 EN respect to neutral condition, while a shallower OMZ (∼50 m) was observed during LN in the main coastal upwelling centers of the PCUS. Using a global model, Yang et al. (2017) found that the denitrifying SETP OMZ core volume expanded during LN (e.g., ∼+35% in 1999-2000 and shrinked during EN (e.g., ∼-75% in 1982-1983/1997-1998. Note that in these recent works, few in situ oxygen measurements were used to validate the simulations during ENSO. In the present work, we make use of a unique interannual DO in situ data set and of a regional physical-biogeochemical coupled model to study the physical and biogeochemical processes involved in the OMZ variability during EN and LN phases. Unlike previous studies, we simulated a large number of ENSO events from 1958 to 2008, which allows a comparison of the impact of extreme and moderate events. After carefully evaluating the model with in situ data, we used the model to investigate the role of the CTWs, the impact of the equatorial circulation and the ventilating eddy fluxes due to the variable mesoscale activity during EN and LN events. Finally, we performed sensitivity experiments forced by contrasted atmospheric and open boundary forcing to investigate the impact of the local (e.g., wind and heat fluxes) and remote forcing (e.g., coastal waves) on the DO changes during EN and LN.

MATERIALS AND METHODS
The ROMS-PISCES Coupled Physical-Biogeochemical Model The Regional Oceanic Modeling System model (Shchepetkin and McWilliams, 2005) was used to simulate the ocean dynamics. The ROMS-AGRIF code (version 3.1) was used (Penven et al., 2006). ROMS resolves the Primitive Equations, based on the Boussinesq approximation and hydrostatic vertical momentum balance (Shchepetkin and McWilliams, 1998). ROMS is coupled to the Pelagic Interaction Scheme for Carbon and Ecosystem Studies (PISCES) biogeochemical model to simulate the marine biological productivity and the biogeochemical cycles of carbon and main nutrients (P, N, Si, Fe; Aumont et al., 2015) as well as DO (e.g., Resplandy et al., 2012). PISCES has three non-living compartments which are the semi-labile dissolved organic matter, small sinking particles and large sinking particles, and four living compartments represented by two size classes of phytoplankton (nanophytoplankton and diatoms) and two size classes of zooplankton (microzooplankton and mesozooplankton).
This coupled model has been used to study the climatological (Echevin et al., 2008;Albert et al., 2010), intraseasonal , and interannual (Espinoza-Morriberón et al., 2017) variability of the surface productivity in the PCUS.

The Oxygen Cycle in the PISCES Model
The oxygen evolution in PISCES is computed taking into account the dynamical transport, biogeochemical processes (sources and sinks) and air-sea fluxes, as follows: Biogeochemical processes + J flux air−ocean exchange u H and w represent the horizontal (zonal and meridional) and vertical currents, respectively. K z represents the vertical diffusivity. The last term on the right hand side represents the air-sea oxygen flux, which depend on the partial pressure airsea difference, solubility and velocity transfer of oxygen. The atmospheric oxygen concentration is constant throughout the simulation. Besides, the biogeochemical processes which produce and consume oxygen are parameterized as follows (see also Resplandy et al., 2012): Diatoms (D) and nanophytoplankton (P) produce oxygen during photosynthesis through new and regenerated production. The corresponding production rates are µ nh4 and µ no3 . Oxygen consumption is due to remineralization of Dissolved Organic Carbon (DOC), respiration of mesozooplankton (G Z Z) and microzooplankton (G M M), and nitrification. λ * DOC represents the remineralization rate. (O 2 ) varies between 0 (in oxic conditions) and 1 (anoxia), mimicking the suppression of remineralization in anoxic conditions. R 1 O C (131/122) represents the change in oxygen relative to carbon when ammonium is converted to organic matter during new production and when DOC is being respired. R 2 O C (32/122) represents the change in oxygen relative to carbon during nitrification (Nit). More detail about the model structure and parameterizations can be found in Aumont et al. (2015).

Model Configuration
The model domain spanned from 15 • N to 40 • S and from 100 • W to 70 • W (Figure 1). Due to the extension of the model area and computational requirements, the horizontal resolution of the grid was 1/6 • (∼18 km). It allowed to reproduce mesoscale structures off Peru owing to the proximity of the region to the equator (see section Ventilation by Eddy Fluxes). The bottom topography from ETOPO2 (Smith and Sandwell, 1997) was used, and the vertical grid had 32 sigma levels.
Open boundary conditions (OBC) for physical variables came from an interannual SODA model solution (version 2.1.6; Carton and Giese, 2008) over the period 1958-2008. As an interannual global simulation of the biogeochemical conditions was not available during the period of study , climatological OBC from CARS2009 climatology (Ridgway et al., 2002) were used for nutrients (nitrate, silicate, phosphate) and oxygen and from World Ocean Atlas climatology (WOA2005; Conkright et al., 2002) for DOC, dissolved inorganic carbon and total alkalinity. As a gridded climatology of iron measurements does not exist, iron OBC from a climatology of a NEMO-PISCES global simulation were used (Aumont et al., 2015).
Wind stress, fresh water, sensible and latent heat fluxes were computed using the bulk parameterization of Liu et al. (1979). The surface wind fields were obtained by summing statisticallydownscaled NCEP daily wind anomalies  and the SCOW monthly climatology (Risien and Chelton, 2008). NCEP daily anomalies and COADS monthly climatology (Da Silva et al., 1994) were summed to obtain surface air parameters. Climatologies of short (COADS) and downward long wave (NCEP) heat fluxes were used. The forcing characteristics of the different simulations are listed in Table 1.
The model was run from 1958 to 2008, the period over which the SODA solution that we use (v.2.1.6) is available. Note that over the last 10 years, there was no occurrence of any strong, or extreme, LN or EN events off Peru. Thus, we are confident that the mean state of ENSO phases evaluated in this work are robust.

Model Experiments
The interannual simulation which is mainly described and studied in the following sections is named the control run (CR). It was forced by interannual atmospheric forcing and boundary conditions. The circulation, productivity and nutrient fluxes during ENSO were evaluated in Espinoza-Morriberón et al. (2017).
Two additional simulations were performed and analyzed ( Table 1): -The Wclim simulation, forced by interannual OBC and monthly climatological atmospheric forcing, was used to evaluate the role of the oceanic remote forcing, e.g., equatorial subsurface currents, equatorial Kelvin waves and coastal trapped waves. -The Kclim simulation, forced by monthly climatological OBC and interannual atmospheric forcing, allowed to study the role of the wind forcing.
Due to computational requirements Kclim and Wclim were run from 1979 to 2008, a period which included two extreme EN events (1982-1983/1997-1998) and one strong LN event (1999)(2000). Five day-averaged outputs were stored for each simulation.

Lagrangian Analysis
The ROMS-offline tracking module (Capet et al., 2004) was used to calculate the trajectories of virtual floats simulating water parcels. The floats trajectories were computed using the 5-day averaged velocity fields. Oxygen concentration were registered along the trajectories. In order to describe the source waters (SW) characteristics reaching the nearshore OMZ, 2000 floats were released over the shelf and slope, between the surface and 350 m depth within ∼150 km from the coast on the first day of each month of spring (October, November, December), during each LN and EN events and during neutral periods. The floats were released from a zonal section at 9 • S, but other sections (e.g., 12 • S) were tested. The floats were then tracked backwards in time during 2 years. The floats reaching the so-called offshore equatorial region (defined by a meridional section at 88 • W, 2 • N−10 • S) were used to compute statistics. The characteristics The Control Run (CR) uses interannual forcings over 1958-2008. Wclim simulation used climatological atmospheric forcing and interannual remote oceanic forcing. Kclim simulation used interannual atmospheric forcing and climatological remote oceanic forcing. Both runs were run over 1979-2008. of these waters (DO concentration and depth in the equatorial region, duration of the transit from the equatorial section to the cross-shore coastal section (at 9 • S) were mapped onto the crossshore section where the water parcels were released (see section Changes in the Origin of the Shelf and Slope Waters During ENSO).

Eddy Fluxes
The oxygen vertical eddy flux was computed as where w and O 2 are the 5 day-averaged vertical velocity and oxygen concentration, respectively. . denotes a moving filter with a 60 days' window. The mean profiles of oxygen mean vertical flux and eddy vertical flux were obtained by averaging values horizontally (horizontal averaging is marked by an overbar) in an offshore box (between 100 and 500 km from the coast and between 6 • S and 16 • S). We chose to average profiles in an offshore box in order to discard the effect of nearshore upwelling and downwelling velocities associated with the CTWs and to retain the role of eddies and filaments.
The meridional (v ′ .O ′ 2 ) and zonal (u ′ .O ′ 2 ) eddy oxygen fluxes were computed using the same time filter. Maps of the intensity and direction of the horizontal eddy fluxes 2 ) at 100 and 400 m during neutral, LN and EN events are described in section Eddy Oxygen Fluxes.

ENSO Index and Statistical Tests
ENSO phases are defined when the 3-month-running mean SST anomaly in the "Niño1+2" region (0-10 • S, 90 • W−80 • W) is higher or < ± 0.5 • C for at least 5 consecutive months for EN and LN, respectively, (as in Espinoza-Morriberón et al., 2017). The ERSST.v4 SST product of Huang et al. (2015) was used. Neutral periods correspond to non-EN and non-LN periods. For EN events, if the SST anomalies are > + 1.6 • C during at least 3 months the event is categorized as "extreme, " otherwise it is a "moderate" EN (see Table in Supplementary Material). The correlation between the Niño1+2 index computed from the model SST and observations was 0.9, thus the majority of EN and LN events were represented by the model.
The non-parametric Wilcoxon test was used to statistically compare the means of some characteristics (e.g., OMZ depth) for EN, LN and neutral period, and to test if these means differ with a 95% confidence level. To compute this test, we used the function "pairwise.wilcox.test" from R software.

Oxygen in situ Observations
Approximately ∼20,000 DO vertical profiles collected by the Peruvian Marine Institute (IMARPE) between 1960 and 2008 Ledesma et al., 2011) were used to evaluate the modeled OMZ. DO was determined by the Winkler method (Carrit and Carpenter, 1966). Observations (Nansen and Niskin bottles) were collected during surveys and at fixed stations along the coast. The data was gridded at the same resolution as the model (1/6 • ). Composite alongshore-averaged cross-shore sections for neutral and ENSO phases between 6 • S and 16 • S and 100 km from the coast (hereafter coastal region) were computed. The depth of the upper limit of the OMZ (i.e., where DO is equal to 22 µmol kg −1 ) was computed from each IMARPE profile (linearly interpolated on a vertical grid with a 1 meter resolution). These values were then averaged horizontally in the coastal region to produce an index (ZO2) characterizing oxygenation in the nearshore region.
In addition, CARS2009 DO observations (Ridgway et al., 2002) were used to evaluate the spatial distribution of the simulated oxygen over the entire model domain. The CARS 1/2 • gridded data was interpolated onto the 1/6 • model grid.

Mean State
The mean modeled and observed OMZ are shown in Figure 2. ZO2 presented a marked cross-shore gradient, with shallower depths (<25 m) near the coast than offshore (∼100 m at 100 km from the coast off the central shelf). The nearshore ZO2 was deeper in the model (∼50-70 m between 6 • S and 14 • S), than in the observations (∼30 m in IMARPE and ∼ 50 m in CARS). ZO2 reached 150 m depth at 86 • W between 6 • S and 14 • S (Figures 2a-c), highlighting the offshore deepening and westward extent of the OMZ.
Figures 2d,e display a meridional section of the OMZ at 86 • W, more than 1,500 km east of the model western boundary. The modeled and observed OMZ boundaries were rather similar between 6 • S and 14 • S. Although ZO2 was relatively similar in the model and CARS observations, the maximum OMZ thickness (found near 8 • S−10 • S at 86 • W) was larger in the model than CARS due to the OMZ deeper (by ∼100 m) lower limit in the model. In the core of the OMZ between 8 • S−11 • S and 200-500 m depth, modeled DO concentrations were lower than 5 µmol kg −1 , lower than CARS (∼10 µmol kg −1 ). Between 2 • S and 2 • N the modeled OMZ vanished, while a thin (∼100 m) OMZ remained between 300 and 450 m depth in CARS.

Interannual Variability
Due to the relative scarcity of in situ data, it was impossible to compute cross-shore sections of DO at distinct latitudes, therefore we present alongshore-averaged cross-shore sections between 6 • S and 16 • S (see section Oxygen in situ Observations) representative of averaged DO conditions near the Peru shelf and slope (e.g., Espinoza-Morriberón et al., 2017). The model and data sections are shown for neutral, EN and LN phases (Figure 3). Near-surface slanted isolines of oxygen were evidenced in the model and IMARPE sections during neutral, EN and LN events, indicating the occurrence of coastal upwelling during each phase (Figures 3A-C,F-H). During neutral periods, between the coast and 200 km offshore, the observed ZO2 was located at ∼180 m depth and the modeled ZO2 at ∼150 m depth (Figures 3A,F). In contrast, during EN (LN) a deeper (slightly shallower) ZO2 was observed in the model (Figures 3B,C). These features were well reproduced by the model (Figures 3G-H). The main differences between model and observations were found in the OMZ core and near its lower limit. In the model, DO values <5 µmol kg −1 were observed in the OMZ core (between 300 and 450 m depth) during neutral and LN periods. During EN, the thickness of lowest DO (<5 µmol kg −1 ) layer was reduced. In contrast, the observations presented a slightly more oxygenated core than in the model during neutral (∼8 µmol kg −1 ), EN (∼15 µmol kg −1 ), and LN (∼12 µmol kg −1 ) periods. Moreover, the modeled OMZ lower limit was deeper than that of the observed OMZ, in particular during LN (Figures 3C-H). During EN, both observations and model presented a thinner OMZ than during neutral periods, owing to a deepening of its upper limit and a shoaling of its lower limit (Figures 3B,G). In contrast, the model was not able to simulate the observed OMZ thinning during LN. Indeed, the modeled OMZ was thicker during LN with respect to neutral conditions due to little change of its lower limit, whereas the shoaling of the lower limit was strong in the observations.
In terms of magnitude, the DO changes in the surface layer (0-150 m) during ENSO phases were qualitatively reproduced by the model. The highest DO anomalies were observed between 10 and 150 m (Figures 3D-I) both in the model (+40 µmol kg −1 ) and in observations (+25 µmol kg −1 ) during EN. During LN, negative anomalies of similar amplitude in the model and observations (∼-10 µmol kg −1 ; Figures 3E,J) were found in the same layer between the coast and 150 km offshore.
Interannual variations of the nearshore modeled and observed ZO2 were highly correlated over 50 years (r = 0.72, p < 0.05, Figure 4). During each EN event, ZO2 increased strongly. The maximum ZO2 (∼300 m) was attained during the peak of the extreme 1982-83 EN for model and observations. In the model and observation, the ZO2 during Neutral and LN periods did not present statistical differences (p > 0.05); however, both presented significant differences respect to EN event (p < 0.05). During neutral periods, the observed and modeled ZO2 were ∼72 m depth on average, while during LN, the oxycline was found at ∼60 m for both. In conclusion, in spite of the discrepancies in FIGURE 4 | ZO2 time series (in meters). Data was averaged each semester for IMARPE (blue line) and model output (black line) in a coastal box (6 • S−16 • S and 100 km from the coast). Error bars represent the standard deviation from IMARPE data. Black bars (at the bottom of the figure) represent the coverage ratio (number of sampled points/ total number of points). The modeled ZO2 was computed using the same sampling as IMARPE observations (gaps mean there is no data). Gray and blue shading represent EN and LN period, respectively. DO concentration within the OMZ (Figures 3-5), the model represented realistically the oxygen changes during ENSO phases.

Coastal Trapped Waves
The impacts of the downwelling and upwelling CTWs during an extreme EN (1997EN ( -1998) and a strong LN (1999LN ( -2000 are highlighted in Figure 5. During the 1997-1998 EN the poleward propagation of two downwelling CTWs was evidenced on the Sea Surface Height (SSH) signal, which showed positive anomalies. The thermocline depth (defined by the position of the 15 • C isotherm, Z15) deepened during the first (∼+80 m) and second (∼+150 m) peak (Figures 5a,b; Espinoza-Morriberón et al., 2017). The slightly tilted isolines indicate poleward propagation of the signals. ZO2 deepened strongly during EN (Figure 5c). The deepest anomalies (+200 m) were observed during the second peak of the event.
During the strong 1999-2000 LN, two strong upwelling IEKW crossed the Central Pacific provoking the shoaling of the 20 • C isotherm (figure not shown). They triggered two CTWs which impacted the alongshore SSH, thermocline depth (ZT15) and ZO2 during the summer-autumn of 1999 and 2000 (Figures 5d-f), with magnitudes weaker than those of the EN downwelling CTWs. The SSH anomalies presented negative values (Figure 5d) and a shallower ZT15 was observed, associated with negative anomalies of ∼50 m north of 6 • S, and ∼30-40 m between 10 • S and 15 • S (Figure 5e). A shoaling ZO2 was also observed during the passage of both upwelling CTWs. Negative anomalies of ∼45 m were found between 8 • S and 14 • S during the passage of the CTW (Figure 5f). Note that the upwelling CTW during LN occurred during the same season (summer) in 1999 and 2000.  1997-1998 (top) and LN 1999 (bottom) events. Model values were averaged within 100 km to the coast. All variables were filtered in time (60 days moving average) and space (100 km alongshore). Empty data in (b,c,e,f) means that ZT15 and ZO2 were not detected.

Changes in the Origin of the Shelf and Slope Waters During ENSO
In this section, we investigate the modification of the pathways of equatorial water masses (hereafter source waters, SW) reaching the PCUS during ENSO phases, which could have consequences on the OMZ ventilation. We characterize the properties of the offshore (88 • W) waters that reached the Peruvian shelf and slope (at 9 • S) during EN and LN. Note that modifying the latitude of the cross-shore section (e.g., 12 • S) did not modify substantially the results presented below.
Expectedly, the DO nearshore section displayed a deeper and shallower ZO2 during EN and LN, respectively (Figures 6a,f,k,p). The equatorial SW were impacted during ENSO phases. During neutral periods the SW DO concentration was ∼60 µmol kg −1 in average (Figure 6g), while SW were slightly less oxygenated (∼40 µmol kg −1 ) during LN (Figure 6b). In contrast, SW were more oxygenated (∼80-120 µmol kg −1 ) during EN depending on the intensity of the events (Figures 6l,q).
The changes in offshore oxygen content were related to changes in the location of the SW, associated with the modified SSCCs fluxes during ENSO. Within 100 km from the coast, during neutral periods, the SW originated from 3.5 • S to 4 • S (Figure 6h) and 110-140 m depth (Figure 6j). SW were then transported toward the coast for ∼12 months (Figure 6i). During LN, SW were located south of 4.5 • S (Figure 6c) and at greater depths (∼160-180 m; Figure 6e), thus were less oxygenated (Figure 6b). Moreover, they were older since it took them more time (∼15 months; Figure 6d) to reach the coastal section. In contrast, during EN, the SW were initially closer to the equator (3 • S−1 • S; Figures 6m,r), shallower (∼100-120 m depth; Figures 6o,t) and thus more oxygenated (Figures 6l,q) than during LN and neutral periods. The SW transit time to the coastal region was much shorter (5-10 months; Figures 6n,s).

Eddy Oxygen Fluxes
As eddy activity is enhanced during ENSO (Chaigneau et al., 2008;Espinoza-Morriberón et al., 2017), we investigate the impact of EN and LN events on the oxygen vertical and horizontal fluxes.
The mean oxygen vertical flux ( w . O 2 ) presented positive values with a peak at ∼40 m depth (Figure 7a). As the input of oxygen due to the mean flux is equal to-∂ z ( w . O 2 ), it adds oxygen above ∼40 m (mainly due to the decrease of mean vertical velocity near the surface, Figure 7d) and takes out oxygen in the subsurface layer between 40 and 150 m depth (due to the upwelling of less-oxygenated deep water from the OMZ core; e.g., Figures 3F-H). At the maximum flux depth (∼40 m) the mean fluxes (∼0.3 µmol m −2 s −1 ) did not present significant differences between neutral and LN periods (Wilcoxon test, p > 0.05), while a decrease was observed (∼0.24 µmol m −2 s −1 , ∼20% reduction) during EN (Wilcoxon test, p < 0.05). This flux reduction was partly driven by the reduction of the vertical velocity (Figure 7d), which counteracted the impact of the subsurface oxygen increase (Figure 7c). Note that the modeled vertical velocity decrease (∼40% at 50 m depth, Figure 7d) cannot be fully attributed to a decreased Ekman pumping during EN (e.g., Halpern, 2002;Chamorro et al., 2018, for the 1997-1998, as it was only ∼5% in our simulation (Figure not shown).
The vertical eddy flux also presented positive values with a maximum at ∼60 m depth (Figure 7b), which corresponds to the depth of maximum oxygen vertical gradient (∂ z O 2 ) (Figure 7c). In the stratified interior, mesoscale eddy fluxes tend to flatten the tilted upper thermocline (Colas et al., 2013) and hence the oxycline. Consequently, the vertical eddy flux brings oxygen into the layer above ∼60 m, while below it acts as an oxygen sink and reinforces the OMZ. As for the mean flux, Neutral period presented a value of ∼0.2 10 −1 µmol m −2 s −1 at 60 m, and was not significantly different to LN (Wilcoxon test, p > 0.05), while lower values (∼0.06 10 −1 µmol m −2 s −1 , ∼70% decrease) and significant differences (Wilcoxon test, p < 0.05) were observed during EN between 50 and 70 m depth. This was associated to a reduction of the oxygen mean vertical gradient during EN (Figure 7c). As a consequence, the vertical gradient of the vertical eddy flux [-∂ z (w ′ O 2 ′ )] was also much weaker during EN than during LN and neutral periods. This highlights the much weaker role of the eddies as sink/source of oxygen below/above the oxycline during EN periods. At depths >100 m, the vertical eddy flux reached values close to zero during neutral and LN periods, whereas slightly positive values (∼0.02 10 −1 µmol m −2 s −1 ) were found during EN. Figure 8 shows the horizontal oxygen eddy fluxes at 100 and 400 m during ENSO phases and neutral periods. Intense eddy fluxes are found in the vicinity of the OMZ boundaries, mainly along the southern limit of the OMZ at 100 m depth and in the equatorial region (4 • N−4 • S) at 400 m depth. Hence this is logically the case at the lateral boundaries of the OMZ, especially at its southern limit where the horizontal gradient of DO is intense at 100 m depth. Eddy advection tends to inject oxygen into the OMZ, and it is more intense at shallower depths. During EN, at 100 m depth, the horizontal oxygen eddy fluxes increased (∼+0.5 µmol m −2 s −1 ; Figure 8C) along the OMZ southern boundary and alongshore (∼5 • S−18 • S) from the coast to 250 km offshore ( Figure 8D). This indicates that eddy fluxes contributed to oxygenate the coastal region. At 400 m, fluxes were much weaker than at 100 m but also enhanced during EN periods (+∼0.25 10 −1 µmol m −2 s −1 ; Figure 8G), both at the southern and mainly at the northern limits of the OMZ and along the northern shelf ( Figure 8H). Eddy fluxes were quite similar during LN and neutral periods (Figures 8A,E). Overall, the OMZ ventilation at its lateral boundaries by mesoscale eddy advection increased during EN, whereas ventilation due to vertical eddy fluxes reduced.

Biologeochemical Drivers of O 2 Changes During ENSO
The OMZ changes due to biogeochemical processes is investigated. The vertical profiles of biogeochemical source and sinks of oxygen (see Equation 2 in section The Oxygen Cycle in the PISCES Model) were averaged horizontally within 200 km to the coast (Figure 9). Most of the modeled biogeochemical processes related to the oxygen cycle take place in the near-surface layers and impact the OMZ indirectly through subduction of oxygen-enriched/depleted waters into the oxycline (e.g., Thomsen et al., 2016). The annual mean profiles of DO production and consumption had an exponential decrease with depth, with the exception of nitrification. The new and regenerated production (through photosynthesis) provided oxygen to the upper ocean above the ZO2 (∼100 m depth). Remineralization was the main process consuming oxygen in surface and subsurface water, above 300 m. Microzooplankton respiration consumed oxygen above the upper limit of the OMZ, while mesozooplankton respiration was also found at depth. Oxygen consumption by nitrification peaked at 30 m in oxygenrepleted waters, as it is inhibited by light (Yoshioka and Saijo, 1984) near the surface (0-10 m depth). Nitrification decreased in the subsurface layer due to oxygen limitation (see Equation  2 in section The Oxygen Cycle in the PISCES Model). Within the OMZ (below 120 m depth), only oxygen consumption was observed (due to remineralization, mesozooplankton respiration, and nitrification).
The impacts of EN and LN on the oxygen sources and sinks were diagnosed in the surface layer (above 100 m depth) and within ∼60 km from the coast (Figure 10). Note that the sinks (i.e., negative terms in Equation 1, section The Oxygen Cycle in the PISCES Model) are shown here with a positive sign, meaning that positive (resp. negative) sink anomalies (e.g., Figures 10C-F,

I-L) indicate an increase in oxygen consumption.
During LN, a positive anomaly (>2 10 −3 µmol m −3 s −1 ) produced by the increase in new and regenerated production was found in a very thin surface layer (above ∼10 m) and ∼20 km offshore. In contrast, just below (10-80 m depth) and nearshore (within 20 km from the coast), new and regenerated production were weaker by ∼-0.45 and ∼-0.8 10 −3 µmol m −3 s −1 , respectively, (Figures 10A,B).
During EN, the production of oxygen decreased (<-6.5 10 −3 µmol m −3 s −1 ), while new (∼+0.46 10 −3 µmol m −3 s −1 ), and regenerated (∼+0.73 10 −3 µmol m −3 s −1 ) productions were enhanced in the subsurface layer between 10 and 80 m depth nearshore and within 60 km from the coast (Figures 10G,H). During ENSO phases, oxygen consumption was also strongly modified (Figures 10C-F,I-L) with anomalies of DOC remineralization and zooplankton respiration presenting marked cross-shore patterns. Nearshore, negative and positive anomalies predominated during LN and EN, respectively, while ∼30 km offshore and above 15 m, opposite sign patterns were found. Nitrification anomalies presented a pronounced vertical gradient. Below ∼30 m depth, the nitrification decreased and increased during LN and EN, respectively, due to the presence/absence of oxygen.
In contrast, remineralization decreased in a thin surface layer (0 -∼10 m) but strongly increased below, down to 100 m depth during EN. Positive anomalies of +1.6 10 −3 µmol m −3 s −1 were found nearshore (within ∼80 km from the coast and between 10 and 80 m, Figure 10I). Expectedly, meso and microzooplankton respiration reduced in a thin surface layer and increased at depth above the shelf, producing positive anomalies of +0.55 and +0.31 10 −3 µmol m −3 s −1 , respectively, between 10 and 80 m and from 0 to 80 km to the coast (Figures 10J,K). Last, nitrification was less intense between the surface and 40 m depth and increased near the bottom over the shelf (+0.21 10 −3 µmol m −3 s −1 between 40-100 m and 0-80 km from the coast; Figure 10L). Nitrification, which was switched off in this depth range during neutral and LN conditions due to the lack of oxygen, took place during EN as shelf waters became oxygenated. Table 1) were performed to investigate the impact of remote (Kclim) and local forcing (Wclim) on the OMZ during ENSO. Figure 11 shows the oxygen changes (in %), in Kclim and Wclim with respect to CR, in an alongshore-averaged (between 6 • S and 16 • S) cross-shore section (see also Figure 3).

Sensitivity experiments (see section Model Experiments and
During LN, the oxygen content increased slightly (∼10%, +10 µmol kg −1 ) in Kclim with respect to CR in the surface layer (between 0 and 50 m depth and within 200 km from the coast, Figure 11A). In this experiment, the ZO2 shoaling forced by the LN upwelling CTW was absent. In contrast, a moderate DO decrease (10-20%, −2 µmol kg −1 ) was found below the upper limit of the OMZ. When the interannual wind variability is suppressed (Wclim), a moderate oxygen decrease (∼3%, ∼-8 µmol kg −1 ) was found above the ZO2 (Figure 11B). Below the ZO2, the DO concentration increased offshore by 10% (∼+2 µmol kg −1 ).
During EN, a strong DO reduction was found over the entire water column in Kclim. It was strongest between 90 and 150 m depth, where DO was ∼60 % (∼-55 µmol kg −1 ) lower than in CR ( Figure 11C). In contrast, more oxygen was found (12%, +10 µmol kg −1 ) in Wclim than in CR between 40-130 m depth and 50-250 km offshore ( Figure 11D).

Bias of the Modeled OMZ
The interannual variability of the main physical parameters, surface productivity and nutrients simulated by our model were validated in Espinoza-Morriberón et al. (2017). The model reproduced well the changes of SSH, temperature, nutrients and chlorophyll-a (Chl) in surface and subsurface waters during both EN and LN phases. Here, we show that the model reproduced deoxygenation and oxygenation during LN and EN, respectively, in agreement with previous studies based on in situ data analysis (Graco et al., 2007(Graco et al., , 2017Gutiérrez et al., 2008). In addition, in annual mean, the shoaling of ZO2 toward the coast (∼50 m depth) was also observed between 14 • S and 16 • S by Graco et al. (2007), likely associated with coastal upwelling. Nevertheless, some discrepancies were observed. First, between 2 • N and 2 • S, observational data evidenced a thin and coastal OMZ (e.g., Karstensen et al., 2008;Czeschel et al., 2015;Llanillo et al., 2018), that is not reproduced by our simulation. This may be due to the intense EUC simulated by the SODA model (used as the physical OBC forcing in our simulation) which ventilates the equatorial region (Montes et al., 2014). Off Peru, overly low DO (<5 µmol kg −1 ) was found within the modeled OMZ core (Figures 2, 3) in line with Karstensen et al. (2008). However, this bias with respect to the CARS climatology and IMARPE observations could be partly due to the numerical artifacts during the construction of the 3D fields and/or to the inclusion of records from old instruments which were not sufficiently accurate at low DO values (Bianchi et al., 2012). Underestimation of OMZ volumes slightly more oxygenated than the core was also reported i.e., Fuenzalida et al. (2009). Blending vertical highresolution CTDO profiles with quality-controlled bottle casts from WOD01, Fuenzalida et al. (2009) found that the extent of hypoxic areas (<60 µmol kg −1 ) in the TSEP was underestimated in global databases such as WOA01. Note that the more recent WOA 2013 DO climatology (Garcia et al., 2013) did not display lower values than CARS (Figure not shown), so for that reason we chose the latter.
Crude parameterization of the oxygen cycle in the PISCES model, due to the fixed-in-time values of biogeochemical parameters (e.g., remineralization rate of DOC and nitrification rate) throughout the interannual simulation may also impact the modeled OMZ. For example, the remineralization rate of DOC only takes into account a proxy of bacterial concentration and DOC concentration in PISCES (Aumont et al., 2015), even though temperature changes (which are particularly strong during EN) impact the bacterial metabolic rates involved in the oxygen cycle (López-Urrutia et al., 2006), an effect seldom included in biogeochemical models (Segschneider and Bendtsen, 2013). Note also that previous regional models of the OMZ using the nitrogen-based BioEBUS model (Gutknecht et al., 2013) also produced a more intense OMZ than the observations (Montes et al., 2014;Vergara et al., 2016;Mogollón and Calil, 2017).
Last, we cannot exclude the possibility that our boundary forcing (SODA) may have partly induced the very low OMZ concentrations simulated by our model. Indeed, Montes et al. (2014) showed that using different sets of physical OBCs, which differed in the SSCCs vertical structure intensity, could modify the shape and intensity of the OMZ.

Impact of Eastward Fluxes From the Equatorial Region
The equatorial circulation was impacted during ENSO, which modified the eastward oxygen fluxes in the equatorial region offshore of the PCUS. To illustrate it, the intensity of the zonal currents (EUC, pSSCC, and sSSCC) and associated oxygen fluxes at 88 • W between 50 and 200 m were computed ( Table 2). The EUC was weaker during EN and LN than during neutral conditions. Stramma et al. (2016) also found a strong weakening of the EUC velocity during the onset (October) of the 2015-2016 EN. In contrast, the SSCCs intensity increased during EN with respect to LN and neutral conditions. These results are in line with Montes et al. (2011), who simulated the circulation during the moderate 2002-2003EN and 1999-2000 Even though stronger SSCCs were observed during EN, the nearshore water column (e.g., at 9 • S, Figure 6) was mainly fed by waters initially transported by the EUC. The water masses transit time from 88 • W to the nearshore region ranged from ∼15 months during LN to ∼6 months during EN (Espinoza-Morriberón et al., 2017). During LN, the longer transit may partly explain the nearshore oxygen deficit due to microbial respiration in the transported water mass. Furthermore, note that the fast pathway observed in the model during EN, from the equatorial regions to the coastal range, was also found by Montes et al. (2011) for the 2002-2003. This shows that not all particles released in the equatorial zone reached the nearshore regions during the same ENSO phase. In the equatorial region, in our simulation, the oxygen fluxes variability is driven by the equatorial currents variability, and not by oxygen variability as oxygen western boundary conditions are climatological. Thus, the oxygen fluxes from the SSCCs and the EUC were less intense during LN than during EN (an exception is the EUC flux during moderate EN; see Table 2). During EN, enhanced SSCCs oxygen contributed to reduce the OMZ volume and displace its offshore boundaries eastward (Montes et al., 2014).

Ventilation by Eddy Fluxes
The horizontal eddy fluxes play an important role in structuring the OMZ lateral extent (Vergara et al., 2016) where strong oxygen horizontal gradients are found. It is well-known that eddy tracer fluxes are downgradient (and pronounced) when there is a strong sink in the tracer distribution (e.g., Wilson and Williams, 2006). In a recent modeling study, Vergara et al. (2016), showed that horizontal eddy fluxes through the OMZ coastal boundary attains its maximum in winter in association with a maximum seasonal Eddy Kinetic Energy (EKE). Similarly, the EKE increase during EN induces a large interannual horizontal eddy flux which ventilates the nearshore OMZ. The magnitude of this enhanced interannual ventilation by horizontal eddy fluxes could be underestimated due to the relatively low horizontal resolution (1/6 • ) of our model (e.g., Capet et al., 2008;Colas et al., 2013). However, a resolution of 1/6 • is considered eddy-resolving off Peru: due to its proximity to the equator, the typical length scale of mesoscale structures would be covered by ∼10 grid points (Belmadani et al., 2012). Indeed, in spite of this caveat, our model represented correctly the spatial and interannual variability of the EKE (Espinoza-Morriberón et al., 2017). Nevertheless, a more complete modeling study with an increased spatial resolution would be needed to better understand and quantify the role of submesoscale variability on the OMZ (Thomsen et al., 2016).

Biological Processes
Remineralization is the major oxygen sink in the water column (Paulmier and Ruiz-Pino, 2009) and is mainly produced by bacterial respiration off Peru (Kalvelage et al., 2015). This is reproduced in our modeling results (Figure 9). Furthermore, zooplankton respiration has been little studied in the PCUS in spite of its high oxygen demand in the water column. Global estimates of mesozooplankton respiration (in fraction of body carbon respired daily) indicated that it represents ∼17-32% of the photosynthetic carbon produced in the open ocean (underestimated in previous works e.g., del Giorgio and Duarte, 2002;Hernández-León and Ikeda, 2005). In our configuration, the zooplankton respiration represented ∼33% of the oxygen produced by photosynthesis (between 6 • S−16 • S and 0-200 m depth), close to global estimates. High biomass of zooplankton has been found in the PCUS, showing that zooplankton can adapt to low-oxygen conditions at great depths (Ayón et al., 2008;Ballón et al., 2011). This implies a high oxygen consumption at greater depths due to the mesozooplankton respiration, which is more tolerant to suboxic condition than microzooplankton (this is explicitly included in PISCES). However, due to a lack of data, the zooplankton biomass in PISCES (see Figure 14 in Espinoza-Morriberón et al., 2017) could not be evaluated, which is needed to estimate reliable respiration rates.
Our model setting has also some limitations, as the nitrogen cycle is crudely parameterized in PISCES (Aumont et al., 2015). High concentrations of nitrous oxide (N 2 O), a potent greenhouse gas produced during nitrification (Codispoti and Christensen, 1985), have been observed off the Peruvian coasts (Kock et al., 2016), especially during pulses of oxygenation. Although not simulated by our model, production of nitrous oxide is likely to occur during EN due to enhanced nitrification in subsurface layers ( Figure 10L). Using the nitrogen-cycle-based ROMS-BioEBUS coupled model, Mogollón and Calil (2017) found a decrease (80 %) of nitrification between the peaks of 1997-1998EN and 1999-2000 LN and production of nitrous oxide during EN. In line with this study, our model simulated a nearshore nitrification decrease of ∼60% between EN and LN (within 20 km from the coast and between 20 and 80 m depth). Furthermore, under suboxic conditions (<5-10 µmol kg −1 ), both anammox (Lam et al., 2009), and denitrification could produce N loss (Dalsgaard et al., 2012). The percentage of N loss contribution by both processes depends on ENSO phases. During 1997-1998EN, anammox, and denitrification contribute to N loss of 40 and 60%, respectively, while during 1999-2000 LN the contribution of anammox increase of 70% (Mogollón and Calil, 2017). In PISCES, anammox is not included, which could explain the overestimation of nitrate in our simulation (see Figure 3 in Espinoza-Morriberón et al., 2017).
In PISCES, denitrification occurred in regions of DO <6 µmol kg −1 (hereafter ODZ) and was simply parameterized using DOC and a remineralization rate (Aumont et al., 2015): Where R NO3 represents the N/C stoichiometric ratio. Our model simulated an increase of the denitrification in nearshore and shallower waters during LN (e.g., during the 1999-2000 LN, Figure 12a), due to the shoaling of the ODZ, despite of the unchanged ODZ volume. This is likely due to a shoaling of the OMZ (Figure 3J), which allowed denitrification to occur at shallower depths, and enhanced remineralization ( Figure 10C). In contrast, during EN, the ventilated water column inhibited denitrification above ∼200 m, and a deeper offshore ODZ was observed (e.g., during the 1997-1998 EN, Figure 12b). Overall, denitrification within the ODZ (defined between 6 • S−16 • S, 0-500 km, and 0-500 m depth) increased by ∼25% during LN and decreased by ∼60% during EN (moderate or extreme), with respect to neutral periods (Figure 12c). Using a global biogeochemical model, Yang et al. (2017) found a denitrification increase of ∼70% between the peak of 1997-1998EN and 1999-2000 which is close to our estimate (∼85%). On the other hand, DO changes during LN and EN were influenced not only by physical processes but also by changes in DO production related to primary productivity changes above the upper limit of the OMZ, between 0 and ∼20 m. The nearshore DO changes induced an increased (decrease) oxygen consumption by remineralization and respiration during EN  (LN). During LN, more (less) oxygen production was found at surface (sub-), due to the thinner mixed layer during LN, in which primary production was confined (Espinoza-Morriberón et al., 2017). In contrast, during EN, the nutricline depth increased and photosynthesis decreased in the surface layer (Espinoza-Morriberón et al., 2017), producing less oxygen.
Oxygen consumption decreased and increased during LN and EN, respectively (mainly at subsurface), due to DO availability. The relationship between the oxygen consumption and DO availability was also evidenced offshore (86 • W) using in situ data during the 2009 LN (Llanillo et al., 2013). The authors found an increase of the oxygen consumption by respiration and The EUC, pSSCC, and sSSCC DO and velocity were averaged between 2 • N−2 • S, 2 • S−6 • S and 6 • S−10 • S, respectively, and from 50 to 200 m depth at 88 • W. remineralization within the OMZ due to the presence of the oxygen-rich Antarctic Intermediate Waters (AAIW).

Remote and Local Drivers of DO Interannual Variability
The equatorial remote forcing was the main driver of changes in the ZO2 and OMZ core during both EN and LN. The passage of downwelling and upwelling CTWs trigger the oxygenation (in the surface and subsurface layers) and deoxygenation (only in the surface layer) during EN and LN, respectively. At subsurface during LN, remotely-forced deep coastal waves (e.g., Pietri et al., 2014) may enhance the equatorward transport of the oxygenated AAIW, which ventilates the deep OMZ between ∼150-350 m and 6 • S−14 • S (e.g., LN 2009; Llanillo et al., 2013). The strong relationship between the variability of ZT15 and ZO2 is well known locally at specific coastal sites off Peru (e.g., 12 • S, Gutiérrez et al., 2008;Graco et al., 2017), and in other regions where CTWs play an important role, e.g., Vallivattathillam et al. (2017) for the western coast of India. These authors mentioned that the effects of biogeochemical processes (e.g., remineralization and respiration) and of vertical mixing were smaller than those due to advection (e.g., CTW). This strong relationship is highlighted in Figure 13A, which shows the spatial pattern of correlation between ZT15 and ZO2 monthly anomalies over the entire PCUS. The correlation was highest (∼0.8) nearshore (between the coast and 50 km offshore) and dropped offshore and south of 18 • S, where the OMZ vanishes.
The high correlation pattern extended westward until 350 km offshore along 10 • S. North of 10 • S, the correlation decreased due to a ZO2 deepening larger than that of ZT15. Figures 13B-C show cross-shore sections of the correlation between ZO2 and ZT15 and the depth difference at 7 • S and further south, at 10 • S. At 7 • S (Figure 13B), the OMZ was more impacted by ENSO (Mogollón and Calil, 2017) and by the eastward SSCCs transporting oxygen-rich waters (Montes et al., 2014), leading to the increasing distance between ZO2 and ZT15 and correlation decrease west of 81 • W. The ZO2 was shallower at 10 • S than at lower latitudes (Figures 2b,c) as it was less impacted by the oxygen-rich waters transported eastward by the EUC and SSCCs. Note that the correlation was very high along the coasts of Colombia and Central America. This indicates that the nearshore interannual variability of ZO2 in the northern hemisphere is also controlled by the shoaling and deepening of isotherms during the poleward propagation of CTWs.
On the other hand, local wind fluctuations played a minor role in oxygenation/deoxygenation processes at interannual time scales. During EN, when the interannual wind variability is suppressed (in Wclim), a relatively weak subsurface oxygen increase is observed. This can be explained as follows: as the alongshore wind increases during EN (see Figure 9 in Enfield, 1981;Espinoza-Morriberón et al., 2017;Chamorro et al., 2018), the coastal upwelling intensifies. This partly compensates the ZO2 deepening, bringing deep oxygen-poor waters in this depth range. As the EN enhanced wind-driven upwelling was absent in Wclim, the ZO2 was deeper and DO higher than in the CR. During LN, an oxygen decrease was observed above 50 m, in the absence of interannual wind variability (Wclim). Primary production was reduced with respect to CR (between ∼4 and 40 m), thus less oxygen was produced by photosynthesis above 30 m within 200 km to the coast (figure not shown). However, the DO increase offshore and below the ZO2 is possibly related to a reduction (∼20%) of the (positive) Ekman pumping during LN in Wclim, which results in less upward transport of oxygen-poor waters than in CR.

Impact of ENSO Diversity
It is important to mention that our analysis did not take into account the diversity of EN events. Indeed, a Central-Pacific (CP) EN (Modoki; e.g., Takahashi et al., 2011) is likely to impact the OMZ differently than a canonical Eastern Pacific (EP) EN. CP EN events are frequently associated to upwelling CTWs during fall and winter, in contrast with the more frequent downwelling CTWs observed during EP EN . Thus, during EN Modoki poor oxygen waters could be more present than during a canonical EN, due to a shoaling of the ZO2 triggered by the passage of the upwelling CTWs (e.g., Figure 5f). Besides, in early 2017 a "coastal El Niño" occurred (Takahashi and Martínez, 2017;Garreaud, 2018), which also had a different impact on the DO. In fact, during canonical EN (e.g., 2015-2016 EN) a deepening of the ZO2 occurred, while during the summer of 2017 (Coastal EN) the impact was greater in the surface layer but the ZO2 stayed shallower (ENFEN, 2017).

CONCLUSIONS AND PERSPECTIVES
In the present work, the physical and biological processes that drive the interannual variability of the TSEP OMZ were studied using the physical-biogeochemical ROMS-PISCES coupled model. First, an evaluation of our interannual simulation over the period 1958-2008 showed that, despite a reasonable bias in the representation of the modeled OMZ, the phases of oxygenation and deoxygenation during EN and LN were qualitatively well reproduced. During EN, the nearshore ZO2 deepened due to the passing of remotely-forced intense "downwelling" CTWs. During LN, the ZO2 shoaled due to the less intense "upwelling" CTWs, enabling oxygen-depleted waters to accumulate on the shelf. The remote equatorial forcing (IEKWs) was responsible for the tight correlation between the thermocline and ZO2 close to the coasts.
Furthermore, the characteristics of the equatorial SW which were transported toward the coast and fueled the nearshore water column were investigated. SW were mainly transported eastward by the SSCCs (located south of the equator) and EUC during LN and EN, respectively, before reaching the Peru coast. During EN the SW were more oxygenated as they originated from the north of the OMZ and transited faster to the coastal region than during LN.
The vertical and horizontal oxygen eddy fluxes were also evaluated. On average, the offshore (between 100 and 500 km from the coast and 6 • S−14 • S) vertical eddy flux brings oxygen to the surface layer (above ∼60 m), while it removed oxygen below. This flux was not modified during LN, but during EN a very strong decrease was observed, due to the reduced oxygen vertical gradient resulting from the deepening of the ZO2. Horizontal eddy fluxes tended to inject oxygen into the OMZ from its lateral boundaries, especially at its southern limit where the horizontal gradient of O 2 is intense. This was not modified during LN with respect to neutral period. However, during EN, the horizontal oxygen eddy fluxes strongly increased due to the increase of eddy activity, intensifying eddy-driven horizontal ventilation.
Regarding the biogeochemical processes, remineralization was the main process consuming oxygen in the water column, followed by zooplankton respiration and nitrification. These oxygen sinks were strongly altered during ENSO phases, depending on oxygen availability in the water column.
Sensitivity experiments performed to evaluate the influence of the equatorial remote (i.e., CTWs) and local forcing (e.g., wind), demonstrated that the equatorial forcing had a greater impact than local winds on the DO changes, in particular during EN.
Future studies will address the mechanisms which drive the DO changes off Peru over the last decades. While a deoxygenation trend has been observed in the open ocean Breitburg et al., 2018), a nearshore shoaling of the TSEP OMZ has also been described in the late 2000s . In contrast, a relative oxygenatedwater column was observed in the last decade on the Peru central shelf, likely due to more frequent downwelling CTW (Graco et al., 2017). Future studies will also address the impact of climate change in the PCUS. As downscaled projections of climate scenarios in the PCUS simulated a strong nearsurface stratification increase Oerder et al., 2015), this could hamper ventilation and intensify the OMZ. However, further deep and poleward proxy data and instrumental observations suggest ventilation of the OMZ core over the last decades (Cardich et al., in review).

AUTHOR CONTRIBUTIONS
DE-M proposed the original idea of the study, produced most of the material and wrote the manuscript. VE, FC, DG, JT, MG, and CQ-C discussed the results and corrected the manuscript. JL participated to the in situ data collection. Numerical simulations were performed on the ADA computer at IDRIS (project i2015011140). The model results are available upon request to despinoza@imarpe.gob.pe and vincent.echevin@locean-ipsl.upmc.fr. This work is a contribution from the Cooperation agreement between the IMARPE and the Institut de Recherche pour le Développement (IRD), through the LMI DISCOH.