Temporal dynamics of surface ocean carbonate chemistry in response to natural and simulated upwelling events during the 2017 coastal El Niño near Callao, Peru

Oxygen minimum zones (OMZs) are characterized by enhanced carbon dioxide (CO₂) levels and low pH and are being further acidified by uptake of anthropogenic atmospheric CO₂. With ongoing intensification and expansion of OMZs due to global warming, carbonate chemistry conditions may become more variable and extreme, particularly in the eastern boundary upwelling systems. In austral summer (February–April) 2017, a large-scale mesocosm experiment was conducted in the coastal upwelling area off Callao (Peru) to investigate the impacts of ongoing ocean deoxygenation on biogeochemical processes, coinciding with a rare coastal El Niño event. Here we report on the temporal dynamics of carbonate chemistry in the mesocosms and surrounding Pacific waters over a continuous period of 50 d with high-temporal-resolution observations (every second day). The mesocosm experiment simulated an upwelling event in the mesocosms by addition of nitrogen (N)-deficient and CO₂-enriched OMZ water. Surface water in the mesocosms was acidified by the OMZ water addition, with pH_T lowered by 0.1–0.2 and pCO₂ elevated to above 900 µatm. Thereafter, surface pCO₂ quickly dropped to near or below the atmospheric level (405.22 µatm in 2017; Dlugokencky and Tans, 2021; NOAA/Global Monitoring Laboratory (GML)) mainly due to enhanced phytoplankton production with rapid CO₂ consumption. Further observations revealed that the dominance of the dinoflagellate Akashiwo sanguinea and contamination of bird excrements played important roles in the dynamics of carbonate chemistry in the mesocosms. Compared to the simulated upwelling, natural upwelling events in the surrounding Pacific waters occurred more frequently with sea-to-air CO₂ fluxes of 4.2–14.0 mmol C m⁻² d⁻¹. The positive CO₂ fluxes indicated our site was a local CO₂ source during our study, which may have been impacted by the coastal El Niño. However, our observations of dissolved inorganic carbon (DIC) drawdown in the mesocosms suggest that CO₂ fluxes to the atmosphere can be largely dampened by biological processes. Overall, our study characterized carbonate chemistry in nearshore Pacific waters that are rarely sampled in such a temporal resolution and hence provided unique insights into the CO₂ dynamics during a rare coastal El Niño event.

1 Introduction

One of the most extensive oxygen minimum zones (OMZs) in the global ocean can be found off central/northern Peru (4–16^∘ S; Chavez and Messié, 2009). High biological productivity is stimulated by permanent upwelling of cold, nutrient-rich water to the surface, supporting remarkable fish production off Peru (Chavez et al., 2008; Montecino and Lange, 2009; Albert et al., 2010). The high primary production also leads to enhanced remineralization of sinking organic matter in subsurface waters which depletes dissolved oxygen (O₂) and creates an intense and shallow OMZ (Chavez et al., 2008). The depletion of O₂ in OMZs plays an important role in the global nitrogen (N) cycle, accounting for 20 %–40 % of N loss in the ocean (Lam et al., 2009; Paulmier and Ruiz-Pino, 2009). Denitrification and anammox processes that occur in O₂-depleted waters remove biologically available N from the ocean and produce an N deficit and hence phosphorus (P) excess with respect to the Redfield ratio ($\mathrm{C}:\mathrm{N}:\mathrm{P}$ = $\mathrm{106}:\mathrm{16}:\mathrm{1}$) in the water column (Redfield, 1963; Deutsch et al., 2001, 2007; Hamersley et al., 2007; Galán et al., 2009; Lam et al., 2009). Upwelling of this N-deficient water has been found to control the surface water nutrient stoichiometry and thus influence phytoplankton growth and community compositions (Franz et al., 2012; Hauss et al., 2012).

Apart from being N-deficient, the OMZ waters are also characterized by enhanced carbon dioxide (CO₂) concentrations and low pH from respiratory processes and are further acidified by increasing anthropogenic atmospheric CO₂ (Feely et al., 2008; Friederich et al., 2008; Paulmier et al., 2008, 2011). Accordingly, surface water carbonate chemistry is influenced by upwelling of CO₂-enriched OMZ water (Van Geen et al., 2000; Capone and Hutchins, 2013). The upwelled CO₂-enriched OMZ water can give rise to surface CO₂ levels > 1000 µatm, pH values as low as 7.6, and under-saturation for the calcium carbonate mineral aragonite (Feely et al., 2008; Hauri et al., 2009). As a result, there is a significant flux of CO₂ from the ocean to the atmosphere off Peru, which is further facilitated by surface ocean warming, making the Peruvian upwelling region a year-round CO₂ source to the atmosphere (Friederich et al., 2008). In contrast, rapid utilization of upwelled CO₂ and nutrients by phytoplankton can occasionally deplete surface CO₂ below atmospheric equilibrium and dampen the CO₂ outgassing (Van Geen et al., 2000; Friederich et al., 2008; Loucaides et al., 2012). The enhanced primary production in turn contributes to increasing export of organic matter, enhanced bacterial respiration, O₂ consumption, and CO₂ production at depth. Such a positive feedback may determine the intensity of the underlying OMZ and promote carbon (C) preservation in marine sediments (Dale et al., 2015).

In response to reduced O₂ solubility and enhanced stratification induced by global warming, OMZs have been intensifying and expanding over the past decades (Stramma et al., 2008, 2010; Fuenzalida et al., 2009). Based on regional observations and model projections, a decline in dissolved O₂ concentrations has been reported for most regions of the global ocean (Matear et al., 2000; Matear and Hirst, 2003; Whitney et al., 2007; Stramma et al., 2008; Keeling et al., 2009; Bopp et al., 2013; Schmidtko et al., 2017; Oschlies et al., 2018). The vertical expansion of OMZs represents shoaling of CO₂-enriched seawater, which has become further enriched by oceanic uptake of anthropogenic CO₂ (Doney et al., 2012; Gilly et al., 2013; Schulz et al., 2019). Since biogeochemical processes in OMZs are directly linked to the C cycle and control surface nutrient stoichiometry, with ongoing ocean warming and acidification, the deoxygenation may have cascading effects on plankton productivity and composition, C uptake, and food web functioning (Keeling et al., 2009; Gruber, 2011; Doney et al., 2012; Gilly et al., 2013; Levin and Breitburg, 2015). Therefore, it is important to monitor the changes in CO₂ when investigating the effects of deoxygenation on marine ecosystems.

To investigate the potential impacts of upwelling on pelagic biogeochemistry and natural plankton communities in the Peruvian OMZ, a large-scale in situ mesocosm study was carried out in the coastal upwelling area off Peru. An upwelling event was simulated in the mesocosms by addition of OMZ waters collected from two different locations where the OMZ was considered to contain different nutrient concentrations and N:P ratios. The ecological and biogeochemical responses in the mesocosms were monitored and compared with those influenced by natural upwelling events in the ambient coastal water surrounding the mesocosms. As part of this collaborative research project, questions specific to the present paper were as follows: (1) how does surface water carbonate chemistry respond to an upwelling event? And (2) how does upwelled OMZ water with different chemical signatures modulate surface water carbonate chemistry? The current study will mainly focus on the temporal changes in surface water carbonate chemistry within the individual mesocosms, including observations made in the ambient Pacific water and a local estimate of air–sea CO₂ exchange, together with the influence of a rare coastal El Niño event (Garreaud, 2018). This provides first insights into how inorganic C cycling links to chemical signatures of OMZ waters in a natural plankton community and its implications for ongoing environmental changes.

