A theoretical study of the effect of subsurface oceanic bubbles on the enhanced aerosol optical depth band over the southern oceans as detected from MODIS and MISR

Submerged oceanic bubbles, which have a much longer life span than whitecaps or bubble rafts, have been hypothesized to increase the water-leaving radiance and thus affect satellite-based estimates of water-leaving radiance to non-trivial levels. This study explores this effect further to determine whether such bubbles are of sufficient magnitude to impact satellite aerosol optical depth (AOD) retrievals through perturbation of the lower boundary conditions. There has been significant discussion in the community regarding the high positive biases in retrieved AODs in many remote ocean regions. In this study, for the first time, the effects of oceanic bubbles on satellite retrievals of AOD are studied by using a linked Second Simulation of a Satellite Signal in the Solar Spectrum (6S) atmospheric and HydroLight oceanic radiative transfer models. The results suggest an insignificant impact on AOD retrievals in regions with near-surface wind speeds of less than 12 m s. However, the impact of bubbles on aerosol retrievals could be on the order of 0.02–0.04 for higher wind conditions within the scope of our simulations (e.g., winds < 20 m s). This bias is propagated to global scales using 1 year of Moderate Resolution Imaging Spectroradiometer (MODIS) and Advanced Microwave Scanning Radiometer EOS (AMSR-E) data to investigate the possible impacts of oceanic bubbles on an enhanced AOD belt observed over the high-latitude southern oceans (also called the enhanced southern oceans anomaly, or ESOA) by some passive satellite sensors. Ultimately, this study is supportive of the null hypothesis: submerged bubbles are not the major contributor to the ESOA feature. This said, as retrievals progress to higher and higher resolutions, such as from airborne platforms, the uniform bubble correction in clean marine conditions should probably be separately accounted for against individual bright whitecaps and bubble rafts.


