Air–Sea Fluxes of Greenhouse Gases and Oxygen in the Northern Benguela Current Region During Upwelling Events

. Ground-based atmospheric observations of CO 2 , δ (O 2 /N 2 ), N 2 O, and CH 4 were used to make estimates of the air– sea ﬂuxes of these species from the Lüderitz and Walvis Bay upwelling cells in the northern Benguela region, during upwelling events. Average ﬂux densities ( ± 1 σ ) were 0.65 ± 0.4 µ mol m − 2 s − 1 for CO 2 , − 5 . 1 ± 2.5 µ mol m − 2 s − 1 for O 2 (as APO), 0.61 ± 0.5 nmol m − 2 s − 1 for N 2 O, and 4.8 ± 6.3 nmol m − 2 s − 1 for CH 4 . A comparison of our top-down ( i.e. , inferred from atmospheric anomalies) ﬂux estimates with shipboard-based measurements showed that the two approaches agreed within 5 ± 55% on average. During the study, upwelling events were sources of CO 2 , N 2 O, and CH 4 to the atmosphere. N 2 O ﬂuxes were fairly low, in accordance with previous work suggesting that the evasion of this gas from the Benguela is smaller than for other Eastern Boundary Upwelling Systems (EBUS). Conversely, CH 4 release was quite high for the marine environment, a result that supports studies that indicated a large sedimentary source of CH 4 in the Walvis Bay area. These results demonstrate the suitability of atmospheric time series for characterizing the temporal variability of upwelling events and their inﬂuence on the overall emissions from the

and Training Centre. A full description of the measurement system is given in Morgan et al. (2015). In brief, observations were made at 21 m height above ground level with a Picarro ESP-1000 cavity ringdown spectrometer (CRDS) for CO 2 and CH 4 , a Los Gatos N 2 O/CO-23d cavity-enhanced absorption spectrometer for N 2 O and CO, and an Oxzilla FC-II dual absolute and differential oxygen analyzer for δ(O 2 /N 2 ).
While the Picarro and Los Gatos analyzers measure continuously (data are recorded approximately 0.5 and 1 Hz, respec-5 tively), the Oxzilla analyzer measures the difference between sample air and a reference tank, at a data rate of 0.01 Hz. Calibration of the instruments was done through four working secondary standards and instrument performance was periodically checked with "target" cylinders (i.e., tanks of known mole fraction which are regularly remeasured). Reference gases were comprised of dry ambient air and stored in 50 L aluminum cylinders. Calibrations were run every 123 hours for the Picarro and Los Gatos analyzers, and every 71 hours for the Oxzilla. All measurements were tied to primary standards on the following The average uncertainty for each species during the study period was 0.028 ppm for CO 2 , 6.5 per meg for δ(O 2 /N 2 ), 0.21 ppb for N 2 O, 0.17 ppb for CH 4 , and 0.15 ppb for CO. CO 2 , N 2 O, CH 4 , and CO are all expressed as dry air mole fractions, in ppb or ppm (1 ppm = 1 µmol mol −1 ; 1 ppb = 1 nmol mol −1 ). 15 By convention, changes atmospheric oxygen were quantified as the O 2 /N 2 ratio relative to a standard, in per meg units (Keeling and Shertz, 1992): In order to isolate the influence of air-sea exchanges on O 2 /N 2 , we have employed the use of a data-derived tracer, known as atmospheric potential oxygen (APO), which masks variations of O 2 /N 2 that are due to terrestrial biosphere exchange (Stephens 20 et al., 1998). Variations in APO are thus primarily due to fossil fuel burning and air-sea gas exchange of O 2 . APO is defined as: (TMI) daily optimally interpolated SST product was selected. The major advantage of this instrument is its ability to measure SST through clouds, which are considerable over the coast, as the TMI measures frequencies in the microwave region (4-11 GHz). The drawback of this dataset is that there is no data within 100 km of the coast. Large upwelling features, however, extend much farther out than this and are readily seen by TMI (Goubanova et al., 2013). Even with this loss of near-shore data, the data coverage is still superior to that of optical sensors like the Moderate-Resolution Imaging Spectroradiometer (MODIS).

