Low-pressure storms drive nitrous oxide emissions in the Southern Ocean

Domain: neuroscience
Published: 2026-01-26
Editorial outcome: major_revision

Abstract
Nitrous oxide is a potent greenhouse gas and the primary ozone-depleting agent of the 21st century, but marine emissions of nitrous oxide remain difficult to constrain due to their spatiotemporal variability. In the Southern Ocean, where extratropical cyclones create conditions conducive to air-sea gas flux, shipboard measurements are unlikely to capture the full extent of nitrous oxide emissions due to the impracticality of sampling said storms. Here, we use machine learning to derive nitrous oxide observations from biogeochemical Argo floats, revealing that low-pressure storms amplify air-sea gradients and create hotspots of emissions. Taking these low-pressure storms into account, rather than assuming 1 atmosphere (the standard condition outside of storms), increases the net annual Southern Ocean nitrous oxide flux by 88%. Our results imply that the Southern Ocean plays a significant role in the global nitrous oxide cycle, and may be a weaker overall sink of greenhouse gases than previously thought.

Introduction
Nitrous oxide (N2O) is the leading ozone-depleting agent of the 21st century1–3 and a potent greenhouse gas, with a global warming potential 273 times that of carbon dioxide (CO2) on a 100-year timescale1. The ocean, a net source, plays a critical role in the global N2O budget, yet the magnitude of marine N2O emissions remains highly uncertain, particularly in dynamic, historically undersampled regions such as the Southern Ocean4–6. This uncertainty limits our ability to assess the impact of marine N2O emissions on climate, because even trace amounts of N2O emitted from the ocean may play a major role in offsetting marine CO2 uptake in terms of radiative forcing7,8. Furthermore, this counteracting effect may be directly stimulated by enhanced biological activity, such as that triggered by natural or intentional ocean iron fertilization7,8.

The Southern Ocean is a major contributor to global marine N2O emissions, yet estimates of its contribution have varied widely. Initial studies suggested that the Southern Ocean accounted for 24–43% of total marine N2O emissions4,5,9, but recent work revised the Southern Ocean flux down to only 19% of global emissions6. These inconsistencies leave a wide range of estimates of the Southern Ocean N2O source (0.8–1.6 Tg N/year; Table 1), making it difficult to assess its significance as a greenhouse gas source to the atmosphere. Furthermore, previous estimates relied on sparse sampling, spatial averaging, and/or climatological monthly means of air-sea N2O disequilibrium, and are thus unlikely to capture the sub-seasonal scales of variability that are important to Southern Ocean biogeochemistry and air-sea gas exchange10–12.

Here, we apply machine learning models to high space- and time-resolution Biogeochemical Argo (BGC-Argo) in situ sensor data to investigate a key driver of sub-seasonal variability: the passage of strong storms through latitudes 50–65°S every 4–8 days13–15. The barometric pressure in the centers of extratropical cyclones can drop to as much as 8% below 1 atm (Fig. S1–S12; S2, 4), creating conditions favorable for N2O outgassing. Recent work on storm-driven CO2 fluxes demonstrated that storms significantly enhance air-sea gas exchange, driving high CO2 fluxes through turbulent entrainment of the deeper inorganic carbon reservoir16–18. Under-sampling storm-driven CO2 fluxes, which occur on timescales of days, can distort seasonal and annual air-sea CO2 flux estimates11,12. Previously impossible due to the lack of observational coverage, we can now quantify the effect of storms on Southern Ocean N2O fluxes using BGC-Argo data, and thus better constrain the role of the Southern Ocean in greenhouse gas exchange with the atmosphere.

Results
Float-based observations of N2O flux
To project sparse N2O measurements onto distributed sensor data, we trained machine learning models on the Global Ocean Ship-based Hydrographic Investigations Program (GO-SHIP) N2O dataset (n=1968 profiles; Fig. S13), then applied these models to predict N2O concentrations from parameters measured by profiling BGC-Argo floats (n=28,112 profiles; Fig. S14). BGC-Argo floats provide year-round, distributed measurements of ocean biogeochemistry, and, along with gliders, have transformed our understanding of Southern Ocean CO2 emissions16,17,19. Because BGC-Argo floats cannot measure N2O directly, we turned to supervised learning algorithms, which have been used to map chemical distributions of iodide20, methane21, pCO222, and surface N2O saturation6. Supervised learning algorithms are a natural framework for modeling the sharp gradients in marine N2O concentrations because, unlike multiple linear regression, they are robust to nonlinearity and perform well on the smaller datasets inherent to labor-intensive ocean sampling23. We found that dissolved N2O is highly predictable from dissolved oxygen, conservative temperature, absolute salinity, and nitrate concentration (Fig. S15). Based on a test dataset comprised of profiles withheld from model training, our models predicted the seawater partial pressure of N2O (pN2O sw) with R²≥0.91, and root mean squared error ≤ 21 natm in the Southern Ocean (Fig. S15).

Storms drive hotspots of N2O flux
Applying these models to BGC-Argo float profile data revealed that high sea-to-air N2O fluxes tend to be collocated with coherent low-pressure centers with high wind speeds (Fig. S1–S12). The air-sea partial pressure disequilibrium of N2O is defined as:

ΔpN2O = pN2O sw − pN2O atm, (1)

where pN2O sw is the partial pressure of N2O dissolved in surface seawater, pN2O atm is the atmospheric partial pressure of N2O, P is the barometric pressure, and χN2O atm is the dry mixing ratio of N2O in air. Low barometric pressures in the centers of storms reduce pN2O atm, widening ΔpN2O (Fig. 1b) and driving greater flux out of the water (Fig. 1c). In other words, the low barometric pressure in storms acts to suck gas out of the water; we term this “the Hoover Effect.”

