Dynamics of the Carbonate System Across the Peruvian Oxygen Minimum Zone

The oxygen minimum zone (OMZ) of Peru is recognized as a source of CO 2 to the atmosphere due to upwelling that brings water with high concentrations of dissolved inorganic carbon (DIC) to the surface. However, the inﬂuence of OMZ dynamics on the carbonate system remains poorly understood given a lack of direct observations. This study examines the inﬂuence of a coastal Eastern South Paciﬁc OMZ on carbonate system dynamics based on a multidisciplinary cruise that took place in 2014. During the cruise, onboard DIC and pH measurements were used to estimate pCO 2 and to calculate the calcium carbonate saturation state ( (cid:127) aragonite and calcite). South of Chimbote (9 ◦ S), water stratiﬁcation decreased and both the oxycline and carbocline moved from 150 m depth to 20–50 m below the surface. The aragonite saturation depth was observed to be close to 50 m. However, values < 1.2 were detected close to 20 m along with low pH (minimum of 7.5), high pCO 2 (maximum 1,250 µ atm), and high DIC concentrations (maximum 2,300 µ mol kg − 1 ). These chemical characteristics are shown to be associated with Equatorial Subsurface Water (ESSW). Large spatial variability in surface values was also found. Part of this variability can be attributed to the inﬂuence of mesoscale eddies, which can modify the distribution of biogeochemical variables, such as the aragonite saturation horizon, in response to shallower (cyclonic eddies) or deeper (anticyclonic eddies) thermoclines. The analysis of a 21-year (1993–2014) data set of mean sea surface level anomalies (SSHa) derived from altimetry data indicated that a large variance associated with interannual timescales was present near the coast. However, 2014 was characterized by weak Kelvin activity, and physical forcing was more associated with eddy activity. Mesoscale activity modulates the position of the upper boundary of ESSW, which is associated with high DIC and inﬂuences the carbocline and aragonite saturation depths. Weighing the relative importance of each individual signal results in a better understanding of the biogeochemical processes present in the area.

The oxygen minimum zone (OMZ) of Peru is recognized as a source of CO 2 to the atmosphere due to upwelling that brings water with high concentrations of dissolved inorganic carbon (DIC) to the surface. However, the influence of OMZ dynamics on the carbonate system remains poorly understood given a lack of direct observations. This study examines the influence of a coastal Eastern South Pacific OMZ on carbonate system dynamics based on a multidisciplinary cruise that took place in 2014. During the cruise, onboard DIC and pH measurements were used to estimate pCO 2 and to calculate the calcium carbonate saturation state ( aragonite and calcite). South of Chimbote (9 • S), water stratification decreased and both the oxycline and carbocline moved from 150 m depth to 20-50 m below the surface. The aragonite saturation depth was observed to be close to 50 m. However, values <1.2 were detected close to 20 m along with low pH (minimum of 7.5), high pCO 2 (maximum 1,250 µatm), and high DIC concentrations (maximum 2,300 µmol kg −1 ). These chemical characteristics are shown to be associated with Equatorial Subsurface Water (ESSW). Large spatial variability in surface values was also found. Part of this variability can be attributed to the influence of mesoscale eddies, which can modify the distribution of biogeochemical variables, such as the aragonite saturation horizon, in response to shallower (cyclonic eddies) or deeper (anticyclonic eddies) thermoclines. The analysis of a 21-year (1993-2014) data set of mean sea surface level anomalies (SSHa) derived from altimetry data indicated that a large variance associated with interannual timescales was present near the coast. However, 2014 was characterized by weak Kelvin activity, and physical forcing was more associated with eddy activity. Mesoscale activity modulates the position of the upper boundary of ESSW, which is associated with high DIC and influences the carbocline and aragonite saturation depths. Weighing the relative importance of each individual signal results in a better understanding of the biogeochemical processes present in the area.

