Quantifying air-sea gas exchange using noble gases in a coastal upwelling zone

The diffusive and bubble-mediated components of air-sea gas exchange can be quantified separately using time-series measurements of a suite of dissolved inert gases. We have evaluated the performance of four published air-sea gas exchange parameterizations using a five-day time-series of dissolved He, Ne, Ar, Kr, and Xe concentration in Monterey Bay, CA. We constructed a vertical model including surface air-sea gas exchange and vertical diffusion. Diffusivity was measured throughout the cruise from profiles of turbulent microstructure. We corrected the mixed layer gas concentrations for an upwelling event that occurred partway through the cruise. All tested parameterizations gave similar results for Ar, Kr, and Xe; their air-sea fluxes were dominated by diffusive gas exchange during our study. For He and Ne, which are less soluble, and therefore more sensitive to differences in the treatment of bubble-mediated exchange, the parameterizations gave widely different results with respect to the net gas exchange flux and the bubble flux. This study demonstrates the value of using a suite of inert gases, especially the lower solubility ones, to parameterize air-sea gas exchange.


Introduction
Noble gases dissolved in seawater are biologically and chemically inert, making them excellent tracers of numerous physical processes that control gas saturation states in the ocean (e.g., bubble-mediated and diffusive gas exchange, temperature change, atmospheric pressure change, ice melting, diapycnal mixing, deepwater formation and ventilation) [1][2][3]. By simultaneously measuring several inert gases with a range of physical properties (e.g., diffusivity, solubility, and dependence of solubility on temperature), we can separately quantify many of these physical processes [4,5]. For example, the solubility of the noble gases decreases with decreasing atomic mass, and the lower solubility noble gases (He and Ne) have a larger portion of their total gas exchange flux driven by bubbles. The dependence of bubble flux on solubility can be explained as follows: for lower solubility gases the atmospheric concentration is high relative to the water concentration, and thus when air bubbles dissolve in the water, the bubbles will generate a larger percent increase in the gas concentration, compared to a higher solubility gas. Additionally, the saturation state of the higher solubility (heavier) noble gases changes more dramatically in response to surface heating/cooling, and mixing of water masses with different temperatures, due to the stronger dependence of the solubility on temperature compared to the lower solubility (lighter) noble gases. Parameterizations of physical processes from inert gases can be applied to bioactive gases, to obtain more accurate estimates of processes including biological productivity [3,6] and denitrification [7,8].
Here, we use a quasi-Lagrangian time-series of the five stable noble gases (He, Ne, Ar, Kr, and Xe) in Monterey Bay, CA, along with a vertical model, to evaluate the performance of four gas exchange parameterizations that differ in their treatment of diffusive and bubble-mediated air-sea gas exchange. These data are some of the few published measurements of a suite of noble gases in a coastal upwelling zone, where environmental factors such as increased levels of surfactants and reduced fetch for wind wave generation may cause the relationship between wind speed and gas exchange rates to differ from the relationship in the open ocean [9,10]. In a forthcoming publication, we will apply the noble gasbased parameterizations of physical processes to O 2 concentration and isotope measurements, in order to quantify net community production and gross primary production during the time-series.

Cruise description
We participated in a seven-day quasi-Lagrangian cruise in Monterey Bay, CA, USA, on the R/V Western Flyer, Sept 27-Oct 3, 2014. During the cruise, rosette casts to ∼180 m were conducted four times per day at roughly 06:00, 12:00, 18:00, and 00:00 local time. Immediately prior to nearly every cast, vertical profiles of microscale turbulence to ∼70 m were obtained using a vertical microstructure profiler (VMP-200, Rockland Scientific). Rates of turbulent kinetic energy dissipation were calculated using two perpendicular air-foil type shear probes on the VMP-200. The Nasmyth spectrum was used to recursively estimate dissipation and diapycnal diffusivity (K z ) [11]. Roughly three high-quality turbulence profiles were obtained every 6 hr. All CTD casts were performed and all noble gas samples were collected in the vicinity of mooring M1 (21 km west of Moss Landing, CA), which is maintained by the Monterey Bay Aquarium Research Institute (figure 1). Starting on Sept 28, 12:35 and throughout the rest of the cruise, we selected where to collect CTD and turbulence profiles by following an autonomous underwater vehicle (AUV) drifting at σ t ≈ 25.2 kg m −3 (25-45 m depth), near the base of the thermocline. Thus, the cruise was quasi-Lagrangian, but with respect to the base of the thermocline, rather than the mixed layer. We also collected noble gas samples from one offshore cast (37 km west of Moss Landing, CA), which indicated that there is minimal spatial variability in noble gas distributions, with respect to both saturation anomaly and concentration.