5
The TRMM SST data presented in this work is deseasonalized by subtracting a second harmonic fit to the data, as it showed a strong seasonal cycle which masked some of the intraseasonal variability when plotted as a time series.
Wind speed data for the South Atlantic was also derived from the TMI instrument on the TRMM satellite. This dataset is a level 3 product which gives the 10 m wind speed over marine areas within the sensor's field of view. The 18.7 GHz channel data product was selected. Like the SST data, a major drawback of this dataset is the absence of data within 100 km of land. 10 Two surface chlorophyll products were used, both level 3 binned products that combined data from multiple satellites, accessed through ESA GlobColour website (http://www.globcolour.info/). The first is denoted CHL1-GSM; this dataset is a merged product of two different sensors (during the time period considered), MODIS and the Visible Infrared Imaging Radiometer Suite (VIIRS). The data is merged using the Garver, Siegel, Maritorena (GSM) model, which blends the normalized water-leaving radiances instead of the end product (chl-a concentrations) (Maritorena and Siegel, 2005). The second product 15 is denoted CHL1-AVW; these data are merged using a weighted average method (AVW). Like CHL1-GSM, it combines data from MODIS and VIIRS for the time frame considered.

Atmospheric Back-Trajectories
Back-trajectories, which trace the path of a particle from a receptor point backwards in time, provide a history of recent atmospheric transport. Back-trajectories were simulated with the the HYbrid Single Particle Lagrangian Integrated Trajectory 20 (HYSPLIT) 4 model (Draxler andHess, 1997, 1998) from NDAO for 120 hours, for the entire time series. A new trajectory was calculated every six hours, starting at 0:00 UTC, i.e., at 1:00, 7:00, 13:00, and 19:00 local time. The model was run with a spatial resolution of 1 • × 1 • and a temporal resolution of 1 hour. A vertical cut-off of 10 km a.g.l. was used. Meteorological fields were obtained from the National Center for Environmental Prediction (NCEP) Global Data Assimilation System (GDAS).

25
A subset of the coastal region was selected to represent the Lüderitz and Walvis Bay upwelling cells.  (Lutjeharms and Meeuwis, 1987;30 Veitch et al., 2009), and where upwelling was downwind of the station during upwelling events. These criteria were considered desirable because they would provide the best opportunities for relating atmospheric anomalies to upwelling events. These determinations were based on analysis of our SST dataset, and HYSPLIT back-trajectories.
Upwelling events were identified based on SST anomalies and 10 m wind speed anomalies. An event was determined to occur if the average deseasonalized SST of the domain was 0.5 • C or lower than a smoothed, second-degree polynomial fit to the entire time series, and the average 10 m wind speed of the study area was 2.5 m s −1 above a smoothed, second-degree polynomial fit to the wind data. These thresholds were determined through visual inspection of maps and time series of SST and wind speed data, and are specific to the domain chosen, since the data considered were averages of the entire area. One 5 standard deviation of the SST anomaly was 0.92 • C, and one standard deviation of the wind speed anomaly was 3.1 m s −1 .
Since the resolution of the SST and wind speed time series was daily, all higher resolution data falling within a day during which an upwelling event occurred was similarly flagged.
Atmospheric anomalies associated with upwelling events were quantified against a baseline determined through a secondharmonic fit that was generated iteratively to all data, excluding all points that lay outside 1 standard deviation from the curve 10 for the subsequent iteration. Small adjustments were made to the anomaly to raise or lower the curve in cases where it did not intersect with background values within ± 5 days of the event. The NDAO data was filtered by wind (wind speeds greater than 2 m s −1 and wind direction within the NNW-SSW sector), as well as by back-trajectory to exclude anomalies which were not associated with marine air masses. For the latter, trajectories could not reside for more than 36 hours of the total 120 hours over land, and could not travel more than 50 km inland past NDAO. A final data selection step was to exclude days that had 15 CO values greater than 15 ppb above baseline, to remove any time periods that may have been influenced by biomass burning.
The filtered time series, with terrestrial influences removed, is shown in Figure 2.