Introduction
The remote sensing community has long noticed anomalously high aerosol optical depth (AOD) retrievals over high wind belts of the southern oceans, North Pacific, and North Atlantic (e.g., Myhre et al., 2005;Zhang andReid, 2006, 2010;Shi et al., 2011a, b;Smirnov et al., 2011;Toth et al., 2013;Kalashnikova et al., 2013;Chin et al., 2014).Some passive retrievals of AOD from satellites observe a belt of high AOD over the southern oceans known as the enhanced southern oceans anomaly (ESOA) that is especially biased when compared with ship-based measurements of AOD.These anomalously high values are thought to be in part due to a combination of cloud contamination and enhanced radiance from the ocean surface from whitecaps and bubble rafts.Given the size of the oceans, even small but consistent biases can have a significant influence in the overall estimated topof-atmosphere (TOA) radiative forcing by aerosol particles.
The University of North Dakota and the Naval Research Laboratory have been systematically investigating persistent Figure 1.An underwater image of bubbles generated by plunging waves.The picture was taken using an underwater bubble camera system designed to measure the number density of bubbles over a size range of 40-800 µm (Zhang et al., 2004).For reference, the metal wire has a width of 200 µm.oceanic biases in satellite AOD estimates.Early studies first verified that the high oceanic AOD belts were in fact highly biased (Smirnov et al., 2011;Toth et al., 2013).This was then followed by the most logical factor, cloud contamination.Indeed, a series of studies suggests that most of the high bias is related to clouds.However, there is a clear lower boundary condition signal as well, with increasing positive AOD bias with wind speed (e.g., Zhang and Reid, 2006;Shi et al., 2011a).Given that sea salt aerosol production, specular reflection (sun glint), and whitecapping all covary with wind speed, AOD retrievals are a potentially confounded system.Some Level 3 products (e.g., Zhang and Reid, 2008;Shi et al., 2011a) include an empirical correction for wind-speedrelated bias to retrieved AOD.Some Level 2 satellite retrievals (e.g., Sayer et al., 2010Sayer et al., , 2012;;Jackson et al., 2013;Levy et al., 2013;Limbacher and Kahn, 2014) also incorporate wind speed data into the radiative transfer calculations using parameterizations of wind effects on whitecaps and bubble rafts.The current study uses a unique combination of data sets to further investigate the mechanics of the ocean lower boundary condition.
Whitecaps and the resulting bubble rafts form an easily identifiable perturbation to the ocean surface reflectivity.However, there is another consideration: subsurface bubbles (e.g., Fig. 1).While whitecaps last for only seconds, subsurface bubbles can have a much longer lifetime (e.g., Johnson and Cooke, 1981).Theoretically, an air bubble in pure water would either rise to the surface under buoyancy (Harper, 1972) or dissolve under surface tension pressure (Epstein and Plesset, 1950).In open ocean environments, bubbles are found to be coated with organic and surfactant materials (Fox and Herzfeld, 1954;Yount, 1979).The coat-ing process prevents gas diffusion and stabilizes the bubbles against buoyancy (Fox and Herzfeld, 1954;Yount, 1979).While rising bubbles burst at the surface and form whitecaps and bubble rafts, stabilized bubbles can stay in water for hundreds to thousands of seconds (Johnson and Cooke, 1981).Under moderate wind conditions (> 3 m s −1 ) most bubbles near the ocean surface are generated by breaking waves (Thorpe and Humphries, 1980;Thorpe, 1982;Thorpe and Hall, 1983;Lamarre and Melville, 1991).Under mid-to high-wind conditions (wind speed > 7 m s −1 ), observationbased studies have shown a horizontal layer of subsurface bubbles that is formed and maintained by a constant supply of bubbles through breaking waves and turbulence (Crawford and Farmer, 1987;Monahan and Lu, 1990;Thorpe, 1982Thorpe, , 1986)).While whitecaps are clear and obvious, they cover < 10 % of the ocean surface for winds as high as 20 m s −1 (e.g., Monahan et al., 1983).However, there is a broad uniform enhancement of the dark ocean surface albedo due to persistent subsurface bubbles.One need only consider the example of a ship wake, which can exist for longer than 3 min, to see the impact of stable submerged bubble populations on water-leaving radiance (Zhang et al., 2004).However, even under low-wind conditions there can be an enhancement in the bubble population due to rain drops (Pumphrey and Elmore, 1990), melting snow (Blanchard and Woodcock, 1957), biological processes (Medwin, 1970), outgassing of sediments (Mulhearn, 1981), growth from stable cavitation nuclei due to gas supersaturation (Johnson and Cooke, 1981), and supersonic pressure (Messinó et al., 1963).
The goal of this study is to evaluate whether subsurface bubbles pose a lower boundary condition problem for aerosol remote sensing.Already these stabilized bubbles have been recognized as a complicating term in retrievals of waterleaving radiance (Zhang, 2001;Zhang and Lewis, 2002;Flatau et al., 2000).While previous efforts to empirically correct the lower boundary condition of aerosol products for wind implicitly incorporates the entirety of the lower boundary condition-specular reflection, whitecaps, bubble rafts, and submerged bubbles, such methods are neither applicable to joint retrievals of ocean and atmospheric products together nor acceptable for higher resolution retrievals such as those performed by aircraft-mounted sensors.
While water-leaving reflectance is a subcomponent of the ocean surface reflectance (e.g., Koepke, 1984;Vermote et al., 1997), the contribution of subsurface bubbles to waterleaving radiance in relation to other ocean features has yet to be explored in the context of aerosol retrievals.Within a pixel of satellite observation, whitecaps are sporadic and scattered whereas bubbles in water form a more or less uniform layer that could exist over regions that are free from whitecap contamination (e.g., Monahan and Lu, 1990).While whitecaps serve as a diffuse reflector, reflecting solar radiation directly at the surface (Frouin et al., 1996;Whitlock et al., 1982), bubbles interact with light below the surface, enhancing water-leaving radiance (Stramski and Tegowski, 2001;Terrill et al., 2001;Zhang et al., 1998).Studies have shown that the contributions of bubbles to water-leaving radiance are rather significant.For example, Zhang et al. (1998) found that organic coatings on bubbles will significantly enhance the backscattering and proposed that bubbles could be the strongest contributor to the light coming out of ocean.Stramski and Tegowski (2001) illustrated the temporal variation of the light field caused by bubble injection under a wave breaking event (wind speed = 10 m s −1 ) recorded by an acoustic backscatter in coastal waters and found a > 2-fold increase in reflectance (400-700 nm) over time periods on the order of several minutes or less.Clearly, ocean surface reflectance patterns can be altered by immersed bubbles under windy conditions, which may further affect satellite aerosol retrievals from passive sensors.
In this study, through a theoretical approach, the impacts of subsurface oceanic bubbles on satellite aerosol-retrieved AODs are studied, especially over the ESOA region.The effects of oceanic bubbles on satellite AODs are examined theoretically, using a linked oceanic and atmospheric radiative transfer model (RTM).The HydroLight oceanic RTM (Mobley et al., 2012) is used to estimate the bubble-induced perturbations in surface reflectance as a function of near-surface wind speed.The HydroLight simulated bubble concentration and surface reflectance relationship is further incorporated into the Second Simulation of a Satellite Signal in the Solar Spectrum (6S) atmospheric RTM (Vermote et al., 2006) for estimating the impact of bubbles to the TOA radiation.Note that in the blue and green parts of the visible spectrum, it is difficult to separate the contributions of bubbles from the total background reflectance due to multiple scattering (e.g., Zhang, 2001).In the red/infrared spectral ranges with strong absorption due to water molecules multiple scattering is less significant, thus the bubble contributions can be identified (Zhang et al., 2002).Accordingly, in this study, the effects of oceanic bubbles on atmospheric aerosol retrievals are studied at the Moderate Resolution Imaging Spectroradiometer (MODIS) 0.66 µm channel.Next, the modeled TOA radiances from the linked RTMs are validated against radiances observed from MODIS.Lastly, the effects of oceanic bubbles on satellite-derived AOD are evaluated theoretically, using simulations from the linked oceanic and atmospheric RTMs.