Noble gas data
In this paper, we present noble gas data (He, Ne, Ar, Kr, and Xe concentrations) from discrete samples, which were collected from the ship's underway seawater line and from Niskin bottles. The samples were collected in copper tubes, sealed with a cold pressure welder, and extracted in the Isotope Geochemistry Facility at Woods Hole Oceanographic Institution [12]. Noble gas abundances were measured on a pulse counting quadrupole mass spectrometer using the system described in   [13]. The dissolved concentration of each gas (in cm 3 ST P g −1 or mol kg −1 ) was determined, using selected ion monitoring [13] and including isotope dilution analysis for Kr and Xe. The estimated error (combined precision and accuracy) is ±0.27, 0.27, 0.24, 0.25, and 0.27 %RSD (relative standard deviation) for He, Ne, Ar, Kr, and Xe, respectively. As we did not analyze samples in replicate from this dataset, the errors are based on a different dataset of duplicate samples that were collected and analyzed by the same methods immediately following this dataset. Three samples with evidence of air contamination (based on anomalies in the relative excess concentration of the five noble gases in the sample, compared to the other samples at a similar depth) were eliminated from the dataset. Because the rate of gas exchange is dependent on the deviation of the gas concentration from equilibrium, accurate solubility functions are critical for accurate parameterization of air-sea gas exchange and other physical processes. Hamme and Emerson (2004) have published high-quality solubility data for Ne and Ar (errors of ±0.30 and 0.13%, respectively) based on equilibration of water with air at a range of temperatures and salinities [14]. For the other gases in seawater, researchers have typically used the solubility data of Weiss (1971) for He, Weiss and Kyser (1978) for Kr, and Wood and Caputi (1966) for Xe [15][16][17]. Errors of ∼2% may be present in these solubility functions because they were determined for 1 atm of each pure gas, and these data must be extrapolated over several orders of magnitude to calculate the solubility with respect to air. For example, Weiss (1971) determined the solubility of Ne in addition to He, and his reported seawater Ne solubility was ∼1.5% lower than Hamme and Emerson (2004) [14,15]. Additionally, the published Kr and Xe solubilities in seawater are subject to uncertainties in the atmospheric mole fraction of each gas [18]. Therefore, for He, Kr, and Xe, we use currently unpublished solubility data from Lott and Jenkins (personal communication, 2015), who determined the solubility of He, Kr, and Xe from 0 to 36.6 PSS, and from <1 to 35 • C by equilibration of water with air. An additional advantage of the Lott and Jenkins dataset is that the solubilities were all measured on the same samples, using the same instrument. Thus the He, Kr, and Xe solubilities of Lott and Jenkins should be more consistent internally compared to data compiled from three different papers (two different laboratories). At the typical sea surface conditions during our study (S = 34.4 PSS, T = 16 • C), the solubilities of Lott and Jenkins are 2.2% greater for He, 1.3% greater for Kr, and 0.1% less for Xe, compared to the published solubilities.
In the main paper, we show the data and model results calculated with the Hamme and Emerson (2004) solubilities for Ne and Ar [14], and the Lott and Jenkins solubilities for He, Kr, and Xe. In the Supporting Information, we show the results using the published solubilities for He, Kr, and Xe, and we provide the noble gas concentration data along with ancillary data (salinity, temperature, atmospheric pressure, etc.) so that it can be used by future investigators with the most accurate solubility functions available. The choice of solubility functions, however, does not affect our overall conclusions; all models simulate the He data better when the Lott and Jenkins solubility is used instead of Weiss (1971) [15].
For the gas molecular diffusivity, which also enters into the gas exchange parameterizations, we used the freshwater data of Jähne (1987) [19] adjusted by -0.138% per ppt salinity (-4.75% adjustment for 34.4 PSS, the typical salinity in our study) for He, Ne, Kr, and Xe, and extrapolated values for Ar [20].

