Marine Nitrous Oxide Emissions From Three Eastern Boundary Upwelling Systems Inferred From Atmospheric Observations

Eastern Boundary Upwelling Systems (EBUSs) are coastal hotspots of the potent greenhouse gas nitrous oxide (N2O). However, estimates of their emissions suffer from large uncertainties due to their significant spatial and temporal heterogeneity. Here, we derive the first multiyear, monthly resolution N2O emissions from three of the four major EBUSs using high‐frequency coastal atmospheric measurements and an inverse method. We find average combined N2O emissions from the northern California, Benguela, and southern Canary upwelling systems to be 57.7 (51.4–63.9) Gg‐N yr−1. We also find an offshore region near the Benguela EBUS that exhibits large pulses of emissions with emissions that reach 677 Gg‐N yr−1 in 1 month. Our findings highlight that atmospheric measurements coupled with inverse modeling can capture the large variability in EBUS emissions by quantifying emissions over large spatial distances and over long time periods compared to previous methods using traditional oceanographic measurements.


Introduction
Nitrous oxide (N 2 O) is a potent greenhouse gas and a major ozone depleting substance (Myhre et al., 2013;Ravishankara et al., 2009). Estimates of emissions from the ocean exhibit significant spread (Battaglia & Joos, 2018;Ciais et al., 2013) due to the challenge in simulating complex biogeochemical pathways, capturing large spatial and temporal variability, and due to sparse measurements. High concentrations of N 2 O in the ocean are found in regions known as Eastern Boundary Upwelling Systems (EBUSs), where high productivity rates due to upwelling lead to low oxygen conditions, favoring N 2 O production. EBUSs are often associated with ocean oxygen minimum zones (OMZs) (Capone & Hutchins, 2013;Oschlies et al., 2018). Strong upwelling in these regions also provides an efficient pathway for release of N 2 O into the atmosphere (Nevison et al., 2004). The four major EBUSs are associated with the California (eastern North Pacific), Benguela (eastern South Atlantic), Canary (eastern tropical North Atlantic), and Humboldt (eastern tropical South Pacific) Current Systems (Chavez & Messié, 2009).
Previous studies have shown that coastal areas can emit disproportionally large amounts of N 2 O compared to their fraction of global area (e.g., Naqvi et al. (2010)). However, previous estimates are based on methods that are not able to capture the significant spatial and temporal heterogeneity in coastal upwelling and thus may be inaccurate. Coastal upwelling events are episodic, occurring on the timescale of hours to days (Nevison et al., 2004) and with spatial extent that can vary seasonally and year-to-year. Previous methods to quantify EBUS N 2 O emissions have used sparse ship-based measurements of seawater N 2 O concentration (e.g., Arévalo-Martínez et al., 2015;Capelle & Tortell, 2016;Fenwick & Tortell, 2018;Frame et al., 2014;Kock et al., 2008;Wittke et al., 2010) or models employing climatological concentration fields and estimates of air-sea exchange (e.g., Buitenhuis et al., 2018;Nevison et al., 2004). Both methods suffer the challenge of capturing variability by being snapshots or by being based on composite N 2 O concentration fields, which have combined sparse measurements over decades. These limitations have resulted in large uncertainties in estimates of EBUS emissions.
Here, we present the first time series of spatially resolved estimates of EBUS N 2 O emissions derived from multiyear records of high-frequency atmospheric measurements and an inverse method, thus capturing variability in time and space. We used data from three coastal stations near the northern California, Benguela, and southern Canary (Mauritanian) EBUSs (no suitable atmospheric measurements are available near the Peruvian EBUS). The data set was composed of 15 years of atmospheric dry air mole fraction measurements from Trinidad Head, California (THD, Prinn et al., 2018), 2 years from the Namib Desert Atmospheric Observatory (NDAO, Morgan et al., 2015) and 4 years from the Cape Verde Atmospheric Observatory (CVAO). Measurements were coupled with the atmospheric transport model NAME (Numerical Atmospheric Modeling Environment) (Jones et al., 2006) at 3-hourly and up to 12 km spatial resolution and a hierarchical Bayesian inverse method (Ganesan et al., 2014;Lunt et al., 2016). We used a global ocean model, ECCO2-Darwin (described in Manizza et al. (2019) and updated to include N 2 O fluxes using Manizza et al. (2012); Nevison et al. (2003); Wanninkhof (1992) as described in the supporting information) at approximately 18 km and monthly resolution to provide a priori ocean N 2 O fluxes for the inversion.
The near-continuous nature of atmospheric measurements means that if wind directions are favorable, many episodic events can be captured in the measurement record and assessed over time. Previous studies (Lueker, 2004;Morgan et al., 2019;Nevison et al., 2004) that have used atmospheric measurements to estimate coastal emissions have had to attribute emissions to predefined source regions and relied on measurements sampling the ocean to not be conflated with terrestrial sources. We show that interpreting atmospheric measurements without a spatially and temporally resolved atmospheric transport model can lead to emissions magnitudes and their spatial distribution to be incorrectly derived.