Data sets and methodology
In this study, winds derived from Advanced Microwave Scanning Radiometer EOS (AMSR-E), ship-based AOD data from Maritime Aerosol Network (MAN), and MODIS radiances and AOD retrievals were collected and collocated.A 6-year 2004-2009 study period is used.Seven years (2002)(2003)(2004)(2005)(2006)(2007)(2008) of collocated AErosol RObotic NETwork (AERONET) and Aqua MODIS DT AOD data are used to aid the analysis.Two radiative transfer models, the 6S atmo-spheric RTM and HydroLight oceanic RTM, are also applied for studying the effects of oceanic bubbles on aerosol retrievals.The ground-based and satellite observations, as well as both RTMs, are discussed in detail in this section.

Observational data
The MAN data, derived from ship-borne measurements of direct solar attenuation by aerosol scattering and absorption over oceans, include retrieved AOD at five wavelengths ranging from 0.34 to 1.02 µm (Smirnov et al., 2011).The reported uncertainty in MAN AOD is ± 0.02 for all five channels (Smirnov et al., 2011).In this study, all MAN AOD data from 2004 to 2009 are used as the ground truth.
MODIS is on board both the Terra (passes over the equator at 10.30 a.m.local standard time) and Aqua (equator overpass at 1.30 p.m. local standard time) platforms.The MODIS instrument measures TOA radiation at 36 spectral channels, with spatial resolutions ranging between 250 and 1000 m at nadir and a wide swath of 2330 km.For this study, the Collection 5 (C5) Aqua MODIS dark target (DT) aerosol products are used (Remer et al., 2005).Note that the Collection 6 (C6) MODIS DT products were released recently, and this product includes a whitecap and ocean-foam effect in the radiative transfer lookup tables (LUTs) used to estimate aerosol properties (Levy et al., 2013).No accounting is made for subsurface ocean bubbles in any existing MODIS product.For this study, the C5 MODIS DT products are chosen to be consistent with the analysis done in Toth et al. (2013).Validated against ground-based observations, Remer et al. (2005) suggests that the uncertainty in over-ocean AODs is on the order of ±(0.03 + 0.05 • AOD).In addition, Shi et al. (2011a) shows that the root-mean-square errors (RMSEs) of C5 Aqua MODIS DT AOD data, estimated based on AERONET data, have a noise floor of 0.06 for AERONET AOD values less than 0.3.Also, the diagnostic estimates of RMSEs reach 0.3 for AERONET AOD values of 1.0.Lastly, for comparison purposes, some analyses based on 1 year (2009) of C6 Aqua MODIS DT data are included in the latter part of the study.Still, unless specifically mentioned, the MODIS data refer to the C5 MODIS data hereafter.
In addition to the Aqua MODIS DT aerosol products, the Collection 5 Level 1b Aqua MODIS radiance data are also collocated with wind speed data from AMSR-E and AOD data from MAN for evaluating the RTM simulations.MODIS channel 1 (0.66 micron) TOA radiance (at a 1 km resolution), along with latitude, longitude, and viewing geometry data, is extracted from the Level 1b MODIS data.The Collection 5 MODIS cloud mask data (MYD35) are also used for excluding cloud-contaminated pixels.
The near-surface wind speed values are obtained from the collocated version 7 AMSR-E data (Wentz and Meissner, 2000).On board the Aqua satellite, the AMSR-E is a conically scanning passive microwave radiometer with sensors in 12 microwave channels.The spatial resolutions of the AMSR-E data range from 5 to 56 km depending on the frequency (e.g., Wentz and Meissner, 2000).Retrieved surface parameters from AMSR-E include precipitation, sea surface temperatures, water vapor, wind speed, and other ancillary data.The AMSR-E data used in this study are formatted as hourly gridded binary files with a spatial resolution of 0.25 • × 0.25 • .As the first step, MAN data (2004MAN data ( -2009) ) are spatiotemporally collocated with the MODIS DT aerosol products.The temporal and spatial thresholds are set to ±30 min and 0.3 • latitude/longitude, respectively, following a study by Shi et al. (2011b).Also, only MAN AODs of less than 0.2 (at 0.55 µm) are used in order to minimize the impact of non-oceanic aerosols.To find the wind speed values for the collocated MAN and MODIS DT data pairs, AMSR-E data from the 0.25 • × 0.25 • grids, where the collocated MAN and MODIS DT data are located, are used.A temporal constraint is also applied to ensure the time difference between the AMSR-E and MAN observations is less than 1 h.As AMSR-E is on Aqua, timing between Aqua MODIS and AMSR-E-derived winds are functionally instantaneous.A total of 141 collocated MAN, MODIS DT, and AMSR-E data pairs are found.To find the TOA radiance at the MAN locations, a 0.1 • latitude/longitude box is centered on each of the 141 collocated MAN and AMSR-E data pairs.All cloudfree Level 1b MODIS data points within the 0.1 • × 0.1 • box are recorded, which is referred to as the MAN-DT-AMSR-E data set.To eliminate cloud contamination, the MODIS cloud mask (MYD35) is applied to this data set and pixels that are labeled as "confidently clear" are chosen.It is this data set that is used to compare MODIS-measured TOA radiance with 6S simulated TOA radiance.
Lastly, 7 years (2002)(2003)(2004)(2005)(2006)(2007)(2008) of Aqua MODIS DT and Level 2 quality-assured AERONET data (Smirnov et al., 2000) are used to study the differences in Aqua MODIS DT and AERONET AOD ( AOD DT−AERONET ) in relation to ocean surface wind speed.Such a study is not new, as the relationship between AOD DT−AERONET and ocean surface wind speed has been reported by several studies (e.g., Zhang and Reid, 2006;Shi et al., 2011a;Levy et al., 2013).This study extends the previous analysis to further evaluate the impacts of bubbles relative to the uncertainties in MODIS DT aerosol retrievals under clean marine conditions.For this purpose, over-ocean Aqua MODIS DT aerosol retrievals are collocated with AERONET data from coastal and island sites.The temporal and spatial differences are set to ± 30 min and ± 0.3 • (latitude/longitude), respectively, for the collocation process, following steps illustrated in Shi et al. (2011a).For each collocated Aqua MODIS DT and AERONET data pair, the surface wind speed value is then obtained from the collocated Navy Operational Global Atmospheric Prediction System (NOGAPS) data.The Level 2 AERONET data are considered the benchmark for validating satellite-based AOD retrievals, and the uncertainty of the Level 2 AERONET data is estimated to be ∼ 0.01-0.015(Holben et al., 1998).