Top-Down Air-Sea Flux Estimates
In order to estimate the surface flux associated with atmospheric anomalies due to upwelling events, the approach of Lueker et al. (2003) was adopted. A simple model was employed to describe the change in the concentration of a species within a 20 well-mixed column of air as it moves over a source region (Jacob, 1999): Here ∆C is the concentration of the species of interest, in mol m −3 , expressed as an anomaly against the background. ∆C is a function of x, which is the distance to the observation point from the flux region. L is the point at which the column (with height h, in m) leaves this region, characterized by a constant flux, F , in mol m −2 hr −1 , and a constant wind speed, U , in m 25 hr −1 . After the column leaves the flux region (x ≥ L, i.e. when the coast is reached), the loss of ∆C from its peak at L (∆C L ) is governed by dilution due to mixing of background air. This requires the dilution rate constant, q, in hr −1 , to be known.
We solved the equation for F by taking the other variables as follows: ∆C was determined from our atmospheric record (see Section 2.4), wind speeds (U ) were obtained from satellite data (Section 2.2), and h was taken as the average height of the planetary boundary layer (PBL) for the Lüderitz/Walvis Bay domain over the course of any given event. PBL data was acquired 30 from the European Centre for Medium-Range Weather Forecasting's (ECMWF) ERA-Interim dataset (Dee et al., 1979). x was taken as the length traveled along a back-trajectory from NDAO to the area affected by upwelling.
The dilution rate constant, q, was estimated by comparing measurements of CO 2 and CH 4 made during the RV Meteor cruise M99 (see Section 2.7). Back-trajectories from NDAO were matched to the closest ship location at the appropriate time.
Any back-trajectory points that were within 100 km of the ship-both horizontally and vertically-within the space of 1 hour 5 were identified, and a dilution rate constant was calculated for both CO 2 and CH 4 , as (Price et al., 2004): Where C M 99 and C N DAO are the mole fractions of CO 2 or CH 4 , as measured in situ on the Meteor or at NDAO, and C b is taken from the fit to the baseline of the NDAO time series for either species, as described in Section 2.4. t is the travel time along the back-trajectory, in hr. While this assumes that there is no flux the gas between the vessel and the station, we excluded 10 instances when the anomaly was heavily influenced by surface fluxes by filtering for poor agreement between CO 2 and CH 4 .
Fluxes of these gases are not well correlated, and can go in opposite directions for summertime air-sea fluxes, so we assume that it is unlikely that both values of q would be biased by the same amount. From 32 values of q we excluded all but 13 for poor agreement (excluding determinations with |q CO2 − q CH4 | > 0.01). The average (± 1σ) was taken to arrive at a single value, 0.011 ± 0.006 hr −1 . Since this is still a rough approximation, we also present estimates of the estimated flux, F , with a 15 value of q tuned to provide the best agreement with an in situ flux estimate (see Section 2.8).

Estimated Top-Down Flux Density Uncertainty
Uncertainties were propagated in quadrature, e.g.: for operations involving addition or subtraction (x = a + b), and Atmospheric measurements of CO 2 and CH 4 were made with a CRDS analyzer (model G1301, Picarro Inc, Santa Clara, 5 CA, USA) located in the atmospheric chemistry laboratory. The instrument's internal pump was used to draw air through a 7 m length of 1/4" SERTOflex tubing, at a flow rate 150 mL min −1 . Inlets identical to those used at NDAO  were placed on the starboard railing of the 6 th superstructure deck, just above the atmospheric chemistry lab, at a total height of ∼21 m above sea level. A second-order, instrument-specific water correction was performed in lieu of physical or chemical drying, identical to the method described in Morgan et al. (2015). As the instrument's pressure control seemed to 10 be affected by strong vessel motion, measurements were excluded if the cavity pressure deviated by more than 0.04 torr from the setpoint of 140 torr. Calibrations were conducted on average every three days and target measurements were made once per day. Reference gases were calibrated at MPI-BGC GASLAB. The uncertainty, derived from the target measurements as at NDAO, was determined to be ± 0.03 ppm for CO 2 and ± 0.43 ppb for CH 4 . The dataset was filtered for contamination by the ship's exhaust using the relative wind direction data from the ship's meteorological instrumentation. 15 Flasks for δ(O 2 /N 2 ) were taken in triplicate and connected in series downstream of a pump. The pump body and valve plates were aluminum, and the structured diaphragms were made of PTFE. When in use the flow rate (3.2 L min −1 ) was higher than the in situ analyzer flow rates (100-200 mL min −1 ). Air was dried with a magnesium perchlorate. During sampling, the line was flushed for 5 minutes before any air was directed to the flasks, then a bypass valve was opened and the flasks were flushed for an additional 15 minutes before they were sealed again. After closure, the pressure of the flask was about 1.6 bars. Flasks 20 were analyzed with a mass spectrometer at MPI-BGC (Brand, 2005). Storage times were less than two months.
Dissolved gas measurements were carried out by means of an autonomous setup for along-track measurements of CO 2 , N 2 O and CO, which combined the analytical approaches from Pierrot et al. (2009)  The estimated uncertainty of the dissolved CO 2 measurements was ± 2 µatm; of dissolved O 2 measurements, ± 4 µmol The uncertainty of the atmospheric measurements of N 2 O was ± 0.9 ppb.