INTRODUCTION
The upwelling region of Peru is one of the most important upwelling systems in the world due to its high productivity that supports abundant fisheries, particularly the Peruvian anchovy fishery Espinoza-Morriberon et al., 2017). Upwelling events mainly occur in spring and summer (Sobarzo and Djurfeldt, 2004;Franco et al., 2018). However, large and productive areas with important fisheries like Peru are considered to be zones that currently experience or that are projected to experience ocean acidification (Shen et al., 2017). The intense biological activity present in these areas produces a large quantity of organic matter (OM), some of which sinks and is degraded by catabolic processes (Bretagnon et al., 2018). Therefore, subsurface OM degradation contributes to the consumption of oxygen (O 2 ) and in conjunction with poor water mass ventilation, leads to the formation of an oxygen minimum zone (OMZ; Paulmier et al., 2008). The Peruvian OMZ presents a poorly oxygenated core (O 2 down to <1 µmol kg −1 : Revsbech et al., 2009), with significant anaerobic activity (e.g., Hamersley et al., 2007) that occurs in the thickest core (340 m in average; Paulmier et al., 2008).
Various studies have shown the existence of intense biogeochemical anomalies (e.g., in O 2 , CO 2 , pH, alkalinity, nitrate, nitrite, and N 2 O) associated with microorganism (e.g., bacteria) as well as mesoorganism and microorganism (e.g., zooplankton) processes that take place near the oxycline. The seasonal cycle of O 2 /oxycline depth is very weak when compared with that of interannual timescales. The OMZ off Peru is characterized by notable interannual variability that is particularly evident with regard to depth, an example of which is the 100 m deepening of the oxycline that is associated with trapped Kelvin waves and which is particularly apparent during El Niño events (Gutiérrez et al., 2011). Furthermore, primary production is 2-to 3-fold higher in summer than in winter (Graco et al., 2007). The Peruvian upwelling region is also considered to be a source of CO 2 to the atmosphere, with an estimated annual flux of 5.7 mol m −2 y −1 (Friederich et al., 2008).
In the Peruvian region, the OMZ is associated with high dissolved inorganic carbon (DIC) concentrations (mean of 2,330 ± 60 µmol kg −1 ) and is defined as a carbon maximum zone (CMZ). High DIC concentrations (2,225-2,350 µmol kg −1 ) have been reported over the entire OMZ, indicating that all studied OMZs may be classified as CMZs (Paulmier et al., 2011). In addition, the existence of the CMZ suggests that it is being formed by the same dynamics (low ventilation and upwelling) and biogeochemical (remineralization that produces DIC and consumes O 2 ) mechanisms (Paulmier et al., 2011) as an OMZ. Surface pCO 2 values in coastal regions can reach 1,500 µatm in subsurface waters, as can be observed in the California Current Takahashi et al., 2009). However, in the Peruvian upwelling system, pCO 2 values that exceed those of other upwelling systems can be found in surface water (Friederich et al., 2008;Shen et al., 2017). When the ocean absorbs more CO 2 , extremely high pCO 2 values (e.g., >1,000 µatm) imply a decrease in pH. "Hotspots" of acidification, like the Peruvian regions, are thus predicted to occur in major fishery zones by mid-century when atmospheric CO 2 is projected to reach 650 µatm (McNeil and Sasse, 2016;Shen et al., 2017). In addition, ocean acidification modeling studies based on the oceanic uptake of anthropogenic CO 2 indicate that the upper water of the Humboldt Current is likely to become more corrosive with regard to mineral CaCO 3 (Franco et al., 2018).
Biomineral saturation, such as omega calcite ( calc ) and omega aragonite ( arag ), is a function of the concentration of CO 2− 3 , Ca 2+ , and temperature and may be expressed with the pressure-dependent stoichiometric solubility product K * sp ( = [Ca 2+ ] [CO 2− 3 ]/K * sp ; Mucci, 1983). The aragonite and calcite saturation states respond directly to changes in the availability of the CO 2− 3 ion. If the ocean absorbs more CO 2 , pH decreases and arag and calc also decrease. For example, when = 1, seawater is saturated with respect to CO 2− 3 . Conditions favor precipitation or the preservation of carbonate minerals when > 1, while conditions favor dissolution when < 1. If the aragonite and calcite saturation states decrease, greater physiological challenges for calcifying organisms are expected to be present Guinotte and Fabry, 2008). Often, large fisheries are present in regions in which ocean acidification has already occurred or is expected to occur (Shen et al., 2017). For example, the Peruvian anchovy (Engraulis ringens) is found in a region where the pCO 2 values can be higher than 1,000 ppm (Friederich et al., 2008). Furthermore, eggs from Peruvian anchovies that spawned in waters with high pCO 2 were found to have survival rate (Shen et al., 2017).
Coastal zones are connected to ocean areas in which OMZs are expanding (Stramma et al., 2008). Therefore, it is important to generate a better understanding of the different components of the carbon system that are involved in the formation and maintenance of OMZs. Very little is known about ocean acidification in Peruvian regions. Moreover, the mechanisms behind the existence of such an intense CO 2 maximum and its relationship with physical and biochemical processes must be elucidated. We evaluated the relevant physical processes that influence the spatial variability of carbonate chemistry and the aragonite saturation state of the Peruvian coastal region. We describe the patterns of relevant water masses with low pH and oxygen concentrations. In particular, we describe shoaling due to mesoscale activity that modulates the position of the upper boundary of the water masses. This shoaling is associated with high DIC and influences the carbocline and aragonite saturation depths. Our oceanographic and biogeochemical data show the effects of physical processes on the spatial and vertical distribution of the CO 2 variables. Due to the shallow OMZ and CMZ in the Peruvian upwelling system, this region can be considered a natural laboratory for investigations that are focused on the CO 2 system.