Oceanic and atmospheric radiative transfer models
The 6S RTM is used to simulate TOA radiation measured by MODIS at the 0.66 µm spectral channel (Vermote et al., 2006).As mentioned in the Introduction, the 0.66 µm spectral channel is chosen as water-leaving radiance is almost negligible.It is for this reason that most aerosol retrievals rely on red and near-infrared wavelength radiances as the core data source.Furthermore, the multiple scattering is less in the red/infrared spectral channels compared with shorter wavelengths, making it much easier to separate the ocean bubble reflectance from the background (Zhang, 2001).For 6S simulations, a standard atmospheric profile (US standard 62) and the default maritime aerosol model included in the 6S model are used (Vermote et al., 2006).
In the 6S model, the ocean surface reflectance is computed based on Koepke (1984) by considering whitecap, glint, and water-leaving reflectance as shown in Eq. ( 1).
where Ref tot is the total surface reflectance, Ref wc is the whitecap reflectance, W is the whitecap coverage area (which is a function of near-surface wind speed), Ref G is the glint reflectance, and Ref wb is the water-leaving reflectance.
HydroLight is designed to simulate radiative transfer processes in oceans (Mobley et al., 2012).In order to estimate surface reflectance changes from ocean bubbles, the parameters of viewing geometry, wind speed, the ocean bubble phase function, and the ocean bubble concentration are needed for the HydroLight model runs.
The ocean bubble phase function (or the general shape of angular scattering) is adopted from Zhang et al. (2002), in which the bubble phase functions are computed from coated spheres based on Mie scattering theory and are intercompared with laboratory (using a volume scattering meter) and field observations.The coating represents surfactant material that adheres onto the bubble surface almost instantaneously after bubble genesis occurs in nature.Including the coating in optical computation is critical for remote sensing applications, because it can increase backscattered light by bubbles by a factor of up to 4 (Zhang et al., 1998).Little is known about the composition and thickness of the bubble coatings.However, a recent study shows that a coating of proteinaceous type (D'Arrigo, 1983;D'Arrigo et al., 1984) with a refractive index of 1.18 relative to water and a thickness of 0.01 µm provides the best match between optical and acoustical observations of bubbles (Czerski et al., 2011).In simulation, the bubble phase function is assumed to remain constant regardless of wind speeds.
Bubbles are frequently formed by breaking waves (Thorpe and Humphries, 1980;Lamarre and Melville, 1991).Because of the rapid rising of bubbles, the density distribution of ocean bubbles decreases exponentially with depth, while the overall concentrations increase with increasing wind speed following a power law.Ocean bubble concentrations in this experiment are obtained from Zhang (2001) and Zhang and Lewis (2002), in which the bubble concentrations at different layers are modeled as a function of wind speed based on field observations.