Data input and model setup
We developed a 1D vertical model for the time-series including diffusion and gas exchange for each noble gas. We ran the model using four different gas exchange parameterizations and evaluated the ability of each parameterization to accurately simulate the mixed layer saturation anomalies and concentrations measured throughout the cruise. The model was initialized with CTD profiles of temperature and salinity, and idealized noble gas profiles, described in more detail below. The model run proceeded as follows: at each time step (20 s), for each gas, the air-sea gas exchange flux was calculated and then the gas concentration in the surface box (upper 1 m) was adjusted, based on this flux. Then the diffusive flux at the edge of each box (1 m spacing) was calculated, which mixed the change in gas concentration from gas exchange throughout the mixed layer and deeper into the water column. We used a Crank-Nicolson diffusion scheme to calculate the flux between each box, based on the vertical diffusivity rates measured every 6 hr. The vertical diffusion rate was variable with depth and with time, and the model did not explicitly prescribe a mixed layer depth. Measured diffusivity was generally orders of magnitude higher near the surface where gases are mixed much more vigorously, compared to deeper waters below the mixed layer. Our vertical diffusivity determined from microstructure profiles is more accurate compared to other published studies where the diffusivity was estimated for the base of the mixed layer only. We ran the model from Sept 28, 00:00-Oct 3, 02:00 local time. We linearly interpolated all data to the model time step of 20 s. CTD profiles of temperature and salinity were binned to the depth grid spacing of 1 m, from 0.5-69.5 m, and then linearly interpolated between casts. Temperature and salinity profiles from each CTD downcast (every 6 hr) were used in the model for calculating the simulated gas Measurements of wind speed, atmospheric pressure, and relative humidity were obtained at 10-min intervals from sensors on mooring M1 (figure 1). There were two wind sensors on the mooring: a Vaisala ultrasonic anemometer and a RM Young propeller anemometer. For this study we used data from the ultrasonic anemometer because it was specified to have higher accuracy and a lower detection limit. The wind sensors were both at 4 m height above the sea surface and the wind speeds were extrapolated to u 10 , the wind speed at 10 m height, using the equation u b = u a (b/a) 0.11 with a and b the sensor and extrapolated heights, respectively [21]. During the time period included in the model (excluding the upwelling event, discussed below), the root mean square deviation (RMSD) of the two sensors was 0.38 m s −1 (9.8%). The mean difference between the two sensors was 0.29 m s −1 , with the anemometer giving higher wind speed 86% of the time. The model results were very similar regardless of which anemometer dataset was used. For example, the net gas exchange flux typically differed by 13% or less when the model was run with data from either anemometer.
The gas saturation anomaly was used in all the gas exchange parameterizations, and is defined in ∆ notation as with C w and C eq the water and saturation equilibrium gas concentrations, respectively. The equilibrium concentration of each gas at each time step was calculated with respect to local sea level pressure and humidity, specifically where C eq,re f is the reference equilibrium concentration at the measured salinity and temperature for 1 atm total pressure of air with water vapor at 100% humidity. The pressure terms are all expressed in atm: P re f = 1 atm is the total reference pressure of air, P sl p is the local sea level pressure, and P H2O and P H2O,eq are the water vapor pressures in situ and at equilibrium (100% relative humidity), respectively, calculated from the salinity, temperature, and relative humidity [22]. The model was initialized with an idealized profile for each gas (figure 2), because we did not have a depth profile at the start of the cruise. The concentration was fit to be similar to the other profiles at depth. The initial surface concentration of each gas was determined by fitting the model to either the concentration or saturation anomaly of the first sample, which was collected on Sept 28, 06:15. The fit was determined by minimizing the RMSD between the sample and model results for the four parameterizations. Fitting the models to the initial sample concentration versus saturation anomaly gave slightly different results because the saturation anomaly of each discrete sample was calculated based on the temperature and salinity of the water sampled, rather than the modeled temperature and salinity, which was interpolated from the CTD downcasts. There were sometimes differences between the realtime underway and interpolated CTD profiles (for underway samples) and offsets in the temperature and salinity profiles between the upcast and downcast (for Niskin samples). Therefore, the model and sample equilibrium concentrations were sometimes different. Below, we separately report the performance of the models optimized to concentration and saturation anomaly (section 4). Our diffusion-based model was not able to reproduce the upwelling event. Therefore, we manually reset the gas profiles at Sept 30, 06:16, the time of the first cast after upwelling. The saturation anomaly of each gas throughout the upper 25 m was reset to the same value for all models, and    . We fit the model to either the mean concentration or mean saturation anomaly of these two samples. This approach was identical to the approach used to determine the initial surface concentrations. We excluded samples collected during the upwelling event ( Sept 29, 19:33 and Sept 30, 01:30) when assessing the model performance.