METHODS
Oceanographic sampling was carried out aboard the RV L'Atalante (Figure 1). The AMOP cruise was planned as a round trip from Callao that departed on January 26 and arrived on February 22, 2014. The first sampling stations were located near the coast. By following the 12 • S transect, the cruise continued northward and sampled the water column in typical offshore stations (0-2,000 m depth). After which, the vessel returned to sample coastal stations following the main gradient of the bathymetry. The ship then headed southward to sample stations along the shelf (typical maximum depths of ∼150 m) until the vessel reached the Pisco sector (15 • S). The last offshore transect was then sampled up to station 28 (14 • 34 ′ S−77 • 16 ′ W), which represents the southernmost point of the cruise (Maes et al., 2014). Discrete samples for DIC and pH analysis were collected over depth gradients at coastal and offshore sites in the Peruvian OMZ area. A total of 31 station profiles of ∼500 meters depth were generated. At each station, vertical profiles were produced using standard oceanographic equipment (CTD-rosette at most stations). Seawater samples for DIC and pH samples were collected by a rosette casts covering the low-oxygen water with an average of 12 depths targeting the key water column features of the oxic photic zone, oxycline (upper), OMZ interface, and anoxic OMZ core.
Water samples for DIC analysis were collected in 250-ml sodium borosilicate bottles and preserved with 50 µl of HgCl 2 saturated solution. For pH analysis, the samples were collected directly from the Niskin bottles using 60-ml syringes. All samples were analyzed aboard the ship. For DIC analysis, a LI-7000 gas analyzer (CO 2 /H 2 O, LICOR, Lincoln, NE, USA) was used. The certificate reference material for DIC analysis was provided by the laboratory of Dr. Andrew Dickson of Scripps Institution of Oceanography (Dickson et al., 2003). The relative difference averaged ±2 µmol kg −1 (0.2% error).
The pH was measured using a glass electrode at 25 • C on the seawater scale (pHsw) as described by Chapa-Balcorta et al. (2015). We use the program CO2SYS (Lewis and Wallace, 1998) and DIC-pH sw to calculate pCO 2 , in situ pHsw, and the aragonite saturation state ( Arg ). We use the constants proposed by Mehrbach et al. (1973) for these calculations. The uncertainty obtained for in situ pH sw , Arg , and pCO 2 was ± 0.04 ± 0.2, and 56 µatm, respectively. The pCO 2 variability range observed in Peruvian waters was 25-fold greater than the error.

Wind Data
Wind speed data were obtained from the SeaWinds database provided by NOAA (Zhang et al., 2006a,b). Wind stress was computed using the bulk equation (Schwing et al., 1996). Monthly means of the wind speed field were also derived from the Copernicus level 4 multi-sensor blended wind product.

Vertical Velocities Due to Wind Stress Curl and Ekman Transport
The curl-driven vertical velocities (w c ) and coastal upwelling associated with Ekman offshore transport (w e ) at the limit of the mixed layer were estimated following Rykaczewski et al. (2015) using Equations (1) and (2): where ρ w = 1,024 kg m 3 is the seawater density, f is the coriolis parameter, τ is the wind stress vector field, τ a is the wind stress parallel to the coastline, and Rd is the local Rossby radius of deformation. Note that using Rd to convert Ekman transport in upwelling rate yields an estimate that is in the low range values of the actual amplitude of vertical velocity since the upwelling scale is smaller than Rd [cf. Marchesiello and Estrade (2010)].