Estimating bubble surface reflectance
The bubble-induced perturbations in surface reflectance are estimated at the MODIS 0.66 µm spectral channel by taking the differences in the HydroLight runs with and without (background reflectance) consideration of oceanic bubbles.The background reflectance is obtained by running Hy-droLight with the bubble concentrations set to 0. Following Zhang (2001), the difference in surface reflectance ( R/π) with and without consideration of oceanic bubbles are plotted as functions of wavelength and near-surface wind speed as shown in Fig. 2.Here the ocean surface reflectance (or R) is defined as π × L w /E d , assuming a Lambertian surface, where L w and E d are in-air water-leaving radiance and surface downward flux, respectively.Note that Fig. 2 is a recreation of a figure from Zhang (2001).Based on Fig. 2, for a given wavelength, the changes in R/π as a function of near-surface wind speed are thus estimated.Following Zhang (2001), the linkage between wind speed and R/π is estimated using Eq. ( 2): where "wspd" is the ocean surface wind speed and L and J are coefficients.For the MODIS 0.66 µm channel, using the bubble concentrations (denoted default bubble concentration) as obtained from Zhang (2001), the L and J values are estimated to be 1.57× 10 −9 and 4.54, respectively.Thus, by assuming that submerged bubbles are uniformly distributed across a study domain, Eq. ( 1) can be modified to account for the effect of ocean bubbles on surface reflectance over whitecap-free regions as shown in Eq. ( 3): For further testing, an experiment is done using upper and lower boundaries of the default bubble concentrations (default bubble) that are estimated from Zhang (2001).The upper boundary is made by doubling the default bubble concentrations (double bubble), while the lower boundary is represented by half of the default concentrations (half bubble).The use of upper and lower boundaries allows analysis of two extreme conditions when compared to a normal set of conditions.Following the steps as mentioned above, the L and J values for the half, default, and double bubble concentrations are estimated and are shown in Table 1 for selected wavelengths.
To compare the relative contributions from whitecaps and submerged bubbles, Fig. 3 shows the magnitude of whitecap  reflectance with that of ocean bubbles for all three bubble concentrations, based on the 6S simulations.Clearly, the reflectance from submerged bubbles is insignificant for wind speed of less than 10 m s −1 .However, at high wind speeds (e.g., 15 m s −1 ), the reflectance from submerged bubbles are 20-30 % of the reflectance from whitecaps.Note that there is significant uncertainty in quantifying the relationship between wind and bubbles and whitecaps.For example, the concentration of bubbles is often related to the wind speeds through a power law.The exponent of the power law, mainly determined empirically, varies from 3.3 to 4.7 for different waters and conditions.For a wind speed of 10 m s −1 , the variation of bubble concentration would be over 1 order of magnitude.In our study, we limited this variability within a very constrained range of a factor of 2 to establish a conserved baseline.It is possible, however, that the effect could be greater or lesser.

Evaluating the linked 6S and HydroLight model
As the first step, the modeled MODIS channel 1 radiances from the linked 6S and HydroLight models are inter-compared with MODIS observations.The MAN-DT-AMSR-E data set is constructed to build an observationbased group of data that includes AOD from MAN, wind speed from AMSR-E, and satellite-measured radiance from MODIS.The MODIS cloud mask data are used to minimize cloud contamination, and only confidently clear pixels are chosen for further analysis.
Figure 4 shows a comparison between the 6S simulated radiance values at the MODIS 0.66 µm spectral channel and the MODIS observed radiance for the default bubble case using the MAN-DT-AMSR-E data set.The standard deviation is also plotted as the horizontal bar, which represents the variance in MODIS radiance within a 0.1 • latitude/longitude box.Data points which have the computed standard deviation values of larger than 4 w/m 2 /sr/µm are removed due to potential cloud contamination.The overall correlation between the modeled and observed radiances is 0.90, indicating that the linked 6S and HydroLight models are performing as expected.Also, similar results are found for both the half and double bubble cases (not shown).All three cases compare reasonably well with the observed radiance, although the effects of ocean bubbles are not readily visible.