Shipboard Air-Sea Flux Density Estimates
In situ oceanographic and meteorological data were taken from the Meteor's instrumentation. In order to determine the total dissolved inorganic carbon (DIC) content of surface waters, total alkalinity was estimated from temperature, salinity, and dissolved O 2 data, using the locally interpolated alkalinity regression (LIAR) version 2 (Carter et al., 2018). The dissociation 30 constants of carbonic acid were also determined from temperature and salinity using the formulations of Millero et al. (2006).
The total DIC content was then estimated from the total alkalinity and f CO 2 . Meteorological data (air temperature, barometric pressure, wind speed, etc.) was observed at a height of 37 m above sea level. The absolute wind speed measured on the Meteor was converted to U 10 through the relationship (Justus and Mikhail, 1976): where U meas is the wind speed in m s −1 , measured at some height z meas , in m.

5
Marine surface flux densities of CO 2 , O 2 , and N 2 O were estimated for the vessel location throughout M99 from shipboard measurements of atmospheric dry air mole fractions and dissolved aqueous concentrations, according to Equation 9. While flux densities of CO were also determined, they are not discussed in this manuscript as the shipboard-measured flux was too small to be detected as an atmospheric anomaly at NDAO (contrary to what was reported in Morgan (2015)). In the case of O 2 , the atmospheric concentration was measured only sporadically with flask samples, so it was taken as a static value, viz. the 10 mole fraction of O 2 in standard dry air, 0.2093 (Aoki et al., 2019). The in situ aqueous solubility of O 2 was calculated using the equations of García and Gordon (1992), of N 2 O using those in Weiss and Price (1980), and of CO 2 using Weiss (1974).
Sea-to-air fluxes (net evasion) are positive.
The flux density (F , in units of mol m −2 s −1 ) is typically determined according to (Garbe et al., 2014): 15 where k w is the gas transfer (or piston) velocity, in m s −1 , C w is the dissolved concentration in the water phase (mol m −3 ), and C a is the concentration of the species in the air in the same units. The formulation can also be altered to accommodate units of partial pressure in both phases. The expression αC a gives the dissolved concentration in the water phase directly at the surface; α is the Ostwald solubility coefficient: the reciprocal of the dimensionless air-water partition coefficient (K AW ) for some temperature, T , and salinity, S (Mackay and Shiu, 1981). 20 As there is no definitive k w -U 10 parameterization, fluxes were computed with two different parameterizations of k w : that of In these equations, U 10 is the wind speed at 10 m's height, and Sc is the Schmidt number of a particular gas at in situ conditions 25 (Jähne et al., 1987;Wanninkhof, 1992). The Schmidt number is dimensionless, and U 10 is in units of m s −1 , producing k w in units of cm hr −1 . The Schmidt number is scaled to the reference conditions of the parameterization; Schmidt numbers were calculated at in situ conditions by dividing the kinematic viscosity of seawater by the diffusivity of a given gas, using an Eyring-style equation (Eyring, 1936;Jähne et al., 1987): where R is the ideal gas constant, T is temperature, and E a is the activation energy for diffusion in water. A and E a are 5 determined through fits to experimental data. For O 2 these were taken from Ferrell and Himmelblau (1967), for CO 2 from Jähne et al. (1987), and for N 2 O from Bange et al. (2001). A salinity correction of 4.9% per 35.5 psu was then applied, this number being the average decrease in diffusivity seen by Jähne et al. (1987) for He and H 2 in an experiment involving artificial seawater. The kinematic viscosity of seawater at in situ conditions was determined by calculating the dynamic viscosity after Laliberté (2007) and the density after Millero and Huang (2009). 10 3 Results and Discussion