Results
The first spatially resolved time series spanning multiple years of EBUS N 2 O emissions from three of the four major EBUSs are presented. Mean emission maps for the northern California, Benguela, and southern Canary EBUSs are shown in Figures 1 and S1 in the supporting information. Area-integrated emissions for each month for the coastal (0-150 km from coast) areas defined in Figure 1 are shown in Figure 2. Estimates from previous studies are provided in Table S1. Due to the episodic and variable nature of EBUS emissions, only a few direct comparisons to previous studies are possible.

Northern California (41°-50°N)
The northern California EBUS exhibits seasonality in wind patterns. The influence of oceanic emissions is observed in the atmospheric mole fraction record of N 2 O at THD each year from April through September. Outside of these months, wind directions are generally not favorable for observing an ocean source with sufficient sensitivity. Furthermore, the northerly winds needed to induce upwelling are found during these months, with winter months exhibiting reduced upwelling or downwelling conditions (Huyer, 1983). We therefore present emissions only for April through September through the inversion period.
Mean (95% confidence interval) area-integrated coastal ocean emissions from 2003-2017 are 19.2 (18.5-19.9) Gg-N yr −1 (Figure 2). If it is assumed that emissions outside of April through September are negligible, a lower-bound of mean annual emissions is 9.6 (9.3-10.0) Gg-N yr −1 . We find that maximum monthly emissions reach 80.4 Gg-N yr −1 , a value that is almost five times larger than the mean, highlighting the large variability in emissions that could be missed with sporadic sampling.
The increase in emissions relative to the ECCO2-Darwin model occurs primarily off the coast of Oregon and Washington, United States, and British Colombia, Canada ( Figure 1). Fluxes off the coast of Vancouver derived from ship-based measurements during the time period of this analysis (Capelle & Tortell, 2016;Fenwick & Tortell, 2018) show Vancouver to be a region of strong upwelling N 2 O emissions with mean emissions similar to those derived here (Table S1). Our mean per-area emissions of 2.4 ng-N m −2 s −1 is however 80% larger than estimated in (Nevison et al., 2004) using composite ΔpN 2 O fields for the same latitude band (Table S1), indicating that using climatological fields to derive emissions could be inaccurate for capturing variable sources.
Two previous studies (Lueker, 2004;Nevison et al., 2004Nevison et al., ) used 2000Nevison et al., -2002 atmospheric data from THD to derive N 2 O emissions using an atmospheric box model. There are several limitations with these previous approaches. First, estimated emissions were highly sensitive to atmospheric model inputs, which included atmospheric dilution, wind speed, planetary boundary layer height (PBLH), and the spatial footprint assigned to the mole fraction enhancement. Atmospheric dilution is described in Nevison et al. (2004) to be the least quantified parameter and a range of values were assigned. Wind speed and PBLH inputs were averages over the spatial area, and the two studies attributed emissions derived from the same data to different predefined spatial areas (i.e., 35°-50°N and 41°-50°N). Third, these estimates used depletion in atmospheric potential oxygen (APO), which is derived from measurements of the O 2 /N 2 ratio and carbon dioxide, to determine times of upwelling (Lueker et al., 2003). Corresponding enhancements in N 2 O mole fraction were then used to infer oceanic N 2 O emissions. As shown in Figure S2, N 2 O mole fraction enhancements during upwelling events at THD could also overlap with enhancements from terrestrial (natural soil Saikawa et al., 2013) and anthropogenic (Janssens-Maenhout et al., 2019) sources and therefore N 2 O enhancements should not be solely attributed to marine emissions. The atmospheric transport model used in this study uses three-dimensional meteorological fields at high spatial and temporal resolution to quantify the footprint of atmospheric enhancements and the inverse method solves for both land and marine emissions to minimize misattribution of enhancements.
Marine N 2 O emissions are governed by physical and biogeochemical drivers, which can vary with climatic conditions (Capone & Hutchins, 2013). The Pacific Decadal Oscillation (PDO) is the leading mode of variability in sea surface temperatures (SSTs) in the North Pacific ( Figure 3a). Variability in SSTs can drive changes in N 2 O emissions through changes in solubility, controlling the amount of N 2 O outgassing to the atmosphere, as well as through changes in upwelling and ventilation (Manizza et al., 2012). Warm and cool PDO phases correspond to higher and lower coastal SSTs, respectively, and higher and lower coastal upwelling strength, respectively ( Figure 3b).
As we have derived estimates of N 2 O emissions from the northern California EBUS spanning nearly two decades, we are for the first time able to correlate patterns in N 2 O emissions with climatic drivers. We long-term quantification not possible from sporadic ocean sampling, but future studies that also employ ocean biogeochemical models and ocean measurements would leverage a powerful combination of tools to diagnose the drivers of emissions.
We used model output from the ECCO2-Darwin model over 2009-2013 to investigate the contributions of thermally driven N 2 O fluxes and those driven by ventilation using a tracer which had biogeochemical processes suppressed but with the same solubility as N 2 O (Manizza et al., 2012). In the northern California EBUS, the model predicts ventilation fluxes from April to September to be around 2-4 times larger than thermal fluxes. Our findings, which show that N 2 O emissions correlate with both SST and upwelling, suggest that both processes could be important. However, further model investigations run over longer time periods that span PDO phases and include important biogeochemical drivers such as pH (Breider et al., 2019) are required to determine the relative contributions of the two processes on decadal timescales.