Adjustment for upwelling
We did not add a continuous upwelling flux throughout the cruise (e.g., using a published upwelling index for the region). Upwelling events are episodic in this region during the fall [23], and our temperature and mooring data indicated that only one upwelling event occurred during our time-series, overnight Sept 29-30. Additionally, we were able to simulate the remainder of the time-series without a continuous upwelling flux.
We omitted advection from the mass balance because we did not have sufficient spatial coverage to resolve advective fluxes. This omission is likely not a serious limitation because our measurements at an offshore station suggested that the dissolved gas concentrations and saturation anomalies were very similar offshore to near M1 (figure 2). Furthermore, we evaluate the performance of the models using only He and Ne, which should be the least affected by advective fluxes and mixing between different water masses. These gases will have the lowest lateral variability in concentration and saturation anomaly, due to the very small temperature dependence of their solubilities (0.2 and 0.7% • C −1 for He and Ne, respectively) [2,5,24]. Additionally, our data show that the surface He concentration was barely altered by the upwelling event, a clearly non-Lagrangian process.  [2,3,25,26]. N11 and S09 were developed from inverse models of noble gas data, with S09 using a three-year timeseries of monthly profiles of He, Ne, Ar, Kr, and Xe concentration at the Bermuda Atlantic Time-series Study (BATS) site in tandem with a 1D oceanic mixed layer model, and N11 using a global dataset of Ne, Ar, N 2 /Ar, and Kr/Ar measurements in tandem with an ocean general circulation model. L13 used a large eddy simulation (LES) model coupled to a bubble population model to parameterize the bubblemediated exchange and did not use oceanic gas measurements to tune the model; however, the L13 parameterization reproduced oceanic data well. With all three parameterizations, the total gas exchange flux is calculated as the sum of diffusive and bubble-mediated exchange:

Choice of gas exchange parameterizations
with F t the total flux, F d the diffusive flux, and F b the bubble flux. Each parameterization further separates the bubble-mediated flux into two components: where F c is the complete trapping bubble flux from typically smaller bubbles that completely dissolve and F p is the partial trapping bubble flux from typically larger bubbles that partially dissolve. The exact equations for F d , F c , and F p differ between authors, as we discuss below. Finally, Sw07 used a global dataset of dissolved inorganic radiocarbon data to quantify the rate of uptake of anthropogenic (bomb) radiocarbon into the ocean. Sw07 does not explicitly include the bubblemediated component of air-sea gas exchange, i.e., it includes F d but not F b . The Sw07 parameterization is within the range of other gas exchange parameterizations derived from oceanic measurements that do not include a separate term for bubble-mediated exchange [27][28][29][30], and is used as F d in N11.
All of the parameterizations use u 10 as the only environmental forcing variable. In L13, u 10 is converted to u * using an empirical formula for the drag coefficient [31,32]. Therefore, none of the models explicitly include variables such as surfactants, precipitation, and/or fetch. S09, Sw07, and N11 all parameterize F d using an equation of the form Here A is an empirical constant and Sc is the Schmidt number of the gas at the water temperature and salinity. Sc is defined as the ratio of the kinematic viscosity of seawater to the gas diffusivity. N11 and Sw07 both use A = 0.27 cm hr −1 , whereas S09 uses A = 0.30 cm hr −1 , and will thus yield a diffusive flux that is 11% greater for a given wind speed and gas concentration. The comparison to F d calculated from L13 is more complex because u * is a nonlinear function of u 10 . At u 10 = 5 m s −1 , F d of L13 is approximately equal to F d of S09 and N11. At lower wind speeds, F d of L13 exceeds S09 and N11 (by a factor of five at u 10 = 1 m s −1 ), and at higher wind speeds F d of L13 is less than S09 and N11 (by a factor of two at u 10 = 10 m s −1 ). The differences in F d for L13 (scaled to u * ) versus the other parameterizations (scaled to u 10 ) also have a slight dependence on the gas Schmidt number. For the bubble-mediated fluxes F c and F p , the differences between parameterizations are not simple to characterize because they are functions of additional variables such as the gas saturation anomaly, diffusivity, solubility, and atmospheric mole fraction, as well as u 10 . The dependence on these factors can differ between models. For example, F p is scaled to the gas diffusivity, D, as F p ∝ D 0.5 in N11 and F p ∝ D 2/3 in S09. L13 sets F p ∝ (Sc/660) −2/3 . All three parameterizations assume that F c is independent of D.
We report the differences between parameterizations for F d , F c , and F p in our model below. We refer the reader to the authors' original papers for more details on each parameterization [2,3,25,26]. In the Supporting Information, we provide MATLAB functions for calculating gas exchange fluxes using each of the four tested parameterizations, as well as several additional parameterizations that only include F d [20].