Upwelling Events
Out of 741 days, 219 days met the SST criteria, but only 173 of these met both the SST and WS criteria, representing 23% of the two-year study period. From these 173 days with upwelling events, 102 had atmospheric conditions that allowed for the detection of an anomaly in the atmospheric time series. 21 of excluded events were filtered based on HYSPLIT back-15 trajectories; the rest were excluded based on local meteorology or elevated CO. This was a conservative approach, since the filtering criteria in part relied on HYSPLIT model results, which could misrepresent the transport pathway and lead to the exclusion of an upwelling event that in fact was favorable for carrying an air mass influenced by upwelling to the station. Despite the greater prevalence of equatorward winds during austral summer, the distribution of events displayed little seasonality (Figure 3), reflecting the fact that upwelling is a short-term, intraseasonal phenomenon, forced by specific atmospheric conditions 20 (Risien et al., 2004;Goubanova et al., 2013). While upwelling in this area is perennial, seasonality is seen in the intensity of upwelling due to the annual migration of the South Atlantic Anticyclone, with a minimum in austral winter (Hagen et al., 2001;Hardman-Mountford et al., 2003;Peard, 2007;Veitch et al., 2009;Hutchings et al., 2009).
An example of an upwelling event is given in Figures 1, 4, and 5. On November 27 th , 2013, high winds resulted in the creation of a very large pool of colder water on the surface that persisted for four days, until winds relaxed and upwelling 25 temporarily ceased until the 4 th of December, when SST dropped again. During these two low SST pulses, chl-a values were higher. A change in the background values of atmospheric potential oxygen (APO), N 2 O, and CH 4 was seen, with a smaller anomaly for CO 2 . The largest anomalies for each species came when high winds were from the coastal sector on November 28 th .
If the area of high flux was close to the coast, anomalies could arrive within a few hours at NDAO, with the development of 30 the afternoon sea breeze. If the region of flux was closer to Lüderitz, the arrival time could be delayed by as much as 50 hours, depending on the wind speed and the degree of meandering of the air mass. Back-trajectories implied that despite the high wind speeds usually seen in this coastal zone, significant travel time (1 to 2 days) could be expected for most air masses of interest.
As a result, the marine surface flux associated with an atmospheric anomaly could only be said to have taken place with 50 hours of its detection. The magnitudes of the average atmospheric anomalies and their corresponding flux density estimates are given in Table 1.  (1977) and from Monteiro et al. (2006) to estimate flux densities of 0.03-8.7 nmol m −2 s −1 , the upper end-member of which compares favorably with the rates found in this study for the same ocean region. Though the latter study observed concentrations at a mooring near Walvis Bay as high as 10 µM, these data were from an uncalibrated probe and are not usable for a direct comparison to other campaigns. In other EBUS, reported flux densities are much lower.

15
Only CO 2 and CH 4 were measured continuously in the atmosphere on M99. During the first days of the cruise, before the main upwelling event, a synoptic event brought elevated mixing ratios of CO 2 and CH 4 offshore ( Figure 6). The shipboard measurements of CO 2 showed sporadic enhancements relative to the smoothed station background. These enhancements coincided with the regions of higher flux closer to the coast encountered under upwelling conditions (Figure 6 and 7). In contrast, apart from the initial synoptic event, CH 4 was consistently at background levels, usually below the value seen at NDAO. The top-down estimates agreed reasonably with the corresponding mean shipboard estimates for a selected 7.5 hour period that coincided with the largest flux densities ( Figure 8 and Table 2), but overestimates the flux density relative to the mean shipboard estimate for the flux event. Agreement was best for CO 2 , but was as poor as a factor of 3 for N 2 O, though given the 15 simplicity of the top-down estimate, this seems fairly reasonable.
The choice of gas transfer velocity parameterization (k w ) was significant in determining the magnitude of the agreement, and the degree of this correspondence between the top-down and shipboard gas transfer velocity-specific estimate varied between species. The cubic relationship of k McG01 produced higher fluxes during the M99 event, and was in better agreement with the top-down estimate. While this provides a measure of confidence in the top-down flux density estimates, it should be noted 20 that the neither of the estimated uncertainties for the top-down or bottom-up approaches account for errors incurred by the simplifying assumptions within their formulations, and these comparisons should not be considered a "calibration" of the flux determination or k w parameterization. We suspect that our top-down estimates are overestimates. Cubic relationships with wind speed, derived mostly from eddy covariance data, have not been disproved per se, but after reevaluation of field data (Landwehr et al., 2014), the weight of evidence seems to favor a quadratic relationship (Roobaert et al., 2018), and hence we 25 prefer k W14 for comparison.