To quantify the integrated effect of low barometric pressure in storms on air-sea N2O fluxes, we calculated air-sea N2O fluxes with the global mean barometric pressure (1 atm) in place of the observed barometric pressures (with wind and all other inputs being equal). We compare against fluxes calculated with 1 atm, as opposed to the mean Southern Ocean barometric pressure, because the characteristically low mean Southern Ocean pressure is itself a product of frequent storm activity13. In particular, storms tend to converge at 60° S18,24, creating the upward motion that dynamically maintains the mean low pressure of the Subpolar Low25–27.

The difference between the fluxes based on true barometric pressure and 1 atm is of a similar order of magnitude to the flux itself (Fig. S16), indicating that barometric pressure is a primary driver of high fluxes in the centers of these storms. This result implies that low-pressure storms create conditions that promote hotspots of N2O flux, and conversely, that estimates based on smoothed barometric pressure fields or 1 atm barometric pressure may be missing a substantial proportion of the net annual Southern Ocean N2O flux. When we take these low barometric pressures, or the Hoover Effect, into account, we find that the net annual Southern Ocean N2O source is 1.6 ± 0.3 Tg N2O-N/year (Fig. 2). In contrast, when we calculate the net annual Southern Ocean N2O source with 1 atm barometric pressure, the net annual flux drops to 0.9 ± 0.3 Tg N/year, or by a margin of 50 ± 20% (Fig. 2). This observation suggests that roughly half of the total flux is driven by deviations in barometric pressure from the global mean.

Furthermore, low barometric pressures and high wind speeds act synergistically in storms, creating a nonlinear response. The mixed partial derivative term of the Taylor expansion of FW14 (equation S4) with respect to wind speed and barometric pressure, centered on the median wind speed Umed and 1 atm barometric pressure P1atm, has the form:

∂2FW14 / ∂P∂U ≈ 2 * Umed * 0.251 * Sc−1/2 * (1 − χN2O atm) * U10 − Umed * P − P1atm (2)

Where U10 is the wind speed at 10 meters, Umed is the median Southern Ocean wind speed, P1atm is a barometric pressure of 1 atm, s is the Henry’s Law constant, and Sc is the Schmidt number. Given that, in storms U10 − Umed > 0 and P − P1atm < 0, the sign of this second-order mixed partial derivative is positive, indicating that the combination of high winds and low pressure in storms produces more flux than would be expected from their individual effects alone. The Taylor expansion of FL13 (eqns. S21–S23) is more complex but yields a similar result (Supplementary Text S2).

To obtain a discrete approximation of the interaction effect, we can compare the observed flux to 1) a “baseline” case calculated with 1 atm barometric pressure and median piston velocities (k); 2) a “wind effect” case calculated with 1 atm barometric pressure and observed piston velocities; and 3) a “pressure effect” case calculated with observed barometric pressures and median piston velocities (Supplementary Text S2). This empirically calculated interaction effect accounts for a significant proportion of the overall flux (Fig. 3).

The importance of storms is also evident in the small number of profiles that contributed disproportionately to the net annual flux. For example, only 5% of the float profiles comprise 34 ± 9% of the net annual flux; 10 ± 3% of profiles make up 50% of the total flux (Fig. S17). Together, these results imply that brief, intense flux events in the middle of Southern Ocean storms account for a disproportionate amount of the net annual Southern Ocean N2O emissions. Previous estimates5,6 do not capture these brief, intense flux events in storms.

Biotic vs. abiotic drivers
It is a surprising result that a difference of only up to 8% in barometric pressure can drive an 88% increase in the net annual Southern Ocean N2O flux. Carbon studies show that air-sea CO2 flux is strongly mediated by the imprint of biology, with outgassing associated with upwelling/entrainment of remineralized carbon. In contrast, for N2O, we hypothesized that storms modulate sea-to-air flux primarily via lowering the barometric pressure, as opposed to entraining the deep N2O reservoir, due to the relatively low, uniform pN2O sw in the Southern Ocean (Fig. S18). To test this hypothesis, we performed a sensitivity test holding pN2O sw constant; this test produced a much smaller change than holding the barometric pressure constant at 1 atm (Fig. S19). The importance of barometric pressure in dominating flux also means that our results are somewhat agnostic to uncertainties in the machine learning prediction: because the effect is driven by low atmospheric pressure, not high surface pN2O sw, the result doesn’t depend on perfectly capturing high pN2O sw in the centers of storms.

The large impact of variations in barometric pressure on N2O flux arises from the relatively small N2O partial pressure disequilibrium in much of the Southern Ocean (Fig. 1b): the median pN2O sw in the Southern Ocean is 334 natm, while the median partial pressure disequilibrium is only 11 natm. This means that small changes in barometric pressure represent a large percentage change in pressure disequilibrium, and thus large changes in air-sea flux. In the Antarctic Southern Zone, the surface water pN2O is close enough to equilibrium with the atmosphere that the pressure effect accounts for 100%—or even reverses the direction—of air-sea flux (Figs. 5, S20).