Results
We evaluated the performance of each parameterization by comparing the model results to near-surface measurements (2-4 m depth, 15 samples total). We show the gas time-series with respect to saturation anomaly as well as with respect to concentration (figure 4), and we calculated the error using both parameters. For the most soluble gases, Ar, Kr, and Xe, all models gave very similar results and simulated the surface measurements well. The four models diverged for He and Ne, which have higher diffusivity and lower solubility, making them more sensitive to differences in the treatment of bubble-mediated gas exchange. The concentration plots for Ar, Kr, and Xe clearly show the impact of the upwelling event overnight on Sept 29-30: near-surface gas concentrations following the upwelling were 2-4% higher than before the upwelling (figure 4c-e). The concentration increase for Ar, Kr, and Xe following the upwelling was driven by the replacement of warmer water by colder water with higher equilibrium gas concentrations (but similar saturation anomalies for each gas). In contrast, the concentrations of He and Ne changed by <1 % due to the upwelling because these gases' solubility has a lower dependence on temperature and therefore the equilibrium gas concentration was similar before and after upwelling. Additionally, because the upwelling event coincided with high winds, some of the concentration increase for He and Ne during this period can be attributed to bubble-driven supersaturation rather than temperature change.
There was one sample following the upwelling, shaded grey in figure 4, that fell below the model results for Ar, Kr, and Xe concentration, but not saturation anomaly. This sample was collected in warmer water (16.4 • C) compared to the four surface samples taken within 24 hr of that sample (15.2-15.7 • C) and therefore it had a lower equilibrium gas concentration (C eq ). It is possible that this sample may have been from a different water mass compared to the rest of the time-series, but it was included in all of our error estimates because we had no specific reason to discount it.
For each gas, the initial surface concentration and surface concentration after upwelling were set by fitting the models to the measured concentrations of the first sample and the first two samples after upwelling, respectively. Due to the measurement uncertainty for these initialization values, we performed a Monte Carlo error analysis where the initial concentration and concentration after upwelling were simultaneously varied randomly 100 times with a normal distribution, with the optimized concentrations (figure 4) set as the mean and the measurement error set as the standard deviation. We evaluated the performance of the different parameterizations by calculating the RMSD between all measured nearsurface sample concentrations and the coincident modeled concentrations, for each noble gas and each parameterization. We then repeated the error analysis, but instead fit the models to the measured saturation anomaly of the first sample, and the first two samples after upwelling, and calculated the RMSD with respect to saturation anomaly. We report the results for each parameterization as the mean ± 09/28 09/ 30 (table 1). Unfortunately, the fact that the highest wind speeds in our time-series coincided with the upwelling event made it more difficult to discriminate between parameterizations. If the models are run through the full cruise duration, without resetting the concentrations after the upwelling, the four gas exchange parameterizations predict substantially different trajectories for He and Ne during and following the high wind event (figure 5). As wind speed increases, the magnitudes of F d , F p , and F c increase and so do the differences between models. Below, we discuss the advantages and disadvantages of including upwelling versus not including upwelling in our model, and the insights gained from both approaches.  Figure 5. Model results without upwelling adjustment. Near-surface gas concentrations (a-c) and saturation anomalies (d-f) from samples and model results that were optimized to give the best fit to the first sample. All plotted samples were included in the error analysis.