Variability of Sea Surface Height Anomalies
In order to calculate the relative importance of the different scales of variability for sea surface height anomalies (SSHa) in the study area, the seasonal, interannual, and mesoscale signals were derived following the methods of Godínez et al. (2010) and León-Chávez et al. (2015). A 21-year database (1993-2014) of mean SSHa distributed by AVISO was used (http://www.aviso. altimetry.fr). The seasonal, mesoscale, and interannual signals were separated. The first was obtained using a harmonic fit, while the second and third were obtained by means of orthogonal empirical functions (Equations 3 and 4). A a and A s are the annual and semi-annual amplitudes of the harmonic fit, respectively, while x is the spatial vector, ϕ a and ϕ s are the annual and semi-annual phases, ω is the annual radian frequency, t is the time referred to the beginning of the year, and n is the number of EOF. The first two terms on the right of Equation (3) are the seasonal components (SSHa season is the first term in Equation 4), which were calculated using a least squares fitting. Fitting errors and uncertainties were calculated as in previous studies (Prentice et al., 2001;Ripa, 2002). F 1 ( x) and f 1 (t) are the spatial and temporal series of the first EOF mode, which in the South Equatorial Tropical Pacific (SETP) typically contain the interannual signal (SSHa inter in Equation 4) that is characterized by ENSO-induced variability (Macharé and Ortlieb, 1993). In the last term of Equation (3), F n ( x) and f n (t) are the spatial and temporal series of the nth EOF mode. The sum of these (2 to N) mainly contains mesoscale variability (SSHa meso , third term of Equation 4), although it may include some interannual variability not captured by the first EOF. Analyzing this information enables us to identify the dominant processes in the region in general as well as during the time of our cruise, allowing us to weigh the significance of the measurements according to the predominant signal type. The local variance (LV) of each variability component (seasonal, interannual, and mesoscale) was calculated as follows: where LV season indicates local variance, σ 2 SSHaseason is the variance of the seasonal component of SSHa, and σ 2 SSHa is the variance of the original SSHa. LV was calculated per pixel for each time series. The same procedure was used for the interannual and mesoscale components. The complete series was used for this calculation . The signal decomposition during the cruise was obtained by reconstructing SSHa data of the separated variability scales for the dates corresponding to the sampling periods. The seasonal data were reconstructed from the harmonic analysis. The interannual data were reconstructed from the first EOF. The mesoscale data were constructed from the sum of the EOF from 2 to N.

Inferred Sea Surface Height Anomalies
With a deep level of no motion in the offshore transect (1,000 db), we used T and S to calculate density and thus inferred sea surface height anomalies, which we compared to the AVISO satellite product. The dynamic height calculation was performed using the functions of the TEOS-10 library (McDougall and Barker, 2011). The measured salinity and temperature data were interpolated vertically with an interval of 1 m using the objective interpolation of Barnes (1964).

RESULTS
We provide a description of the circulation patterns that drive the physical and biogeochemical patterns observed in the study region during the cruise, and we use this as a frame of reference to define the importance of each water mass in the region. We also provide a description of the oceanographic conditions during the cruise, which is interpreted considering the analysis of altimeter data to estimate anomalous conditions in seasonal and interannual timescales.

Wind, Wind Stress Curl, and Ekman Transport Conditions
At the beginning of the sampling period in January 2014, the wind speed was <5 m s −1 throughout the coastal area spanning from Pisco to Malabrigo. The wind speed was stronger offshore and presented values of 6-7 m s −1 (Figures 2A,B). The most intense winds 10 m s −1 were observed offshore in February, although winds of 7 m s −1 were observed in the coastal region near Pisco. During the January sampling period, the winds oscillated from ∼3 to 5 m s −1 and produced relatively weak upwelling compared to the remainder of the period ( Figure 2C). Stronger yet opposite wind stress can be observed in February (Figure 2D). The combination of calm (Figure 2A) and subsequently strong winds ( Figure 2B) produced different oceanographic conditions; namely, regions with high vertical mixing and areas with coastal upwelling (Figure 3). The mean vertical velocities associated with Ekman transport and the monthly vertical transport means due to winds stress curl for January and February are presented in Figure 3. Of note, February presented the highest Ekman offshore transport and Frontiers in Marine Science | www.frontiersin.org 6 October 2019 | Volume 6 | Article 617 FIGURE 6 | Vertical distribution of the potential density anomaly (σ θ , kg/m 3 ), dissolved oxygen (DO, µmol kg −1 ), and omega aragonite ( arag ) along the transect.
vertical transport values. However, the vertical transport near to Chimbote did not occur because of the presence of an anticyclone eddy (Figure 4).

Variability in Sea Surface Height Anomalies
The SSHa in the study area for January and February 2014 are shown in Figure 4. The cruise track variability that was attributed to mesoscale eddies during February was present during the anticyclonic eddy from Malabrigo to the middle of Huacho for both odd shore and coastal stations. In front of Huacho and to the north of Callao, cyclonic eddy influence may be observed in both January and February (Figures 4A,B). While the effect of each of the conditions on the chemistry of the CO 2 system will be described later, the spatial distribution of the CO 2 variables was directly associated with the various oceanographic conditions.

Water Masses, Water Types, and Vertical Structure
The distribution of water masses was examined to analyze the potential effects on the carbonate chemistry in the study area. We also found it very useful to study the water-type distribution because of the broad T-S variability. The T-S plot in Figure 5 incorporates several characteristics for comparison on the Z axis, such as latitudinal distribution, oxygen concentration, pH SW , pCO 2 , and omega aragonite ( Arg ). The T-S diagrams for all sections in the study area indicate the presence of four water masses in the first 500 m based on the classifications of Silva and Neshyba (1976), Silva et al. (2009). In addition, three water types were observed as defined by the conditions of Swartzman et al. In surface waters, MESC and SSW were detected in the first 25 m ( Figure 5A). MESC is less salty and it was located north of 8 • N, while SSW had salinities >35 and was observed south of 12 • N ( Figure 5B). The corresponding MESC values following 24 kg m −3 isopycnal (located in <10 m) for oxygen, pH, pCO 2 , and arag were 250 µmol Kg −1 , ∼8.1, <500 µatm, and >3.5 units, respectively. The values for SSW between 10 and 60 m for the 25 kg m −3 isopycnal were >200 µmol Kg −1 , >7.8, ∼550 µatm, and >2.0 units (Figure 5). CCW located below MESC was also found at 8 • N between 20 and 60 m depth (Figures 5A,B). When following the 25 kg m −3 isopycnal in <20 m depth, the average values of CCW characteristics indicated that this water type presented less oxygen, lower pH, higher pCO 2 , and lower arag (close to zero to 50 µmol Kg −1 , ∼7.7, <1,000 µatm, and <2.0, respectively) when compared to that of MESC. The water type CAW was observed between 55 and 80 m. This water mass presented the lowest oxygen concentrations, low pH, higher pCO 2 , and low arag (0 and 20 µmol kg −1 , ∼7.7, >1,300 µatm, and <2.0, respectively). In subsurface waters, the main water mass of the OMZ and CMZ in the study area considering the observed characteristics (low  oxygen and pH and high pCO 2 and DIC) was ESSW. The upper ESSW limit presented a potential density anomaly of 26 kg m −3 (Figure 5) in the offshore transect, and this isopycnal was found at ∼70 m depth (Figure 6). However, the isopycnal was found at ∼30 m near the coast at 14 • S. The oxygen concentration declined to <20 µmol kg −1 at this depth but also presented pH values of 7.7 and <1.5 units for arag (Figure 6). Figure 7 presents the spatial distribution at 30 m for temperature, density, oxygen, pCO 2 , pH, and omega aragonite. With the exception of a few stations located in the offshore region between 8 and 11 • S, the area was influenced by subsurface waters and was therefore characterized by low concentrations of oxygen (<10 µmol kg −1 ), subsaturation values of omega aragonite (∼1 near the coast Callao-Pisco), low pH (<7.5), and high pCO 2 (values from ∼ 500 to 1,200 µatm). The highest values of pCO 2 were observed in the coastal region between Callao and Pisco at the southern end of the transect.

Vertical Mixing
Two different oceanographic scenarios were found during the cruise and consisted of the presence or absence of upwelling. The wind stress near the coast was low during January (<0.03 N m −2 ) prior to sampling (Figures 2, 3) and presented wind speed values <3 ms −1 . However, in February, the monthly wind speed mean was 7 ms −1 while the wind stress was ∼0.07 N m −2 . In February we found the highest values for both Ekman offshore transport and vertical transport (Figure 3). During February, the upper ESSW limit was found near the surface between 20 and 50 m depth in the coastal area between 12 and 14 • S. The lowest surface temperatures (∼17 • C), lower oxygen concentrations (<100 µmol kg −1 ), the highest surface DIC concentrations (∼2,300 µmol kg −1 ), and low omega aragonite values (<1.5) were present at ∼20 m depth (Figure 8).
During the cruise, weakly stratified regions were identified based on temperature profiles with a N-S gradient (Figure 8). The stratification parameter was calculated according to Simpson (1981) for the upper 300 m of the water column. This value expressed in J m −3 and represents the amount of work per volume that is needed to mix the water column to a specified depth. The lowest stratification values were detected (<200 J m −3 ) between 11 and 15 • S near the upwelling area found between Callao and Pisco (see vertical profiles in Figure 8). However, higher stratification values (>400 J m −3 were found between 8 and 11 • S. For all sampled areas, the calculated stratification average was 350 J m −3 . This combination of upwelling and stratification plays an important role in the spatial and vertical distribution of DIC, omega aragonite, oxygen, and temperature.

Circulation Based on SSHa Variability Scales
The decomposition of SSHa variability for the Peruvian area based on 21 years of AVISO data shows that the interannual signal is dominant (∼30-80%), followed by the mesoscale signal (20-60%), and finally the seasonal signal (∼7 and 30%; Figure 9). However, the cruise took place during a normal year. The seasonal cycle presented weak variance, which is consistent with previous studies in the eastern tropical Pacific .
Seasonality between 6 and 16 • S had a greater effect on SSHa between 80 and 84 • W in offshore areas while a small influence was observed near the coast between 8 and 13 • S, which coincides with the location of the sampling network ( Figure 9A). The interannual variability was mainly the result of ENSO as shown in Figure 9B where the temporal component Frontiers in Marine Science | www.frontiersin.org 9 October 2019 | Volume 6 | Article 617 of the first mode (EOF1) can be seen to be correlated with the El Niño multivariate (MEI) and ONI indices. It should be noted that the region is partially located within the area referred to as El Niño region I (Trenberth and, 2019). Figure 9B also shows that the ENSO effect extended south of the study area near the coast with a local variance between 57-77%. Mesoscale eddy activity is generated by variations in density that arise from coastal currents in their own vortices. The greatest effect of this activity was detected with SSHa in the southern region between 13 and 16 • S, although mesoscale eddy activity presented latitudinal variation and minimum values were present in the northern region ( Figure 9C). This is likely because most eddies are generated near the coastal region between 10-16 • S and reach their greatest diameters between 13-16 • S.
During neutral El Niño conditions, as in February 2014, the ENSO signal may be overshadowed by high mesoscale variability. The signal decomposition for the sampling dates during the cruise is shown in Figure 10. January and February were neutral ENSO periods (ONI = −0.4), and the interannual signal was weak (0.5-1.8 cm). The strongest driver of SSHa variability was the mesoscale signal (−5.9 to 5.1 cm), which was mainly due to the presence of eddies, with fluctuations of −1.3 to 3.0 cm generated by seasonality.

Inferred Sea Surface Height Anomalies, Satellite Product, and Cruise Track Variability Attributed to Mesoscale Eddies
To evaluate if mesoscale eddy activity played a large role in carbon system variability within the inner and outer transects, two approaches were considered. Inferred SSHa by calculate density and thus to compare with the satellite product ( Figures 11A,C,D). We found that both estimations showed similar trends. However, a linear regression analysis was carried out to determine if the inferred SSHa were associated with the satellite dynamic height product. The R² value of this analysis was 0.78. Likewise, a statistically significant Pearson linear correlation coefficient of r p = 0.88 was obtained, indicating that increases in sea surface height were associated with SSHa (confidence level of 95%).
After SSHa were estimated, they were used to geo-reference each carbonate system profile to determine the degree of variability in the carbonate system that may be attributed to mesoscale eddies (Figures 11B,E). For example, the large horizonal and vertical variability associated with DIC values of 2,250 µmol kg −1 was due to eddy presence. The linear regression was performed considering SSHa and DIC depth (Z). An association was found between Z DIC= 2250 (depth at which DIC = 2,250 µmol kg −1 ) and SSHa for the deepest portions of the main offshore transect. The model obtained has x 0 = −29.11, x 1 = 8.08, R 2 = 0.46, and r p = −0.68. We found that mesoscale processes change the upper depth limit of ESSW and the depth of this water mass, while affecting DIC, the aragonite saturation depth, and the upper limit of the OMZ.

DISCUSSION
This study describes the Peruvian upwelling conditions observed during January and February 2014 and the associated spatial variation of the CO 2 system. During upwelling events, the aragonite saturation depth was observed to be close to 50 m, but values of <1.2 were detected close to 20 m and were accompanied by low pH (<7.5), high pCO 2 (>1,250 µatm), and high DIC (>2,250 µmol kg −1 ). These chemical characteristics were linked with the portion of the water-column that spanned from the upper limit of Equatorial Subsurface Water (ESSW) to the upper limit of the OMZ (<5 µmol kg −1 ) and the carbocline (∼2,300 µmol kg −1 ). Our analysis of the decomposition of SSHa variability also shows that the coastal-oceanic signal was dominated by interannual variability. However, during our sampling period, mesoscale processes were key drivers of surface water dynamics. In particular, the upper depth limit of ESSW may be altered by physical phenomena, such as eddies or interannual events. This results in changes in the depth of this water mass and affects DIC, the aragonite saturation depth, and the upper limit of the OMZ.

Vertical Mixing, Water Masses, and Water Types
Large spatial variability in sea surface anomalies were found (Figure 4). Part of this variability may be attributed to the influence of mesoscale eddies, which modify biogeochemistry. For example, shallower saturation horizons and high DIC values may be due to cyclonic eddy activity, while low DIC concentrations will be deeper due to anticyclonic eddy activity. The sampling area showed eddy production in the oceanic region, which may account for a large percentage of the variability that was observed in the offshore stations located south of 13 • S (>40%) as may be seen in Figure 9C. Additional information regarding the effect of eddies on the biogeochemistry inside an OMZ has been reported by Bettencourt et al. (2015) who indicate that mesoscale structures have relevant dual roles at depths between 380 and 600 m. Firstly, their mean positions and paths delimit and maintain the boundaries of oxygen minimum zones. Secondly, their high frequency fluctuations entrain oxygen across these boundaries. As a result, these eddy fluxes contribute to the ventilation of the oxygen minimum zone. Eddy structures can introduce water types with different pH and arag values (Figure 5). For example, MESC and SSW at the surface were affected by anticyclonic and cyclonic eddies (Figure 6), respectively, and MESC had higher pH and oxygen values than SSW (Figure 5).
The offshore surface waters were saturated with regard to aragonite ( arag >1.5 to >4.5) in the first 50 m of the water column (Figures 6, 7). However, for depths below 150 m, the arag values were <1. These values are similar to those reported for the open ocean (300 km from the coast) in the Eastern South Pacific (Feely et al., 2012;Jiang et al., 2015). The upper limit of ESSW can be identified by the 26 kg m −3 isopycnal (Figure 5). The depth of this isopycnal surface was closely related to the thermocline and with the aragonite saturation horizon (Figure 6). Therefore, most of the chemical features associated with this isopycnal were similar regardless of depth. In most cases, the sinking of this density surface directly resulted in changes in the depth of the aragonite saturation horizon. Near Pisco, the upper limit of ESSW was as shallow as ∼20 m depth.
During the cruise, the low surface arag values (<1.2) that were observed near the coastal stations were due to the influence of upwelled water. However, a cyclonic eddy was present during sampling along the coast between 11 and 12 • S (Figure 4; February). This eddy resulted in an upward vertical transport of aragonite-undersaturated ESSW, which reduced the surface values of arag (Figures 6, 7). Furthermore, an anticyclone was observed in most of the oceanic region (8-12 • S and 80-81 • W) that resulted in a depression of the aragonite saturation horizon. The highest arag values (near 3) that were found in the Peruvian area were observed offshore in the northern transects. Nevertheless, the high offshore pH values (>8.1) and low pCO 2 values (<600 µatm) that were found in this region suggest that the aragonite saturation state is driven by two processes: outgassing of CO 2 and the biological consumption of dissolved inorganic carbon. However, submesoscale processes are another physical driver of CO 2 outgassing. Köhn et al. (2018) carried out a CO 2 variability study across an upwelling front near Peru using high-resolution underway measurements and provided evidence of the complex submesoscale distribution of surface CO 2 in the Peruvian upwelling system. Thus, physical processes may also influence the distribution of pCO 2 on different spatial and time scales.
Based on the analysis of 1993-2014 SSHa data, we found that interannual variability plays the most important role (up to 60% of local SSHa variance) in determining the circulation of the coastal areas and in the southern region of our study area, followed by mesoscale variability. This result is different from what was observed by Godínez et al. (2010) for the OMZ located in the Eastern Tropical Pacific (ETP) off the coast of Mexico (15-24 • N). In their study, seasonality was found to be the dominant signal, which was followed by the mesoscale signal. The mesoscale contribution to SSHa variability was four times larger than what was observed for the ETP (10-20% of the local SSHa variance). In contrast with our study area, interannual variability was three times larger than mesoscale variability near the coast. The waters off Peru are intimately associated with the dynamics of the equatorial Pacific and are thus subject to large temporal variations. The ENSO cycle, with its warm El Niño phase and cold La Niña phase, is comprised of a natural interannual oscillation frequency of the ocean-atmosphere system in the tropical Pacific . Over several months, this cycle substantially alters the functioning of ecosystems associated with the coastal outcrop (Grados et al., 2018).
In our study, mesoscale variability was less important than interannual variability, but we also found that mesoscale eddy activity notably influenced carbon system variability within the inner and outer transects (Figure 11). We found that the large spatial variability observed in our study area was attributed to the influence of mesoscale eddies, which can modify the distribution of biochemical variables, such as the aragonite saturation horizon, in response to shallower (cyclonic eddies) or deeper (anticyclonic eddies) thermoclines. According to Willett et al. (2006), eddies are the main contributors of mesoscale variability in the ETP. Eddies affect the spatial distribution of dissolved chemicals by Ekman pumping and by exporting dissolved and particulate materials into the open ocean, fostering horizontal advection (Chen et al., 2007;Samuelsen and O'Brien, 2008;Gaube et al., 2014;Nagai et al., 2015). In the study area, mesoscale processes influenced the upper depth limit of ESSW, changing the depth of this water mass and affecting DIC, aragonite saturation, and the OMZ depth. Mesoscale eddies, whether they are driven by wind forcing or otherwise, can persist for weeks (cyclonic) or months (anticyclonic; Gonzalez-Silvera et al., 2004;Palacios and Bograd, 2005). Although mesoscale processes were not as important as interannual processes, they were still very important in the region during the cruise and affected the spatio-temporal variability of the biogeochemical patterns in the region and are thus key factors to consider when studying the carbonate system in the ocean offshore of Peru.

Variability of Dissolved Inorganic Carbon
The DIC concentrations detected in the study area were within the data ranges that have been reported for SSW in the ETP (Paulmier et al., 2011;Franco et al., 2014;Chapa-Balcorta et al., 2015). The DIC variability observed during sampling indicates that upwelling and mesoscale variability affect DIC values. The corresponding changes in pCO 2 and arag affected both regional biogeochemical conditions and habitat quality for calcifying or pH-sensitive species.
The impact of climate change on the dynamics of the Eastern Boundary upwelling system (EBUS) has been debated and a consensus has not been reached so far, mostly because of competing effects due to the increase in upwelling favorable winds induced by the expansion of the Hadley cell (Bakun et al., 2015;Rykaczewski et al., 2015) and the increase in surface warming that increases vertical stratification and tends to reduce upwelling (e.g., Echevin et al., 2012 for the coast of Peru). There is also uncertainty with regard to whether upwelling favorable winds will actually increase in the mid-latitudes  given that increasing trends in alongshore favorable winds are barely detected in existing data sets (e.g., Figure 1 of Belmadani et al., 2014). It has been postulated that along the Peruvian coast, global warming has increased the heat difference between the continent and the adjacent sea. This change has resulted in an intensification of the winds, which has resulted in the coastal outcrop. However, recent studies suggest a notable cooling of the marine-coastal strip between central Peru and northern Chile over the last 35 years, although records of coastal winds are scarce and the causes behind this trend remain unclear (Yáñez et al., 2018). The first simulations of climate change effects in the waters off Peru predicted a decrease of the coastal outcrop, an intensification of coastal retention processes, and a reduction in marine productivity, in which thermal stratification plays an important role (Gutiérrez et al., 2011). These changes in upwelling intensification will also be reflected in the carbon biogeochemistry of coastal areas.

Aragonite Saturation Horizon Maintenance
A lower oxygen concentration at the aragonite saturation depth was found in the southern stations near the coastline (Figures 6, 7). The remineralization of organic matter produced at the surface may intensify the consumption of oxygen in the subsurface layer due to low ventilation driven by strong stratification in the area (Fiedler and Talley, 2006;Fiedler et al., 2013;Franco et al., 2014). The latter is consistent with what was reported by Fiedler and Talley (2006), who concluded that 90% of organic matter was remineralized in the first 200 m of the water column. In turn, Paulmier et al. (2006) showed that organic matter remineralization in the oxycline of the Humboldt system was 3-fold more intense than in the core of the OMZ, which favors the maintenance of these zones. In highly productive coastal zones, eutrophication due to runoff from terrestrial nutrient addition may have a greater effect on the decrease of aragonite saturation than on ocean acidification driven by anthropogenic CO 2 absorption (Borges and Gypens, 2010). Water with low oxygen and carbonate ion concentrations may have different biological responses, as this situation may lead to habitat increases or decreases for the species that reside in OMZs, particularly because these chemical environments affect their metabolic functions.
Like biological processes, physical drivers play a very important role in the vertical distribution of several CO 2 variables in the ocean. The results reported in this study clearly demonstrate that DIC concentrations, pCO 2 , and arag depend not only on the upwelling season but also on the presence, absence, and intensity of mesoscale eddies. Therefore, proper evaluations of physical drivers should always include an analysis of the indicators of these eddies, such as dynamic height, geostrophic currents, or altimetry in models that are able predict which areas are likely to experience ocean acidification.
The upwelling system in Peru is a natural laboratory for the study of regional interactive processes between the land, ocean, and atmosphere while providing a unique opportunity to understand carbon variability in different scales. For example, the Peruvian anchovy Engraulis ringens lives and spawns in a region where pCO 2 values may be higher than 1,000 ppm (Friederich et al., 2008), and low survival rates have been attributed to these high values (Shen et al., 2017). These authors also highlight the Frontiers in Marine Science | www.frontiersin.org 13 October 2019 | Volume 6 | Article 617 need to understand the effect of pCO 2 levels on spawning and mortality rates to predict the effects of ocean acidification on fisheries management. In response to future warming, changes in timing, intensity, and spatial heterogeneity are expected to occur in most EBUS (Wang et al., 2015).

CONCLUSIONS
During the upwelling events, the aragonite saturation depth was observed to be close 50 m, but values <1.2 were detected close to 20 m along with low pH (<7.5), high pCO 2 (>1,250 µatm), and high DIC (>2,250 µmol kg −1 ) values. These chemical characteristics were associated with Equatorial Subsurface Water (ESSW). In addition, the aragonite saturation depth became gradually shallower from north to south, beginning at 140 m and decreasing to <40 m with pCO 2 values close to 1,250 µatm. Large spatial variability in surface values was also observed. Part of this variability can be attributed to the influence of mesoscale eddies, which can modify the distribution of biogeochemical variables. The results reported here clearly demonstrate that DIC concentrations, pCO 2 , and arag depend not only on upwelling but also on the presence, absence, and intensity of mesoscale eddies. Mesoscale eddies changed the upper depth limit of ESSW, changing the depth of this water mass and affecting DIC, aragonite saturation, and OMZ depth.
Our analysis of the decomposition of SSHa variability that was based on a 21-year data set, indicated that the interannual signal was the dominant coastal-oceanic signal. However, during our sampling period, mesoscale processes were key drivers of surface water dynamics. The 2014 cruise was characterized by weak Kelvin activity, and physical forcing was more associated with eddy activity. Large spatial variability at the surface was also found. Part of this variability was attributed to the influence of mesoscale eddies which modify biochemistry. For example, the aragonite saturation horizon was affected by shallower or deeper thermoclines, which were the result of cyclonic or anticyclonic eddies, respectively. Near the coast and south of 12 • S, pH T values from 8 to 7.6 were observed in upwelling zones but were also distributed in several shallow water types. Values <1.2 were detected close to 20 m along with low pH (minimum of 7.5), high pCO 2 (maximum of 1,250 µatm), and high DIC concentrations (maximum of 2,300 µmol kg −1 ). Weighing the relative importance of each individual signal leads to a better understanding of the biogeochemical processes in the area.

AUTHOR CONTRIBUTIONS
The conceptualization and arrangement of the original manuscript was under the charge of JH-A, AP, and VG. The conceptual design and financial support of the oceanographic expeditions were under the direction of VG, AP, BD, and CM. All CTD data were facilitated through JS. Preprocessing of all data was carried out by CC-B, JS, IM, GD, and MB. Laboratory analyses were conducted by JH-A. All authors contributed equally to the original manuscript.