Seasonal Variability of the Southern Tip of the Oxygen Minimum Zone in the Eastern South Paci ﬁ c (30° ‐ 38°S): A Modeling Study

We investigate the seasonal variability of the southern tip (30° – 38°S) of the eastern South Paci ﬁ c oxygen minimum zone (OMZ) based on a high horizontal resolution (1/12°) regional coupled physical ‐ biogeochemical model simulation. The simulation is validated by available in situ observations and the OMZ seasonal variability is documented. The model OMZ, bounded by the contour of 45 μ M, occupies a large volume (4.5x10 4 km 3 ) during the beginning of austral winter and a minimum (3.5x10 4 km 3 ) at the end of spring, just 1 and 2 months after the southward transport of the Peru ‐ Chile Undercurrent (PCUC) is maximum and minimum, respectively. We showed that the PCUC signi ﬁ cantly impacts the alongshore advection of dissolved oxygen (DO) modulating the OMZ seasonal variability. However, zonal transport of DO by meridionally alternating zonal jets and mesoscale eddy ﬂ uxes play also a major role in the seasonal and spatial variability of the OMZ. Consistently, a DO budget analysis reveals a signi ﬁ cant contribution of advection terms to the rate of change of DO and the prominence of mesoscale variability within the seasonal cycle of these terms. Biogeochemical processes and horizontal and vertical mixing, associated with subgrid scale processes, play only a secondary role in the OMZ seasonal cycle. Overall, our study illustrates the interplay of mean and (mesoscale) eddy ‐ induced transports of DO in shaping the OMZ and its seasonal cycle off Central Chile.