Theoretical estimations
As the next step, an approach based on a LUT is adopted to estimate the impact of bubbles on AOD retrievals.Using the linked 6S and HydroLight model, LUTs of simulated MODIS TOA reflectance values are constructed as functions of solar zenith angle and viewing zenith angle (SZA and VZA; each varied from 0 to 60 • at intervals of 10 • ), relative azimuth angle (AZM; from 0 to 180 • at intervals of 30 • ), AOD (from 0.01 to 0.4 at intervals of 0.01), and wind speed (varied from 1.0 to 20 m s −1 ).The LUTs are constructed for two scenarios: with and without the inclusion of ocean bubbles.
For any given observing condition, the simulated reflectance values (both with and without bubble cases) are used to compute the errors in the retrieved AOD without considering bubbles.For example, for a given bubble concentration with a given wind speed and fixed viewing geometry, a relationship between AOD (AOD bub ) and simulated TOA reflectance is established.For the same observing conditions, the simulated reflectance is then used to search for the corresponding AOD value for bubble-free conditions (AOD no_bub ).The difference between AOD bub and AOD no_bub ( AOD) represents the simulated retrieval error without considering ocean bubbles in the retrieval process.Note that this process is done for default, double, and half bubble concentration cases.
To illustrate the concept, Fig. 5 shows the averaged AOD and near-surface wind speed relationships for three selected SZA values (0, 30, and 60 • ) and for the default, half, and double bubble concentration cases.There are numerous data entries in the LUTs.Thus, for illustration purposes, to construct Fig. 5 for a given set of SZA, wind speed, and bubble concentration, AOD values from all VZA and AZM entries are averaged into one value.The standard deviations of the variations, although not shown here, are on the order of 10-20 % of the mean values for wind speeds larger than 10 m s −1 for the SZA = 0 • case and are on the order of 30-40 and 50-70 % for SZAs of 30 and 60 • , respectively.Evident in Fig. 5 is that for wind speeds less than 10 m s −1 , AOD is almost negligible for all three bubble concentrations and all SZAs.However, AOD can be significant for wind speeds 12 m s −1 or higher.For comparison, the averaged AOD DT−AERONET (0.66 µm) and NO-GAPS surface wind speed relationship are overplotted on Fig. 5, represented by the triangle symbols.Each triangle symbol is computed by averaging AOD DT−AERONET values within a given 2 m s −1 NOGAPS wind speed bin using 7 years (2002)(2003)(2004)(2005)(2006)(2007)(2008) of collocated over-ocean Aqua MODIS DT and Level 2 AERONET data.Again, only Aqua MODIS DT and AERONET data pairs that have AERONET AOD values of less than 0.2 are used, to avoid influences from significant aerosol episodes.
Also, MODIS DT AOD retrievals with cloud fraction larger than 80 % (e.g., Shi et al., 2011a) and bad retrievals as identified by the quality assurance (QA) flags (QA = 0) included in the Aqua MODIS DT data are excluded.The thick black line is the linear fit through the averaged AOD DT−AERONET values.Clearly, the AOD DT−AERONET values increase as NOGAPS wind speed increases, indicating larger uncertainties in MODIS DT AOD retrievals at high ocean surface wind conditions.Note that such a finding is not new and is consistent with what has been reported in several previous studies (e.g., Zhang and Reid, 2006;Shi et al., 2011a;Levy et al., 2013).At a wind speed of 15 m s −1 , the averaged AOD DT−AERONET value is found to be around 0.05.In comparison, the AOD is approximately 0.01 (0.005, 0.02) for the default (half, double) bubble concentration case for a wind speed of 15 m s −1 and SZA of 0 • .Thus, approximately 20 % of the wind-related uncertainty for the MODIS DT products over remote oceans can be attributed to subsurface bubbles.Also, AOD values vary as a function of bubble concentration and the greatest impacts are observed for the double bubble case.Note that this conclusion is derived without considering uncertainty in the Koepke (1984) whitecap model (e.g., the area of whitecap coverage).We leave the latter issue for a future study.
Based on Fig. 5, from a theoretical analysis perspective, at low wind speeds there is no significant impact from ocean bubbles.However, as wind speeds increase, so does AOD.The impact becomes more significant once wind speeds are above 12 m s −1 , as AOD values increase exponentially with wind speed.While surface wind speeds of this magnitude do not occur in broad spatial or long-term averages (average global wind speed is around 6-7 m s −1 ), for oceanic regions with high near-surface wind speed, the impacts of ocean bubbles on the satellite-retrieved AOD values can be significant.For example, over the high-latitude southern oceans (30-70 • S), average MAN AOD is around 0.03-0.07(at 0.55 µm ;Smirnov et al., 2011;Toth et al., 2013).Even a change in satellite-retrieved AOD of 0.01 can be considered significant.
Using the MAN-DT-AMSR-E data set, the impact of bubbles on the difference between MODIS DT and MAN AOD ( AOD MODIS−MAN ) and wind speed is also studied (with and without the correction of the bubble effect) as shown in Fig. 6.The red cross signs show the original collocated data pairs.The blue-filled circles are the averaged raw AOD MODIS−MAN values for every 2 m s −1 wind speed bin.The blue line is the linear fit to the blue-filled circles.Again, the difference between MODIS DT and ground-based AODs and wind speed is not new as it have been reported using the Collection 4 MODIS DT and AERONET data (Zhang and Reid, 2006), as well as the Collection 5 MODIS DT and AERONET data (Shi et al., 2011a).Figure 6 confirms that such a relationship can also be observed using the collocated MAN and Collection 5 MODIS DT data.The green-filled circles are the averaged AOD MODIS−MAN values after applying a bubble correction (using the default bubble concentration) based on the constructed LUTs.The bubble correction is implemented to each MAN and MODIS DT pair by subtracting bubble-induced AOD to the AOD MODIS−MAN of the pair.The green line is a linear fit through the green-filled circles.Again, the averaged AOD corrections are small for the wind speed of 12 m s −1 or less.The AOD correction is noticeable (∼ 0.01), however, for the averaged wind speed bin of 12-14 m s −1 .