Benguela (23°-35°S) and Offshore South Atlantic
The location and wind patterns of NDAO make it useful for estimating ocean trace gas fluxes from atmospheric data. We investigated the degree to which land and ocean emission contributions are separated in mole fraction measurements at NDAO. As shown in Figure S3, N 2 O enhancements coincide both with depleted APO as well as with low carbon monoxide (CO) mole fractions. This implies that at many of the times when the ocean upwelling source is picked up at NDAO, there is little terrestrial and anthropogenic influence.
We averaged emissions from the Benguela EBUS over all months in the period Aug 2013-Sep 2015 because wind patterns at NDAO indicate that marine emissions from coastal and offshore regions can be picked up in the atmospheric record year-round. In addition, upwelling in the northern Benguela EBUS has been shown not to exhibit significant seasonality (Chavez & Messié, 2009).
Mean (95% confidence interval) area-integrated coastal emissions are found to be 28.4 (26.2-30.7) Gg-N yr −1 (Figure 2). Upwelling filaments, which can transport N 2 O offshore, have been found to occur 150-400 km from coast and these emissions should be considered as part of the upwelling system (Arévalo-Martínez et al., 2019). Offshore emissions 150-400 km from the coast are found to be 7.1 (5.5-9.0) Gg-N yr −1 (Figure 4).
Our mean per-area emissions of 3.5 ng-N m −2 s −1 are around 10 times larger than the climatological fluxes derived from composite ΔpN 2 O fields in (Nevison et al., 1995) (Table S1). However, mean emissions from August 2013 are consistent with those derived from ship-based measurements in the same month (Morgan et al., 2019). Atmospheric measurements from NDAO were also used in Morgan et al. (2019) to estimate upwelling emissions from the Walvis Bay and Lüderitz cells for August 2013, but estimated emissions were larger than those derived in this study and in the ship-based estimates (Table S1). We propose that the main reason for the lower coastal emissions estimated here using the same data set is that the atmospheric transport model and inverse method attributes some of the NDAO N 2 O mole fraction enhancements to offshore regions rather than to the coastal margin ( Figure S4). This finding highlights that studies that have predefined source regions to interpret atmospheric measurements could inaccurately quantify emissions.
Our emissions exhibit a similar spatial pattern to those derived in Arévalo-Martínez et al. (2019), with large emissions at 23°S (Walvis Bay) and reduced emissions between 23°and 27°S in the Lüderitz upwelling cell (Figure 1). Although measurements from other studies are not available for comparison, we find emissions south of 27°S to be of similar magnitude to those north of 23°S.
The southern boundary of the Benguela EBUS interacts with the very energetic Agulhas current, where filaments and large eddies can form offshore (Hutchings et al., 2009). Offshore emissions from this boundary region could be related to the Benguela EBUS, but because of the vague boundary definition, we quantified them separately. We estimate mean area-integrated emissions from an open ocean region to the south-west of the station (Figure 2) to be 56.8 (44.7-71.8) Gg-N yr −1 with these emissions reaching over 100 Gg-N yr −1 in several months and a maximum value of 677 Gg-N yr −1 in 1 month (Figure 4). One possible explanation is that these emissions could be associated with mesoscale eddies, which were present in the NDAO