Introduction
Despite the large primary productivity that takes place along the coast off central Chile during the upwelling season, surface waters are commonly undersaturated with oxygen (e.g. Letelier et al., 2009). Near the coast subsurface waters with very low dissolved oxygen (DO) are transported to the surface by the upwelling circulation cells. In fact, the presence of an intense and relatively shallow oxygen minimum zone (OMZ) is one of the most striking oceanographic features along most of the eastern boundary of the Eastern South Pacific (ESP). A relative high rate of microbial decomposition of organic matter and weak subsurface circulation combine to generate these oxygen-poor environments (e.g. Helly & Levin, 2004;Karstensen et al., 2008;Wyrtki, 1962). The OMZ of the ESP is quite intense, there, regions of the open ocean may frequently reach values of DO lower than 20 μM, in the suboxia range (Fuenzalida et al., 2009;Naqvi et al., 2010;Paulmier et al., 2006;Paulmier & Ruiz-Pino, 2009), or even anoxic zones may also be found (e.g. Schunck et al., 2013;Ulloa et al., 2012). Such a low DO concentration largely impacts marine communities and biogeochemical cycling (Helly & Levin, 2004;Niemeyer et al., 2017;Stramma et al., 2011).
During the last two decades several studies have focused on the physical variability, the biogeochemical processes and bacteria communities present in this OMZ of the ESP (e.g. Frenger et al., 2018;Levin et al., 2002;Montes et al., 2014;Morales et al., 1999;Paulmier et al., 2006;Stramma et al., 2010;Ulloa et al., 2012;Ulloa & Pantoja, 2009;Vergara, Dewitte et al., 2016), which have led to a better understanding of the natural and forced variability of DO within OMZs as a whole. However, most of the studies have been conducted in the tropical part of the ESP, off Peru and northern Chile, where the OMZ extends deeper and further offshore and experiences only relatively small seasonal variability. In this work we focus on the extratropical region, near the southern tip of the OMZ. In contrast to the tropical ESP, where zonal flow originating from the equatorial region dominates the surface and subsurface circulation feeding the OMZ (e.g., Montes et al., 2010;Stramma et al., 2010), close to the coast of Chile the subsurface circulation is governed by a poleward undercurrent, namely, the Peru-Chile Undercurrent (PCUC) (e.g., Silva & Neshyba, 1979;Strub et al., 1998). Since the pioneer work of Gunther (1936), the PCUC has been associated with the presence of oxygen depleted waters at intermediate depths along northern and central Chile. Subsequently, many other works have confirmed this relationship (Wooster & Gilmartin, 1961;Silva & Neshyba, 1979;Codispoti 1989;Hormazábal et al., 2006;Silva et al., 2009) and have related the variability of the OMZ to changes in the PCUC transport. This undercurrent transports oxygen depleted Equatorial Subsurface Water (ESSW) southward along the continental slope and shelf edge off Chile, which can be traced as far south as 48°S (Silva & Neshyba, 1979). The PCUC undergoes large variability at a wide range of timescales, standing out the intraseasonal one (from about 20 to 90 days) and the interannual that characterize the El Niño-La Niña cycles (e.g. Shaffer et al., 1997Shaffer et al., , 1999Pizarro et al., 2001Pizarro et al., , 2002Dewitte et al., 2012;. To date, the paucity of DO data has prevented a thorough assessment of the temporal variability of the OMZ off Chile associated with temporal and spatial changes of the PCUC. Direct estimations of the PCUC off Peru, based on data from current meter moorings and hydrographic stations, showed that the PCUC transports about 1 Sv at 10°S (Huyer et al., 1991). Nevertheless, recent shipboard ADCP observations conducted between 4°S and 18°S off Peru -between 2008 and 2012 by Instituto del Mar del Perú (IMARPE)-showed a much larger southward transport related to the PCUC (Chaigneau et al., 2013). According to those estimates the mean PCUC transport increases from about 1.8 Sv at 5°S to a maximum of 5.2 Sv at 15°S. South of this last latitude the transport reduces slightly (Chaigneau et al., 2013). PCUC transport estimation at 30°S -based on current-meter and hydrographic observations-showed a mean value of about 1.0 Sv (Shaffer et al., 1999). All the above studies - Huyer et al. (1991), Shaffer et al. (1999), Chaigneau et al. (2013)-emphasize that the seasonal variability of the PCUC is rather small compared to the large fluctuations observed at intraseasonal and interannual periods. Table 1 shows a comparison of the mean and seasonal transport of the PCUC off Peru and Chile obtained from different observational and model studies.
In contrast to that observed off Peru and northern Chile, off central Chile (~30°S-40°S) seasonal variability may be relatively more important due to the large increase of the seasonal fluctuation of the alongshore wind-stress and the wind-stress curl (e.g. Aguirre et al., 2012;Bakun & Nelson, 1991;Shaffer et al., 1999) associated to the seasonal migration of the energetic atmospheric low-level jet (Muñoz & Garreaud, 2005;Renault et al., 2009) and the dissipation of perturbations of equatorial origin. The relatively weak seasonal variability of the PCUC associated to marked latitudinal fluctuations questions the extent to which the OMZ characteristics are linked to the PCUC transport at seasonal timescale. This issue is hard to tackle with from low to medium resolution global models owing to their limitations in simulating realistically coastal upwelling dynamics and parametrisations in the biogeochemical models not necessarily tuned for eastern boundary current systems. High-resolution modeling has been proven more useful to address many aspects of the oceanic dynamics and circulation of the ESP-OMZ (Aguirre et al., 2012(Aguirre et al., , 2014Combes et al., 2015;Leth & Shaffer, 2001 The realism of the regional oceanic models allows in particular addressing quantitatively aspects of the physical-biogeochemical interaction and OMZ variability as shown by recent studies off Peru (Bettencourt et al., 2015;Montes et al., 2014;Vergara, Dewitte, et al., 2016). These latter studies have in particular shown a realistic representation of both spatiotemporal variability of the OMZ and the PCUC, and allows addressing the effect of the eddy-induced circulation on the OMZ variability. Thus, here we take advantage of these recent progresses in coupled modeling to investigate the seasonal variability of the southern tip of the eastern South Pacific Oxygen Minimum Zone between 30°S and 38°S. Our study focused on the role of the meridional changes of the PCUC transport and its seasonal cycle on the OMZ. We question the extent to which the PCUC transport controls the OMZ variability taking also into account other important aspects of the OMZ variability associated to mesoscale activity (Bettencourt et al., 2015;Vergara, Dewitte, et al., 2016). The rest of the paper is organized as follows: In section 2 we described the data used for model validation, followed by a brief description of the coupled physical-biogeochemical model. An extensive model assessment is included in the Appendix. In section 3 we present (1) a description of different metrics that characterize the southern tip of the ESP OMZ and its annual variability, (2) the seasonal variability of the PCUC and its relationships with the OMZ, and the role of the zonal transport of the DO in the OMZ, (3) an integrated seasonal budget of DO in the OMZ and the role of mesoscale variability, and (4) an estimation of the oxygen eddy fluxes. In section 4 we discuss our key findings and the main strengths and limitations of the presented results, and finally, a summary of our main finding and conclusion are presented in section 5.
2. Data, model description and methods

Data
Satellite-derived sea surface temperature (SST) and chlorophyll-a (Chl-a) concentration data used in this study were the level-3 products of the Moderate Resolution Imaging Spectro-radiometer-MODIS-Aqua mission available from the National Aeronautics and Space Administration Jet Propulsion Laboratory PO.DAAC (https:// oceancolor.gsfc.nasa.gov). These products have been pre-processed by PO.DAAC to obtain a daily resolution and a spatial resolution of~4 km (NASA Goddard Space Flight Center, O. E. L and O. B. P. G., 2014). Sea level height and geostrophic currents were obtained from the Archiving, Validation, and Interpretation of Satellite Oceanographic (AVISO; currently processed and distributed by Copernicus Marine and Environment Monitoring Service (CMEMS) with no changes in the scientific content; see http://marine.copernicus.eu). These data were used to evaluate the simulated velocities near the surface. A comparison of model and satellite derived eddy kinetic energy (EKE) can be found in Vergara et al. (2017) (their Figure ).
In situ data of temperature, salinity and dissolved oxygen concentration from different oceanographic cruises conducted off south-central Chile (35.5°-40°S) were also used to assess the model performance. The cruises were supported by the "Fondo de Investigación Pesquera" (FIP; MOBIO-BIO project) and were carried out during Nov-2004, Dec-2005, Oct-2006and Mar-2008. A total of 358 CTD profiles between 0 and 600 m depth were collected during those cruises and processed by the Physical Oceanography Group at the University of Concepcion. A monthly ship-based time series from the University of Concepcion at 36°30'S over the continental shelf (~90 m depth) for the period 2002-2008 was used to evaluate the seasonal variability near the coast at this latitude [see Sobarzo et al. (2007) and Escribano et al. (2012) for more details].
Climatological information of temperature, salinity, dissolved oxygen and nutrients for the study region were obtained from the Commonwealth Scientific and Industrial Research Organisation (CSIRO) Atlas of Regional Seas (CARS) gridded compilation (CARS 2009; www.marine.csiro.au/~dunn/cars2009/), it comprises gridded (0.5°latitude by 0.5°longitude every one month) fields of the different variables recorded over the period of modern ocean measurements. The seasonal cycle is based on an annual and semi-annual harmonic over the first km of the ocean. We use this atlas, along with other data sets, for the model assessment described in the Appendix.
Ocean current observations from a current-meter mooring located over the continental slope at 30°S, were used to assess the time variability of the model PCUC. The mooring was equipped with an upward-looking 300 kHz acoustic Doppler current profiler (ADCP) at~120m depth -set up with a bin size of 4 m-and 4 Aanderaa current meters (RCM 8) at 220, 330, 480 and 750 m depth. The sampling interval for all the equipments changed during the measuring period between 0.5 h and 1 h. These current time series started originally in November 1991 and was maintained until 2010 (the ADCP was added to the mooring in May 2003). Here comparisons between model and observed currents were performed during the period 2003 and 2006.

Coupled physical/biogeochemical model
The regional high-resolution coupled physical-biogeochemical model used for this study takes into account the main processes linked with the eastern boundary upwelling systems (EBUS) and associated OMZ, it extends up to the equatorial region (12°N) in order to grasp the connection with the equatorial variability through coastal wave dynamics and transport by the PCUC. However, the region analyzed in this study is limited to the southern tip of the OMZ, comprised from 30°to 38°S (Figure 1). The model encompasses the period 2000-2008.
The hydrodynamics was modeled by the Regional Ocean Model System (ROMS), an ocean model widely used for regional studies. It is a free-surface, terrain-following coordinate model with split-explicit time stepping and with Boussinesq and hydrostatic approximations that allows achieving adequate resolution to resolve mesoscale dynamics (see Shchepetkin & McWilliams, 2005, for a full description of the model). The ROMS simulation extend from 1958 to 2008 and has been validated using mean sea surface temperature data, EKE (based in geostrophic surface currents estimated from satellite altimetry) and vertical current structure (Dewitte et al., 2012;Vergara et al., 2017). Only the last nine years were used here. The complete domain extends from 12°N to 40°S, and from the coast to 95°W having a horizontal resolution of 1/12°(~8 km) and 37 sigma levels in the vertical. Most of the sigma levels (23) are distributed in the first 500 m depth, the remaining levels are distributed down to 4500 m depth, which allows a good representation of the OMZ, even in the deep ocean. The atmospheric forcing consists in a statistical downscaling product from NCEP-NCAR (2.5°x 2.5°) data that provides wind stress and wind speed (see Goubanova et al., 2011), this method refines the NCEP-NCAR resolution to 0.5°x 0.5°and correct for the biases of the NCEP reanalyses near the Peru-Chile coast. The latter is used in the bulk formula combined with COADS monthly climatology (1°x 1°) (da Silva et al., 1994), air temperature and humidity to estimate latent heat fluxes. The statistical model is constructed from the QuikSCAT data so that the atmospheric forcing is almost equivalent to the original data over the period 2000-2008, which also motivates to focus on this period. Other flux terms also are from the COADS. The lateral open boundary conditions for temperature, salinity and horizontal velocity were obtained from SODA 1.4.2 reanalysis (Smith et al., 1992). The SODA 1.4.2 horizontal resolution is 0.25°( lat) x 0.4°(lon) and with 40 vertical levels with 10 m spacing near the surface. For further description of the hydrodynamic model simulation and validation see Dewitte et al. (2012) and Vergara et al. (2017).
The ROMS model was coupled to a biogeochemical model named BioEBUS following similar methodology than Montes et al. (2014) and Vergara, Dewitte, et al. (2016). Hereafter to this physical-biogeochemical coupled we will refer as ROMS/BioEBUS. The BioEBUS model is a nitrogen-based model developed for EBUS Gutknecht, Le vu, Marchesiello, et al., 2013) derived from N 2 P 2 Z 2 D 2 model (Koné et al., 2005). This model consists of 12 compartments having two classes of phytoplankton ("small" and "large", the small class is representative of small-flagellates named PS, while the large is representative of diatoms named PL), two classes of zooplankton (representing small-ciliates ZS and largecopepods ZL) and two classes of detritus (small DS, large DL). Dissolved organic nitrogen (DON) is represented in one compartment, while dissolved inorganic nutrients are represented by nitrate, nitrite and ammonium. The biogeochemical initial and open boundary conditions for nitrate and DO were obtained from CSIRO Atlas of Regional Seas (CARS, 2006). The nitrous oxide concentrations were estimated from DO using the parameterization given by Suntharalingam et al. (2000Suntharalingam et al. ( , 2012. The phytoplankton biomass is initially derived from satellite information (climatological data based on SeaWiFS observations) and vertically extrapolated using the Morel and Berthon (1989) methodology. Finally, initial and lateral boundary conditions for other biogeochemical tracers (NO 2 -, NH 4 + and DOM) are established using constant vertical profiles based on Koné et al. (2005), similar to those used by Montes et al. (2014) and Vergara, Dewitte, et al. (2016). In general, the biogeochemical parameters used in our model simulation were the same as those used by Montes et al. (2014) and Vergara, Dewitte, et al. (2016) for Peru and northern Chile. See appendix in Montes et al. (2014) for a quantitative description of the parameters used for BioEBUS in this study.
The evolution of any tracer in BioEBUS is determined by the advection-diffusion equation . In particular, for DO we can write: The first term on the right-hand side of (1) is the advection, where u is the fluid velocity, the second term is horizontal subgrid-scale diffusivity, where K h is the horizontal eddy diffusion coefficient -equal to 100 m 2 s -1 in this version of the model-, and the third term is the vertical mixing. (with a turbulent diffusion coefficient K z calculated using the K-profile parameterization mixing scheme; Large et al., 1994). Note that the model also has numerical diffusion associated with inherent spurious diapycnal mixing of the numerical scheme, so that K h has to be empirically adjusted. The last term, SMS (sources minus sinks), includes all biogeochemical processes considered by the model that act as sources and sinks, in this case, for DO. Note that the tracer equation (1) was evaluated online at each time step. The total physical term (named PHYS below) is the summed-up contribution of the advective term (ADV = -∇ · (uDO)) plus the horizontal diffusion (K h ∇ 2 h DO) and vertical ( ∂ ∂z K z ∂DO ∂z ) mixing, hereafter called Hmix and Vmix, respectively. The ADV can also be represented separately through its contributions related to the different components as: the zonal (Xadv = -u∂(DO)/∂x), meridional (Yadv = -v∂(DO/∂y) and vertical (Zadv = -w∂(DO)/∂z) advection e.g. , where u, v and w are the zonal, meridional and vertical velocity components, respectively. Note that ADV contains also diffusion inherent to the numerical advection scheme. Hmix and Vmix are explicit source of mixing, however, they are not the only term contributing to mixing in the model. The term mixing will refer to the integrated effect of all processes contributing to mixing directly or indirectly. Besides horizontal diffusion and vertical mixing, mixing can be also caused by non-linear advection. e.g. eddy fluxes (cf. Vergara, Dewitte, et al., 2016). Here non-linear advection (of DO) refers to the cross products of the anomalous velocity field (u') and the anomalous DO (DO') -i.e. in the context of Reynolds decomposition, this corresponds to the term (u'DO') where prime stands for an anomaly with regards to any preferred baseline (either mean climatology or 3-month running mean in this study). Note that equation (1) would involve the divergence of this term. Linear advection refers to either the mean advection of anomalous gradient in DO (u∂DO'=∂z) or anomalous advection of mean gradients in DO (u∂DO =∂z). Here also mean and anomalies are defined by the a priori Reynolds decomposition that is considered. In this context the mean quantities are "constant" while the anomalies are variable in time.

Characterizing the extratropical OMZ
At the present time, there is still no full agreement about the best threshold in oxygen that defines the OMZs.
The first global study of hypoxic waters considered a DO threshold of <8 μM (Kamykowski & Zentara, 1990). However, Karstensen et al. (2008) used three thresholds to analyze the characteristics of the OMZ in the eastern tropical Pacific and Atlantic: DO~4.5 μM in the suboxic range, DO~45 μM, and a more relaxed level of 90 μM in order to include in their analysis the less intense OMZ of the eastern Tropical Atlantic. Paulmier and Ruiz-Pino (2009), in contrast, excluded the Atlantic OMZ when they used (~20 μM) to describe the different OMZs using the WOA2005 climatology. Helly and Levin (2004) used also a threshold of~20 μM in their global study of the OMZs on continental margin seafloors. In contrast, Stramma et al. (2012) used a much larger value (~150 μM) to analyze the expansion of the OMZ in the tropical Atlantic.
Here we define the OMZs as the region where DO levels are ≤ 45 μM, because this threshold value is in the range of values characterizing hypoxia (e.g. Naqvi et al., 2010). For more practical reason, we also note that waters with DO~45 μM are present most of the time inside our domain of interest. Also, as the model overestimates somewhat DO in the southern part of the domain (see Appendix), a lower value, like DO <20 μM, would narrow the extent of the OMZ within our domain compared to the observations. Considering that low value in DO concentration (i.e DO< 20 μM) is required to stimulate denitrification, this model bias needs to be kept in mind when interpreting the results (Section 3).
Finally, it is worth noting that in our analysis we use the volume of the OMZ as a metric to characterize its seasonal variability, but for analyzing the budget of DO we use a fixed volume defined as the region with 9year (2000-2008) mean DO ≤ 45 μM.

Estimation of DO eddy fluxes
To evaluate the contribution of the eddy-induced circulation on the DO seasonal variability, we consider a Reynolds decomposition of the DO and the three velocity fields, where the departure from the "mean state" accounts for fluctuations having a typical time scale lasting from weeks to a few months so that eddy variability refers here mostly to the mesoscale turbulence. This writes as follows: where ∅ represents DO or any velocity component, e ∅ is the 3-month running mean and ∅′ is the deviations from e ∅. Time series with a temporal resolution 5-day mean model output are used for this calculation. To estimate the mesoscale eddy fluxes of DO we compute the covariance between horizontal velocity fluctuations (u', v') and DO', i.e. ⟨u′DO′⟩ 3-m and ⟨v′DO′⟩ 3-m , where "⟨⟩ 3-m " stands for a 3-month moving average. Hereafter, for simplicity, we refer to eddy fluxes for theses quantities (instead of mesoscale eddy fluxes). The methods to derive the eddy flux is similar to that used in Vergara, Dewitte, et al. (2016).
We thus obtained covariance time series every 5 days, from which we estimate the climatology. In this paper, we focus on the seasonal cycle of the horizontal eddy flux, ⟨u′DO′⟩ 3-m and ⟨v′DO′⟩ 3-m , that consists in the climatology of the eddy flux (⟨u′DO′⟩, ⟨v′DO′⟩) where the braket indicates here the mean over the entire period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008). The latter will be compared to the mean seasonal flux that corresponds to the flux associated

Seasonal and spatial variability of the southern tip of the OMZ
As a first step, we document the main characteristics of the southern tip of the OMZ and its seasonal variability, including dynamical and morphological features that are used for the interpretation of the subsequent diagnostics. Along the Chilean coast the zonal extension and thickness of the OMZ (represented here by the region with DO ≤ 45 μM or~1 mL L -1 , i.e. near the suboxia range, e.g. Naqvi et al., 2010) rapidly decrease southward, but near the coast it remains at intermediate depth as far south as 38°S ( Figure 1).
Observations show that the OMZ may extend even farther south, exceeding the southern border of our study region (e.g., Silva et al., 2009). The tongue of low-oxygen waters is rather well associated with a tongue of high salinity (salinity > 34.55) and both are largely reduced southward by mixing with low-salinity and well ventilated AAIW and SAAW (Figures 1b and c, and Figure A4; cf. Silva et al., 2009 andLlanillo et al., 2012). Although on average the offshore extension of the OMZ decreases southward , its offshore boundary zigzags with maximum ( The seasonal cycle of salinity and DO are inversely correlated inside the whole OMZ. To illustrate this, we present results for 3 cases, in all of them both DO and salinity time series are obtained by spatiallyaveraging the corresponding variable in the whole OMZ region. The 3 cases are the following: 1) Just the 12-month time series that conform the seasonal cycle (r = -0.87; this represent the correlation between the time series shown in the lower panels of Figures 4d and 4e); 2) 9-year (2000-2008) monthly time series of DO an Salinity (r=-0.97), and 3) similar to case 2, but only the region between 35°S and 38°S was considered (r=-0.83). These last 2 correlation values are significantly different from zero at the 95% confidence level (using an integral time scale of 3 months; Emery & Thompson, 2004). South of~36°S, the seasonal cycle of DO increases and it differs from that observed in the northern part of the study region ( Figure 4). For instance, the minimum DO found during fall in most of the domain is not present south of 36.5°S ( Figure 4d). In this region the relative contributions of the different water masses present inside the OMZ display an abrupt change. ESSW decreases from about 64% at 36°S to~40% at 37°S, while AAIW increases from~30% to 40% between those latitudes. The increase of SAAW in this region accounts for the other 14% ( Figure 1c).
The above results suggest that the seasonal and spatial variability of DO inside the model OMZ are related to changes in the proportion of ESSW more than to biological processes. The OMZ is less intense and it volume is smaller during spring ( Figure 2 and Figure 3; in this context a less intense OMZ means that the concentration of DO inside the OMZ, i.e. inside the volume for which DO < 45 μM, is larger). It is worth noting that the PCUC is weaker in austral winter (August-September), just a couple of months before the decrease in intensity of the OMZ is less intense (Figure 4). In the next subsection we analyze in more detail the variability of this poleward flow and it relationship with the OMZ.

Seasonal variability of the PCUC and the OMZ
The meridional component of the velocity within the OMZ is mostly southward year-round, but around August-September this poleward flow is largely reduced and may be slightly northward. Conversely, maximum poleward currents are observed near late April-May ( Figure 4f). North of 36.5°S this flow shows an important semi-annual component, with a secondary maximum in December, and a secondary minimum during the end of the austral summer. This semi-annual component in the poleward flow has been well documented based on direct current observations at 30°S and its forcing has equatorial and local origins (cf. Shaffer et al., 1999;Pizarro et al., 2002;Ramos et al., 2006; see section 4). Note that the OMZ is located in a region dominated by the PCUC (see white contours on Figure 3) and the weakness of this flow -centered around August (Figure 4f)-is followed by an increase of the DO in the OMZ about 2 months later ( Figure 4d). Conversely, the maximum poleward flow centered around April-May coincides with a reduction of the mean DO.
In our model the mean position of the poleward undercurrent does not always coincide with the mean position of the core of the OMZ (Figure 3). The undercurrent commonly extends deeper than the 45 μM isopleth, and this isopleth may extend farther offshore than the undercurrent, particularly in the northern part of the 10.1029/2019JC015201 Journal of Geophysical Research: Oceans study region. We estimate the transport of the PCUC by spatially averaging the southward flow between 80 and 800 m depth, and between the coast and 200 km offshore (only values smaller than -0.02 m s -1 were considered in the calculation). The mean transport of the PCUC decreases southward from about 1.2 Sv at 30°S to 0.5 Sv at 37°S (Figure 5a). The former value is slightly larger than the~1.0 Sv estimated using direct currentmeter observations and hydrographic data by Shaffer et al. (1997Shaffer et al. ( , 1999) over a different period. The mean southward transport of the PCUC in the study region was~0.7 Sv, consistent with the mean value estimated by Aguirre et al. (2012) and Vergara, Echevin, et al. (2016), both estimates are also based on numerical simulations using the same source of data for the atmospheric forcing (i.e QuikSCAT data). These studies also showed that the PCUC transport decreases southward and has a minimum value in winter off central Chile. Unfortunately, there are not enough published data available to contrast these model estimations southward of 30°S (see Table 1). The variability of the PCUC transport in our model showed an important semi-annual component, particularly north of 36°S, consistent with the semi-annual variability described above and also observed in the alongshore flow over the slope at 30°S (see Figure 5a). South of 36°S the transport of the PCUC is drastically reduced (Figure 5a). A similar, relatively abrupt, change was observed in the meridional velocity and DO concentration inside the OMZ (Figures 4d and  4f), but somewhat farther south. The seasonal cycle of DO averaged inside the whole OMZ showed a maximum of~32 μM in austral spring (i.e. October-November, see bottom panels of Figures 4d), about 2 month later than the PCUC velocity reached a minimum (see bottom panels of Figures 4f). This result is consistent also with the estimated transport of the PCUC and with the DO averaged inside the region dominated by this current ( Figure 5).
To further illustrate the PCUC control on the OMZ variability we show the variability of DO inside the core of the PCUC (Figure 5b). The mean DO concentration within PCUC core varies from~26 μM at 30°S to > 60 μM at 37°S. This concentration tends to be smaller when the southward transport is weaker, that is in part, because when the PCUC is weaker its core tends to be restricted to the core of the OMZ (Figure 3). On the other hand, the minimum poleward transport of the PCUC in August is related to a rapid increase in DO concentration particularly south of 34°S (Figure 5b).

Seasonal variability of the zonal flow and the OMZ
The zonal components of the velocity, averaged inside the OMZ region (which is near cross-shore in our study region), shows small seasonal variability (~0.5 cm s -1 ) with an important semi-annual component (Figure 4c). Given the large cross-shore gradient of DO (Figures 3 and 4), small changes in the zonal speed may contribute to modulate the DO seasonal cycle in the OMZ through anomalous zonal advection of mean oxygen gradient. The most striking feature of the zonal speed is its large spatial variability, with the presence of alternating zonal jets of positive (eastward) and negative (westward) currents. For instance, bands of eastward speeds are centered around~31°S, 34°S and~37°S, although their latitudes may change somewhat during the year. These zonal jets are closely related with the offshore extension of the lateral boundaries of the OMZ (contrast Figure 2 with 4c), that is, where the cross-shore velocity is predominantely eastward, the oxygen isopleths that bound the OMZ approaches the coast, while at those latitudes where the velocity is predominantly westward, the lateral OMZ boundary extends further offshore.
From figure 4c we can distinguish 3 bands of zonal velocity (here we refer to the zonal velocity averaged inside the OMZ at each latitude: u OMZ (Lat,t), where Lat is latitude and t is time) dominated by positive (onshore) flow. The first band is centered around 31°S and has a maximum (onshore) flow in winter (u OMZ~1 .5 cm s -1 , this value is the average of the velocity in the band from 30.3°S to 31.6°S), the second band is located between~33°S and~35°S with a small maximum also in winter (u OMZ~0 .5 cm s -1 ) and the third one, laying around 36°S and 37°S, with a small maximum in summer (u OMZ~0 .8 cm s -1 ) ( Figure 4c). Thus, the onshore flow is somewhat more intense in winter north of 36°S, when the poleward flow is minimum, contributing to the increased DO -and to the reduction of the OMZ offshore extension. The increase in the onshore transport of DO and the weakening of the PCUC result in a well-defined maximum of DO inside the OMZ and a reduction of its total volume during spring. Despite the fact that the amplitudes of the seasonal cycles of u OMZ averaged inside these bands are small, they are significant at the 90% confidence level. To calculate the significance of the amplitude of the seasonal cycle we use a Monte Carlo method. It involves the generation of an artificial time series by randomly sampling the original one (with replacement), and to fit an annual harmonic to the artificial time series. This process was repeated 1000 times. Then, based on the 1000 artificial amplitudes a P% percentile was calculated. If the observed amplitude was larger than this percentile the amplitude was considered to be significant at (100-P)% of confidence level.
To summarize, our analysis suggests that the mean OMZ (i.e. zonal and meridional extent) is heavily constrained by the mean circulation as evidenced by the existence of a "high salinity tongue" matching the limits of the OMZ and the meridional variability of the off-shore limit following approximately latitudinally alternating bands of weak and stronger zonal "jets". The seasonal variability of the OMZ seems also tightly linked to the variability of the PCUC with a stronger PCUC in austral spring associated to the intensification and southward extension of the OMZ. There is also a significant contribution of the semi-annual cycle to the seasonality in both the PCUC and OMZ structure. While the PCUC variability appears as a key driver of the OMZ seasonal change, our description of the OMZ mean state and seasonal variability also suggest that other processes are occurring simultaneously that are documented below in the light of the results of a DO budget.

DO budget 3.2.1. Integrated seasonal budget of DO in the OMZ
To analyze the contribution of the different drivers to the DO seasonal cycle we consider the seasonal variability of the RHS terms in the advection-diffusion equation (1). The tendency terms were averaged inside the volume delimited for the mean position of the 45 μM isoline (i.e. inside the fixed volume that encompass the annual mean of the OMZ shown in Figure 1a). We present the climatological variation of the different terms including their mean value to assess their relative magnitude ( Figure 6). The zonal advection (Xadv), SMS and the horizontal mixing (Hmix) are negative, while meridional advection (Yadv), vertical advection (Zadv) and vertical mixing (Vmix) are positive. Xadv shows the largest magnitude, while Yadv is the second largest. In contrast, the mixing terms (Vmix and Hmix) are, in the total average, 2 orders of magnitude smaller than the other terms conforming the RHS of (1) with vertical mixing (positive) larger in absolute value than horizontal mixing (negative). The amplitude of the annual harmonic of DO (fitted to the dashed black curve shown in Figure 6a) is 4.3 μM and represents about~14% of the mean DO (30.7 μM) concentration inside the whole OMZ control volume. This harmonic explains~82% of the seasonal cycle of DO and has 10.1029/2019JC015201 Journal of Geophysical Research: Oceans a maximum in austral spring. The seasonal cycle of DO within the OMZ volume is also represented in Figure 4d (bottom panel), but in that case the amplitude of the variability is much smaller because the volume changes according to the seasonal variability of the 45 μM isopleth, so we are only considering waters with DO lower than 45 μM. The rate of change of the DO (∂DO/∂t, red line in Figure 6a) peaks in July, about 2 months before the DO peaks (Figure 6a, red line), when the DO is increasing most rapidly. There is only a qualitative consistency when the averaged seasonal cycles of DO and ∂DO/∂t are compared, revealing the relative importance of non-linear terms (i.e non-linear advection and mixing). The seasonal cycle of ∂DO/∂t has a secondary maximum in April, which accounts for the presence of semi-annual variability. The main peak in this series is in August and may be associated with the weakening of the southward transport of DO (cf. Figure 4f).
The summed-up contribution of the physical tendency terms (PHYS formed by ADV + mixing) is the largest contributor to the seasonal cycle of ∂DO/∂t within the OMZ (Figure 6b). Consistently, PHYS and ∂DO/∂t time series show similar seasonal variability, which emphasizes the prominent role of physical processes in controlling the averaged seasonal variability of DO inside the OMZ. Furthermore, the ADV terms, and in particular, the zonal advection of DO (Xadv; Figure 6b) dominates the seasonal variability of PHYS, while the mixing terms (Hmix, Vmix) are much smaller ( Figure 6c) and their mean contribution to the seasonal variability of DO (averaged inside the OMZ) is only significant during the transition phase when ∂DO/∂t changes sign or remains weak (like in May-June). It is worth noting that ADV may also generate mixing (i.e. eddy-induced advection; see section 3.3) associated with mesoscale turbulence, so Hmix and Vmix are mainly related here to small-scale (subgrid-scale processes) mixing.
The biogeochemical term (SMS in Eq. (1)) has small seasonal variability compared with PHYS and it permanently behaves as a sink of DO inside the OMZ with a single maximum in austral winter. The yearly means of PHYS and SMS have similar magnitude (1.8 and 1.4 x10 -6 μM s -1 , respectively), but opposite signs. On average, the meridional (Yadv) and vertical (Zadv) advection behaves as oxygen sources, while the zonal advection (Xadv) acts as a sink for the whole volume. The amplitude of the annual and semi-annual harmonics adjusted to the different terms showed in Figure 6 is provided in Table 2.
The sum of all averaged terms shown in Figure 6 does not necessarily have to be zero, even if ∂DO/∂t was zero, due to the important influence of non-linear terms that make up PHYS. In fact, the magnitude of PHYS exhibits large mesoscale variability (see below), which suggests a relevant contribution of non-linear eddy fluxes -related to mesoscale processes-to the variability of the OMZ. Thus, the spatially-averaged terms inside the whole volume may not represent properly the characteristics of the seasonal cycle of the DO budget. On the other hand, the mean of ∂DO/∂t inside the whole volume for the period 2000-2008 was of the order of 10 -8 μM s -1 . This small trend can be associated with decadal variability present in the study region (note that, for the full simulation period, 1958-2008, the model is stable and the mean ∂DO/∂t is practically zero). Table 3 shows the total mean values inside the OMZ for the period 2000-2008 and the climatological variability (standard deviation) of ∂DO/∂t and ADV, including the Xadv, Yadv and Zadv components (all values are in μM s -1 ×10 -6 ). The explained variance of ∂DO/∂t in term of the advection (R 2 , in percentage) is also shown. As the different nonlinear advection components are correlated among them (not shown), they cannot be regarded as independent contribution to ∂DO/∂t. In general, all advective terms contribute to ∂DO/∂t in the OMZ, but Xadv seems to play a major contribution as also was suggested by Figure 6b.

The role of the mesoscale activity
While the previous analysis indicates seasonality in the processes, there are also indications of high frequency fluctuations which may result from the intricate circulation in relation with mesoscale features. To visualize the spatial structure of the seasonal variability associated with the different terms that control the oxygen budget inside the OMZ we use an EOF analysis applied to the seasonal fields of PHYS (particularly Xadv, Yadv and Zadv) and SMS. These fields were previously averaged between the isopycnal surfaces of 26.4 and 26.8 kg m -3 , where the core of the OMZ is located (Figure 3). Note that the seasonal time series of each field are monthly i. e each time series have only 12 values. The analysis was restricted to the region delimited by the mean position of the 45 μM isopleth, consistently with the above analysis (i.e. we remove from the analysis the white area showed in Figure 7. To evaluate the sensitivity of the EOF results to the selected region we also considered a wider region, including the region outside the OMZ. The structure of the EOF of the different terms was rather similar with most of the variance near the coast, i.e. where the OMZ is observed (not shown).
The first EOF mode for PHYS explains 16% of the total variance and its spatial structure shows large variability with length-scale typical of mesoscale processes (Figure 7a). The explained variance for the first EOF mode of different ADV terms was around 20% (Figures 7b to d). The EOF mode 1 of Xadv and Yadv are spatially well anti-correlated (r =-0.7; Figure 7b and c). Eddy flows and their rectification on the mean flows seem to be relevant, contributing to this large spatial variability. Time series (EOF principal components) of Xadv and Yadv show an important semi-annual component with maxima in March and July-August, coincident with the periods of weaker poleward undercurrent (Figure 4f). In contrast to PHYS, the first EOF of SMS showed a uniform (in phase) structure in the whole study region and its time series showed a well-defined annual variation with a maximum in July. This EOF mode explains 72% of the total variance. In this density range (26.4 and 26.8 kg m -3 ) SMS was always a sink for DO, so the EOF time series indicates that DO consumption was minimum during austral winter.
In order to gain further insight into the spatial variability of ADV inside the OMZ, we perform an EOF analysis in three cross-shore coastal sections: 30°S, 34°S and 37°S, delimited by the 45 μM and the first 150 km (Figure 8). The first EOF mode of all advective terms showed rather similar structures, with maxima near the continental slope, just inside the core of the PCUC (dashed black contour in Figure 8). Nodal lines, centered around 250 m depth for Xadv and Zadv, and 150 m depth for Yadv, are observed in the different sections ( Figure 8). The principal components of Xadv in the 3 locations showed a rather consistent semi-annual fluctuation ( Figure 8m) similar to that observed in the horizontal EOF analysis. Conversely, Yadv and Zadv time series were quite variable (Figure 8n, o). The first EOF mode of SMS explain a large fraction of the total variability and was vertically in phase in all sections, with larger amplitude toward the upper limit of the OMZ. Their time series are similar to those observed in the horizontal sections with a maximum in July, but the time series from 30°S, which showed a maximum centered around April.

Journal of Geophysical Research: Oceans
To summarize, the three advective terms (Xadv, Yadv and Zadv) are the main contributors to the seasonal variability of the DO budget in the OMZ, these terms have large and correlated mesoscale variability and vertical structure with a node typically centered around 200 m depth. The typical annual change is not dominated by one harmonic, it shows large semiannual and shorter variability. These results, along with previous studies from the northern region of the eastern South Pacific OMZ [cf. Montes et al., 2014;Bettencourt et al., 2015;Vergara, Dewitte, et al., 2016), as well as from other OMZ (Resplandy et al., 2012), suggest that the DO budget in our study region is impacted by mesoscale processes. The formalism of the oxygen budget does not allow an estimate of the separate contribution of mesoscale variability to the rate of change in DO since eddies are embedded into the total circulation from which advection terms are derived. To get insights on the role of mesoscale variability in shaping the seasonal cycle of the OMZ, the DO eddy flux can be estimated (see section 2.4) and compared to the mean seasonal flux (i.e. the DO flux associated with the seasonal variability of the "mean" circulation) to infer its role on the OMZ seasonal cycle. : First we note that the mean eddy flux, although smaller than the mean seasonal flux, is of the same order of magnitude, with regions (e.g. north of 34°S) where the amplitude of the mean eddy flux is comparable to that of the mean seasonal flux. MEFs showed large spatial variability with lower values (<0.3 m s -1 μM) in the OMZ. In the offshore region outside the OMZ the fluxes were larger (~1 m s -1 μM) and predominantly eastward. The regions with larger eastward MEFs, like those centered around 32°S and 34.5°S, are related to a smaller offshore extension of the OMZ, as can be observed following the 90 μM contour (Figure 9a). Thus, the MEFs seem to contribute to ventilate the OMZ in some regions, particularly during summer and spring (Figure 9b to e), and to shape the offshore extension of the OMZ. The mean horizontal DO advective fluxes (e u g DO and e v g DO) are, in general, larger than the MEFs and they have also large spatial variability with eddies and meanders in the offshore region. Near the coast these fluxes are mainly southward, consistent with the predominance of the PCUC in this density range (i.e. 26.4 and 26.8 kg m -3 ). The seasonal anomalies of these advective fluxes during the different seasons are also consistent with the shape of the 45, 60 and 90 μM contours showed in Figures 9g to 9j.
To better visualize the relative contributions of the horizontal fluxes of DO to the oxygen budget of the OMZ and its seasonal variability we show the mean flux across the 60 μM surface between 26.4 and 26.8 kg m -3 isopycnals ( Figure 10). The mean position of this surface is shown in Figure 9 as a black contour (surfaces of 90 μM and 45 μM, are also shown there). We use the 60 μM surface because it better represents the boundary of the OMZ than the 90 μM one, and it is present in the complete domain most of the year, while the 45 μM surface is not present in the southern part of the domain during some periods or is restricted to a small range of depths and does not cover the entire chosen density range. Results showed that the MEFs of DO

10.1029/2019JC015201
Journal of Geophysical Research: Oceans predominantly transport DO toward the OMZ, particularly in the region located northward of 34°S. These fluxes have small annual variability (about 20% of the total variability of the DO fluxes) and they do not have a clear seasonal signal (Figure 10a and c). Conversely, the mean advective fluxes of DO are predominantly negative (Figure 10d), except in August, September and October (Figure 10f, bottom panel). The seasonal cycle of this flow has also a semi-annual component that is coherent with the seasonal variability of the zonal current shown above (contrast Figures 4c and  10f). Thus, in average for the whole study region, while zonal MEF contribute to transport DO toward the OMZ, the mean advective flows tend to transport DO in the opposite direction and both fluxes have the same order of magnitude, but the second one showed a larger seasonal variability, with minimum in fall and spring consistently with the seasonal variability observed in different metrics of the OMZ described above (Figure 4) and with the total DO budget ( Figure 6). These results also show the relevance of the MEFs in ventilating the OMZ.

Discussion
The lack of observations of both dissolved oxygen and currents off central Chile has limited the ability to evaluate the relationship between the PCUC and the OMZ and to explore the role of other concurrent processes that may also contribute to modulate the OMZ in this region. Here, we used a coupled physical/biogeochemical model to tackle this problem.
Although the model may have important limitations associated to uncertainties in the parameter values of the model parametrisations or to inaccurate representation of key biogeochemical processes taking place in the OMZ , it provides a powerfull tool to explore quantitatively the diversity of forcing mechanisms associated to seasonal variability of the OMZ in a region where observations are crucially lacking.
The model assessment showed that the model has biases in both mean state and seasonality (Appendix). The extent to which the model bias could influence our results would require undertaking sensitivity experiments to a number of processes, which is beyond the scope of the present work. The improvement of the mean state could be achieved through using more realistic wind forcing in particular near the coastal region where along-shore winds of most products including scatterometer data are stronger than observed (cf. Astudillo et al., 2017) producing a cool bias near the coast. The realism of the model simulations could be also improved through tuning the biogeochemical parameters to the study region since the model parameters were originally best fitted for simulating the OMZ off Peru (see Montes et al., 2014; Table A1) (e.g. maximum growth rate of phytoplankton, hydrolysis rate of detritus, nitrification and anamox processes among others). On the other hand, an improvement of the seasonal variability would require to include river run-off and consider air-sea interacction at mesoscale (e.g. Oerder et al., 2018). Sediment processes, that can be very relevant for the dynamics of DO near the coast (Bianucci et al., 2012;Siedlecki et al., 2015), could be also included to improve the realism of the model in term of DO near the coast. Despite of these limitations, we consider that our regional ocean simulation is realistic enough for carrying out process studies focused on seasonal timescale.
Dissolved oxygen and salinity in the model simulation were spatially well correlated inside the OMZ and both co-varied simultaneously throughout the year. Most of the DO variability was associated with

10.1029/2019JC015201
Journal of Geophysical Research: Oceans changes in water mass composition rather than with changes driven by biogeochemical processes. Nevertheless, according to our results those changes are not only driven by changes in the PCUC, as was firstly hypothesized, but by zonal advection and by eddy fluxes associated to mesoscale variability wich contribute to ventilate the OMZ. In the following we discuss in detail this and the rest of our main findings.
The seasonal cycles of the PCUC showed to be well correlated with all the metrics used here to characterize the OMZ (Figures 4 and 5) and also with several terms that are relevant in the DO budget of the OMZ (Figures 6a, 6b, 7f, 7h and 10f). Thus, our results support the interpretation that the variability of the PCUC is an important driver of the OMZ seasonal cycle. A striking feature of the seasonal cycle of the PCUC is its semi-annual variability. In fact, north of 36°S both the annual and semiannual harmonics of the undercurrent have about the same amplitude. This semi-annual variability was also evidenced in the metrics of the OMZ (Figures 4 and 5) and in the dominant advective terms in the DO budget of the OMZ (Figures 6a, 6b, 7f, 7h and 10f). For instance, the OMZ volume (DO concentration inside the OMZ) has two maxima (minima). This volume reaches the absolute minimum -and the DO concentration inside the OMZ reaches an absolute maximum-about one or two months after the PCUC southward transport reaches its absolute minimum. According to Pizarro et al. (2002) the semi-annual component of the PCUC seasonal cycle is of equatorial origin, while the annual harmonics results from a combination of equatorial and local wind forcing (see also Shaffer et al., 1999). Nevertheless, the seasonal variability of the alongshore wind varies in phase and amplitude along the west coast of South America, which may also introduce semi-annual variability in the PCUC. Off Peru these winds are maximum in winter, while off Chile they are more intense during spring and summer. Off northern Chile (between~18°S and 25°S) the seasonal cycle is weaker and increases toward the central coast off Chile associated to the presence of an atmospheric coastal jet that develops in spring and summer (Muñoz & Garreaud, 2005). This jet migrates meridionally during the calendar year and thus may also induces semi-annual alongshore variability of central Chile (cf. Dewitte et al., 2008, their Figure 7a). Note. The values of R 2 (explained variance in percentage) among ∂DO/∂t and the different advection terms are also shown.
Abbreviations: ADV = advection term; Xadv = zonal advection; Yadv = meridional advection; Zadv = vertical advection.  (1) Integrated in the Omz, as Showed in Figure 6 Annual harmonic Semiannual harmonic Note. The physical term (PHYS) is the summed-up contribution of all advection and mixing terms. The advection terms are the zonal (Xadv), meridional (Yadv), and vertical (Zadv) advections. The biogeochemical fluxes (SMS) represents the "source minus sink" contribution to the rate of change of DO (∂DO/∂t) due to biogeochemical processes. The subgrid mixing terms are the summed-up contribution of the horizontal diffusion (Hmix) and vertical diffusivity (Vmix). The DO amplitude is in μM and the terms of the equation (1)

Journal of Geophysical Research: Oceans
In contrast to the semi-annual variability observed in the different metrics representing the core of the OMZ, over the continental shelf the seasonal cycle of the oxycline, and oxygen concentration in the whole water column, is dominated by a single annual harmonic. There, the minimum depth of the oxycline coincides with the season of maximum upwelling-favorable winds (cf. Figure 5a). This is also observed off northern Chile between around 20°S and 24°S (e.g. Reyes et al., 2007, see their Figure 158). Thus, seasonal changes in DO over the continental shelf seem to be more tightly related to the annual upwelling cell than to the drivers that control DO in the OMZ core. Nevertheless, different ship-based, monthly, time series from northern and central Chile show also large month to month and interannual variability of the oxycline not well related to the local upwelling favorable winds, but to coastally trapped waves of equatorial origin, which may largely disturb the pycnocline, as well as the oxycline and the PCUC (Shaffer et al., 1997, Morales et al., 1999Ulloa et al., 2001;. Thus, the dynamics of the OMZ over the continental shelf seems to be different from that described here for the core of the OMZ. The variability of DO over the continental shelf and its relationship with the coastal upwelling and the offshore OMZ was not specifically addressed in our study, because it would require a model with better spatial resolution to correctly represent the continental shelf and the biogeochemical processes taking place therein. e.g. a sediment compartment which our model did not incorporate. Along with the alongshore advection, the cross-shore transport of DO plays also a key role in the seasonal cycle of the OMZ and so in the DO budget in our study region. This cross-shore or near zonal transport in the oceanic region was related to the presence of inserted zonal bands of positive and negative striations (Figure 4c). These striations have been observed in different regions of the global ocean through satellite altimetry (e.g. Maximenko et al., 2005), Argo float (Cravatte & Marin, 2012) and cruise data (Brandt et al., 2010), and have been reproduced in global and regional models (Belmadani et al., 2017;Nakano & Hasumi, 2005;Sinha & Richards, 1999). The presence of eastward and westward jets seems to result from the eddy field generated in the presence of planetary vorticity gradient and they efficiently develop near eastern boundaries related to particular coastal features (Belmadani et al., 2017;Davis et al., 2014). They may play also a role in ventilating the OMZs located at the eastern ocean boundaries as was shown for the tropical North Atlantic, where these meridionally alternating jets contribute to modulate the zonal advection of DO along the isopycnals surfaces at interannual and decadal timescales (e.g. Brandt et al., 2010;Hahn et al., 2017).
Our results suggest that this kind of zonal jets contribute as much as the PCUC to the transport of DO toward the OMZ and both play a major role controlling its seasonal cycle. These jets also shape the offshore boundary of the OMZ modulating its extension. Regions that are dominated by westward (eastward) velocities extend (reduce) the offshore distance of the OMZ boundary. Thus, these meridional jets add a dimension to the paradigm of PCUC-driven OMZ variability, which was the starting point of this study. On average,

10.1029/2019JC015201
Journal of Geophysical Research: Oceans the zonal currents transport low oxygen waters offshore extending the OMZ, except during end of winter and beginning of spring (Figure 10f), when the zonal advection of DO is toward the OMZ, contributing along with the weakening of the PCUC to reduce the total volume occupied by the OMZ in the study region and to increase the DO concentration inside it. In contrast to the mean westward transport of DO due to the zonal advection, mesoscale eddies fluxes contribute, on average, to transport DO eastward, i.e. toward the OMZ, during all the year (Figure 10c) playing also a key role in the dynamics of the OMZ. A relevant uncertainty related to eddy fluxes results from the role of subsurface mesoscale eddies. Those eddies are ubiquitous in the study region (Colas et al., 2012;Combes et al., 2015;Hormazábal et al., 2013), and they are not well enough captured by satellite data nor in situ data which is too sparse to have a reliable quantification of them. The role of these subsurface eddies in the offshore transport of oxygen depleted waters and in the overall dynamics of the OMZ and on its biogeochemistry is poorly understood (e.g., Frenger et al., 2018). Recent studies from the west coast of South America have shown that the PCUC plays also a key role in the generation of subsurface mesoscale eddies (e.g. Thomsen et al., 2016). Thus, the PCUC contributes directly to the southward transport of oxygen depleted ESSW and indirectly to its offshore transport by promoting the generation of subsurface mesoscale eddies.
The above discussion calls for further investigation and sensitivity experiment with the model considering that the representation of mesoscale activity in the model is critically dependent on a number of parameters, including horizontal resolution and the way the air-sea interaction at mesoscale is taken into account. This may have a significant influence on the net primary production and thereby on oxygen content through export (see Renault et al., 2016). While sensitivity experiments would be important to investigate further Journal of Geophysical Research: Oceans the role of eddy activity in shaping the OMZ off Central Chile, as a preliminary step, it may be necessary to consider other processes important in our study region. In particular, the biogeochemical model used here tends to overestimate the DO near the coast, which has been also observed off Peru for this particular model (Montes et al., 2014;Vergara, Dewitte, et al., 2016). Such a bias may be detrimental for our conclusions on the role of the PCUC on the OMZ variability since, in our study region, DO concentration is sensitive to the amount of oxygen advected southward along the slope by the PCUC. Sensitivity experiments to the northern boundary conditions of the domain considered here (i.e. 30°-38°S) may thus be necessary to evaluate the extent to which the amount of DO is critical for explaining the OMZ seasonal fluctuations. Another approach would be to tune the biogeochemical model for the study region considering that the model configuration used in our study was calibrated for the tropical part of the OMZ (Montes et al., 2014), and no specific calibration was performed to improve comparison among model outputs and observations off central Chile. In particular, in different OMZs the biogeochemical BioEBUS model was shown to be highly sensitive to the rate of organic matter decomposition and its vertical distribution (as well as to the sinking rate of particles; . Finally, it is worth noting that sediment processes, that can be relevant for the dynamics of DO near the coast (Bianucci et al., 2012;Siedlecki et al., 2015), are not included here, although sediment may act as an important sink of DO near the coast. Future work will need to consider these improvements in the model setup in order to refine our estimate of the physical control of the OMZ by the PCUC and eddy activity.

Summary and conclusions
A high resolution coupled physical-biogeochemical model was used to assess the seasonal variability of the southern tip of the Southeastern Pacific OMZ between 30°S and 38°S. The model was able to capture relatively realistically the mean structure of the OMZ and its spatial and temporal variability in the study region, although DO is somewhat overestimated near the coast and the very southern tip. We used different metrics to characterize the seasonal variability of the OMZ: (1) The amplitude of the typical seasonal variability of DO concentration inside the OMZ itself (i.e. inside the volume enclosed by the isosurface of DO = 45 μM) is about 2 μM (about 7% of the mean DO), with larger values in austral spring.
(2) The mean volume of the OMZ between 30°S and 38°S was~4×10 4 km 3 changing annually by~25%, from 4.5 ×10 4 km 3 during early austral winter to 3.5 ×10 4 km 3 in austral spring. Thus, we can conclude that the amplitude of the seasonal cycle of the OMZ is relatively weak. However, south of 35°S the volume of the OMZ is reduced, DO increases rapidly and seasonality becomes more relevant. The OMZ volume and its seasonal variability is highly correlated to subsurface waters with high salinity (salinity >34.5) in the study region consistent with the common idea that the OMZ is closely related to the warm, salty and nutrient rich ESSW.
Our analyses suggest that there is a close relationship between the PCUC transport variability and the seasonal cycle of OMZ. Intense poleward transport by the PCUC is related to a larger OMZ volume off southcentral Chile. However, this modulation would not be homogeneous alongshore due to other simultaneous processes such as disturbances in the transport by the meridionally alternating jets and eddy activity. These processes play also a role in shaping the mean characteristics of the OMZ, including changes observed in the offshore extension of the OMZ at different latitudes and the alongshore rate of decreasing of its volume. In addition, by evaluating the DO budget, it was determined that advective terms (ADV) are the main contributors (~90%) to the rate of change of DO, while subgrid mixing and biogeochemical (SMS) terms play a minor role in the seasonal variability of the OMZ off Chile. The seasonal cycle of these advective terms showed large mesoscale variability in the OMZ region, consistent with the significant contribution of non-linear advection, and accordingly, it was shown that mesoscale fluxes of DO play also a significant role in the seasonal variability of the OMZ.

Appendix: Model assessment
In this section we assess the model's realism in terms of the mean state and seasonal cycle of the OMZ. Many aspects of the circulation have been previously validated using remote and in situ observations in Dewitte et al. (2012) with a main focus on mean state and interannual variability (e.g., SST, EKE, the Peru Chile Undercurrent) off Peru. Vergara et al. (2017) complemented this assessment by validating the thermocline depth, EKE, coastal currents and sea level anomalies over the whole model domain. The coupled Journal of Geophysical Research: Oceans ROMS/BioEBUS simulation has been evaluated by Montes et al. (2014) and Vergara, Dewitte, et al. (2016), also with an emphasis in the OMZ off Peru. Here we focused our assessment off central Chile from 30°S to 38°S, which comprises our study region. Figure A1 shows the study region and the positions of several oceanographic stations and a currentmeter mooring located over the slope (at 1000 m total depth) near 30°S used for the model assessment.

Sea Surface Temperature and Chlorophyll-a
Monthly SST and Chl-a data from MODIS-Aqua satellite were used to evaluate the model skills in reproducing the observed mean fields and seasonal cycles of these variables ( Figure A2). The model reproduces the main features of the observed SST and Chl-a mean fields. Near the coast and north of 34°S the model showed larger mean values of Chl-a than the observations, while south of 35°S the model values were a bit underestimated. The coastal SST was about 0.8°C to 1°C colder in the model, which is largely due to an overestimation of the coastal winds used in the forcing (Astudillo et al., 2019). In the offshore region the differences between model and observed mean fields were, in general, much smaller than the differences observed near the coast ( Figure A2c and g). The seasonal cycle of the SST near the coast (first 80 km) showed also important differences with the observations. Consistently with the coastal cold bias associated to a stronger upwelling rate than that observed, the model overestimates Chl-a during most of the year. It does not reproduce an observed increase of Chl-a in austral spring as shown by observations ( Figure A2d). In contrast to the coastal region, the model and observed seasonal cycles of Chl-a and SST were rather similar offshore. The larger Chl-a and the lower SST observed in the model simulation near the coast may be associated with an overestimation of upwelling, which may be,

10.1029/2019JC015201
Journal of Geophysical Research: Oceans in turn, related to an overestimation of the wind stress near the coast when satellite-based winds are used. Experiments carried out by one of the authors showed that coastal SST is more realistic when a shoreward decrease of the wind stress near the coast is considered. This may also have an impact on net primary production and thereby on the export production and oxygen inventory (see Renault et al., 2016). Presently, we are working with a new coupled ROMS/BioEBUS simulation to evaluate the response of Chl-a when a realistic wind drop-off near the coast is considered.

Vertical structure of the variables and water masses
To evaluate the model ability to represent the mean vertical structure of temperature, salinity, DO and nitrate we use the CARS climatology as a benchmark. Noteworthy, oxygen data are rather scarce in the study region and the oxygen climatology from CARS may have also limitations, which should be kept in mind. The Figure A3 shows zonal sections of DO at 30°S, 34°S and 37°S. The model exhibits marked biases consisting mostly in a sea-sow pattern resulting from an overestimation of DO in the core and an underestimation below the OMZ ( Figures A3c, A3f, A3i). Despite this mean bias, the model is able to simulate realistically its mean depth and the positions of the upper and lower oxyclines, but the offshore extension of the OMZ (bounded by the isopleth of 45 μM in our study) is much smaller in the model than in CARS. The simulated vertical gradients of DO in the upper and lower oxyclines are also smaller than the observed ones, particularly in the lower part of the OMZ. The overestimation of DO is more pronounced in the southern part of the study region and may exceed~70 μM close to the coast. This misfit is analyzed further by comparing the model outputs and the data at a coastal site located at 36°30'S. Spatial correlation between model and CARS (using zonal sections that expand from the surface and 700 m depth and from the coat and 80°W) at different latitudes indicates a general good agreement for salinity, temperature and nitrate (figures not shown). Nevertheless, nitrate showed a large bias (5-10 μM) near the coast, between surface and 200 m depth. For temperature, salinity and nitrate the spatial correlations between simulated and CARS variables were larger than~0.85 at all latitudes ( Figure A3j), while for DO, the correlations tend to decrease southward reaching the lowest value (r = 0.79) at 37°S. Model T-S diagrams were also compared with those from CARS in the first 200 km and from the surface to 1500 m depth. The core of the OMZ in the model coincides with a salinity maximum centered around the 26.5 kg m -3 isopycnal surface and is consistent with CARS based T-S diagrams ( Figure A4). Despite the model underestimating the salinity maximum related to ESSW and overestimating the salinity minimum related to AAIW (i.e ESSW is less salty and AAIW is saltier in the model) it reproduces reasonably well the subtropical water masses structure around the OMZ, namely: Subantarctic Water, Equatorial Subsurface Water and Antarctic Intermediate Water ( Figure A4).

Seasonal cycle of DO at 36°30'S
We used cruise data at 36°30'S to evaluate the seasonal cycle of the model DO in the southern tip of the OMZ near the coast ( Figure A5). Since September 2002 a monthly, ship-based time series has been maintained over the continental shelf (~90 m depth) at 36°30'S by the University of Concepcion (e.g. Escribano et al., 2012;Escribano & Schneider, 2007). The data were depth-averaged over the whole water column (between 0 and 90 m depth) and a seasonal cycle was estimated over the period 2003-2015 for the observations and 2000-2008 for the model data. The comparison indicates that the model has skill in simulating the seasonal variability of DO, both amplitude and phase. Negative DO anomalies are observed during the upwelling season in austral spring and summer and positive ones during fall and winter. While both seasonal cycles were well correlated (r = 0.85), the amplitude of the simulated annual harmonic was somewhat smaller (32.5μM) than that from the observations (44.5μM).

Comparison of model and observed in situ DO
Cruise data offers the opportunity to compare some mesoscale features that are relevant for our study and that are not present in a climatology like CARS. Firstly, we compare the simulated vertical structure of DO and alongshore velocities with in-situ data (and geostrophic velocity) in a zonal transect at 35.5°S carried out during November 21, 2004 ( Figure A6). During that month the model and the observations showed a mesoscale eddy at that latitude (our simulation did not assimilate observations, so it is expected that mesoscale eddies do not coincide in space and time). The model reproduced fairly well both the magnitude and vertical structure of DO related to mesoscale eddies, and the alongshore velocities also were consistent with observations, although concentrations of DO were overestimated by~20-30 μM at the core of the OMZ.
Below 500 m depth DO was underestimated by~50 μM. As shown above (A.3) the simulated AAIW has relatively high oxygen and low salinity below~500 m depth. The comparison of the T-S diagrams near the coast (first 200 km) corroborates that observed in Section A.2 using CARS data and the whole period (compare Figures A6 and A4).
We use also cruise data to evaluate the vertical oxygen structure near the continental slope, where the oxycline is more intense. We considered DO profiles between 0 and 600m depth from two different sites (35.5°S; 73.2°W and 37.5°S; 74.0°W). The observations (and model profiles) were obtained in November-2004, December-2005, October-2006and March-2008. The Figure A7 shows a comparison between the different profiles at both sites. As mentioned before, the model reproduces the OMZ, but it shows, in general, larger values of DO than the observations, while below the OMZ the model underestimates DO by~16 μM at 35.5°S Figure A4. T-S-DO diagrams for ROMS/BioEBUS model (a) and CARS climatology (b). Only water encompassing the region within the first 200 km from the coast and from the surface to 1500 m depth in whole domain were considered. DO (in μM) is represented by color. The Surface Subantarctic Water (SAAW, 11.5°C; 33.8), Equatorial Subsurface Water (ESSW, 12.5°C; 34.9) and Antarctic Intermediate Water (AAIW, 3°C; 34) are shown. The corresponding values of temperature and salinity to each water mass is indicated between parenthesis. These values were obtained from Silva et al. (2009). and 19 μM at 37.5°S. The upper oxycline is relatively well reproduced, except during 2005 at 37.5°, when a very sharp and unusually deep upper oxycline was observed. In contrast, the lower oxycline was, in general, smoothed resulting in a large bias between modeled and observed DO.

The PCUC near 30°S
To assess the ability of the model in reproducing the vertical structure and seasonal variability of the PCUC, we used data from a current-meter mooring located over the slope at~1000 m total depth near   Figure A1). The data used here expand only from 2003 to 2006 to include a period where an ADCP was present in the mooring. The ADCP (an RDI workhorse 300 kHz) was located at 120 m depth (upward-looking) with a vertical resolution of 4 m. In addition, the mooring had 4 punctual Aanderaa RCM 8 currentmeters located at 220 m, 330 m, 480 m and 750 m depth. As we are interested here in the seasonal cycle, the original hourly data were monthly averaged and then used to estimate a mean vertical profile and typical annual variability ( Figure A8). The seasonal cycle of the current at 220m (near the PCUC core) was estimated by averaging monthly values from 2003 to 2006. In general, the PCUC shows very large intraseasonal variability at this site [cf. Pizarro et al., 2002, their Figure 3] and consistent semi-annual and annual fluctuations. The mean profile of the model PCUC in the same location is somewhat deeper and weaker than the observed one ( Figure  A8a), resulting in a small vertical shear near the upper 100 m. The variability, based on monthly means, is also smaller in the model PCUC (shaded areas in Figure A8b). Here, we compare also the annual and semi-annual harmonics of the current at 220m ( Figure A8b). The amplitudes of both harmonics are well reproduced by the model, but the harmonics are lagged by 1-month. Despite, the fact that the model slightly underestimates a bit the magnitude of the mean PCUC in this location, the total transport at this latitude is rather similar to the transport estimated using in situ currents and hydrography (see main text).