Model Uncertainty and Limitations
The model makes many simplifying assumptions, which we have attempted to incorporate into our uncertainty assumptions.
The assumption of a constant boundary layer height, h, is of course not realistic, but we have attempted to account for natural variations in this parameter with our uncertainty range. The uncertainty in the distance traveled depends on the accuracy of the 30 back-trajectories, although major errors in transport (e.g. a continental air mass causing the anomaly) are likely to be excluded due to the filtering criteria of the atmospheric record. We also note that the model is sensitive to the dilution rate constant, q, and that this is the input we know the least well. Previous studies suggest that our in situ determination is reasonable; (Price et al., 2004) used multiple species and a model approach to arrive at a mean q of 0.010 ± 0.004 hr −1 , which is nearly identical to our value of 0.011 ± 0.006 hr −1 , although the spatial and temporal scale they considered was larger. Dillon et al. (2002) observed rates that were higher for a Sacremento pollution plume, ca. 0.2 hr −1 .
While we cannot determine for specific events when the simplifying assumptions of the model are problematic, we assume that variations in the dilution rate, the boundary layer height, and errors in transport generally average out, and that the estimated 5 average flux from upwelling events is representative. If we "tune" the model to match the M99 event to the k W 14 fluxes for each species, we get a lower value of 0.0058 ± 0.007 hr −1 . Using this value of q we obtain the tuned fluxes reported in Table   2. Though these fluxes have a higher uncertainty associated with them, due to the higher uncertainty in q, we recommend these values, as the top-down and bottom-up comparison suggests that the model overestimates the fluxes.