10.1029/2020GL087822
Geophysical Research Letters measurement footprint during times of N 2 O mole fraction enhancements ( Figure S5) but may not be well-resolved in global ocean biogeochemistry models. Mesoscale eddies have been shown to have different biogeochemical properties and N 2 O emissions from the surrounding ocean (Grundle et al., 2017). However, to investigate this hypothesis, future work should directly sample eddies in this region, as has been done for other EBUSs (Arévalo-Martínez et al., 2016;Grundle et al., 2017). Coupling these measurements with a longer time series of data from NDAO would allow for greater process-level information to be inferred.
Our results highlight the large variation in emissions from this open ocean region, which would not likely be captured by sparse measurements. The only previous ship-based measurements from this region, which occurred prior to the beginning of the measurement record used in this study, show that very large ocean N 2 O mole fraction enhancements can exist (Weiss et al., 1992). These measurements show that enhancements are episodic with enhanced N 2 O only found in one leg of two ship transects separated by a period of 1 month. Global ocean estimates using composite ΔpN 2 O maps derived from these shipbased measurements could therefore substantially overestimate or underestimate fluxes for such strongly variable regions.

Southern Canary (16°-23°N)
Upwelling in the southern Canary EBUS is semicontinuous (Chavez & Messié, 2009). However, our estimates are constrained by measurements only in a subset of months as discussed in section 2.4. Mean (95% confidence interval) area-integrated emissions for the months during 2013-2017 that are constrained by measurements are 12.7 (10.4-15.8) Gg-N yr −1 (Figure 2). This corresponds to mean per-area emissions of 2.7 ng-N m −2 s −1 , which is four times larger than the climatological flux shown in (Nevison et al., 1995) (

Sensitivity Tests
Three sensitivity inversions were performed to test the robustness of our results: (i) Atmospheric measurements reflect the net effect of all sources and sinks of atmospheric trace gases upwind of a receptor. To demonstrate that the estimation of non-ocean sources did not significantly impact the estimation of the ocean source, we derived emissions using a subset of data that were filtered to exclude data points that were heavily influenced by anthropogenic or soil sources. The procedure for data filtering is discussed in the supporting information. Emissions derived from the filtered data set ( Figure S6) show that results for the three EBUSs and the offshore emissions from the Benguela are consistent with the unfiltered estimates. This finding suggests that the inversion is able to partition land and ocean sources because when some contribution of the land-based sources is removed, the inversion still estimates similar ocean emissions. (ii) We tested for the influence of the a priori emissions on derived emissions by scaling the total a priori (ocean and land) emissions by 2-10 times the original value, keeping the remainder of the methodology the same.
In the northern California and Benguela regions, similar emissions were derived for all months, confirming the atmospheric constraint on the ocean source. In the southern Canary, results are consistent for a majority of months but those that are not have been excluded from the analysis ( Figure S7). Because the total a priori emissions were scaled, this resulted in a large perturbation to emissions from the land sector, particularly for the northern California region where there are more significant land sources than in the Benguela or southern Canary regions. The consistent emissions derived for the northern California EBUS provide confidence in the ability of the inversion to separate ocean and land emissions, for if this separation were dependent on the a priori emissions, a large perturbation to the land emissions would affect the derived ocean emissions. (iii) We assessed the effect of the a priori boundary condition field on derived emissions by using two global model fields, MOZART (Emmons et al., 2010;Palmer et al., 2018) and CAMS LMDZ (Thompson, Chevallier, et al., 2014), which are discussed in the supporting information. We show in Figure S8 the prior and posterior boundary conditions at each site from the two models as well as the emissions estimated using each of these boundary conditions. These results show that emissions are consistent within uncertainties for the different boundary conditions and when offsets to the model boundary conditions are also estimated in the inversion that a single site can constrain both boundary conditions and emissions.
We also carried out tests to investigate whether the representation of the coast in the model could have strongly impacted our results. Underlying model processes important for resolving coastal features, such as land-sea breezes, are strongly dependent on model resolution. The spatial resolution of the meteorological fields driving NAME increased from 60 km in 2003 to 12 km in 2017. The fact that emissions are being derived with similar magnitude throughout the period provides confidence in the ability of the model to partition emissions along the coastal boundary. As we aggregated our emission maps into total coastal emissions using a land-sea border that is defined at the resolution of the modeloutput (0.352º x 0.234º), we also aggregated these emissions using different border definitions (i.e., by moving the border one or two grid cells or approximately 30-60 km inland or offshore), keeping the coastal definition the same (0-150 km from the border) ( Figure S9). The main impact on aggregated emissions using different coastal boundaries is in northern California, where land sources are more significant than at the other two EBUSs. Mean April-September emissions in the northern California EBUS could range between 14.4 and 24.6 Gg-N yr −1 , compared to our result of 19.2 (18.5-19.9) Gg-N yr −1 , depending on whether the border is moved one grid cell offshore or one grid cell inland. This represents an extreme perturbation (i.e., a~30 km change to the coastal boundary) but suggests that while differences lie outside of the confidence interval of the main results, that a similar magnitude of emissions is being derived. Differences in the other EBUSs are minimal. This experiment indicates that transport model uncertainty at the coast (e.g., in representing land-sea breezes) and uncertainty in coastal definition, which if significant would result in large changes when aggregating emissions using different borders, do not substantially alter the conclusions of this study.
Despite robustness of the results to the above tests, there could still be systematic uncertainties that are not accounted for. These could be due to, for example, vertical mixing in NAME, the representation of the planetary boundary layer, or other structures in the inversion framework. Quantifying these uncertainties would require models to be assessed regularly through for example, independent tracer release campaigns. Model intercomparison studies (e.g., Bergamaschi et al., 2015; and observing system simulation experiments (e.g., Wells et al., 2015) have attempted to quantify some of the uncertainties in current inverse modeling capability for N 2 O.

Conclusions and Discussion
The average combined coastal N 2 O emissions for the three EBUSs are 50.6 (45.6-56.1) Gg-N yr −1 which increases to 57.7 (51.4-63.9) Gg-N yr −1 when including the 150-400 km band of emissions from the South Atlantic. Mean emissions from each of the three EBUSs are of similar magnitude; however, the largest pulses of emissions occur from the Benguela EBUS. Both the northern California and Benguela EBUSs have maximum monthly emissions that are 4 and 5 times greater than the mean, while there is little variability from the southern Canary. The time series of emissions that we derive from atmospheric measurements make it possible for the first time to quantify this variability.
Significant offshore emissions were only found in the South Atlantic near the Benguela EBUS and were not present in the eastern North Pacific or eastern tropical North Atlantic. In the South Atlantic, we identified several months with very large pulses of emissions (>600 Gg-N yr −1 ) from an open ocean area where there could be interaction between the Benguela EBUS and the Agulhas current. These pulses are an order of magnitude larger than annual average emissions, highlighting the significant variability in these sources.
Previous studies have estimated these three EBUSs to contribute 42 Gg-N yr −1 for larger latitude extents than used here (Nevison et al., 2004). Our estimates do not cover more southerly extents of the California EBUS, regions north of NDAO in the Benguela, or the Humboldt EBUS, and there could be important additional contributions from these areas. A study based on oceanographic measurements shows that emissions off Peru could be large (200-900 Gg-N yr −1 ) (Arévalo-Martínez et al., 2015).
If atmospheric measurements could be implemented in the EBUSs not captured by this study as well as in open ocean regions that are influenced by the EBUSs, the method used here could quantify the magnitude and variability in these emissions over time and provide a more complete account of global ocean N 2 O emissions. Measurement stations should be situated near upwelling regions and ideally, far from other sources. The primary limitation of this approach is that land-based sources need to be robustly accounted for, and in some regions, these emissions may be much larger than coastal ocean emissions. Including measurements of APO and anthropogenic tracers such as CO would help to diagnose any such influences. In addition, more frequent campaigns of simultaneous ocean and atmospheric measurements would allow for regular assessment and comparison of the fluxes derived from the two methods.
Recent studies have shown that over the previous decades, ocean warming and its reduced ventilation have caused deoxygenation and expansion of OMZs, including in the EBUSs (Oschlies et al., 2018). While responses depend on time-scales and regions, model studies predict significant changes in N 2 O production and emissions in the future (Battaglia & Joos, 2018). Coupled with intensified coastal upwelling (Wang et al., 2015), increased production could lead to greater emissions to the atmosphere, reenforcing the positive feedback between ocean biogeochemical processes and climate warming. As we have shown here, atmospheric measurements coupled with high-resolution transport modeling and an inverse method could provide us with the means to quantify any such long-term changes in the EBUSs.