Seasonal cycles in N2O flux
While the magnitude of air-sea flux is driven by barometric pressure, the seasonality of flux is driven by the seasonal cycle in surface seawater dissolved N2O (pN2O sw), and, to a lesser extent, wind speed. The overall seasonal cycle suggests an autumnal peak in Southern Ocean integrated flux (Fig. 4), matched by April-June flux maxima in four out of five zones (Fig. 5). Drivers of this autumnal maximum are evident from sensitivity tests performed by calculating fluxes with 1 atm barometric pressure, median piston velocities (k), and median pN2O sw (Fig. 4). Fluxes calculated with 1 atm barometric pressure have a lower magnitude but the same seasonal cycle as the observed fluxes (Figs. 4a, 5); in contrast, fluxes calculated with median k values display a muted autumnal maximum (Fig. 4b), and fluxes calculated with median pN2O sw display a muted autumnal maximum (Fig. 4c). This effect is driven by the seasonal cycle in pN2O sw in each zone (Fig. S21). In the Subantarctic Zone and Polar Frontal Zone, there are particularly strong peaks in pN2O sw in April through June (Fig. S21f, i), which correspond to the autumnal entrainment maximum in these zones29. Future work should investigate the role of circulation in bringing N2O produced at depth to the surface.

There is no seasonal cycle in barometric pressure in our dataset, but rather a large degree of spread throughout the year, reflecting the fact that floats encounter low-pressure centers associated with storms and intervening high-pressure ridges (Figs. S21, S22). The distribution shift towards lower pressures in the southernmost zones (Fig. S21) reflects the increasing storm density from north to south, which has been previously observed in the Southern Ocean18. Averaging over these points produces a general decline in barometric pressure moving from north to south, with the lowest barometric pressures in the Antarctic Southern Zone and Southern Ice Zone (Fig. S21j, m).

Sea ice causes fluxes to decrease over the course of austral winter in the southern ice zone; in the absence of sea ice limiting exchange, N2O fluxes in the southern ice zone could reach a September maximum of >0.1 Tg N/year, more than twice as high as the peak flux in any other zone (Fig. 5j). Estimating the net annual flux based on peak months (April-June) vs. trough months (October-December) leads to statistically significant differences in each zone (Fig. S20), indicating the importance of sampling the full seasonal range of sea-to-air N2O fluxes in the Southern Ocean. Thus, extrapolating from summertime shipboard measurements may lead to an underestimate of the total annual flux, especially if those measurements take place during the early part of austral summer.

A higher estimate of the annual Southern Ocean N2O flux
Our float-based estimate of the net annual Southern Ocean N2O source (1.6 ± 0.3 Tg N/year) represents a significant upward revision from earlier work6 (0.8 ± 0.2 Tg N/year; Table 1). The discrepancy may be explained by the pressure effect: when using 1 atm barometric pressure, our estimate (0.9 ± 0.3 Tg N/year) is within error of previous work6 (Table 1). Previous work6 focused on predicting monthly mean N2O disequilibria (ΔχN2O) with monthly mean World Ocean Atlas30 data, which may have smoothed over low-pressure centers in the centers of storms. In contrast, the BGC-Argo in situ observations used here allowed us to predict pN2O, and, thus, calculate fluxes on sub-daily timescales. Furthermore, our models were trained on both surface and subsurface data, which may have allowed them to better capture Southern Ocean winter conditions (i.e., source waters).

Our improved quantification of N2O flux suggests that the Southern Ocean is a more important component of the global N2O cycle than previously thought. Given a global N2O flux of 4.2 ± 1 Tg N/year6, our estimate indicates that the Southern Ocean represents 38 ± 12% of all marine N2O emissions, which represents a significant upward revision. An important caveat is that this estimate of the global source (4.2 ± 1 Tg N/year)6 is subject to the same uncertainties and potential aliasing as the corresponding earlier estimates of Southern Ocean N2O emissions. Our estimate also implies that N2O flux from the Southern Ocean offsets a greater proportion of CO2 uptake by the same region: using a N2O global warming potential of 2731, our updated flux predicts that Southern Ocean N2O outgassing offsets 218 ± 41 Tg CO2 uptake per year, twice the 109 ± 27 Tg CO2/year suggested previously6. Assuming an annual net Southern Ocean CO2 uptake of 0.87 Pg C/year31, this means that N2O outgassing from the Southern Ocean offsets 7% of the annual CO2 uptake by the same region.

Limitations of this analysis and avenues for future study
There are several key limitations of this analysis that could be improved with further work. Firstly, BGC-Argo oxygen measurements tend to be biased low32. A first-order correction for this bias resulted in a ~10% reduction in the net annual N2O flux, which is within our range of uncertainty (Supplementary Text S3). Correcting for float oxygen biases is still an area of active research32, however, so no correction was applied here. The observed reconstruction bias (“Methods”, vide infra) may have resulted in part from float oxygen bias, in addition to spatiotemporal sampling differences. Secondly, given typical ocean decorrelation length and time scales of ~160 km and ~60 days, respectively33,34, float profiles every 10 days may not be entirely independent, leading to an underestimate of the overall uncertainty. Third, the parameterization of N2O air-sea gas exchange is still an active area of research, and the subjectivity in the choice of parameterization is one of the leading sources of uncertainty (Fig. 2). Since CO2 and N2O have similar solubilities, eddy covariance CO2 flux studies35–38 could help provide insight into N2O air-sea gas exchange, following previous work using eddy covariance to better constrain bubble fluxes of CO239,40. Future work should attempt to use the insights from these eddy covariance studies to help narrow the uncertainties due to the subjectivity of parameterization of air-sea gas exchange. Relatedly, the interactions between air-sea gas exchange and sea ice are more complex than the simple fractional correction applied here would suggest41–43. More work is needed on biogeochemical N2O cycling within and around sea ice itself, as well as a more nuanced treatment of how different types of sea ice (pancake ice, leads in sea ice) modulate air-sea gas exchange.