10
Correlation slopes of atmospheric species can provide further confidence and insight into source processes, if there is an underlying biogeochemical relationship. The well-known inverse relationship between N 2 O and O 2 in the ocean, for instance, is a result of organic matter decomposition and nitrification (Cohen and Gordon, 1979;Nevison et al., 2003;Naqvi et al., 2010;Frame et al., 2014). The tight coupling between N 2 O and O 2 seen in surface concentrations during M99 is preserved during air-sea gas exchange, as these gases behave similarly (Figure 9). The observed linear correlation in surface waters between 15 these two species is also an indication that intense denitrification was not occurring in the OMZ at the time of sampling, which if it were taking place, would lead to a breakdown in this relationship at low O 2 (Cohen and Gordon, 1979). The approximate molar ratio of these two species, −0.8 × 10 −4 to −1 × 10 −4 (N 2 O:O 2 ; mol mol −1 ), is the same observed by Lueker et al. (2003) for the Trinidad Head region, and appears to be a globally consistent value (Nevison et al., 2005;Manizza et al., 2012).
This linear regression slope is often expressed in terms of the excess N 2 O (measured N 2 O minus N 2 O at saturation) and 20 apparent oxygen utilization (saturation minus observed), ∆N 2 O-AOU in nmol µmol −1 . Quantified in this way, it has been used as an estimate of the yield of N 2 O as a function of the amount of oxygen consumed . However, the relationship is not strictly linear, since N 2 O production is enhanced at low oxygen levels Naqvi et al., 2010;Trimmer et al., 2016). ∆N 2 O-AOU is also sensitive to mixing, as N 2 O production rates vary widely in the ocean, meaning that the mixing of water masses with different compositions can overwhelm the in situ production signal 25 (Suntharalingam and Sarmiento, 2000;Nevison et al., 2003). The ∆N 2 O-AOU for M99 was 0.088 ± 0.003 nmol µmol −1 , with an intercept of 1.6 ± 0.04 nmol, and an R 2 of 0.69. This is a low value, nearly identical to results from the eastern basin of the sub-tropical North Atlantic, where South Atlantic Central Water is found (Cohen and Gordon, 1979;Suntharalingam and Sarmiento, 2000;Walter et al., 2006), and the eastern equatorial Atlantic (Arévalo-Martínez et al., 2017). Frame et al. (2014) showed that most of the N 2 O in the Benguela Current region is produced in the water column and in the sediment by nitrifier 30 denitrification (Frame et al., 2014). Nevertheless, a substantial portion of this N 2 O remains at depth (with a concentration maximum at 200-400 m) and is advected away from the region, without the chance for atmospheric release (Gutknecht et al., 2013b, a;Frame et al., 2014). Hence, the low ∆N 2 O-AOU value found for the Meteor cruise probably reflects both physical and biogeochemical dynamics.
In the case of variations of O 2 and CO 2 , the stoichiometry of surface waters is not preserved after air-sea exchange, as the majority of carbon is speciated in the carbonate system, and only the portion that remains as dissolved CO 2 is available for air-sea gas exchange. This leads to a change in the ratio, for instance, from 0.58 ± 0.03 in surface waters to −6.53 ± 0.42 in the atmosphere, for the upwelling event encountered during the R/V Meteor cruise M99 (Figure 8). These two species can become decoupled through the influences of changing solubility, which would drive evasion of both gases, and net biological 5 production, which would drive evasion of O 2 and invasion of CO 2 . These complicating influences are the likely reason for the poorer correlation seen between these two species when compared with N 2 O and O 2 .
While the top-down flux density estimates cannot be confirmed with shipboard estimates of flux densities for CH 4 , it is worth considering that concentrations of methane in bottom waters on the Namibian shelf are likely the highest ever measured in an open coastal system. Values as high as 475 µM in the bottom waters and greater than 5,000 µM in sediment porewaters 10 have been observed (Scranton and Farrington, 1977;Monteiro et al., 2006;Brüchert et al., 2009;Naqvi et al., 2010). In the water column, the concentration maxima is usually at the seabed or in bottom water, but it is variable and can even occur at the surface (1 m) (Brüchert et al., 2009). Dissolved methane concentrations are tightly coupled with O 2 and show considerable variability, with elevated concentrations being triggered by episodes of hypoxia (Monteiro et al., 2006;Brüchert et al., 2009).
The pulse-like nature of CH 4 in the Benguela means that the full range of dynamics cannot be captured with a campaign-based 15 sampling approach (Brüchert et al., 2009). What is clear is that there is a tremendous amount of methane production at depth, but that the source is variable in strength (Emeis et al., 2004;Brüchert et al., 2009). In light of the fact that large pockets of free methane gas are contained in the sediment in the Walvis Bay region, as well as the existence of craters and pockmarks on the seafloor, combined with observation of bubble streams from the seabed, suggest a mechanism by which methane produced in sediments can be abruptly transported to the surface and hence, avoid oxidation (Emeis et al., 2004;Brüchert et al., 2006Brüchert et al., , 20 2009). Consequently, the Benguela Current is a suspected source of CH 4 to the atmosphere, but the amount is ill-constrained (Naqvi et al., 2010); Emeis et al. (2018) recently estimated the total annual emission for the entire system to be less than 0.17 Tg CH 4 yr −1 .
Interestingly, methane was not well-correlated with either APO or CO 2 in the atmosphere during all upwelling events, suggesting a spatial decoupling (since a cross-correlation analysis indicated this was not a result of lag/temporal decoupling) 25 between methane and these two species. While background observations of CH 4 were generally well-correlated with CO 2 and O 2 at NDAO, only some upwelling events showed such coupling; it seems there is a general relationship between methane and oxygen, but it is not consistent and is occasionally non-existent. Unfortunately, since there are still very few measurements of water-column CH 4 in the Benguela, a full explanation of the methane source remains elusive. From these atmospheric trends it can only be deduced that there is some separate biogeochemical influence on methane that is not exerted over CO 2 , O 2 , 30 or N 2 O. This observation is arguably consistent with the concept of a dominant sedimentary source of methane that is more localized within the inshore mud belt, where high POC fluxes have created a thick layer of diatomaceous ooze containing free methane gas pockets (Emeis et al., 2004;Brüchert et al., 2006;van der Plas et al., 2007;Brüchert et al., 2009).