Applying ocean bubble correction to the ESOA
It has been shown that, theoretically, ocean bubbles can affect satellite aerosol retrievals in red wavelengths under conditions with near-surface wind speeds greater than 12 m s −1 .As the last step of this analysis, the impacts of oceanic bubbles on the ESOA feature detected from MODIS are evaluated.For this exercise, 1 year (2009) of Aqua MODIS AOD and AMSR-E wind speed data are used.For each day, AMSR-E wind speed and gridded Aqua MODIS AOD, solar zenith, viewing zenith, and relative azimuth angle data are constructed at a 1 × 1 • (latitude/longitude) resolution.
Figure 7a is a map of the average AMSR-E global wind speed for the year 2009 at a resolution of 1 • latitude/longitude.A relatively high wind speed band is found for the observed wind conditions at the default bubble concentration, ocean bubbles do not have a significant impact on satellite-retrieved AOD.Also, almost all areas sustain an annual mean AOD change of less than 10 %, with a majority having less than a 4 % change in AOD.
Figure 7d and e are similar to Fig. 7c but show results for the double (Fig. 7d) and half (Fig. 7e) bubble cases.The double bubble concentration has the largest impact on AOD, yet as evident in Fig. 7d, the majority of areas still experience less than a 10 % change in annual mean AOD.Therefore, these results indicate that ocean bubbles do not have a major impact on the ESOA.This is likely because the average wind speed throughout the ESOA region is not high enough to sustain a large contribution from ocean bubbles, as evident from Fig. 7a.
Recently, the C6 Aqua MODIS DT data were released.For comparison purposes, a portion of Fig. 7 is regenerated using the C6 MODIS DT data.Figure 8a shows the yearly mean C6 MODIS DT AOD in the ESOA region and Fig. 8b shows the AOD values with the use of default bubble concentration.As is evident in Fig. 8, there is no significant change in the ocean bubble correction when using C6 Aqua MODIS DT data.
The long-term means, however, may not represent individual cases.Figure 9 shows the cumulative distribution of instantaneous wind speeds in the ESOA region using AMSR-E data from the ascending orbits for 2009.Similar results are found using the data from the AMSR-E descending orbits and thus are not shown.Although the median wind speed of the ESOA region is under 10 m s −1 , the wind speeds can exceed 12 m s −1 in more than 10 % of cases.Thus, the subsurface ocean bubble effects may need to be considered for applications that use instantaneous MODIS DT retrievals, such as operational aerosol data assimilation (Zhang et al., 2008(Zhang et al., , 2011(Zhang et al., , 2014)).
To complement previous results, variations in the seasonal wind speed over the region as shown in Fig. 7 are studied as shown in Table 2 for the year 2009.The annual mean wind speed is 9.50 m s −1 with a standard deviation of 3.95 m s −1 .The northern hemispheric spring (March-April-May) and fall (September-October-November) have values similar to the annual mean.During the northern hemispheric winter (December-January-February), the mean wind speeds drop about 1.5 m s −1 .The highest mean wind speed of 10.27 m s −1 is found for the northern hemispheric summer (June-July-August).Each seasonal mean wind speed is not high enough for ocean bubbles to have a significant impact based on previous results.However, 1 standard deviation above each mean is either very close or above the 12 m s −1 threshold.This indicates that there are scenarios in which subsurface ocean bubbles could play a major factor in satellite AOD measurements.

Conclusions
In this study, the effects of ocean bubbles on satellite aerosol measurements are studied through a theoretical approach using a linked HydroLight oceanic and 6S atmospheric radiative transfer model.LUTs of the bubble-induced uncertainties in oceanic aerosol optical depth values retrieved from passive sensors ( AOD) are constructed as a function of satellite viewing geometry, near-surface wind speed, bubble concentration, and AOD.The AOD and wind speed relationships are studied for selected collocated MAN, MODIS, and AMSR-E data pairs.The contributions of AOD to the ESOA are also analyzed.This study suggests the following: 1.It is evident that at low wind speeds there is no significant impact from ocean bubbles on AOD retrievals using passive-based remote sensing techniques.However, the impact becomes much more significant at wind speeds above 12 m s −1 .At the wind speed range of around 15 m s −1 , the bubble-induced AOD may account for 20-30 % of the wind-related bias in overocean MODIS DT AOD retrievals.
2. The impacts of oceanic bubbles on the ESOA phenomenon are evaluated using 1 year of MODIS Collection 5 data, Advanced Microwave Scanning Radiometer Earth Observing System (AMSR-E) data, and lookup tables of AOD.It is found that ocean bubbles are not a major contributing factor to the ESOA, as the average wind speeds in that region are not high enough to produce a significant impact.The residual ESOA feature is therefore likely to be caused by other factors such as whitecaps and sub-pixel cloud contamination and thus deserves further study.
3. To avoid multiple scattering, the effect of subsurface bubbles on AOD retrievals is only evaluated at the 0.66 µm channel.However, as implied from Fig. 2, a much stronger effect may be expected at shorter wavelengths.Future studies may be needed to evaluate the impacts of subsurface bubbles on AOD retrievals from passive sensors that are centered on widely used wavelengths like the 0.55 µm spectrum channel.
4. Recently, the Collection 6 Aqua MODIS DT aerosol products have been released.New changes to the C6 MODIS DT aerosol products include a modified cirrus cloud detection scheme as well as the dependency of ocean surface reflectance as a function of wind speed.As a result, the ESOA feature is much reduced (Levy et al., 2013).Still, the submerged bubbles are not considered, and thus most of the discussions in this paper are valid for the C6 MODIS DT aerosol products.
5. There are several limitations in this study.Only theoretical calculations are included in the study for simulating the effects of bubbles on aerosol retrievals.The spatial and temporal variations of submerged bubbles and their optical and physical properties are not considered; this variation may significantly perturb the results from this study.It is likely that oceanic bubble contributions to the ocean surface reflectance are partially accounted for in empirical approaches based on direct estimation of overall wind speed impact (Shi et al., 2011a).Also, this is only a theoretical analysis.For practical applications, the uncertainties in whitecap estimates (fractional coverage and spectral reflectance, e.g., Frouin et al., 1996;Anguelova et al., 2006) need to be fully considered and incorporated into the analysis.The theoretical derivations used to estimate wind speed effects on surface reflectance for MODIS and MISR satellite aerosol retrievals explicitly include whitecaps and bubble rafts but not subsurface bubbles (Levy et al., 2013;Limbacher and Kahn, 2014).For future applications that require accurate estimations of atmospheric aerosol concentrations from satellite observations, oceanic bubble concentration is a factor that needs to be taken into consideration for ocean regions with strong near-surface winds.

Figure 2 .
Figure 2. The reflectance difference ( R/π) as functions of wavelength and near-surface ocean wind speed.The reflectance difference is defined as the difference in ocean surface reflectance (divided by π) between the HydroLight simulations with the use of the default bubble concentrations and a bubble concentration of 0.

Figure 3 .
Figure 3. Reflectance of whitecaps and ocean bubbles (for all three bubble concentrations) as a function of wind speed (estimated from 6S).

Figure 4 .
Figure 4.A comparison of 6S-HydroLight modeled radiance versus Aqua MODIS channel 1 radiance (0.66 µm).The default bubble concentrations are used in the 6S-HydroLight model simulations.The horizontal lines are the standard deviation of the radiance.

Figure 5 .
Figure 5. Averaged AOD as a function of wind speed for three solar zenith angles (SZA) of 0, 30, and 60 • as well as for using the default, half, and double bubble concentrations.The averaged differences (triangles) in Aqua MODIS DT and AERONET AOD (0.66 µm) are also plotted as a function of NOGAPS surface wind speed.The thick black line is the linear fit through the triangle symbols, and the 95 % confident intervals of the means are shown in vertical and horizontal lines.

Figure 6 .
Figure 6.A plot of the difference in MAN and MODIS DT AOD (AOD MODIS − AOD MAN , 0.66 µm) as a function of AMSR-E wind speed.The red crosses are the raw data.The blue-filled circles are the averaged AOD MODIS−MAN values for every 2 m s −2 wind speed bin.The green-filled circles are the same AOD MODIS−MAN values as the blue-filled circles after correcting for the bubble effect.

Figure 7 .
Figure 7. (a) Yearly averaged AMSR-E wind speed for 2009 for the latitude range of −30 to −70 • .(b) Yearly averaged Aqua DT MODIS AOD for the same region.(c) AOD correction for the default bubble concentration.(d) Same as Fig. 7c but for the double bubble concentration.(e) Same as Fig. 7c but for the half bubble concentration.

Figure 8 .
Figure 8. Yearly mean MODIS Collection 6 AOD paired with the AOD correction for the default bubble concentration.

Figure 9 .
Figure 9. Cumulative distribution function (CDF) of wind speed frequency over the ESOA region for 2009 using AMSR-E data from the ascending orbits.

Table 1 .
L and J coefficients for Eq.(2) for the half, default, and double bubble concentration scenarios at several wavelengths.

Table 2 .
Annual and seasonal mean wind speeds and standard deviations for the year 2009 in the ESOA region.