While BGC-Argo data improve spatial and temporal coverage compared to shipboard data, there are still locations where no float profiles exist, and additional data would likely improve our estimate. Although we perform a bootstrap analysis to assess uncertainty due to sparse coverage (Eq. 15), the estimated absolute flux values presented here apply only to the regions and seasons where float data are available. Missing coastal and shelf processes may be an especially important source of potential bias in our N2O fluxes. BGC-Argo provides good coverage of the open ocean but is missing coastal and shelf regions, where both emissions and uncertainties are high44; future work could take advantage of other autonomous technologies, such as gliders and moorings, to better constrain N2O emissions in the near-coastal regions of the Southern Ocean and elsewhere. There is also a need for more high-quality shipboard N2O measurements, especially in coastal regions, to improve the representativeness of the training dataset. Comparing artificially-seeded, autonomously-measured N2O fluxes to a model ground truth45 could help evaluate where additional N2O measurements would be most impactful.

Future work could also situate the air-sea N2O fluxes calculated here in the context of storm geometry. A recent study found weaker storm activity in the Atlantic sector of the Southern Ocean18; the storm masks from that study could be used to quantify the percentage of high-N2O-flux profiles that occur within defined storms, and break this down by sector to quantify whether there are indeed fewer low-pressure points in the Atlantic sector. Visually, the highest N2O fluxes appear to occur on the northern edge of each low-pressure center (Fig. S1–S12); this aligns with the same recent study, which showed that the highest wind speeds tend to occur on the northern edge of the low-pressure center with respect to storm geometry18. An analysis of N2O fluxes with respect to storm geometry could be used to quantitatively evaluate this pattern.

Potential climate feedbacks
While changes in N2O flux under different warming scenarios are difficult to predict, we can perform sensitivity tests to demonstrate the potential response of N2O flux to different climate forcings such as increasing wind stress46–48, decreasing barometric pressure47,49, and increasing storm frequency28. The mean barometric pressure over the Southern Ocean is predicted to change by only 0.002 atm by 210026, not enough of a change to generate a significant increase in the annual net N2O flux, but a more drastic shift in the mean barometric pressure of 0.01 atm increased the annual net N2O flux by 31% to 2.1 ± 0.4 Tg N/year (Table 1). Models tend to disagree on both the amount and distribution of Southern Ocean sea ice loss50,51, so we performed a sensitivity test removing the capping effect of sea ice in its entirety. This increased the overall flux by 25% to 2.0 ± 0.4 Tg N/year (Table 1; Fig. S23), driven primarily by strong outgassing in the Southern Ice Zone (Fig. 5). Increasing wind speeds by 25%, which is the predicted increase in Southern Ocean wind stress under RCP8.525, increased N2O flux by 38% 2.2 ± 0.7 Tg N/year (Table 1; Fig. S23). The combined effects of increasing wind speeds by 25% and removing the sea ice barrier increased the net annual flux by 81% to 2.9 ± 0.8 Tg N/year (Table 1; Fig. S23). Finally, as shown above, low barometric pressures and high wind speeds have an interaction effect in storms, which could lead to a nonlinear response of flux to increasing storm frequency. While these sensitivity tests are not simulations of future climate scenarios with associated feedbacks, they demonstrate that the Southern Ocean N2O source is a dynamic system that is poised to experience a significant ramp-up in N2O emissions as the region changes.

Methods
Datasets
We compiled the Random Forest model training dataset from N2O measurements collected during CLIVAR (Climate and Ocean: Variability, Predictability and Change) and GO-SHIP (Global Ocean Ship-Based Hydrographic Investigations Program). The GO-SHIP dataset provides full water column, quality-controlled N2O data from 2005 to 2023 (Fig. S13; Supplementary Table 3). We restricted our training dataset to the GO-SHIP program due to the high level of quality control applied to the N2O measurements and accompanying hydrographic and nutrient measurements taken from the same Niskin bottle. GO-SHIP data were filtered to retain only samples for which the CTD pressure, temperature, salinity, and oxygen were flagged as good (flags 2, 6, and 7), the nitrate and N2O water samples were flagged as good (flags 2, 6, 7, and 8), and the Niskin bottle was flagged as good (2). This filtering left 35,664 data points from 1968 unique profiles from 24 cruises in the training dataset. GO-SHIP data from the full water column (as opposed to just the surface) were included in model training to better capture the conditions and water masses from which N2O in the Southern Ocean originates.

BGC-Argo is a global, autonomously drifting float array with in situ sensors for temperature, salinity, dissolved oxygen, nitrate, and other biogeochemically relevant parameters52. BGC-Argo floats sample for these parameters from 2000m depth to the surface every 10 days, then transmit the data back to data assembly centers on land52. For this study, BGC-Argo float data were obtained from a snapshot of the BGC-Argo dataset from October 9th, 202453. The float dataset was limited to data points where [O2] and nitrate concentration ([NO3−]) were both available; excluding [NO3−] from the predictor features would increase the number of float profiles available for this analysis, which could be explored in further studies. Float data were filtered to retain only delayed-mode-adjusted pressure, temperature, salinity, [O2], and [NO3−] with quality flags 1 (“good”), 2 (“probably good”), 5 (“value changed”), and 8 (“interpolated, extrapolated, or other estimation”). Float data were further restricted to points south of 30°S. While the full water column was considered in model training, only surface points were retained in the BGC-Argo dataset because the parameter of interest was air-sea flux. To handle the ice avoidance algorithm, which prevents floats from rising shallower than 15–20m when under ice, surface points were extracted by taking the shallowest point in each profile that was also shallower than 30m. This left 28,112 surface data points from an equivalent number of unique profiles from 325 total floats in the Southern Ocean, starting in 2008 (Fig. S14). Floats were assigned to zones based on climatological mean locations of Southern Ocean fronts19,54 (Fig. S14a). The effective decorrelation length scales, calculated on a monthly mean basis, ranged from 160 to 240 km per float, with the best coverage in the Southern Ice Zone (Fig. S14b). Effective decorrelation timescales ranged from 0.04 to 0.06 days per float. While each float samples every 10 days, they sample at random intervals, providing much finer temporal resolution in aggregate55.