Regional Flux Context and Future Steps
To put our estimated flux densities into the context of regional greenhouse gas budgets, it is necessary to have an approximation of the area involved. For the sake of discussion only, we can very roughly estimate the total annual flux of each species from the Walvis Bay/Lüderitz domain by assuming that grid cells with wind speeds above 3.5 m s −1 and a deseasonalized SST below −0.1 • C were grid cells with upwelling, that upwelling extends to the coast, and that each grid cell with upwelling had 5 the same flux density, i.e. the top-down estimated flux density for a given upwelling event. This approach is very crude, and would result in an underestimate, since transport conditions were not always conducive to observing an upwelling event. This would result in annual fluxes of ≥ 206 ± 151 Gmol yr −1 for CO 2 , ≥ −1.6 ± 1.1 Tmol yr −1 for O 2 , ≥ 272 ± 248 Mmol yr −1 for N 2 O, and ≥ 3.3 ± 2.8 Gmol yr −1 for CH 4 from upwelling events.
Such an estimate for CO 2 is substantial when compared to the net flux that Laruelle et al. (2014) estimated of −424.9 10 Gmol yr −1 for the entire Benguela region, or the −141.5 Gmol yr −1 found by Gregor and Monteiro (2013) for the southern Benguela. A high flux of CO 2 is possible due to the higher wind speeds and the more remineralized character of the South Atlantic Central Water that upwells at Lüderitz. Using inverse methods, Gruber et al. (2001) constrained the net flux of oxygen for the temperate South Atlantic (an area of 1.5 × 10 7 km 2 ) to be 15.5 Tmol yr −1 , which suggests that upwelling events in the Lüderitz region could be regionally significant. An annual source of this magnitude for N 2 O is modest when compared to 15 the budgetary calculations of other coastal regions. As the global coastal upwelling source of N 2 O is estimated to be 7,140 Mmol yr −1 (Nevison et al., 2004), the emissions for the Lüderitz/Walvis Bay region would represent 3.2% of the area but only 1.7% of these emissions. For CH 4 this approximation is two to three times higher than the net evasion from the Arabian Sea (Bange et al., 1998), and 10 to 20 times greater than the annual release of CH 4 from the entire Mauritanian upwelling system (Kock et al., 2008;Brown et al., 2014), but modest when compared to the estimate of  for the entire Bengulea A full regional top-down quantification of the annual fluxes-as opposed to the event-based averages we present here-of the Benguela from land-based atmospheric measurements would require additional efforts, such as a Bayesian synthesis inversion of the data, where a prior surface flux field is adjusted to best match the atmospheric observations. This approach relies on an accurate knowledge of atmospheric transport, and ideally a dense network of measurement sites. Currently there are only two 25 such sites in the Benguela region known to the authors, NDAO and Cape Point Observatory (CPT) in South Africa. Until the observational coverage is improved, uncertainties will remain large. Given the cost and manpower involved in measuring fluxes in situ with a research vessel, it seems advantageous to pursue a denser network of regional sites along the southwestern coast of southern Africa. Such data would be of use to ongoing efforts to accurately determine terrestrial and marine greenhouse gas budgets for the region (e.g. Valentini et al., 2014;Bange et al., 2019).

30
In other EBUS, monitoring of some or all of atmospheric CO 2 , CH 4 , N 2 O, and δ(O 2 /N 2 ) is ongoing. Most of these sites are flask sampling locations, which are not as suitable as continuous measurements for capturing the short-lived nature of upwelling events. The California Current System is probably the best sampled, through two Advanced Global Atmospheric Gases Experiment (AGAGE) site, and numerous NOAA/ESRL activities. Long-term measurements are also made at the Cape Verde Atmospheric Observatory off of the Mauritanian upwelling region. While these efforts are extremely valuable, a greater density of continuous atmospheric measurement sites could provide great benefit in quantifying greenhouse gas fluxes from upwelling systems.

Summary and Conclusions
We have shown that coastal atmospheric anomalies of CO 2 , O 2 , N 2 O, and CH 4 can be related to upwelling events in the This upwelling area also functions as a significant source term in the CO 2 budget of the Benguela Current. 15 We have focused here on upwelling events, because they are distinguishable from other sources of intraseasonal variability in the atmospheric record. A full top-down accounting of the greenhouse gas budget of the the Benguela could be accomplished through a Bayesian atmospheric inversion of one or more coastal stations. Based on our results and previous studies (Nevison et al., 2004;Lueker et al., 2003) we suggest to establish a network of high-resolution atmospheric GHG measurements adjacent to major coastal upwelling regions in order to monitor their GHG emissions. Continuous land-based monitoring will help to Emeis, K.-C., Brüchert, V., Currie, B., Endler, R., Ferdelman, T., Kiessling, A., Leipe, T., Noli-Peard, K., Struck, U., and Vogt,           . Comparison of the variability of O2 with respect to N2O and CO2 with respect to O2 at NDAO and in surface water. Displayed are the data corresponding to atmospheric anomalies associated with upwelling events (left), of all marine boundary layer air masses as selected by back-trajectories (center), and dissolved concentrations of CO2, N2O, and O2 during M99 (right). Atmospheric O2 is expressed as APO in ppm equivalents, and dissolved concentrations are expressed as the difference between the measured concentration and the concentration at saturation, i.e., an excess, except for the bottom right panel, where it is shown as a deficit. In that plot, the Redfield ratio of 1.45 is plotted as a dotted red line for reference. Slopes (m) are given at the top of each plot. The black circles correspond with the upwelling event encountered during M99 (see Figure 8).