Modeled surface concentrations and saturation anomalies
When the model is adjusted for upwelling, the L13 model has the most skill in simulating the surface He and Ne data, yielding the lowest RMSDs of the four parameterizations (table 1). L13 is the most accurate parameterization for predicting surface He concentration and saturation 65% and 64% of the time, respectively. Additionally, L13 is the most accurate parameterization for predicting surface Ne concentration and saturation 56% and 58% of the time, respectively. N11 is most accurate 25-33% of the time for He and Ne concentration and saturation. The S09 and Sw07 models were each the most accurate parameterization in 10% or less of the Monte Carlo simulations for He and Ne. In general, the RMSD for He and Ne concentration and saturation is similar for L13, N11, and Sw07, and higher for S09. The modeled concentrations and saturation anomalies for He and Ne, when including upwelling, looked very similar for L13 and Sw07, although as shown in table 1, L13 was usually the most accurate parameterization. However, a large part of the apparent skill of Sw07 at reproducing the observed He data comes from the upwelling adjustment. The gas added by upwelling compensates for the lack of gas added by an explicit bubble-generated flux. For He, when the model is run without adjustment for upwelling, the differences between parameterizations, and the need for an explicit bubble-generated flux to produce the observed He supersaturation become clearer (figure 5). Without the upwelling adjustment, the Sw07 parameterization predicts lower He concentrations and saturation anomalies than the other parameterizations, especially toward the end of the time-series. If the time-series were extended, the He concentration in Sw07 would continue to decay toward equilibrium (He eq ), and Sw07 would therefore underestimate the true He concentration. While we do not have a long enough time-series to say for certain, the fact that all but one surface He sample in our time-series are supersaturated suggests that a model that predicts equilibrium concentrations at steady state, such as Sw07, would be incorrect. An explicit bubble flux into the ocean is needed to generate consistently supersaturated surface waters for gases that are insensitive to temperature change. The solubility of He only changes by 0.2% • C −1 .
Without the upwelling adjustment, L13 is most accurate for predicting surface He concentration and saturation anomaly 55% and 64% of the time, and N11 is most accurate 45% and 36% of the time, respectively. The performance of L13 and N11 for He is similar with and without the upwelling adjustment. However, notably, without the upwelling adjustment, the S09 parameterization appears to overestimate the bubble flux for He during the high winds that coincided with the upwelling.
For Ne, when the model is run without an adjustment for upwelling (figure 5), S09 is the most accurate parameterization for Ne concentration and saturation anomaly 100% of the time and has the lowest RMSDs of the four parameterizations, even though S09 has the highest RMSDs for Ne when the Again, if the time-series were extended, the Ne concentration modeled with Sw07 would decay toward equilibrium, which is not consistent with the observations. All of the surface Ne samples are supersaturated, despite the solubility of Ne having a week dependence on temperature (0.7% • C −1 ), which supports the need for an explicit bubble-mediated flux into the ocean to generate some portion of the observed supersaturation. For Ar, Kr, and Xe, all parameterizations simulate the data well after the correction for upwelling. For these gases, the RMSD of each parameterization is very similar and we conclude that all parameterizations have similar skill in simulating the heavier noble gases for our dataset. With a somewhat longer time-series, and/or without the upwelling event, the differences between the parameterizations would likely be clearer, and we might be able to see the importance of bubble-mediated exchange for somewhat more soluble gases, such as Ar. When the model is run without adjustment for upwelling, all four parameterizations predict similar trajectories for Ar, Kr, and Xe and underestimate the concentrations and saturation anomalies of these gases after the upwelling event (figure 5).
These results demonstrate the different factors that controlled the change in concentration for each gas overnight on Sept 29-30. For He, bubble-mediated exchange generated a small concentration increase (<1%). For Ne, upwelling and bubble-mediated exchange were both important in producing a moderate concentration increase (∼1%). For Ar, Kr, and Xe, upwelling/mixing caused a large (2-3%) concentration increase. Because Ar, Kr, and Xe were barely affected by the high wind event that coincided with the upwelling, their concentrations could not be used to constrain the proportion of the surface water that was replaced with colder water with higher gas concentrations.
Over the whole time-series, the surface concentrations of Ar, Kr, and Xe were primarily controlled by upwelling and diffusive gas exchange. These more soluble gases contrast with He, which was primarily controlled by bubble-mediated and diffusive gas exchange, and Ne, which was controlled by a combination of upwelling, bubble-mediated gas exchange, and diffusive gas exchange.   The lower skill of the S09 parameterization compared to N11 and L13 in simulating He (both with and without the upwelling adjustment), and the improved performance of S09 in simulating Ne when the upwelling adjustment is not performed, both suggest that S09 may be overestimating bubble-mediated exchange in this environment. One potential explanation is that the S09 parameterization was tuned using He solubility data that was ∼1% too low (saturation anomalies that were ∼1% too high), compared to the LJ15 solubility used in this model. Assuming that the S09 model was fit to reproduce He saturation anomalies that were too high, it may overestimate bubble fluxes. As we observe, S09 does predict the largest bubble-mediated fluxes of the four parameterizations. When we run the model using the He solubility of Weiss (1971) [15], which is ∼2% lower than the He solubility of LJ15, all models underestimate the near-surface He concentrations and saturation anomalies (figures 2 and 3 in Supporting Information). The S09 parameterization comes closest to simulating the He data when using the Weiss (1971) solubility [15] because it predicts the largest bubble fluxes.
The N11 and L13 parameterizations were fit to inert gas data not including He, and the L13 parameterization was not tuned using any oceanic gas data. Therefore, N11 and L13 would not have errors related to uncertainties in He solubility. This result underscores the need for accurate gas solubility functions in order to realistically interpret oceanic gas data [14].
The modeled air-sea fluxes help to explain why the parameterizations all give very similar trajectories for Xe, Kr, and Ar in figure 4. For Xe, the Sw07 and N11 models predict nearly the same F d because the diffusive flux greatly exceeds the bubble flux (table 2). The absolute and relative values of F c and F p for Xe vary widely between the different parameterizations. L13 is unique among the parameterizations in simulating a bubble stripping effect in which heating-induced supersaturation exceeds the bubble overpressure effect, resulting in a net removal of gas by partial bubble trapping (F p ). The S09 and N11 parameterizations always give net bubble flux into the ocean for F p , even for supersaturated gases. N11, S09, and Sw07 predict similar F d for Xe and have a small or zero bubble flux (the magnitude of F b is 5% or less of the magnitude of F t ). For L13, the sum of F d and F b for Xe gives a net flux that is very similar to F d for the other parameterizations.
The choice of optimizing the initial conditions to concentration versus saturation anomaly does not significantly affect our interpretations (table 1). For example, the modeled total air-sea flux (F t ) for each gas varies by 7% or less between the two optimization values. Additionally, the Monte Carlo error analysis results are very similar when the wind speed from the propeller anemometer is used rather than the sonic anemometer. The modeled F t for each parameterization is generally within 13% for either anemometer, and both wind speed sensors predict L13 is most accurate 64-77% of the time for He and 56-72% of the time for Ne.