2 Material and methods

2.1 Study site

The experiment was conducted in the framework of the Collaborative Research Center 754 “Climate-Biogeochemistry Interactions in the Tropical Ocean” (https://www.sfb754.de/en, last access: 16 June 2020) and in collaboration with the Instituto del Mar del Perú (IMARPE) in Callao, Peru (Fig. 1a). The coastal area off Callao lies within the Humboldt Current System and is influenced by wind-induced coastal upwelling (Bakun and Weeks, 2008).

2.2 Mesocosm setup

Eight Kiel Off-Shore Mesocosms for Future Ocean Simulations (KOSMOS) units (M1–M8), extending 19 m below the sea surface, were deployed by the research vessel Buque Armada Peruana (BAP) Morales and moored at 12.06^∘ S, 77.23^∘ W in the coastal upwelling area off Callao, Peru (Fig. 1a) on 23 February 2017 (late austral summer). The technical design of these seagoing mesocosms is described by Riebesell et al. (2013). For a more detailed description of the mesocosm deployment and maintenance in this study, please refer to Bach et al. (2020a).

The mesocosm bags were filled with surrounding seawater through the upper and lower openings. Both openings were covered by screens with a mesh size of 3 mm to avoid enclosing larger organisms such as fish. The mesocosm bags were left open below the water surface for 2 d, allowing free exchange with surrounding coastal water. On 25 February, mesocosm bags were closed with the screens removed, tops pulled above the sea surface, and bottoms sealed with 2 m long conical sediment traps (Fig. 1b). The experiment started with the closure of the mesocosms (day 0) and lasted for 50 d. Each mesocosm bag enclosed a seawater volume of ∼ 54 m³. After the bags were closed, daily or every-second-day sampling was performed to monitor the initial conditions of the enclosed water before simulating an upwelling event on day 11 and 12 (see Sect. 2.4 for details).

Figure 1 The study site of the mesocosm experiment (a) created and modified using Ocean Data View (Reiner Schlitzer, Ocean Data View, https://odv.awi.de/, last access: 6 April 2021) and a schematic illustration of a KOSMOS mesocosm unit (b). We acknowledge reprint permission from the American Geophysical Union as parts of this drawing were used for a publication by Bach et al. (2016). The star symbol marks the approximate location of mesocosm deployment.

2.3 Simulated upwelling and salt addition

To simulate an upwelling event in the mesocosms, OMZ-influenced waters were collected from the nearby coastal area and added to the mesocosms. Two OMZ water masses were collected at Station 1 (12^∘01.70^′ S, 77^∘13.41^′ W) at a depth of ∼ 30 m and at Station 3 (12^∘ 02.41^′ S, 77^∘22.50^′ W) at a depth of ∼ 70 m, respectively, using a deep-water collection system as described by Taucher et al. (2017). These two water masses were sampled for chemical and biological variables as done in the mesocosms (see Sect. 2.4). The OMZ water collected from Station 3 had a dissolved inorganic nitrogen (DIN) concentration of 4.3 µmol L⁻¹ (denoted as “Low DIN” in this paper) and was added to M2, M3, M6, and M7. The OMZ water from Station 1 had a DIN of 0.3 µmol L⁻¹ (denoted as “Very low DIN” in this paper) and was added to M1, M4, M5, and M8. Before OMZ water addition, approximately 9 m³ of seawater was removed from 11–12 m of each mesocosm on 5 March (day 8). During the night of 8 March (day 11), ∼ 10 m³ of OMZ water was added to 14–17 m of each mesocosm. On 9 March (day 12), ∼ 10 m³ of seawater was removed from 8–9 m followed by an addition of ∼ 12 m³ OMZ water to 0–9 m of each mesocosm.

To maintain a low-O₂ bottom layer in the mesocosms and avoid convective mixing induced by heat exchange with the surrounding Pacific, 69 L of a concentrated sodium chloride (NaCl) brine solution was added to the bottom of each mesocosm (10–17 m) on day 13, which increased the bottom salinity by ∼ 0.7 units. Following that, turbulent mixing induced by sampling activities continuously interrupted the artificial halocline. Hence, on day 33, 46 L of the NaCl brine solution was added again to the bottom of each mesocosm (12.5–17 m), which increased the bottom salinity by ∼ 0.5 units. At the end of the experiment after the last sampling (day 50), 52 kg of NaCl brine was added again to each mesocosm to calculate the enclosed seawater volume from a measured salinity change by ∼ 0.2 units (see Czerny et al., 2013a, and Schulz et al., 2013, for details). The average final volume for each mesocosm bag was calculated at ∼ 54 m³. With known sampling volumes and deep-water addition volumes during the experiment, the enclosed volumes of each mesocosm on each sampling day could be calculated. The NaCl solution for the halocline establishments had been prepared in Germany by dissolving 300 kg of food industry grade NaCl (free of anti-caking agents) in 1000 L of deionized water (Milli-Q, Millipore) and purified with ion exchange resin (Lewatit^® MonoPlus TP 260, Lanxess, Germany) to minimize potential contaminations with trace metals (Czerny et al., 2013a). The NaCl solution for the volume determination was produced on-site using locally purchased table salt. For a more detailed description of OMZ water and salt additions, please refer to Bach et al. (2020a).

2.4 Sampling procedures and CTD operations

Sampling was carried out in the morning (07:00–11:00 local time) daily or every second day throughout the entire experimental period. Depth-integrated samples were taken from the surface (0–10 m for day 3–28) and bottom layer (10–17 m for day 3–28) of the mesocosms and the surrounding coastal water (named “Pacific”) using a 5 L integrating water sampler (IWS, HYDRO-BIOS, Kiel). Due to the deepening of the oxycline as observed from the CTD profiles, the sampling depth for the surface was adjusted to 0–12.5 m while that for the bottom was changed to 12.5–17 m from day 29 until the end of the experiment (day 50).

For gas-sensitive variables such as pH and dissolved inorganic carbon (DIC), 1.5 L of seawater from each integrated depth in each mesocosm was taken directly from the fully filled 5 L integrating water sampler. Clean polypropylene sampling bottles (rinsed with deionized water in the laboratory; Milli-Q, Millipore) were pre-rinsed with sample water immediately prior to sampling. Bottles were filled from bottom to top using pre-rinsed Tygon tubing with overflow of at least one sampling bottle volume (1.5 L) to minimize the impact of CO₂ air–water gas exchange. Nutrient samples were collected into 250 mL polypropylene bottles using pre-rinsed Tygon tubing (see Bach et al., 2020a, for details). Sample containers were stored in cool boxes for ∼ 3 h and protected from sunlight and heat before being transported to the shore. Once in the lab, sample water was sterile-filtered by gentle pressure using syringe filters (0.2 µm pore size), Tygon tubing, and a peristaltic pump to remove particles that may cause changes to seawater carbonate chemistry (Bockmon and Dickson, 2014). For DIC measurements, the water was filtered from the bottom of the 1.5 L sample bottle into 100 mL glass-stoppered bottles (DURAN) with an overflow of at least 100 mL to minimize contact with air. Once the glass bottle was filled with sufficient overflow, it was immediately sealed without headspace using a round glass stopper. This procedure was repeated to collect a second bottle (100 mL) of filtered water for pH measurements. The leftover seawater was directly filtered into a 500 mL polypropylene bottle for total alkalinity (TA) measurements (non-gas-sensitive). Filtered DIC and pH samples were stored at 4 ^∘C in the dark, and TA samples were at room temperature in the dark until further analysis. Samples were analyzed for DIC and pH on the same day of sampling, while TA was determined overnight (see Sect. 2.5 for analytical procedures).

CTD casts were performed with a multiparameter logging probe (CTD60M, Sea & Sun Technology) in the mesocosms and Pacific on every sampling day. From the CTD casts, profiles of salinity, temperature, pH, dissolved O₂, chlorophyll a (chl a), and photosynthetically active radiation were obtained (see Schulz and Riebesell, 2013, and Bach et al., 2020a, for details).

2.5 Carbonate chemistry and nutrient measurements

Total alkalinity was determined at room temperature (22–32 ^∘C) by a two-stage open-cell potentiometric titration using a Metrohm 862 Compact Titrosampler, Aquatrode Plus (Pt1000), and 907 Titrando unit in the IMARPE laboratory following Dickson et al. (2003). The acid titrant was prepared by preparing a 0.05 mol kg⁻¹ hydrochloric acid (HCl) solution with an ionic strength of ca. 0.7 mol kg⁻¹ (adjusted by NaCl). Approximately 50 g of sample water from each sample was weighed into the titration cell with the exact weight recorded (precision 0.0001 g). After the two-stage titration, the titration data between a pH of ∼ 3.5 and 3 were fitted to a modified non-linear Gran approach described in Dickson et al. (2007) using MATLAB (The MathWorks). The results were calibrated against certified reference material (CRM) batch 142 (Dickson, 2010) measured on each measurement day. In this paper, measured TA values refer to the measured values that have been calibrated against the CRM.

Seawater pH_T (total scale) was determined spectrophotometrically by measuring the absorbance ratios after adding the indicator dye m-cresol purple (mCP) as described in Carter et al. (2013). Before measurements, samples were acclimated to 25.0 ^∘C in a thermostatted bath. The absorbance of samples with mCP was determined on a Varian Cary 100 double-beam spectrophotometer (Varian), scanning between 780 and 380 nm at a 1 nm resolution. During the spectrophotometric measurement, the temperature of the sample was maintained at 25.0 ^∘C by a water bath connected to the thermostatted 10 cm cuvette. The pH_T values were calculated from the baseline-corrected absorbance ratios and corrected for in situ salinity (obtained from CTD casts) and pH change caused by dye addition (using the absorbance at the isosbestic point, i.e., 479 nm) as described in Dickson et al. (2007). To minimize potential CO₂ air–water gas exchange, a syringe pump (Tecan Cavro XLP) was used for sample and dye mixing and cuvette injection (see Schulz et al., 2017, for details). For the dye correction, a batch of sterile filtered seawater of known salinity was prepared. The pH_T was determined once for an addition of 7 µL of dye and once for an addition of 25 µL at five pH levels (raised to 7.95 with NaOH and lowered to 7.74, 7.58, 7.49, and 7.36 with HCl stepwise). The pH change resulting from the dye correction addition was calculated from the change in measured absorbance ratio for each pair of dye additions (see Clayton and Byrne, 1993, and Dickson et al., 2007, for details). The dye-corrected pH_T values measured at 25.0 ^∘C and atmospheric pressure were then re-calculated for in situ temperature and pressure as determined by CTD casts (averaged over 0–10/12.5 m for surface and 10/12.5–17 m for bottom). For carbonate chemistry speciation calculations (see Sect. 2.6), the dye-corrected pH_T values were used as one of the input parameters.

Dissolved inorganic carbon was measured by infrared absorption using a LI-COR LI-7000 on an AIRICA system (MARIANDA, Kiel; see Taucher et al., 2017, and Gafar and Schulz, 2018, for details). The results were calibrated against CRM batch 142 (Dickson, 2010). Unfortunately, due to a malfunctioning of the AIRICA system, we obtained measured DIC data only up to 7 March (day 10). Therefore, measured TA and pH_T were used for calculations of carbonate system parameters at in situ temperature and salinity, but we used DIC measurements from day 3–10 for consistency checks of calculated carbonate chemistry parameters. In this paper, measured DIC values refer to the measured values that have been calibrated against the CRM.

Inorganic nutrients were analyzed colorimetrically (NO${}_{\mathrm{3}}^{-}$, NO${}_{\mathrm{2}}^{-}$, PO${}_{\mathrm{4}}^{\mathrm{3}-}$, and Si(OH)₄) and fluorimetrically (NH${}_{\mathrm{4}}^{+})$ using a continuous-flow analyzer (QuAAtro AutoAnalyzer with integrated photometers, SEAL Analytical) connected to a fluorescence detector (FP-2020, JASCO). All colorimetric methods were conducted according to Murphy and Riley (1962), Mullin and Riley (1955a, b), and Morris and Riley (1963) and corrected following the refractive index method developed by Coverly et al. (2012). For details of the quality control procedures, see Bach et al. (2020a).

2.6 Carbonate chemistry speciation calculations and propagated uncertainties

Calculations of carbonate chemistry parameters (in situ pH_T, DIC, pCO₂, and calcium carbonate saturation state for calcite and aragonite) were performed with the Excel version of CO2SYS (Version 2.1; Pierrot et al., 2006) using K1 and K2 dissociation constants from Mehrbach et al. (1973) which were refitted by Lueker et al. (2000). The dissociation constants for KHSO₄ from Dickson (1990) and for total boron from Uppström (1974) were applied in the calculations (see Orr et al., 2015, for details). The observed pH_T and TA as well as inorganic nutrient concentration (phosphate and silicic acid) were used as input CO₂ system parameters. In situ salinity and temperature were obtained by CTD casts and averaged over surface (0–10 or 0–12.5 m) and bottom (10–17 or 12.5–17 m) waters for each sampling day. In situ pressure was approximated for surface (5 dbar) and bottom (13.5 or 14.75 dbar) waters. For details of calculation procedures and choices of constants, see Lewis et al. (1998) and Orr et al. (2015).

To evaluate the performance of pH_T and TA measurements, quality control procedures were performed. First, standard deviations of pH_T measurements were graphed over time. Measured TA values of a control sample (CRM batch 142; Dickson, 2010) were plotted over time and compared to the warning and control limits calculated from their mean and standard deviation (for details please see Dickson et al., 2007) as well as the certified value of the CRM. To compute a range control chart for the evaluation of measurement repeatability, the absolute difference between duplicate measurements of CRMs on each sampling day was calculated and plotted over time and compared to the warning and control limits calculated from their mean and standard deviation (for details see Dickson et al., 2007).

We used the R package “seacarb” with a Gaussian approach and an input variable pair (pH_T, TA) to calculate uncertainties for calculated CO₂ system parameters (Orr et al., 2018; Gattuso et al., 2020). The contributions of input uncertainties in nutrient concentrations and in situ salinity and temperature to the uncertainties in the CO2SYS-based calculations are often small (< 0.1 %; Orr et al., 2018), so they were not considered in our propagation. The input uncertainties in pH_T and TA were estimated based on our measurements (Table 1). Standard uncertainties include random and systematic errors. For TA, systematic errors were removed by calibrating the measured results using CRMs (see Sect. 2.5). Hence, the random error in TA, estimated by the averaged standard deviation of all the CRM measurements (4.4 µmol kg⁻¹; n=62), was used as the standard uncertainty. For pH_T, an uncertainty of 0.01 was used as the standard uncertainty. Due to the unavailability of CRMs that correct for systematic error in pH measurements, the standard deviations of repeated measurements (0.0012; n=377) only accounted for the random components of standard uncertainties (Orr et al., 2018). Therefore, we used 0.01 in our uncertainty propagation as an approximation of the total standard uncertainty for pH_T, which has been used in previous assessments (Orr et al., 2018).

Table 1 Standard uncertainties in pH_T and TA estimated based on our measurements are denoted by μ(pH_T) and µ(TA). Based on μ(pH_T) and μ(TA), propagated uncertainties were estimated for each data point in R and averaged for each reported variable (μ), with standard deviation (σ), minimum (min), and maximum (max) values presented. The relative percentages (%) of propagated standard uncertainties were calculated by dividing the propagated uncertainty by the corresponding data point and averaged for each reported variable (μ), with σ, min, and max values presented.

The air–sea flux of CO₂ (FCO₂, mmol C m⁻² d⁻¹) in the Pacific was determined based on

$$\begin{array}{}\text{(1)}& {\mathrm{FCO}}_{\mathrm{2}}=k{K}_{\mathrm{0}}\mathrm{\Delta }p{\mathrm{CO}}_{\mathrm{2}},\end{array}$$

where k is the gas transfer velocity parameterized as a function of wind speed, K₀ is the solubility of CO₂ in seawater dependent on in situ salinity and temperature (Weiss, 1974), and ΔpCO₂ is the difference between pCO₂ in the surface water and in the atmosphere (Wanninkhof 2014). Wind data were averaged over 2 sampling days for the sampling location from a satellite-derived gridded dataset (Global Land Data Assimilation System model, near-surface wind speed, $\mathrm{0.25}\times \mathrm{0.25}{}^{\circ }$, 3 h temporal resolution, 12.375 to 11.875^∘ S, 77.375 to 76.875^∘ W), obtained from NASA Giovanni (Rodell et al., 2004; Beaudoing and Rodell, 2020). In situ salinity and temperature were obtained from the CTD casts (see Sect. 2.4). Calculated pCO₂ based on (pH_T, TA) and an estimated atmospheric pCO₂ of 405.22 µatm (referenced to year 2017; Dlugokencky and Tans, 2021; NOAA/GML) were used in the air–sea flux estimation.

3 Results

3.1 Responses of surface layer nutrient concentrations

The OMZ-influenced water masses were collected from two locations and added to the mesocosms to simulate an upwelling event (see Sect. 2.3). The two water masses were named Low DIN and Very low DIN, respectively, based on their DIN concentrations (Table 2). Both water masses shared similar silicic acid (Si) and phosphate (PO${}_{\mathrm{4}}^{\mathrm{3}-})$ concentrations but differed in DIN concentration. The Low DIN water had a DIN concentration of 4.3 µmol L⁻¹, 14 times as high as that of the Very low DIN water (0.3 µmol L⁻¹, Table 2).

Table 2 Inorganic nutrient concentrations of the two collected deep-water masses. Please note that DIN is the sum of nitrate, nitrite, and ammonium. P is phosphate. Si is silicic acid.

On day 10 before OMZ water addition, the average surface DIN concentrations of the two treatment groups were similar (3.4 µmol L⁻¹) but lower than that in the Pacific (9.8 µmol L⁻¹, Table 3). Surface layer DIN concentrations in the mesocosms ranged between 2.0 and 6.0 µmol L⁻¹ before OMZ water addition (Fig. 2a). The addition of OMZ water elevated surface DIN in the Low DIN mesocosms to 3.6–6.4 µmol L⁻¹ but lowered that in the Very low DIN to 0.9–2.0 µmol L⁻¹. The average surface DIN concentration in the Very low DIN decreased to 1.6 µmol L⁻¹ while the Low DIN slightly increased to 4.7 µmol L⁻¹ (Table 3), followed by a sharp depletion on day 16 except for M3. M3 received the highest input of DIN (6.4 µmol L⁻¹) and was not depleted until day 24. Despite several small peaks in M3, M4, M5, and M6 (≤1.6 µmol L⁻¹), surface DIN concentrations in the mesocosms were at around the limits of detection (LODs – NH${}_{\mathrm{4}}^{+}$, 0.063 µmol L⁻¹; NO${}_{\mathrm{2}}^{-}$, 0.054 µmol L⁻¹; NO${}_{\mathrm{3}}^{-}$, 0.123 µmol L⁻¹) most of the time after depletion. A slight rise could be observed from day 44 towards the last sampling day (day 48). In the Pacific, surface layer DIN concentration was mostly greater than 5 µmol L⁻¹ (except on day 16 and 18) and became considerably higher during the second half of the experiment (> 10 µmol L⁻¹ for day 26–44, Fig. 2a).

Table 3 DIN concentration (µmol L⁻¹) in the surface layer of each mesocosm (M1–M8) and the average DIN concentration (µmol L⁻¹) for each treatment (Low DIN and Very low DIN; n=4) before (t10) and after deep-water addition (t13). The DIN concentration in the surface Pacific water is also shown.

Surface layer PO${}_{\mathrm{4}}^{\mathrm{3}-}$ concentrations in the mesocosms initially ranged between 1.1 and 1.5 µmol L⁻¹ and were elevated by OMZ water addition to around 1.9 µmol L⁻¹ (Fig. 2b). Thereafter, PO${}_{\mathrm{4}}^{\mathrm{3}-}$ exhibited a slow but steady decline until the end of the study with a slightly higher decrease in Low DIN mesocosms (blue symbols, Fig. 2b). Throughout the study, PO${}_{\mathrm{4}}^{\mathrm{3}-}$ in the mesocosms was never lower than 1.1 µmol L⁻¹. Surface layer PO${}_{\mathrm{4}}^{\mathrm{3}-}$ in the Pacific was generally higher, fluctuating between 1.4 and 2.9 µmol L⁻¹. In the mesocosms, enhanced chl a concentrations were observed at depths shallower than 5 m and below 15 m before OMZ water addition (Fig. 2c). Following OMZ water addition, a chl a maximum occurred at ∼ 10 m and persisted until day 40, except for M3 and M4 with a 1-week-delayed increase in the former and a lack of bloom in the latter (Fig. 2c). After day 40, chl a concentrations in all mesocosms (except for M4) increased to 12–38 µg L⁻¹ with a bloom occurring in 0–10 m (Fig. 2c). Throughout the study, a chl a maximum was continuously observed above 10 m in the Pacific (Fig. 2c).

Figure 2 Temporal dynamics of depth-integrated surface DIN concentration (a) and PO${}_{\mathrm{4}}^{\mathrm{3}-}$ concentration (b) and vertical distribution of chl a concentration determined by CTD casts (c). The solid black lines on top of the colored contours represent the average values over the entire water column, with the corresponding additional y axes on the right. The vertical white lines represent the day when OMZ water was added to the mesocosms. Color codes and symbols denote the respective mesocosms. Abbreviations: OWA, OMZ water addition; SA, salt addition. The dataset is available at https://doi.pangaea.de/10.1594/PANGAEA.923395 (Bach et al., 2020b).

3.2 Temporal dynamics of carbonate chemistry

Before OMZ water addition, surface layer pH_T in the mesocosms ranged between 7.80–7.94 with a slight decline by ∼ 0.1 over time (Fig. 3a). The initial surface layer TA ranged between 2310 and 2330 µmol kg⁻¹ (Fig. 3b, day 3–12). Surface layer pCO₂ and DIC ranged from 541 to 749 µatm and 2119 to 2180 µmol kg⁻¹, respectively (Fig. 3c, d).

The two collected OMZ water masses shared similar carbonate chemistry properties despite the differences in DIN concentrations. In both water masses, pH_T was ∼ 7.48, DIC was ∼ 2305–2310 µmol kg⁻¹, TA was ∼ 2337 µmol kg⁻¹, and pCO₂ was between 1700 and 1780 µatm (Table 4).

Table 4 The in situ pH_T, TA, DIC, pCO₂, Ω_Ar, and Ω_Ca of the two collected OMZ water masses.

Surface DIC and pCO₂ were elevated from ∼ 2150 µmol kg⁻¹ and ∼ 600 µatm to ∼ 2200 µmol kg⁻¹ and ∼ 900 µatm (except M7), respectively, by OMZ water addition without distinct differences between the two treatments (Mann–Whitney U test, p>0.05; Fig. 3c). Following OMZ water addition, surface pCO₂ in the mesocosms decreased quickly and reached minima at 340–500 µatm (except M3 and M4) on day 24 and 26. These minima corresponded with DIC minima at 2040–2110 µmol kg⁻¹ and pH_T maxima at 7.9–8.1 (except M3 and M4; Fig. 3c, d). After reaching the minima, surface layer pCO₂ exhibited a steady increase to 410–680 µatm from day 24 to day 38 and later declined in M3, M5, and M7 while the rest remained relatively stable until day 42 (Fig. 3c). Interestingly, and unlike the other mesocosms, after OMZ water addition, pCO₂ in M3 steadily declined from 928 to 342 µatm until the end of the experiment while that in M4 remained constantly higher than the other mesocosms (>700 µatm), with a slightly decreasing trend to 645 µatm towards the end of the study (Fig. 3c).

In the Pacific, much lower surface pH_T and higher surface pCO₂ and DIC were observed compared to the mesocosms, with an average of 7.7 (7.6–7.8), 1078 µatm (775–1358 µatm), and 2221 µmol kg⁻¹ (2173–2269 µ mol kg⁻¹; minimum-to-maximum range in parentheses; Fig. 3c, d), respectively. TA in the Pacific was initially similar to that in the mesocosms, fluctuating between 2310 and 2330 µmol kg⁻¹, and later decreased to ∼ 2310 µmol kg⁻¹ for the rest of the study.

Surface waters in the mesocosms and the Pacific were always saturated with respect to calcite and aragonite throughout the entire experimental period, with lower values observed in the Pacific (Fig. 4a, c). Bottom waters in the mesocosms and Pacific were always saturated with respect to calcite during the experiment (Fig. 4b), while bottom waters in the Pacific were under-saturated with respect to aragonite before day 13 (0.88–0.99) and had Ω_Ar values slightly above 1.0 for the rest of the study period (Fig. 4d).

Figure 3 Temporal dynamics of measured depth-integrated surface pH_T (a) and TA (b) and calculated pCO₂ (c) and DIC (d). The error ribbons present measurement and propagated standard uncertainties in the calculations, respectively. Color codes and symbols denote the respective mesocosms. Abbreviations: OWA, OMZ water addition; SA, salt addition.

Figure 4 Temporal dynamics of depth-integrated surface calcite saturation state (a), bottom calcite saturation state (b), surface aragonite saturation state (c), and bottom aragonite saturation state (d) in the mesocosms and the surrounding Pacific. The error ribbons present the propagated standard uncertainties in the calculations. When Ω>1 (above dashed red line), seawater is supersaturated for calcium carbonate. When Ω<1 (below dashed red line), seawater is under-saturated for calcium carbonate. Color codes and symbols denote the respective mesocosms. Abbreviations: OWA, OMZ water addition; SA, salt addition.

3.3 Air–sea CO₂ fluxes in the Pacific

Positive FCO₂ values indicate CO₂ outgassing from the surface waters to the atmosphere, while negative values indicate a CO₂ flux from the atmosphere to the ocean. The air–sea CO₂ flux in the Pacific was constantly positive throughout our study, fluctuating from 4.2 to 14.0 mmol C m⁻² d⁻¹ over time (Fig. 5a). The minima of FCO₂ occurred on day 26 and 30, while the maximum occurred on day 32 when near-surface wind was the highest (2.89 m s⁻¹, Fig. 5b), corresponding to the minima and maxima of surface pCO₂. Co-occurring with a decrease in surface temperature to below 19 ^∘C after day 36 (Fig. 5c), FCO₂ slightly declined from ∼ 10 to ∼ 6 mmol C m⁻² d⁻¹ (Fig. 5a). FCO₂ was positively correlated with near-surface wind speed (r²=0.4). No correlation was found between FCO₂ and temperature (r²=0).

Figure 5 Temporal dynamics of surface air–sea CO₂ flux (a), near-surface wind speed (b), and surface temperature (c) in the Pacific. FCO₂>0 (above dashed red line) indicates CO₂ outgassing from the sea surface to the atmosphere. FCO₂<0 (below dashed red line) indicates a CO₂ flux from the atmosphere to the sea.

4 Discussion

4.1 Quality control and propagated uncertainties

To compare the sensitivity of different calculated variables to uncertainties in the input variables, the propagated uncertainties were averaged for each calculated variable, reported in numerical values and percentages relative to the calculated values of each variable (Table 1). Among the 4 reported variables, Ω_Ca and Ω_Ar were the most sensitive to uncertainties in pH_T and TA with an average uncertainty of 5.1 %. This adds ambiguity to whether the bottom water (10–17 m for day 3–28, 12.5–17 m for day 29–50) in the Pacific was under-saturated with respect to aragonite when Ω_Ar was oscillating near 1 (Fig. 4d). The propagated uncertainty in pCO₂ was slightly lower (3.8 %), while DIC was the least sensitive (0.3 %).

We examined the internal consistency between DIC measurements and calculations. DIC was measured from day 3 until the malfunction of the instrument on day 10. In total, 53 sets of measured DIC and calculated DIC (from measured pH_T and TA) values were obtained from day 3 to day 10 and compared to test their consistency (Fig. 6a). The calculated DIC values were generally in agreement with the measured values (r²=0.985; p<0.005), showing that the calculations made an overall good prediction for the measured DIC values. The average of the residuals (calculated DIC − measured DIC) was −8.27 ± 6.9 µmol kg⁻¹, indicating an underestimation of calculated DIC. This result is consistent with a previous observation of underestimated calculated DIC (pH_T, TA) compared with measured DIC when applying the same set of constants (−6.6 ± 7.9 µmol kg⁻¹; Raimondi et al., 2019). The reasons for such underestimation have not been addressed in previous studies and remain unclear. No significant relationships with input variables pH_T and TA (r²=0.12 for both) and temperature (r²=0.30) were found in the DIC residuals (salinity remained the same from day 3 to day 10). The lack of correlation with pH_T and TA indicated that the underestimation in calculated DIC was not a result from changes in pH_T and TA. Although dissociation constants are known to be salinity- and temperature-dependent, the lack of correlation between DIC residuals and temperature may be attributed to the relatively narrow ranges of temperature in the mesocosms (17.9–20.9 ^∘C from day 3–10). The offsets were typically larger at lower temperatures (e.g., samples from the Arctic; Chen et al., 2015).

To assess the quality of carbonate chemistry measurements in this study, the stability and performance of measurements were evaluated. The standard deviation of triplicate pH_T measurements varied by up to 0.003 with an average of 0.0012 throughout the whole experiment (Fig. 6b). The average standard deviation was in agreement with reported analytical precisions of pH (0.003, Orr et al., 2018; 0.002, Raimondi et al., 2019; Ma et al., 2019).

For TA, triplicate measurements of CRM distributed before and after the sample measurements were carried out on each measuring day to monitor the stability of the measurement process and the performance of the system. Based on the offsets, a correction factor was applied to the measured values of samples on each sampling day to calibrate for instrument drift. As shown in Fig. 6c, 90.5 % of the measured TA values of CRM fell between warning limits (upper warning limit, UWL, and lower warming limit, LWL) with one data point falling outside the control limits (upper control limit, UCL, and lower control limit, LCL), overall suggesting a relatively stable measurement system. The average measured TA was 2209.9 µmol kg⁻¹, which was 17.69 µmol kg⁻¹ lower than the certified concentration of the CRM (2227.59 µmol kg⁻¹), indicating a relatively poor accuracy (compared to the suggested bias of less than 2 µmol kg⁻¹; Dickson et al., 2003, 2007). The poor accuracy could be attributed to the fact that the concentration of the acid titrant was not checked after being prepared, as suggested in the protocol (Dickson et al., 2003). A range control chart was computed based on duplicate measurements of CRM made prior to the sample measurements on each sampling day to evaluate the consistency of the offset between measured and certified TA values over the course of the study (Fig. 6d; Dickson et al., 2007). The absolute difference (range) between the repeated CRM measurements was on average 1.4 µmol kg⁻¹. All the range values fell below the UWL (3.50 µmol kg⁻¹, Fig. 6d), suggesting a relatively good precision of the measurement system.

Figure 6 Comparison of calculated values of DIC (pH_T, TA) and measured values (a). The black line is the regression line, with the corresponding equation and r² shown in the top-left corner. The dashed blue line shows the regression line forced through the origin. Standard deviations are of all the triplicate pH_T measurements on each sampling day over the study period. The dashed orange line shows the average (n=377) of the standard deviations (b). TA values of CRM measurements on each sampling day over the study period. The dashed orange line shows the average (n=62) of the measured values, and the dashed green line indicates the certified value of the CRM (c). The absolute difference in TA values between duplicate CRM measurements (range) on each sampling day over the study period. The dashed orange line shows the average (n=21) of the ranges (d). Abbreviations: UCL, upper control limit; UWL, upper warning limit; LWL, lower warning limit; LCL, lower control limit.

4.2 CO₂ responses to the simulated upwelling event

At the beginning of the experiment, surface pCO₂ levels in the mesocosms were > 500 µatm (Fig. 3c). This suggests that we initially enclosed an upwelled water mass that was enriched with respiratory CO₂. The addition of OMZ water with high concentrations of CO₂ to the mesocosms reduced the surface pH_T by 0.1–0.2 and increased the surface pCO₂ to > 900 µatm (except for M7, which was 819.4 µatm on day 13). The simulated upwelling substantially reduced the variability in CO₂ between mesocosms because OMZ water addition replaced ∼ 20 m³ of seawater in each mesocosm (out of ∼ 54 m³). The enhanced pCO₂ level is comparable with our observations in the ambient Pacific water (> 775 µatm, Fig. 3c). These values also agree with reported observations for our study area in 2013 (> 1200 µatm in the upper 100 m and > 800 µatm at the surface; Bates, 2018).

In the days after OMZ water addition, surface pCO₂ in the mesocosms dropped near or below the atmospheric level (405.22 µatm; Dlugokencky and Tans, 2021; NOAA/GML) with a decline in DIC by ∼ 100 µmol kg⁻¹(except M3 and M4; Fig. 3c, d). The declining pCO₂ could be partially attributed by CO₂ outgassing due to a high CO₂ gradient from the sea surface to the air. Due to a rare coastal El Niño event (Garreaud, 2018), the CO₂ loss process may have been enhanced by a rapid surface warming (19.8–21.0 ^∘C from day 14 to 36, Fig. 5) which reduced surface CO₂ solubility (Zeebe and Wolf-Gladrow, 2001). However, air–sea gas exchange could not explain surface CO₂ under-saturation in relation to the atmosphere, as observed in response to OMZ water addition in some mesocosms (Van Geen et al., 2000; Friederich et al., 2008; Fig. 3c). Biological production typically has impacts on CO₂ drawdown that are 1 to 4 times greater than air–sea gas exchange in the equatorial Pacific where surface waters are exposed to local wind stress (Feely et al., 2002). This interpretation is supported by the continuously high DIC in M4 where photosynthetic biomass buildup was substantially lower (Fig. 3d). Hence, the depletion of nutrients (Fig. 2a, b) and increase in chl a concentration (Fig. 2c; Bach et al., 2020a) strongly suggest that the loss of DIC (except M4) was primarily driven by biological uptake and phytoplankton growth. Nevertheless, it is difficult to dissect how much CO₂ was outgassed and how much was taken up photosynthetically as we did not measure air–sea gas exchange in the mesocosms (please note that equations from Wanninkhof, 2014, are not applicable to mesocosms; Czerny et al., 2013b). In previous mesocosm studies, N₂O addition has been a common practice to monitor air–sea gas exchange in the mesocosms (Czerny et al., 2013b). However, this was not carried out in our study because it would have interfered with ¹⁵N label incubations to determine N loss processes (Schulz et al., 2021). Due to high variability in DOC data and the poorly constrained gas exchange of CO₂, C budgeting often comes with high uncertainties and large errors, even for a relatively simple dataset (Czerny et al., 2013b; Boxhammer et al., 2018). It becomes even more difficult for the current dataset because the water column was not homogeneously mixed like in previous studies. With the inability to estimate CO₂ gas exchange, it was impossible to calculate a reasonable C budget for this study. Before OMZ water addition, dissolved inorganic N:P ratios in the mesocosms ranged from 1.6 to 3.5 (data not shown), indicating N is the limiting nutrient in the water column (Bach et al., 2020a). Not surprisingly, the uptake of DIC was higher in the Low DIN mesocosms, which received more input of DIN from OMZ water addition, with on average 41.0 µmol kg⁻¹ higher drawdown compared to the Very low DIN from day 13 to day 24 (excluding M3 and M4; Mann–Whitney U test, p=0.05; Table S1). This observation agrees with the general expectations that the addition of limiting nutrients to the water column should enhance biological biomass buildup. Such differences in DIC uptake, however, were not reflected in the buildup of particulate organic carbon (POC) in the mesocosms (excluding M3 and M4; Mann–Whitney U test, p>0.1). As mentioned above, the differences in OMZ water DIN between the two treatments were minor, and hence, their potential to trigger treatment differences was small. Also, Akashiwo sanguinea was persistent in the water column in the mesocosms (except M4 where it never bloomed), retaining the biomass in the water column and not sinking out until the very end of the experiment (Bach et al., 2020a). Due to the developing N limitation after the biomass buildup there, much of the consumed DIC could have been channeled to the dissolved organic carbon (DOC) pool. Indeed, we observed a pronounced increase in DOC following OMZ water addition (except for M4; Igarza et al., 2022). The increase in DOC may be attributed to extracellular release by phytoplankton due to nutrient limitation or cellular lysis of phytoplankton cells by bacteria (Myklestad, 2000; Igarza et al., 2022).

After day 24, variability in carbonate chemistry between individual mesocosms increased, with a general trend of recovering from CO₂-under-saturated conditions during the peak of the bloom (except for M3 and M4, Fig. 3c). One factor that may have controlled the differences in CO₂ increase are the mesocosm-specific phytoplankton succession patterns. A shift from a diatom-dominated community to a dominance of dinoflagellates (in particular Akashiwo sanguinea) occurred when DIN was exhausted, which was absent in M3 and M4 (Bach et al., 2020a). The different succession patterns in the plankton community are the most likely explanation of why M3 and M4 behaved differently from the others in terms of surface layer productivity and hence carbonate chemistry. Although the rates of DIN depletion in M3 and M4 were similar to the others, the reduction in pCO₂ in M3 experienced a 1-week delay, which is consistent with the delayed buildup of chl a biomass (Figs. 2c, 3c). On the other hand, the pCO₂ level in M4 remained constantly elevated throughout the experiment, indicative of a lack of a phytoplankton bloom (Figs. 2c, 3c). M4 was the only mesocosm where A. sanguinea remained undetectable, whereas a delayed and reduced contribution by A. sanguinea was observed in M3. This strongly suggests that A. sanguinea was a key factor driving the trend of carbonate chemistry in the mesocosms.

Near the end of the experiment, a slight decline in pCO₂ became apparent in the mesocosms, which co-occurred with a second phytoplankton bloom observed in the uppermost layer of the water column (Figs. 2c, 3c). This bloom was likely fueled by surface eutrophication due to defecating sea birds. During the last part of our experiment, Inca terns (Larosterna inca) were frequently observed to rest on the roofs and the edges of the mesocosms (Bach et al., 2020a). Bird excrements, dropped into the mesocosms, are known to be enriched in inorganic nutrients (Bedard et al., 1980). They generally contain 60 % water, 7.3 % N and, 1.5 % P, and the main form of N is uric acid and ammonium, which makes them slightly acidic (De La Peña-Lastra, 2021). Therefore, the droppings may lower the pH of the surface water. However, this should have been visible in decreasing alkalinity, which was not the case. The excrements may also be high in dissolved organic nitrogen (DON), evidenced by a substantial increase in DON concentrations in the mesocosm surface from day 38 onward (Igarza et al., 2022). The triggered surface eutrophication and phytoplankton blooms were noticeable from an accumulation of chl a biomass above the mixed layer in the mesocosms near the end of the study (Fig. 2c). As a result, another drawdown of DIC could be observed in the mesocosms except for M4, M6, and M8. While the buildup of chl a was comparable with that triggered by OMZ water addition, the drawdown in DIC was less pronounced, potentially counteracted by the release of CO₂ by enhanced respiration and remineralization following the previous bloom. Also, the second bloom occurred in the top 2 meter in the mesocosms (Fig. 2c), where gas exchange can quickly replete the DIC drawdown during photosynthesis and biomass build up.

There are potential complications when monitoring the carbonate chemistry dynamics in an enclosed mesocosm. First, the brine addition to the mesocosms could have influenced the community composition which played an important role in driving the carbonate chemistry. However, the effects of the brine on the enclosed organisms have been discussed in past studies and considered negligible for a salinity increase of less than 1 (∼ 0.7 and ∼ 0.5 increase for both salt additions in our study, respectively; Czerny et al., 2013a). The difference in salinity in the mesocosms from the Pacific was less than 1 throughout the study, so we believe salinity was not a stressor to the system (Bach et al., 2020a). Second, multiple factors may introduce variabilities into air–sea gas exchange in the mesocosms. With the water surface sheltered from direct wind forcing by the 2 m high plastic walls of the mesocosm bag, air–sea gas exchange could be very low (but not zero due to other energy inputs such as thermal convection and surrounding wave movements; Czerny et al., 2013b). On the other hand, the extensive daily sampling that actively perturbs the mesocosm surface may enhance the gas exchange (Czerny et al., 2013b). With a high variability in DOC data and a lack of direct measurements of CO₂ gas exchange, it is impossible to estimate a reasonable C budget for our dataset (Czerny et al., 2013b; Boxhammer et al., 2018). Other complications include heterogenous initial conditions in the mesocosms and light limitation due to self-shading (please see Bach et al., 2020a, for an extensive discussion). It would have been ideal to have had control mesocosms that were treated the same way except in terms of the OMZ water addition to rule out the effects induced by mesocosm manipulations and changes in hydrodynamics. Not having this was a compromise to ensure enough replicate numbers for both treatments despite the enormous cost of mesocosm experimentation. Nevertheless, a previous study has examined impacts of different mixing techniques in outdoor mesocosms and found no effects on phytoplankton biomass and minor effects on phytoplankton and zooplankton community composition (Striebel et al., 2013). In our study, various measures were also taken to minimize the mixing (brine additions and slow casting of CTD).

4.3 Temporal changes in carbonate chemistry in the coastal Pacific near Callao

According to estimations by Takahashi et al. (2009) of global air–sea CO₂ fluxes, our study site in the equatorial Pacific (14^∘ N–14^∘ S) is a major source of CO₂ to the atmosphere. Our near-coastal location showed high pCO₂ levels over the study period (with an average of 1078 µatm), with a sea-to-air CO₂ flux of 4.2–14.0 mmol C m⁻² d⁻¹ (Fig. 5). Compared to the criterion of high CO₂ fluxes (5 mmol C m⁻² d⁻¹ or more) as proposed by Paulmier et al. (2008), our study site was a strong CO₂ source to the atmosphere most of the time. These results of air–sea CO₂ fluxes were slightly higher than observations by Friederich et al. (2008) along the coast of Peru in February 2004–2006 (0.85–4.54 mol C m⁻² yr⁻¹, spatially averaged for 5–15^∘ S along the coast of Peru). This is not surprising because Friederich et al. (2008) averaged the air–sea CO₂ fluxes for 0–200 km from the shore, where much lower pCO₂ levels were observed offshore (< 600 µatm) compared to our nearshore study site. The decline in pCO₂ with increasing distance from shore was driven by biological uptake and outgassing to the atmosphere (Friederich et al., 2008; Loucaides et al., 2012). However, when compared to the magnitude of DIC drawdown triggered by upwelling events in the mesocosms, the flux of CO₂ to the atmosphere was insignificant. Assuming a 10 m mixed layer in the Pacific with a DIC concentration of 2200 µmol kg⁻¹, the DIC content below a 1 m² surface area would be ∼ 22 mol m⁻². With an upper bound outgassing of 14.2 mmol C m⁻² d⁻¹ over 10 d (day 13–24), the loss of CO₂ would only be 0.142 mol m⁻². On the other hand, the average DIC drawdown of 118.2 µmol kg⁻¹ in the Very low DIN and 160.3 µmol kg⁻¹ in the Low DIN mesocosms (M3 and M4 excluded) during this period accounts for 1.18 and 1.60 mol m⁻², respectively, over the same water column. This shows that biological processes drawing down CO₂ are stronger than loss by air–sea gas exchange.

During our study, we experienced a coastal El Niño event which was the strongest on record (compared to those recorded in 1891 and 1925) and induced rapid sea surface warming of ∼ 1.5 ^∘C and enhanced stratification (Garreaud, 2018). Previous investigations showed that the impact of reduced upwelling on CO₂ fluxes is pronounced for upwelling areas (Feely et al., 1999, 2002). A decline in upwelling of CO₂-enriched OMZ water results in a decrease in sea-to-air CO₂ fluxes. For example, during the 1991–1994 El Niño year, a total reduction in CO₂ fluxes to the atmosphere was reported for the equatorial Pacific. They were only 30 %–80 % of that of a non-El Niño year (Feely et al., 1999, 2002). This is likely to be the case for our study location. Most studies have investigated air–sea CO₂ fluxes at larger timescales and regional scales (Feely et al., 1999; Friederich et al., 2008; Takahashi et al., 2009). Therefore, it is difficult to conclude the magnitude of the coastal El Niño influence on the local CO₂ fluxes in our study by comparing our results with previous observations. Nevertheless, our observations can serve as the first evidence of carbonate chemistry dynamics in the coastal Peruvian upwelling system during a coastal El Niño event. Observations of sea surface carbonate chemistry with a high temporal resolution (every second day) in nearshore waters are scarce, as they are rarely covered by typical research expeditions in the open ocean (Takahashi et al., 2009; Franco et al., 2014), especially during such an extremely rare coastal El Niño event. Comparisons of our data with previous or future observations may enhance our understanding of how inorganic carbon cycling interacts with extreme climate events in upwelling systems.

CO₂-enriched OMZ water has been occasionally reported to be under-saturated with respect to aragonite (Feely et al., 2008; Fassbender et al., 2011). In our study, calcite under-saturation did not occur in the mesocosms or in the Pacific (Fig. 4). Aragonite under-saturation, however, was observed below the surface (10–17 m for day 3–28, 12.5–17 m for day 29–50) of the Pacific at the start of the experiment (Fig. 4d), when pCO₂ was the highest (pCO₂>1100 µatm, Fig. 3c). Aragonite under-saturation was also observed in the two deep water masses collected at deeper depths (30 and 70 m) in the Pacific (Table 4). Throughout the study period, the aragonite saturation state fluctuated close to around 1 below the surface (Fig. 4d). Considering the water column we sampled in the Pacific still belonged to the upper surface ocean, we could expect deeper and more CO₂-enriched water in the underlying OMZ to be most likely under-saturated with respect to calcite and aragonite. Hence, our observations of aragonite under-saturation in the Pacific suggest a potential risk of dissolution for marine calcifiers in response to the ongoing intensification and expansion of acidified OMZ water (Comeau et al., 2009; Lischka et al., 2011; Maas et al., 2012).

5 Conclusions

Our observations in the mesocosms revealed that, following the addition of two OMZ water masses with different nutrient signatures, there was a higher drawdown of DIC in response to slightly more DIN input from the OMZ water addition but no difference in the buildup of POC and chl a (Figs. 2a and c, 3d). The timing of the first phytoplankton bloom was consistent with a shift from a diatom-dominated community to A. sanguinea dominance in most mesocosms, indicating that A. sanguinea was a key factor driving the changes in carbonate chemistry under N-limited conditions. A second phytoplankton bloom was triggered by defecation of Inca terns, which eased the N limitation in the mesocosms (Fig. 2c). These findings provide improved insights into the links between upwelling-induced N limitation, phytoplankton community shifts, and carbonate chemistry dynamics in the Peruvian upwelling system.

The surrounding Pacific waters at the study site were characterized by constantly high pCO₂ levels (with an average of 1078.1 µatm). Most CO₂ flux estimates have been conducted in the open ocean, and few studies have surveyed coastal regions (Takahashi et al., 2009; Franco et al., 2014). Our study site was a strong CO₂ source to the atmosphere most of the time (4.2–14.2 mmol C m⁻² d⁻¹), despite a rare coastal El Niño event. However, evidence from our mesocosm experiment suggests biological responses that draw down DIC can quickly turn a CO₂ source into a sink in the upwelling system. The influence of the co-occurring coastal El Niño event on the local CO₂ fluxes remains unclear. Nevertheless, future carbonate chemistry fluctuations are expected to be enhanced by expanding and intensifying ocean deoxygenation, as well as by reducing buffer factors (Schulz et al., 2019). Hence, it is essential to improve our understanding of the mechanisms driving the inorganic carbon cycling in upwelling systems. As a unique dataset that characterized nearshore carbonate chemistry with a high temporal resolution during a rare coastal El Niño event, our study gives important insights into the carbonate chemistry responses to extreme climate events in the Peruvian upwelling system.