The atmospheric mixing ratio of N2O (χN2O atm) was obtained from NOAA flask dataset measurements at the Palmer Station (64.6°S, 64°E) and South Pole (90°S, 0°E) monitoring stations56. Assuming that the atmosphere is zonally well-mixed, the Southern Ocean monthly mean was estimated by weighting each station by the cosine of its latitude (meaning that χN2O atm was primarily driven by the Palmer Station). Monthly mean χN2O atm were linearly interpolated in time to the date and time of each float profile.

Barometric pressure and wind speed were obtained from surface-level ERA557 and NCEP38 reanalysis data, matching each float profile to the nearest 6-hourly ERA5 and NCEP values in space and time (Figs. 1; S4b, S4c). Fractional sea ice cover was obtained from ERA5, because NCEP sea ice cover is only 0 or 138. NCEP winds tended to be higher than ERA5 winds (Fig. S24a), while the barometric pressures from the two datasets were more similar (Fig. S24b).

Derived variables
Absolute salinity (SA), potential temperature (θ), conservative temperature (Θ), density (ρ), potential density anomaly (σθ), oxygen saturation, and spiciness were calculated with the Thermodynamic Equation of Seawater 2010 (TEOS-10) Python implementation58. Apparent oxygen utilization (AOU) was calculated from the dissolved oxygen concentration ([O2]) and solubility of oxygen at the in situ temperature and salinity59. The solubility of N2O in seawater was calculated based on ref. 60 using an updated version of a gas exchange toolbox61. The partial pressure of N2O, the target variable for the machine learning models, was calculated from:

pN2O sw = [N2O] * ρ / kH (3)

where pN2O sw represents the seawater partial pressure of N2O in natm, [N2O] is the concentration of N2O in nmol kg−1, ρ is the in-situ density in kg/m3, and kH is the Henry’s Law constant in mol m−3 atm−1. Additional derived variables included the day of the year, transformed as cos(2πDOY/365.25) to eliminate the discontinuity at the end of each calendar year, and longitude, transformed into two different predictors following refs. 62,63 as cos(π(LONGITUDE – 110°E)/180) (“LON1”) and cos(π(LONGITUDE – 20°E)/180) (“LON2”) to ensure that regions of maximum explanatory power occurred over the oceans.

Random forest model tuning and training
We trained random forest regressor models to predict pN2O sw from variables measured by BGC-Argo floats. Prior to random forest model training, a “test” dataset was generated by randomly withholding 20% of profiles from the full GO-SHIP dataset (Fig. S13). A random 20% of profiles, as opposed to 20% of data points, were withheld to ensure that inherent correlations among the data points from a single profile did not artificially boost the apparent effectiveness of the machine learning algorithm. The test dataset was completely withheld from the model selection, hyperparameter tuning, and training process (vide infra), to be used to evaluate each random forest model’s performance on unseen data after all training and validation steps.

The “validation” dataset was generated by withholding 20% of the data points from the remaining 80% of profiles. The validation dataset was used during the model selection and hyperparameter tuning process to fine-tune the model and prevent overfitting, helping to choose the best model parameters. This left 64% of the training dataset for model training during feature selection and hyperparameter tuning.

Six random forest regressor models, based on six candidate sets of predictor variables (Fig. S25), were tested as potential contributors to the final model ensemble. Instead of testing all possible combinations of predictor features, predictor features were chosen that were likely to provide information about water masses and N2O production mechanisms, and whose range in the training dataset covered the full dynamic range of values in the BGC-Argo dataset. Additionally, a model containing a larger number of predictor features was included in the model selection process as a point of comparison. To assess the uncertainty in each candidate model’s performance, a K-fold cross-validation was performed by training each model five times, withholding five distinct sets of data (folds) of validation data such that each point appeared in the validation dataset exactly once. Models were trained using the RandomForestRegressor from sci-kit-learn64, with 600 decision trees following ref. 62 and a minimum of one sample per leaf to maximize model performance. For each of the five folds, model R2, root mean square error (RMSE), and mean absolute error were obtained by applying the model to the validation dataset.

To account for potential biases in the BGC-Argo data compared to the GO-SHIP data, we calculated the reconstruction bias for each model as the difference between the median predicted pN2O sw and the median Southern Ocean surface pN2O sw in the GO-SHIP dataset. Reconstruction bias was calculated as the difference between these two values:

bias = median(pN2O)BGC − median(pN2O)GO-SHIP (4)

where median(pN2O)BGC is the median Southern Ocean surface pN2O sw predicted from BGC-Argo float data in natm, and median(pN2O)GO-SHIP is the median Southern Ocean surface pN2O in the GO-SHIP dataset. This procedure was repeated five times per candidate model, for each of the five folds, to obtain an uncertainty in that model’s bias.

The K-fold cross-validation revealed that adding more predictor features lowered RMSE but increased reconstruction bias (Fig. S25), likely due to overfitting to the training data, reducing the influence of core oceanographic features, and/or decreasing the effective sample size per feature65,66. Thus, to minimize overfitting and bias, we chose the simplest four models—including only absolute salinity, conservative temperature, either [O2] or AOU, and [NO3−]—to include in the final model ensemble. Because each Random Forest model consists of 600 decision trees, our approach is essentially an ensemble of ensembles.

