Bubble Clouds in Coastal Waters and Their Role in Air-Water Gas Exchange of CO 2

Bubbles generated by breaking waves can drive significant gas exchange between the ocean and atmosphere, but the role of bubble-mediated gas transfer in estuaries is unknown. Here, backscatter data from 41 acoustic Doppler current profiler stations was analyzed to assess subsurface bubble distributions in nine estuaries along the U.S. East and Gulf Coast. Wind speed, wind direction, and current velocity were the dominant controls on bubble entrainment, but the relative importance of these physical drivers depended on local geomorphology. Bubble entrainment in high-current or shallow, long-fetch estuaries began at wind speeds <5 m s−1. In deep or fetch-limited estuaries, bubble entrainment was less frequent and generally began at higher wind speeds. Data observed during several storms suggests that episodic bubble-driven gas exchange may be an important component of annual CO2 fluxes in large, shallow estuaries but would be less significant in other coastal systems.


Introduction
Air-water gas exchange is an important component of global biogeochemical cycling, and extensive research has focused on measuring gas transfer between the ocean and atmosphere [1,2].Much less is known about air-water gas exchange in estuaries and coastal waters, despite their significance to global OPEN ACCESS and regional carbon budgets [3][4][5].The underlying mechanisms of gas transfer are influenced by a range of environmental factors that vary between ocean and coastal waters as well as among different types of coastal systems [6,7].Nevertheless, after correcting for gas solubility, wind speed is the only environmental forcing in many widely-used gas transfer dependencies, e.g., [1,[8][9][10][11][12][13]. Single-variable parameterizations contribute to large uncertainty in air-water flux estimates for CO2 and other trace gasses because they overlook important environmental controls on gas exchange [1,[14][15][16][17][18].
Gas exchange between the ocean and atmosphere is a function of the gas transfer velocity and the gradient across the air-water interface.The gas transfer velocity is regulated by boundary layer interactions which can be separated into three fundamental components: (1) smooth-surface exchange in the absence of waves; (2) rough-surface exchange in the presence of non-breaking waves; and (3) bubble-mediated exchange driven primarily by breaking waves [1,19].Wave breaking increases turbulence while generating bubbles that experience enhanced gas transfer from surface tension and hydrostatic forcing.Bubbles can also drive gas exchange through buoyancy effects, bubble dissolution, and aerosol production upon bursting.The importance of these processes in air-water exchange is well-recognized but is not well understood due to the difficulty of measuring bubble-mediated fluxes [9,[20][21][22][23].
Attempts to quantify the role of bubbles in oceanic gas exchange have yielded widely-variable dependencies on wind speed and wave spectra that span four orders of magnitude [21][22][23][24][25][26][27].Uncertainties among these parameterizations are amplified in coastal waters by several compounding factors.Surface turbulence and wave breaking exhibit large spatial heterogeneity linked to wind speed, fetch, water currents, water-column stratification, coastal topography, water quality, and bathymetry [7,14,25,[28][29][30][31][32].The lack of constraint at high wind speeds is another major obstacle in estimating air-water fluxes [14].For example, the range of wind speeds represented in estuary gas-transfer data is approximately half the range of oceanic data, even though storms are frequent in land-ocean convergence zones [33,34].The influence of local controls on bubble entrainment and implications for gas exchange have yet to be examined at a scale that compares different coastal system types.
Acoustic sensing methods have been frequently used to study subsurface bubble plumes formed by breaking waves in deep and open-ocean waters [35][36][37][38][39][40][41][42][43][44][45][46][47][48].Recently, Vagle et al. [49] and Wang et al. [50] deployed acoustic instruments in continental shelf waters to characterize the static bubble cloud that persists below the water surface due to wave breaking during storms.These studies provided new data on bubble plumes that are key to understanding gas-transfer velocities at high wind speeds.Here, I employ similar methods at a course but broad scale to examine the environmental controls on bubble entrainment in nine coastal systems along the U.S. East and Gulf Coast (Figure 1, Table 1).This study stemmed from prior research that was aimed at quantifying air-water CO2 fluxes during storms, and accordingly, implications for gas exchange are discussed with specific consideration of CO2.

Acoustic Measurement of Bubbles
Breaking waves generate plumes of varying bubble concentration, size distribution, and penetration depth [51][52][53][54].Smaller bubbles have a lower buoyancy and thus longer average lifetime than larger bubbles and are more susceptible to sustained entrainment by subsurface turbulence.At higher wave breaking frequencies, bubble plumes overlap to form a continuous bubble cloud that extends to a depth determined by the interaction of local physical forces.Acoustic methods of observing these plume characteristics are based on the assumption of selective attenuation and sound scattering by resonant bubbles in water.The different scales at which these processes can be studied require specific observational methods.Multi-frequency instruments are required to determine bubble size distribution and void fraction due to size-dependent resonance frequencies, while direct resolution of individual waves and bubble plumes requires high-performance, high-frequency SONAR that are specifically designed for the task [48,49,55].The cost and operational demands of these instruments often restrict the extent and time of deployment.This study employs a simplified, low-fidelity approach to estimate mean bubble cloud depth using a large, existing dataset from acoustic Doppler current profilers (ADCPs).and small systems (right scale).10° lines were projected from all stations (shown here at CHB0302 only) for estimation of fetch.Station-specific geomorphology features are listed in Table 1.
Table 1.Geomorphology features of ADCP stations.

ADCP Observations
Echo-intensity data from 41 bottom-mounted, single-frequency ADCPs (Teledyne RD Instruments, Poway, CA, USA) were obtained by request from the National Oceanic and Atmospheric Administration Center for Operational Oceanographic Products and Services, USA (NOAA CO-OPS) (Figure 1, Table 1).Ancillary data for each station was obtained through the NOAA public data access portal [56].ADCPs were deployed as part of current surveys occurring between 1999 and 2011.In following sections, stations are either referenced individually by the station name defined in the NOAA CO-OPS database, or multiple stations are referenced as a group by listing only the common station letters, e.g., CHB99 would represent both CHB9903 and CHB9905 (Figure 1, Table 1).
The frequency and depth bin sizes for internal averaging of data varied from 300 to 1200 kHz and 0.35 to 2.0 m, respectively, depending upon the depth at the deployment site (Table 2).Echo intensity profiles for all ADCPs were recorded as average ensembles of 345 pings over 6-minute intervals along with temperature and pressure sensor data over each interval.Vertical profiles of relative backscatter in decibels (dB) were determined from raw data following the methods of Gostiaux and van Haren [57] and Rossby et al. [58].
The mean depth of the bubble cloud for each ensemble was defined as the depth at which the backscatter profile dropped below a site-specific noise threshold (Figures 2 and 3).For simplicity, I henceforth refer to the mean bubble cloud depth as "bubble depth".The backscatter threshold was determined based on ambient noise levels in the water column when wind speeds were in the lower 10th percentile of the study period distribution.Only backscatter data from 10% below the surface to 50% of the total depth were used to estimate the threshold value, as this excluded potential surface interference while still representing the higher noise level of the upper water column.The threshold was then set at the mean plus three standard deviations of the low-wind backscatter data.
Side-Lobe Interference ADCP data from the top 6% to 10% of the water column are normally rejected when calculating current velocities because the side-lobe echo off the air-water surface may not be adequately suppressed.The resulting signal appears as a peak in the backscatter intensity at a depth proportional to the transducer angle to vertical (Figure 2).However, it is unnecessary to exclude all near-surface bins if the goal is to resolve the one-dimensional depth of threshold exceedance.Here, the effect of instrument tilt was exploited to suppress the relative influence of bins that were contaminated by the side-lobe signal.The depth of the side lobe signal was determined based on geometric corrections in Equation 1: where is the beam-specific depth of the side-lobe peak, is the distance from the transducer to the surface, and , , and are the transducer configuration angle, the roll angle, and the pitch angle to vertical, respectively.Finally, the four beams were integrated into a single vertical backscatter profile using the MATLAB R2012b spaps() smoothing function with each bin value weighted proportional to the side-lobe peak depth (Equation 2) on a primary pass and proportional to the default trapezoidal rule approximation an a secondary pass.
, , exp ⁄ 0.001 where , is the weighted backscatter, , is the observed backscatter, and is the depth of bin i.The constant term in Equation 2 keeps data at from being entirely removed before the second filter pass.Figure 2 shows examples of smoothed profiles for low and high-frequency ADCPs.An example of the time-series representation of the threshold depth is shown for station CHB0304 in Figure 3. Profiles where the shallowest bin was not above the threshold and profiles that were not continuously decreasing between the surface and threshold depth, i.e., contaminated by zooplankton or sediment, were excluded.Subsamples comprising at least 10% of the data at each station were visually inspected to confirm the absence of systematic error in raw and processed data.Infrequent anomalies (<1%) were observed in single-beam backscatter profiles but these data could not be distinguished as either contamination (e.g., fish or boats) or relevant physical processes (e.g., boils or windrows) and were not removed due to their negligible impact.
Current speed and direction for each ensemble were extracted from two bins below the threshold depth.This method allowed for the closest approximation of surface currents while presumably scaling with the wave height such that no pings would be contaminated by wave troughs over the ensemble period.Side-lobe interference could only be reduced in the backscatter data and this '2-bin' method ensured that current speed data were well below the side-lobe contamination depth for all beams at all stations.Signal processing and statistical analyses were performed in MATLAB R2012b.

Meteorological Observations
Hourly winds at each ADCP location were estimated using data from the nearest National Climatic Data Center meteorological [59] or National Data Buoy Center [60] meteorological station.All wind speeds were scaled to a height of 10 m (U10) following Large and Pond [61].Fetch estimates for each station were determined in ArcGIS 10.1 using the NOAA Medium Resolution Digital Vector Shoreline dataset [62].Lines were projected from each station on 10° intervals to correspond with wind vector resolution, and each line terminated where it intersected with the shoreline (Figure 1).Fetch was defined by the orthodromic length of each line in the respective direction with a maximum value set at 25 km.
Bubble depths and current vectors were averaged hourly to align with wind data and reduce the effect of anomalies that were not removed by prior processing.To prevent unequal weighting, hourly intervals were removed entirely if >20% of the ensembles had been previously flagged.Current velocities relative to the wind direction were derived from the hourly current vector with positive values indicating along-wind currents and negative values indicating opposing currents.Absolute current velocities were all positive to prevent bias resulting from water quality differences between ebb and flood tides.The hourly change in the wind vector, ΔU10, was defined as an additional environmental control related to sea state and shear generated by changing winds.

Statistical Models
Environmental controls on bubble depth were examined using multiple linear regression analysis.The regression model used forward-backward stepwise selection to identify first-order independent variables and interactions that significantly affected bubble depth.Selection was based on the sum of squared errors of prediction and was constrained by the Bonferroni rule.The Bonferroni rule sets the P-value-to-enter <0.05/q, where q is the number of predictor variables considered.This constraint allowed the full dataset to be used and ensured selection of the most parsimonious model [63].Evaluation of model performance was based on the fit parameters of adjusted R 2 and root-mean-squared error (RMSE).Standardized coefficients and the percent variance explained were used to assess parameter sensitivity.
Regression analysis was tested on two scales of application.The 'general model' included all data from all stations, while the 'site model' was run at each station individually using data from only that site.Input parameters were U10, ΔU10, current velocity (CR), wind-relative current velocity (CW), fetch (F), wind direction (WD), and bubble depth (Zb).For the general model, bubble depth was also represented as the fraction of water column depth (Zb/Z0), which ensured that the predicted depth would not exceed the total depth at shallow sites.Model performance was tested for the three scenarios shown in Table 3. CW, F, and WD were included only as the interaction terms U10:CW, U10:F, and U10:WD, which represent the product of each predictor pair.
WD was functionally a site-specific, angular proxy variable that accounted for the effects of bathymetry, stratification, and the surrounding topography that were not resolved by the other variables.The angular-linear transformation of U10:WD was nested as a non-linear regression within the model following the equations defined by SenGupta and Ugwuowo [64], Equation 3: where is the intercept, is the regression coefficient for predictor variable , is the amplitude, is the angular frequency, is the wind direction, is the acrophase, is the parameter of skewness, and is the random error component.For stations where the angular relationship was sharply peaked, e.g., due to the placement next to a man-mad structure, cos was replaced by sin , where is the parameter of kurtosis.
Three descriptive approaches were explored as a means to pre-group stations independent of geographic location: k-means, hierarchical clustering, and binomial classification of mean fetch, depth, and current speed at wind speeds <10 m s −1 .Descriptive statistics of wind speed and bubble depth distributions, including the mean, standard deviation, and relevant confidence intervals were similarly investigated as a comparative index for pre-and post-analysis grouping.

Bubble Depth Summary Statistics
The nine study systems represented a broad range of estuary features with variable magnitudes of wind-driven and current-driven forcing over the deployment periods (Tables 1 and 2).Bubble depth distributions were generally similar among stations in the same system, but feature-dependent classification was otherwise highly sensitive to minor changes in selection criteria and correlated poorly with bubble depth statistics.The only consistent trend among stations was a linear decrease in bubble depth from low to moderate U10 (Figure 4).There were two exceptions to this trend: HUR0503, a deep, low-fetch station that showed no bubble entrainment over the deployment period, and PIR0701, a high-current station located in a broad coastal channel.
Bubble depth distributions were influenced predominantly by site-specific geophysical attributes rather than instrument setup or ambient noise levels.For example, CHB and MOB stations shared similar instrument configurations, deployment depths, and observed wind speeds, but bubble depth at the MOB stations showed a much stronger correlation with U10 (Figure 4).On the other hand, large variability and some of the highest bubble depth ranges were observed at CHB and PIR despite that these stations had the lowest potential for measurement error.Bubble depth and distribution statistics showed low sensitivity to manual adjustment of the site-specific threshold value.This low sensitivity is consistent with prior observations in ocean studies, but background noise at most of the 41 stations was greater than reported at ocean sites [35,48,49].The background noise level varied by as much as 35% between intra-system stations, e.g., within CHB and SAB, due to either instrument configuration or spatial and temporal changes in water quality.This suggests that ADCP backscatter thresholds should be based on site-specific noise levels over the entire deployment period, as done in this study, rather than only at the time of deployment.

Storm Events
Several storm events occurred during ADCP deployments including Hurricane Isabel (2003) at CHB03, Hurricanes Dennis and Floyd (1999) at CHB99, and multiple strong frontal systems at LIS, MOB, and SAB stations (Figure 3).As expected, depth and fetch appeared to a have a large influence on bubble entrainment at high winds.Bubble depths during storms were only a few meters in short-fetch areas like SAB and CHB99 stations.For example, CHB99 experienced a constant U10 between 15 and 20 m s −1 for 7 days during Hurricane Dennis (1999), but bubble depths reached no more than 3 m (Figure 3: Days 242-249, Figure 4).Storm-driven bubble depths were much deeper at long-fetch stations, reaching 10 m in the LIS and reaching the sediment at CHB03 and MOB stations (Figures 3 and 4).
The influence of wind direction on bubble plume depth was not limited to fetch alone in some estuaries.MOB and CHB03 show deep bubble entrainment during strong northeast winds and relatively less bubble entrainment during comparable winds from the south and southwest, despite large fetch in both directions (Figures 1, 3, and 4).South winds decrease density stratification in both systems and increase vertical mixing.North winds enhance stratification and deeper bubble depth would at first seem counterintuitive.However, bubble subduction and duration may be enhanced by shear turbulence from the two-directional flow, which can also influence wave steepness and breaking frequency [65].

Regression Analysis
Model performance data are presented in Tables S1 and S2, and parameter sensitivities are shown in Figure 5 and Tables S3 and S4.Preliminary model runs indicated that ΔU10 was insignificant at all but a few stations and contributed little to output variability.Wave age and near-surface shear generated by changes in the wind vector may influence bubble distributions [66], but ΔU10 as it was represented here was a poor predictor of bubble depth and was excluded from final model runs.
General model: Wind speed alone showed poor correlation with bubble depth under scenario 1 (Tables 3 and S1).The general model fit a larger portion of the response variation under scenario 2, which included current-based and fetch-based variables.However, the RMSE indicated that the added complexity and larger dataset did little to improve overall model prediction (Table S1).This result was not surprising given the geophysical differences among stations and the different forcing conditions observed over deployment periods (Tables 1 and 2).S3).
Site model: Performance of the site model progressively increased from scenario 1 to scenario 2 and scenario 3.In this order, the model better explained the variance and reduced the RMSE at nearly every station (Table S1).Model performance was most improved under scenario 2 at several stations with high currents (FEB, FPI0902, PIR) and stations with asymmetric fetch limitation (CHB03, FEB, LIS1027, LIS1029, SAB0805).However, performance did not improve equally at all stations that met these criteria, and the parameter sensitivity could not be inferred a priori based on current speed or fetch characteristics from Tables 1 and 2. For example, LIS1011 and LIS1012 had some of the highest current speeds but were equally or more sensitive to the wind direction (Table 2, Figure 5).Under scenario 3, model performance was most improved at stations where the GIS-based fetch parameter used in scenario 2 was unable to resolve man-made structures and bathymetric influences.MOB1104, MOB1105, and SAB0803 were located near inlets that were surrounded by shallow sand bars, and several LIS and BOS stations were located near small islands and jetties.The U10:WD parameter of scenario 3 better accounted for these features and reflected large directional sensitivity (Tables S1-S4, Figure 5).To test whether increasing the fetch resolution would improve scenario 2 predictions, fetch at station MOB1106 was manually corrected based on satellite imagery to account for nearby structures.At this higher spatial resolution, model performance under scenarios 2 and 3 were approximately equal.Higher-resolution fetch data were also applied to CHB03 stations, but here the improvement was negligible.CHB03 stations experience frequent stratification that is influenced by wind direction and cannot be explained by fetch alone.The deep, high-fetch stations in LIS showed the lowest overall increase in model performance and the lowest sensitivity of bubble depth to U10 (Tables S1-S4, Figure 5).

Geophysical Controls on Bubble Plumes
Empirical data from diverse system types is key to understanding the physical processes that generate and shape bubble plumes.Prior analyses of in situ bubble distributions relative to wind speed are limited to a handful of studies.Thorpe [35] found a decreasing, near-linear trend for the hourly U10-bubble depth relationship at a deep (34 m), fetch-limited coastal site near Oban, UK.The study used higher resolution sampling methods (0.45 m bins sampled every 10 seconds), but the U10-bubble depth data were still highly scattered and resembled that of station LIS1012.Of the 41 stations, LIS1012 was most similar to the Oban station in terms of mean depth, fetch, and current velocity.U10-bubble depth regressions at shelf and oceanic sites have shown relatively good correlation and a well-defined U10 threshold at which deep bubble entrainment begins [37,49,50].LIS1011 is geomorphology similar to these open-water stations and showed a comparable U10-bubble depth relationship.The remaining 39 stations showed fewer similarities to bubble distributions observed in prior studies.
There was less bubble entrainment at the deep, long-fetch stations in the LIS relative to most other stations.A possible explanation is that wave breaking is less frequent at these stations because the wave field is able to reach a fully developed state.Shoaling of larger waves that occur in systems like LIS may also be linked to the deeper bubble depths observed at bordering shallow-water stations (LIS1029, FEB1108) [53,54].Shallow breaking waves, known as spilling waves, appear to be the dominant wave type at stations with limited fetch, because even at high wind speeds bubble depths were less than ~3 m (CHB99, SAB, LIS1038).These data agree with recent evidence that wave height in short-fetch waters remains constant from moderate to high winds, and only the frequency of wave breaking changes with wind speed [52,67].
When comparing bubble depth distributions among stations, it is important to consider that bubble entrainment is affected by physical processes other than wave breaking alone.Studies of Langmuir circulation in lake and oceanic waters have shown that bubbles generated by a small number of breaking waves can be concentrated in convergence regions and can be subducted to depths of >10 m [39,43,68].The depth and duration of bubble plumes in shallow and stratified systems can also be influenced by shear-driven turbulence that is generated at the seabed or across horizontal density surfaces [38].It is unlikely that spilling waves can generate enough turbulence to sustain bubbles at depths far below the surface [40], and in short-fetch, high-current systems like PIR, bubble depth may closer reflect surface turbulence rather than the presence of large plunging waves.Indeed, model output confirmed that strong currents were correlated with deep entrainment of bubbles (Figures 4 and 5).However, bubble depths were generally shallower in areas where currents were high only over short distances, e.g., due to funneling through inlets or man-made structures (CHB9903, FEB1102, FPI0902, MOB1106, SAB0803).
Deep bubble entrainment was observed in shallow, long-fetch estuaries like MOB and CHB, but was intermittent and mostly during storms (Figures 3 and 4).

The Role of Bubbles in Air-Water CO2 Exchange
Studies of surface-water CO2 in coastal waters have increased substantially in recent decades, but there are still relatively few gas-transfer parameterizations to choose from when estimating air-water CO2 fluxes [3,4].Despite the paucity of gas exchange data, large variation in transfer velocities has been observed between estuaries [6,7,[30][31][32][33][34]69].Accordingly, it is critical to apply a parameterization that best represents the relevant geophysical controls on gas exchange.For slightly-soluble gases like CO2, air-water exchange is limited by liquid-phase transport, and injection of bubbles by breaking waves can be a highly effective means of gas exchange [70][71][72][73].Lab and modeling studies have estimated that the contribution of bubble-mediated gas transfer scales with wind speed on the order of U10 3 to U10 6 , while gas transfer across the unbroken water surface scales linearly with U10 [14,21,22,25,26,70].However, bubble-mediated transport processes are complex and depend on the mechanism of bubble formation, the dynamic behavior of bubbles, and time evolution of bubble plumes.Bubbles can enhance gas exchange by increasing the air-water surface area, but also by asymmetric transfer that can supersaturate surface waters relative to the atmosphere.The main drivers of asymmetric exchange are hydrostatic pressure and surface tension, which increase the partial pressure of gas inside the bubble [70,73].The slower the equilibration time of a bubble relative to its time evolution, the larger the effect of asymmetric transfer.Asymmetric transfer of CO2 by bubbles is much smaller than bubble-driven symmetric transfer, but it can be significant over long time scales when the concentration of CO2 in the water (CO2w) and atmosphere (CO2atm) are near equilibrium [73].
Estimates of mean bubble cloud depth presented here provide only anecdotal evidence related to bubble-mediated gas exchange because of the coarse resolution of ADCP data in time and space.For example, a given backscatter signal over an hourly averaging interval could be due to deep, dense bubble plumes that are generated periodically by plunging waves.On the other hand, the same signal could equally be due to a more consistent, but less dense bubble layer sustained by water turbulence or frequent spilling waves.Bubble-mediated gas exchange under these two scenarios could be quite different, and it would therefore be impossible to quantitatively link backscatter profiles to gas exchange.Despite this, the backscatter signal indicates the presence of bubbles in either case, and model results can be used as a crude indicator of sites where bubbles are likely relevant to gas exchange.

Comparison of Coastal Environments
A recent gas tracer study by Ho et al. [11] presented one of the best spatial representations of estuarine gas transfer velocity and included concurrent ADCP measurements.The study area of Ho et al. [11] was approximately 150 km upstream from stations HUR0401 and HUR0502 in a narrow section of the Hudson River Estuary that is surrounded by steep topography.ADCP data from HUR and from Ho et al. [11] show an absence of wave breaking at wind speeds observed during either study period.Ho et al. [11] found one of the lowest U10-gas transfer dependencies observed in any estuary.Our data from HUR stations suggest that this low dependency may be due in part to the lack of bubble-driven gas exchange.Data from the other 39 stations show that HUR represents the estuary type where bubbles are least significant to CO2 fluxes.
Bubbles appeared to be a more important pathway for gas exchange in high-current or shallow, long-fetch estuaries like BOS, CHB, and PIR, where bubble entrainment was frequent and began at relatively low wind speeds (Figure 4).Statistical analysis confirmed that wind direction and current velocity can have a significant impact on bubble depth, and prior estuary studies have directly linked these parameters to gas exchange under low-to-moderate wave conditions [6,7,[30][31][32][33][34]69].Estuarine flux data that capture bubble-mediated gas exchange are sparse due to the inherent difficulty of sampling in rough conditions.Floating dome and gradient flux equipment that are commonly used in gas exchange studies are limited to an absence or low-frequency of breaking waves [6,30,31,69].A few shallow-water dome studies have noted the presence of spatially-limited wave breaking [74,75].These studies found some of the highest gas-transfer velocities yet measured, but the relative influence of wave breaking and bottom-generated turbulence could not be differentiated.Whitecap coverage has been suggested as a more practical option for estimating bubble-mediated exchange [1,14,27].Callaghan et al. [76] found that whitecap coverage was influenced by wind speed, wind direction, tidal currents, and swell waves at a shallow coastal site exposed to the North Atlantic Ocean.Similar sensitivities were found in this study at stations with large exposure angles to the ocean (BOS1106, BOS1108, BOS1111, FEB1102, FEB1108, PIR0701), and the influence of swell may help explain the higher RMSE at these stations (Table S1).However, it should be noted that assumptions about gas exchange based on broad classification criteria are hazardous.Several methods to pre-group systems and individual stations based on geomorphology and forcing conditions from Tables 1 and 2 were tested, but no significant correlation with the observed bubble distributions was found.Instead, results have been compared in generalized terms like long-and short-fetch estuaries.This type of feature-based classification can be notoriously subjective, e.g., Table 2 [77], and a more robust method is still needed to extrapolate results to other estuaries.

Ambient Water Conditions
Scattered bubble depth distributions and different parameter sensitivities between sites may also be linked to water quality.Estuaries are often characterized by large gradients in water temperature, salinity, and the presence of surfactants, which can influence bubble plume behavior and gas exchange across the bubble surface [14,[70][71][72]78].Surfactants act as an additional surface layer through which gas must diffuse.For individual bubbles, this would reduce symmetric exchange but concurrently enhance asymmetric exchange due to the longer time required for equilibration.Surfactants also affect wave breaking and bubble plume characteristics.Asher et al. [72] showed that the presence of surfactants in a wave-breaking simulation tank decreased CO2 flux in seawater but increased CO2 flux in freshwater.The authors attributed this result to greater bubble formation in surfactant-rich freshwater, but note that the effect of surfactants on gas exchange in natural environments remains uncertain.The role of bubbles in air-water CO2 flux also depends on CO2w.When CO2w and CO2air are near solubility equilibrium, asymmetric processes can dominate.Zhang [73] recently estimated that supersaturation of surface waters from asymmetric CO2 flux accounts for over 20% of the annual oceanic CO2 uptake.However, asymmetric processes are negligible when the difference between CO2w and CO2air is large.Coastal waters are highly productive and show large spatial and temporal variability in CO2w, e.g., [79].Thus, the most likely role of bubble-mediated gas exchange in estuaries would be to ventilate surface waters where CO2w is driven by biological activity.

Storm-Driven Fluxes
Figure 4 illustrates how bubbles can be mixed throughout the entire water column of shallow, long-fetch estuaries during storm events.In frequently-stratified systems like CHB, the mixture of bubbles with high-CO2 bottom water and porewater could release a substantial quantity of CO2 to the atmosphere [80].Hence, episodic perturbations by storms can have a large impact on estuarine gas fluxes relative to the short time scales on which they occur.Anecdotal evidence from prior studies supports this hypothesis.Sediment resuspension and intense vertical mixing have been observed during several storm events in CHB and MOB [81][82][83][84][85][86].Akulichev and Bulanov [55] found that bubbles reached the seafloor during a storm in the shallow (10 m) Sea of Japan and Evans et al. [87] showed that the elevated CO2 in the Columbia River Estuary, USA was almost entirely ventilated during a major storm event.Crosswell et al. [80] estimated that in a shallow, seasonally-stratified estuary, the CO2 efflux due to a single tropical cyclone (Hurricane Irene, 2011) could equal the net air-water CO2 transport over several cyclone-free years.Storm-generated bubble plumes had a smaller impact in deep or fetch-limited systems, e.g., LIS and CHB99 (Figure 4).In large, deep systems like LIS, storms may be less significant to air-water CO2 fluxes compared to the seasonal effect of convective mixing [88].Alternatively, bubble-driven CO2 exchange may be an important component of annual CO2 fluxes in shallow systems where storms are frequent, e.g., MOB.

Conclusions
This study used acoustic methods of bubble detection developed in prior research [35,41,[48][49][50]57] to compile one of the most extensive datasets of its kind.Collectively, these data amount to six years (2160 deployment days) of ADCP measurements that were gathered during two decades of NOAA current surveys.The dominant environmental controls on bubble entrainment were identified using statistical analysis and used to compare a diverse range of coastal environments.The inclusion of easily-measurable parameters into a site-specific regression model significantly improved model performance, whereas a widely-applicable model showed less improvement and poor overall prediction capability.When considering currently available flux data, parameterizations based on estuary features may be a better option than combining data from multiple estuary types into a single parameterization.I caution however, that system classification may not be intuitive.
The influence of wind speed, fetch, current velocity, and wind direction on bubble depth were site-specific, but general trends were observed within estuaries and between similar sites.Differences between stations were most pronounced during storms.For example, within the same U10 range, there was no detectable bubble entrainment as some sites, while bubbles were dispersed throughout the entire water column at others.It has been argued that bubbles and the forces that generate them can substantially increase CO2 exchange.If so, the different bubble distributions between estuaries may imply differences in biogeochemical cycling.The role of episodic gas exchange via bubbles appeared to be minor relative to tidal forcing in high-current estuaries and relative to seasonal forcing in deep estuaries.However, storm-driven bubble entrainment could have a much larger impact on carbon cycling in broad, shallow estuaries.I emphasize that these conclusions are preliminary and highly speculative, as coarse resolution of bubble depth alone can hardly be linked to gas exchange.Additionally, the effect of water quality on bubble exchange processes is not well understood, particularly at the high surfactant concentrations and variable water conditions found in estuaries.Future research efforts focused on high-resolution, multi-frequency acoustic measurements are needed to characterize bubble plume dynamics on a scale at which they occur.Gas exchange studies that capture a broader range of estuary types and forcing conditions are needed to improve estimates of bubble-mediated CO2 flux.

Figure 1 .
Figure 1.Acoustic Doppler current profilers (ADCP) deployment locations showing National Oceanic and Atmospheric Administration Center for Operational Oceanographic Products and Services, USA (NOAA CO-OPS) station names in large systems (left scale) and small systems (right scale).10° lines were projected from all stations (shown here at

Figure 5 .
Figure 5. Sensitivity of bubble depth to wind speed (U10), current velocity (CR), and wind direction (U10:WD), represented by standardized coefficients under scenario 3. Wind-relative current velocity (CW) was excluded, as it was insignificant at nearly all stations (TableS3).

Table 2 .
ADCP configuration and parameter statistics.Current velocity is the average of the daily mean and maximum values over the deployment period.All others averages are based on hourly data.
* ΔU 10 excluded based on preliminary model assessment.