Conclusions and future work
This study is complementary to others that have demonstrated the value of using oceanic inert gas measurements in tandem with models to quantify air-sea gas exchange fluxes [2,3,6]. We demonstrated that short-term, high-frequency measurements of inert gases and diapycnal diffusivity can be used to quantify air-sea gas exchange in coastal regions. We found that accurately parameterizing bubblemediated exchange was necessary to simulate the near-surface measurements of He and Ne collected during our time-series. The tested parameterizations gave a wide range of results for the direction and magnitude of the net bubble flux, and the proportion of partial versus complete bubble trapping, indicating that there are still large uncertainties in models of bubble-mediated gas exchange. Higher wind speed conditions, and/or a longer period of Lagrangian observations would have resulted in a greater divergence between the parameterizations for He and Ne. The parameterizations that displayed the most skill in simulating the He observations were Liang 2013 [25] followed by Nicholson 2011 [2]. This result was observed regardless of whether we ran the model with or without an adjustment for upwelling. For Ne, the parameterization of L13 was most accurate when the model was adjusted for upwelling partway through the cruise. For Ar, Kr, and Xe, all parameterizations gave very similar results and simulated the observations well, after adjustment for upwelling. For these heavier gases and the moderate wind speeds observed during our study, we conclude that diffusive exchange driven by temperature change was more important than bubble-mediated fluxes in controlling the gas concentrations of Ar, Kr, and Xe. Due to the complication of upwelling during our cruise, a longer uninterrupted Lagrangian time-series and/or higher average wind speeds would have enabled clearer discrimination between parameterizations for the heavier gases, and an evaluation of the importance of bubble-mediated exchange for Ar.
The four tested parameterizations performed fairly well, despite not explicitly incorporating factors other than wind speed that may affect gas exchange rates in the coastal ocean (e.g., fetch and surfactant concentration). Evaluating these and other gas exchange parameterizations against in situ measurements in a wide range of environmental conditions, and for longer periods of time, is an important goal of future work.