The hyperparameters for each of the four chosen models were tuned by manually searching ranges of minimum samples per leaf, the number of decision trees, and minimum samples per split to find the combination of hyperparameters for each model that minimized that model’s RMSE in the Southern Ocean. The optimum values for minimum samples per leaf and minimum samples per split were one and two, respectively, so these were used for final model training. The number of estimators had minimal impact on model performance, so we used 600 estimators for final model training to minimize overfitting.

Using these hyperparameters, each of the four chosen models was re-trained on the combined training and validation dataset (e.g., 80% of stations that were not withheld as test data). Model R2, RMSE, and MAE were then evaluated on the test data (e.g., the 20% of stations that were initially withheld). [O2] and AOU were the most important predictor variables; the exclusion of [NO3−] made little difference in the Southern Ocean, indicating that future work could exclude this variable to access a greater number of float profiles (Fig. S15). The models tended to underpredict extremely high values of pN2O sw (Fig. S15b, e, h, k), but this problem does not persist in the Southern Ocean data, which span a much smaller range of pN2O sw (Fig. S15c, f, i, l). To expand this analysis to global scales, more high-quality training data are needed from oxygen-poor upwelling regions where these high values of pN2O are found.

Application of random forest models to BGC-Argo data and associated uncertainties
Each of the four final trained models was applied to the surface BGC-Argo dataset in the Southern Ocean to predict pN2O sw from the temperature, salinity, dissolved oxygen, and nitrate measured by the BGC-Argo floats. This produced four datasets of predicted pN2O sw based on the four random forest models.

To generate a true prediction interval around model predictions for unseen points, it is necessary to consider uncertainty due to training dataset noise (aleatoric uncertainty) in addition to model uncertainty due to limited training data and/or model architecture (epistemic uncertainty). Based on tests with injected artificial noise (Supplementary Text S3), we approximate uncertainty due to training dataset noise with Mean Absolute Error, such that

Δ2 residual = 1/4 * Σ i=1 to 4 (MAE i) 2 (5)

where Δ2 residual is the residual variance of the training data, or the uncertainty due to the inherent “noise” or randomness of the data themselves, and MAE i is the mean absolute error of model i for the test dataset.

Because a Random Forest consists of an ensemble of decision trees, each of which excludes a certain number of training data points, it is designed to handle this sort of uncertainty. Techniques such as Quantile Regression Forests67 take advantage of the Random Forest model architecture by calculating not just the mean output of all the decision trees but also different quantiles. We take a simple version of this approach, taking the standard deviation of all of the decision tree predictions as the uncertainty around a given model’s output:

Δ2 trees = 1/4 * Σ i=1 to 4 (σ trees,i) 2 (6)

where σ trees,i is the standard deviation of the 600 decision tree predictions for a given point by model i:

σ2 trees,i = 1/(600−1) * Σ j=1 to 600 (pN2O j,i − pN2O i) 2 (7)

where pN2O j,i is the “leaf,” or prediction by tree j in model i. pN2O i is the mean predicted pN2O by model i for a given point, equal to the mean of all 600 decision trees:

pN2O i = 1/600 * Σ j=1 to 600 pN2O j,i (8)

We found that initializing the Random Forest model with a different set of predictor features produces different final results. Thus, in our analysis, we use four different Random Forest models initialized with four different sets of predictor variables. Because each Random Forest model consists of 600 decision trees, our approach is essentially an ensemble of ensembles. The disagreement between the mean predictions by the four models tends to exceed the disagreement between the decision trees for an individual model. So, as a final source of uncertainty, we include the uncertainty of the model ensemble:

σ2 ensemble = 1/4 * Σ i=1 to 4 (pN2O i − pN2O ensemble) 2 (9)

where pN2O ensemble = 1/4 * Σ i=1 to 4 pN2O i (10)

Together, the combined uncertainty in the predicted pN2O at a given point is calculated from:

Δ2 pN2O = Δ2 residual + Δ2 trees + σ2 ensemble (11)

Where Δ2 residual is the uncertainty due to noisiness in the training data, Δ2 trees is the uncertainty due to each model’s architecture, and σ2 ensemble is the standard deviation of predictions by the four Random Forest models.

The air-sea flux of N2O was calculated according to two wind speed-based parameterizations of air-sea gas exchange (Supplementary Text S1). The first of these parameterizations68 has a quadratic dependence on wind speed that includes bubble-mediated exchange implicitly in the nonlinear component of the wind speed dependence. In the second parameterization69, diffusive and bubble-mediated gas exchange are modeled explicitly with linear and cubic dependences on wind speed, respectively.

Uncertainty in air-sea flux
Uncertainties in individual air-sea fluxes were propagated from 1) uncertainties in wind speed; 2) uncertainties in χN2O atm and pN2O sw:

ΔF2 i = ΔF2 winds,i + ΔF2 N2O,i (12)

where ΔF2 i is the uncertainty in air-sea N2O flux associated with float profile i. To quantify uncertainties due to wind speed, fluxes were calculated with combinations of two gas exchange parameterizations (Supplementary Text S1) and two wind products (ERA557 and NCEP70). The uncertainty due to wind speed was estimated as the standard deviation of fluxes calculated from these four methods:

ΔF winds,i = σ(Fi,W14,ERA5, Fi,L13,ERA5, Fi,W14,NCEP, Fi,L13,NCEP) (13)

where Fi,W14,ERA5 denotes the flux calculated from ref. 68 and ERA5 winds57 for individual float profile i.

Predicted uncertainty in pN2O sw was calculated using a Monte Carlo simulation (Fig. S26), which yielded ΔF N2O,i; we accounted for covariances between different gas exchange parameterizations and wind speed products. For each profile i, 1000 realizations of pN2O sw were sampled from a normal distribution with μ = pN2O and σ = ΔpN2O and 1000 realizations of χN2O atm were sampled from a normal distribution with μ = the mean atmospheric mixing ratio of N2O and a standard deviation derived from the errors published in the NOAA flask dataset35. Then, from these 1000 values, four flux estimates were calculated with different wind product-parameterization pairs. The combined uncertainty incorporates both the Monte Carlo uncertainties within each method (parameterization-wind product pair) and the covariance between methods:

ΔF2 N2O,i = 1/4 * Σ j=1 to 2 Σ k=1 to 2 Cov(F j,i, F k,i) (14)

where Cov(F j,i, F k,i) is the covariance between methods j and k for profile i. The sum is calculated over the 1000 Monte Carlo iterations of the four parameterization-wind product pairs:

Cov(F j,i, F k,i) = 1/(1000−1) * Σ iter=1 to 1000 (F j,i,iter − F j,i) * (F k,i,iter − F k,i) (15)

where F j,i,iter is the flux for profile i, method j, Monte Carlo iteration iter, and F j,i is the mean flux for profile i, method j. It was assumed that the uncertainties in sea ice cover were negligible compared to the uncertainties in the other parameters.

The uncertainty in the net annual flux from the Southern Ocean (ΔFt) was calculated as the sum in quadrature of 1) the random uncertainties in individual fluxes; 2) the standard deviation of the net annual flux calculated from four different parameterization-wind product pairs; and 3) the uncertainty due to sparse sampling, calculated from bootstrapping:

ΔF2 t = ΔF2 t,N2O + σF2 t,winds + ΔF2 t,bootstrap (16)

where ΔF2 t,N2O is the cumulative sum of the individual uncertainties ΔF N2O,i from the Monte Carlo analysis:

ΔF2 t,N2O = Σ i=1 to 28,112 (ΔF N2O,i * A zone * D month) 2 (17)

Where A zone is the area of the zone, and D month is the number of days in the month, in which a profile occurs.

The term σF t,winds is the standard deviation of the cumulative fluxes calculated from each parameterization-wind product pair:

σF t,winds = σ(Ft,W14,ERA5, Ft,W14,NCEP, Ft,L13,ERA5, Ft,L13,NCEP) (18)

where

Ft,W14,ERA5 = Σ i=1 to 28,112 (Fi,W14,ERA5 * A zone * D month) (19)

where Fi,W14,ERA5 denotes the flux calculated from ref. 68 and ERA5 winds57 for individual float profile i.

To evaluate the uncertainty associated with interpolating over locations where no float profiles exist, we performed a bootstrap analysis. Repeated this 1000 times, we randomly selected a subset of the floats, re-integrated the fluxes from those floats in space and time to account for fewer floats per zone and month, and summed the resulting integrated fluxes. We then obtained an uncertainty from the standard deviation of the 1000 bootstrap sums. As expected, this uncertainty increased with decreasing bootstrap sample size (Fig. S26), which yielded ΔF t,bootstrap = 0.06 Tg N/year for 162 floats (50% of the total).

In our estimate of the overall Southern Ocean N2O flux, σF t,winds or the subjectivity in the choice of gas exchange parameterization and wind product, are the leading sources of uncertainty, with much smaller contributions from ΔF2 t,N2O, and ΔF t,bootstrap (Fig. 2).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Comments to manuscript ‘Low-pressure storms drive nitrous oxide emissions in the Southern Ocean’\nThank you for your eOorts in addressing my previous comments. I can see that the\nmanuscript has been substantially improved, and many of the issues that previously\nconfused me are now much clearer.\nHowever, I remain cautious as to whether my major concerns have been fully\nresolved.\n1. Comment on the primary reason of the low mean sea level pressure in the Southern\nOcean.\nIn both the authors’ reply and the revised manuscript, it is argued that the low air\npressure in the Southern Ocean results from frequent storm activity. I am not entirely\nconvinced by this explanation.\n• Although I am not an atmospheric science expert, as I understand it, the main\ncause of the low mean sea level pressure in the Southern Ocean is the large-scale\natmospheric circulation associated with the Hadley-Ferrel-Polar cell system,\nrather than the influence of individual weather-scale storms. Air tends to sink near\n~30°S and over the South Pole, while it rises around ~60°S. Consequently, the\nregion near ~60°S is inherently a low-pressure belt due to this ascending branch of\nthe circulation. The Ferrel cell (and associated westerlies) indeed contributes to\nthe development of mid-latitude cyclones and anticyclones, but the mean low-\npressure zone itself is primarily a product of the large-scale circulation pattern.\n• This also explains why the air pressure does not vary substantially between\nseasons, even though wind speeds change considerably (see Fig. S9 in your\nmanuscript, attached below). While storms can locally reduce air pressure on\nshort timescales, they are unlikely to cause a persistent large-scale reduction in\nmean sea level pressure.\n• Moreover, Fig. 1A (attached below) shows that air pressure at the float locations\ndoes not exhibit any “hotspot” features, but instead follows a clear latitudinal\ngradient. This pattern indicates that the low pressure in the Southern Ocean is part of\na quasi-stationary, large-scale system rather than a synoptic-scale feature.\nFurthermore, regional storm intensity varies across the Southern Ocean, for example,\nthe Atlantic sector experiences weaker storm activity (see the figure from Turner et\nal., 2025, attached below), yet the air pressure shown in Fig. 1A in your paper remains\ncomparably low across sectors at similar frontal zones. This further supports the\ninterpretation that the low pressure is not primarily driven by storm frequency or\nintensity.\n• To be clear, this does not substantially undermine the core findings of the study. I\nsimply believe that the main signal the authors observe (i.e., the diOerence in fluxes\ncalculated using actual air pressure versus 1 atm) is primarily due to the inherent low-\npressure belt in the Southern Ocean, rather than storm activity although I am not an\natmospheric. And of course, high wind coupling with the low-pressure system will\nenhance the N2O outgassing. Storms undoubtedly cause additional local reductions\nin pressure, but their contribution to the overall low-pressure system and their large-\nscale impact across the Southern Ocean require further quantification. Therefore,\nattributing the entire low-pressure pattern to storm activity seems overstated.\n2. Comment on whether the limited float data can represent the entire Southern Ocean\nI appreciate the authors’ eOorts in conducting an uncertainty analysis, which I find\nmethodologically sound. However, I am not convinced that adding uncertainty\nestimates alone can justify that the derived values are representative of the entire\nSouthern Ocean because the bootstrap analysis can still only reflect the flux\nsensitivity to the regions and times that having data.\nIt is evident that the float data are unevenly distributed across both seasons and\nregions, while storm activity and N2O concentrations also vary significantly with\nlocation and time. Therefore, the current dataset may not adequately capture the\nspatial and temporal heterogeneity of the system.\nI suggest that the authors either:\nInterpolate or model the available data to obtain a basin-wide estimate for the\nSouthern Ocean, or\nMore simply, revise the conclusions to clearly state that the estimated absolute flux\nvalues apply only to the regions and seasons where float data are available.\n3. Comment on the bubble eOect on N2O\nIt is valuable to learn from the authors’ analysis that small bubbles contribute only a\nminor fraction to the total bubble-mediated flux. Since the authors mention that\nfuture eddy covariance (EC) studies could provide further insight, I would encourage\nthem to consider the existing EC datasets for CO2, which has a solubility comparable\nto that of N2O.\nSeveral recent modeling studies (e.g., Deike et al. and Yang et al., 2022) have already\nincorporated EC-derived CO2 flux data to constrain bubble-mediated components. In\nthe Southern Ocean, at least four EC CO2 studies are available to my knowledge\n(Butterworth & Miller, 2016; Dong et al., 2024; Landwehr et al., 2018; Yang et al., 2021).\nTherefore, I suggest that the authors either make use of these existing datasets or, at\nminimum, discuss their relevance and implications for interpreting bubble eOects on\nN2O fluxes in the Southern Ocean.", "reviewer_2": "I reviewed a previous version of this manuscript, originally submitted to Nature Geoscience. As can be seen on my\ncomments to that version (within the response letter provided by the authors), my general assessment was positive. Despite\nsome caveats, I argued that this manuscript represents a solid contribution. Overall, the major strength I see in this study, is\nthat it advances our mechanistic understanding of the complex relationships between N2O gradients at the sea-air interface\nand synoptic events driving large-scale wind forcing. This is especially relevant in a region which remains (and will perhaps\nremain for some time ⎼given the constraints associated to our technical capabilities for conducting direct N2O\nmeasurements) largely under-sampled. Importantly, this manuscript describes and quantifies a novel aspect of the total\nvariability of N2O partial pressures and sea-air fluxes, which so far did not receive attention; namely barometric pressure\nvariability. A positive “side-effect” of this study is its potential application to other climate-relevant trace gases, whose direct\nmeasurements in the region are also scarce.\nMost of the criticism I provided with my previous review was focused on missing information/lack of clarity with regards to\nsome of the approaches followed by the authors. I also inquired about the potential role of sea ice dynamics (in particular the\nseasonal variability of sea ice coverage) on the observed N2O sea-air flux densities distribution. Lastly, I inquired about the\nextent at which near-coastal areas in the Southern Ocean (in a way the “Glacier-Ocean Continuum”) are represented in this\nstudy; particularly because these areas are still considered a large source of uncertainty in global ocean N2O emission\nestimates. In my opinion, with the revised version the authors addressed adequately all my concerns (and those from the\nother reviewers). While N2O fluxes across sea-ice-air interfaces is not a process that could be addressed mechanistically\nwith the approach followed in this study, the current version makes it clear where the limitations are, and where future\nstudies could (hopefully will) continue. An improved quantification of the seasonality in N2O sea-air fluxes from near-coastal\necosystems in the Southern Ocean is still elusive, since both direct and indirect observations (with sufficient resolution) are\nextremely challenging to obtain.\nKind regards,\nDamian L. Arévalo-Martínez"}, "author_response": "We thank the editor and reviewers for their thorough evaluation of the manuscript and\nconstructive comments. We have revised the manuscript in response to each of these\npoints, with special emphasis on the comment about whether it is storms or general low\nair pressure which is the underlying control on emissions.\nAs before, this Response to Reviewers file provides complete documentation of the\nchanges made in response to each Reviewer comment. Reviewer comments are shown in\nbold. Author responses are shown in plain text. Quotations from the revised manuscript\nare shown in italics and indented.\n\n[... Author response text regarding Reviewer 1's points on atmospheric circulation, float data representativeness, and bubble effects ...]\n\nWe are grateful to the reviewer for evaluating this manuscript a second time, and thank the\nreviewer for this positive assessment. We will just quickly highlight a few additions to the\nmanuscript that relate to the reviewer’s comments. Firstly, to further evidence the\nrelationships between N2O gradients at the sea-air interface and synoptic events, we\nadded a series of figures to the supplement (Figs. S1-S12, reproduced above) showing that\nhigh N O fluxes tend to be associated with coherent low-pressure centers that also have\n2\nhigh wind speeds. Second, we added more language to the discussion about how we are\nmissing coastal and shelf processes in both the BGC-Argo and shipboard datasets, and\nhow this may be an important source of bias. Finally — and perhaps most importantly —\nwe highlighted the need to continue to make high-quality N O measurements, especially in\n2\ncoastal regions."}]

Decision letter