A Model for Air Entrainment Rates in Oceanic Whitecaps

Air‐entraining whitecaps provide an important source of bubbles over the global oceans, yet the rate at which the associated air is entrained is not well known. This lack of understanding limits the ability to accurately parameterize bubble‐mediated gas exchange and sea spray aerosol flux. In this paper I present a model to predict the total volume of air entrained by individual whitecaps and extend it to estimate the rate at which air is entrained per unit sea surface area. The model agrees well with existing models and measurements and can be forced by the rate at which energy is dissipated by the wavefield which can be routinely provided by spectral wave models. I then use the model to present the first distributions of the estimated total volume of air entrained by individual whitecaps, as well as their rate of air entrainment and air degassing.


Introduction
The air entrained by breaking gravity waves at the ocean surface plays a critical role in the regulation of the earth's climate by facilitating bubble-mediated gas exchange and sea spray aerosol production (De Leeuw et al., 2011;Deike, 2021;Keeling, 1993;Lewis & Schwartz, 2004;Melville, 1996;Woolf, 1993).The entrained air forms bubbles which results in the broadband scattering of light and, when present in sufficient quantities, can be seen as whitecaps at the ocean surface (Koepke, 1984;Monahan, 1986).These bubbles alter the optical properties of the upper ocean and increase the albedo of the global oceans (Frouin et al., 2001;Zhang et al., 1998).
Despite the importance of entrained air in the upper ocean to a wide range of processes, very little is known about the amount, and rate, of air entrained by either individual, or populations of, breaking waves.This lack of knowledge stems from the difficulties associated with (a) making time and space-resolved measurements of bubble size distributions inside actively breaking ocean waves and (b) characterizing the scale and frequency of occurrence of individual oceanic whitecaps in any given sea state.Consequently, many bubble-mediated processes are indirectly parameterized in terms of wind speed only, despite the fact that these processes are fundamentally wave-driven.
To address the knowledge gap Deike et al. (2017) combined results from direct numerical simulations of air entrainment in a three-dimensional simulation of a unidirectional breaking wave with the Phillips (1985) statistical distribution of breaking waves.The resulting model provides an estimate of the total volume flux of air by a population of oceanic whitecaps in a given sea state.Estimates of the third moment of the Phillips distribution, a measure of the wave spectrum and the wave slope are needed as model input parameters (Deike et al., 2017).The air entrainment model has been used to parameterize the bubble-mediated component of the gas transfer velocity of CO 2 (Deike & Melville, 2018).
In this paper I construct a model of the volume of air entrained in a single whitecap following the framework described in Callaghan et al. (2016Callaghan et al. ( , 2017)).The framework outlines how the surface foam features of individual whitecaps, which can be measured via optical remote sensing of the ocean surface, can be used to estimate the total energy dissipated by a single breaking wave along with the associated average bubble plume injection depth (Callaghan et al., 2024).By explicitly estimating the bubble plume injection depth, the two-dimensional surface whitecap signature can be linked to the evolving sub-surface three-dimensional turbulent two-phase flow which is important when considering air entrainment.Subsequently I construct a model for the air entrainment rate per unit sea surface area by a population of whitecaps following the energy balance arguments developed in Callaghan (2018).The result is a model for air entrainment rate that can be forced by the wavefield energy dissipation rate associated within a given sea-state.
The paper proceeds as follows.Section 2 describes the model and constrains key variables using available field and laboratory data.Section 3 applies the model to a field data set of 508 individual oceanic whitecaps taken at the Air Sea Interaction Tower (ASIT) at the Martha's Vineyard Coastal Observatory (MVCO) in 2008 as part of the SPACE08 campaign.The SPACE08 data set is extensively described in Callaghan et al. (2012); Callaghan et al. (2017Callaghan et al. ( , 2024) ) and Callaghan (2013Callaghan ( , 2018)).Finally, conclusions are presented in Section 4.

The Air Entrainment Model
The air entrainment model presented here is based upon the framework developed and discussed in Callaghan et al. (2016); Callaghan et al. (2017Callaghan et al. ( , 2024) ) which describes, and physically interprets, the time-dependent foam area evolution of individual whitecaps, A(t).Whitecap foam area evolution is characterized by a quasi-linear increase in time which is driven predominantly by air-entrainment during active breaking.The decay phase of foam area is characterized by a quasi-exponential decrease in time, driven by bubble plume degassing and surfactant-driven foam stabilization.The maximum whitecap area, A o , occurs at time t = t A o and separates the whitecap growth phase from its decay phase.Integral timescales of the growth and decay phases are described by dt, respectively.To explicitly separate the effects of bubble plume degassing and surfactant stabilization, τ decay can be written as the sum of τ degas and τ stab .If surfactant effects are negligible, then τ stab ≈ 0. Callaghan et al. (2016) describes how these whitecap area and timescale measurements can be used to estimate the total energy dissipated by a single whitecap and its associated bubble plume injection depth.Callaghan (2018) builds on this and presents an energy dissipation based model for total whitecap coverage, W, and growth phase whitecap coverage, W growth .

Model Development
The total volume of air entrained by a single whitecap can be written as where α(t) and V(t) are the time-evolving air fraction and volume of the two-phase flow beneath the surface whitecap, respectively.The leading dimensionless scaling factor ξ accounts for the fact that some air may be lost during the whitecap growth period and some air may be entrained after the time of maximum foam area, t A o .For deep to intermediate-water breaking waves with a finite duration considered in this study I assume that ξ = 1, but for depth-limited shoaling breaking waves it would be expected that ξ > 1.The value of the dimensionless coefficient c o is dependent on the time history of α(t)V(t).Here I use a value of 2 which assumes that α(t)V(t) increases linearly in time to t = t A o .
Constraining α(t) is key to estimate the total air entrained by a whitecap, however this is beyond the scope of the present work.Instead, I use an effective volume-weighted air fraction, α eff , defined as By writing V(t) = A(t)z p (t), where z p (t) is the time-evolving bubble plume injection depth defined (formerly referred to as a penetration depth in Callaghan et al. (2016)), Equation 1can now be written as The value of z p (t) is assumed to be uniform in x and y beneath the surface whitecap foam patch.By invoking the definitions for a foam area-weighted bubble plume injection depth, ẑp , and an integral whitecap area growth timescale, τ growth , both given in Callaghan et al. (2016), the total volume of air entrained by a single whitecap is given by On the right hand side of Equation 4I have assumed ξ = 1 and defined an average vertical downward entrainment velocity as w ent = ẑp / t A o .It should be noted that because ẑp is an average quantity, its value will be less than the maximum bubble plume injection depth during active breaking.Consequently instantaneous values of the entraiment velocity will be greater than w ent , probably by up to a factor of two or more (Callaghan et al., 2016).
Normalizing Equation 4by A o τ growth gives the average rate of air entrainment per unit whitecap area for a single whitecap, Vwc , which is given as Having established a model for the rate of air entrainment by a single whitecap, I now extend this to quantify the average rate of air entrained per unit sea surface area by a population of whitecaps within a given observational area, A obs , and time period, T obs .This is denoted as Vss and is given by In Equation 6I use c o = 2, as discussed above.
Evaluation of Equation 6 requires a combination of simultaneous above and below-water measurements for all individual whitecaps in a population.This has yet to be achieved in field measurements and to proceed I simply assume that α eff is constant for all whitecaps.The quantity w ent can be inferred from surface whitecap features since the t A o can be measured directly and ẑp can be inferred from τ degas .As quantified in Callaghan et al. (2013); Callaghan et al. (2016Callaghan et al. ( , 2017)), deeper bubble plumes take longer to degas such that there is a physical link between τ degas and ẑp .
I choose to further simplify Equation 6to allow Vss to be determined from a measure of the whitecap coverage of the sea surface.To do so, I define the growth-phase weighted average value of the entrainment velocity to be By adopting the definition for growth-phase whitecap coverage, W growth , in Callaghan (2018), which is given as Equation 6 is now written as Under the limitations of the assumptions adopted above, it can be seen that the rate of air entrained by oceanic whitecaps is linearly related to W growth with an additional dependence on the average entrainment velocity for a given population of whitecaps along with their effective air fraction.Following Callaghan (2018), Vss can also be expressed in terms of total whitecap coverage as where δ depends on how long the decaying foam patches reside on the water surface.This time period is in turn dependent on bubble plume injection depth and surfactant-driven stabilization of whitecap foam.In the case of the former, deeper bubble plumes sustain the surface whitecap for longer, and, for the latter, surfactants typically cause foam to linger on the surface after the bubble plume has degassed.Callaghan (2018) provide further details on how to parameterize δ, and the component of δ that is driven by bubble plume degassing only (See their Section 4.2).Moreover, Callaghan (2018) shows how both W growth and W can be estimated from the rate of energy dissipation of the surface gravity wavefield.

Constraining Model Parameters α eff and w ent
Values of w ent for individual whitecaps have yet to be directly measured in the field, but can be inferred using above-water measurements of foam area evolution because t A o can be measured and ẑp can be inferred from measurements of τ degas (Callaghan et al., 2016(Callaghan et al., , 2017)).Resulting distributions of w ent derived from the MVCO data set are shown in Figure 1 for four observational periods along with the average for the entire data set.As stated above, instantaneous values of w ent should be expected to be larger by a factor of two or more than this mean value and ẑp will typically be smaller than the instantaneous value of the bubble plume depth at the end of active breaking (see Supporting Information in Callaghan et al. (2016)).The distribution of all values of ẑp in Figure 1 follows a lognormal behavior and has an ensemble average value of 6.8 cm s 1 compared to the value of w ent = 6.5 cm s 1 as defined by Equation 7. The associated lognormal fits to ẑp estimates from each observational period are given in Table S1 in Supporting Information S1.Also included in Table S2 in Supporting Information S1 is the probability that the values of w ent from different pairs of observational periods are drawn from the same continuous distribution as evaluated by a two-sample Kolmogorov-Smirnov test.
The value of α eff is now constrained using Equation 5 together with estimates of Vwc from Cipriano and Blanchard (1981) and w ent = 6.5 cm s 1 estimated from the MVCO data presented in Figure 1.As described in Woolf (1993), Cipriano and Blanchard (1981) created a steady-state breaking wave analogue in the laboratory The gray, red, green and blue lines represent the distributions for Periods I-IV respectively.The associated average wind speeds for these periods were 11.4,8.2, 5.7 and 13.7 m s 1 .Period IV directly followed Period III and is representative of unsteady, rapidly changing conditions.
that entrained air at a volumetric rate of 1.25 × 10 4 m 3 s 1 .This was determined for bubble sizes with radii in the range 50 μm to 4 mm which overlaps well with bubble sizes measured in laboratory and oceanic breaking waves (Deane & Stokes, 2002).That air entrainment rate took place through a nominal whitecap area of 0.02 m 2 giving a value Vwc ≈ 6.25 × 10 3 m 3 m 2 s 1 .
Using values of Vwc = 6.25 × 10 3 m 3 m 2 s 1 and w ent = 6.5 cm s 1 in Equation 5gives α eff ≈ 0.1 when c o = 1 is chosen.This latter choice is reasonable given the steady state nature of the breaking wave analogue in Cipriano and Blanchard (1981).Previous measurements of the air fraction inside laboratory breaking waves show enormous spatial and temporal heterogeneity (e.g., Blenkinsopp and Chaplin (2007); Lim et al. (2015)).Indeed, representative average values of α have been reported using a variety of methods including (a) conductivity measurements (e.g., Lamarre and Melville (1991) α ≈ 0.2), (b) fiber optic probes (e.g., Rojas and Loewen (2010) α ≈ 0.012 0.37 for laboratory spilling breaking waves and α ≈ 0.17 0.29 for laboratory plunging breaking waves) (c) impedance probes (e.g., Cox and Shin (2003) α ≈ 0.15-0.2) and (d) digital imagery (e.g., Deane and Stokes (2002) α ≈ 0.065 for oceanic breaking waves and α ≈ 0.045 for laboratory breaking waves).The value of α eff ≈ 0.1 found here is therefore very reasonable given the indirect nature of the calculation using a combination of a simplified model and a collection of diverse data sets.

Application to Field Data
Having developed models for the total volume of air entrained by a single whitecap (Equation 4) and the average rate of air entrainment per unit whitecap area (Equation 5) I now apply these models to the MVCO oceanic whitecap data set.I follow this with estimates of the average rate of air entrainment per unit sea surface area forced using measurements of W reported in Callaghan et al. (2008).

Individual Oceanic Whitecaps
Figure 2 shows the estimated total volume of air entrained (V air ) by, and the rate of air entrainment per unit whitecap area ( Vwc ) of, individual oceanic whitecaps in panel (a) and (b), respectively.Also shown in panel (c) is the estimated air degassing rate per unit whitecap area, Vdegas , which is calculated by normalizing Equation 4 by the quantity A o τ degas .To the best of the knowledge of the author, the distributions in Figures 2a-2c) for winddriven oceanic whitecaps are the first to be presented in the literature.
The estimated values of V tot span over three orders of magnitude which reflects the variation in A o and τ growth values.The former ranges from 0.2 to 26 m 2 and the latter from 0.1 to 2.2 s.The entire data set in Figure 2a is wellapproximated by a lognormal distribution, as indicated by open black circles.All lognormal fit coefficients for the entire data set, and each observational period, are provided in Table S3 in Supporting Information S1.There is a clear shift in the V tot distributions to higher values with an increase in wind speed from 5.7 m s 1 in Period III, to 11.4 m s 1 in Period I. Period IV with the highest average wind speed is the exception because, as discussed in Callaghan et al. (2012Callaghan et al. ( , 2017)), it is representative of unsteady conditions.These conclusions are supported by the statistical analysis in Table S4 in Supporting Information S1.The standard deviation of the fitted lognormal distributions to all four periods is remarkably similar (See Table S3 in Supporting Information S1), indicating that differences in the distributions are driven by different mean values.
In contrast to V tot , the spread of the Vwc distributions in Figure 2b is much narrower, only spanning about one order of magnitude, with a very weak wind speed dependence.The Vwc estimated values are well-approximated by a lognormal distribution, with associated fit coefficients provided in Table S5 in Supporting Information S1.A two-sample Kolmogorov-Smirnov test presented in Table S6 in Supporting Information S1 indicates a large increase in statistical similarity in Vwc between the observational periods compared to that for V tot , with Period III being an exception.
Distributions of the air degassing rate per unit whitecap area ( Vdegas ) are shown in Figure 2c and exhibit the least spread.Corresponding lognormal distribution coefficients are presented in Table S7 in Supporting Information S1. Figure 2c implies that the degassing process may be largely independent of the scales of waves observed here, and is dominated by the dynamics of rising bubbles in a turbulent flow as argued in Blenkinsopp and Chaplin (2011).Indeed, the cross-sectional area of bubbles in size distributions typical of oceanic breaking waves is dominated by bubbles with radii larger than about 0.5 mm, and the rise velocity of bubbles with radii between approximately 0.5 to 10 mm shows relatively little variation (Deane et al., 2013;Deane & Stokes, 2002).It is therefore physically plausible that the average rate at which air degasses from a decaying whitecap shows relatively little variation.

Estimates of the Air Entrainment Rate per Unit Sea Surface Area
Having looked at distributions of various air fluxes for individual oceanic whitecaps, I now estimate the rate at which air is entrained per unit sea surface area as a function of both growth-phase whitecap coverage and wind speed.The former uses Equation 9 and the latter uses Equation 10 in combination with measured W data from Callaghan et al. (2008) and the model presented in Callaghan (2018).
The solid black line in Figure 3a shows the rate of air entrainment per unit sea surface area by actively breaking waves as predicted by Equation 9.The shaded area around this line represents the range of estimated Vss values when α eff is set to vary between 0.045 and 0.37, chosen to encompass the majority of air fraction values given in Section 2.2.There are no available field measurements to compare the model output to, but a comparison can be made to model estimates presented in Deike et al. (2017).It is seen that using Equation 9 with α eff = 0.1 predicts Vss values that lie in the range of values presented in Deike et al. (2017), and that their values lie within the shaded region.Even though there remains a great deal of uncertainty about the most appropriate values of α eff and w ent , and how these may vary with sea state, the level of agreement between the models is encouraging given their different formulations.
Estimated air entrainment rates for the W data set of Callaghan et al. (2008) are presented in Figure 3b.In doing so, I have assumed that the W values do not have any contribution from surfactant-driven foam stabilization and reflect the entrainment and degassing processes only.Consequently, I have used equation 21 in Callaghan (2018) to evaluate δ.This is an oversimplification and not likely to be true, but any surfactant influence cannot be determined in retrospect.Also shown for comparison are model estimates digitised from Figure 3d

Conclusions
I have presented a simple model to quantify air entrainment by whitecaps with key variables constrained with a combination of laboratory data and estimates from field data.I have used the model to provide the first reported distributions of the volume of air entrained, along with its entrainment and degassing rates, for individual oceanic whitecaps.The volume of air entrained varies by several orders of magnitude for individual whitecaps within a sea state and shows a dependence on wind speed.In contrast the rate of air entrainment per unit whitecap area shows a much narrower distribution and no clear wind speed dependence.
When forced with measurements of whitecap coverage from the North Altantic the model produces estimates of the rate of air entrainment per unit sea surface area that compare favorably with the model of Deike et al. (2017) and the limited laboratory measurements of Toba (1961).Given that constant values of air fraction and entrainment velocity have been used here, variability in the estimates presented here is ultimately driven by variability in fractional whitecap coverage.A key point of the previous work by Deike et al. (2017), and widely held in the broader community, is that wind speed alone is insufficient to fully constrain wave breaking and hence air entrainment, and wave information is also needed.To this end, the energy dissipation model for W presented in Callaghan (2018) can be easily used in spectral wave models to provide estimates of air entrainment rates using a combination of wind and wave variables.
To conclude, the physically based model presented here can provide routine estimates of the rate of air entrainment by whitecaps per unit sea surface area.Improvements can be made by better constraining the air fraction and entrainment velocity of bubble plumes in actively breaking waves (i.e., short timescales) through dedicated field campaigns, but the effort to achieve this is not trivial.The success of the model presented here will be evaluated through its ability to provide predictions of bubble-mediated processes such as air-sea gas transfer and sea spray aerosol production that agree with observations.
for the entrainment of air by oceanic whitecaps is presented which agrees well with existing models and measurements • Distributions of the estimated volume of air entrained by individual oceanic whitecaps are presented for the first time • Key uncertainties in the air fraction and entrainment velocities of individual whitecaps remain due to a lack of measurements Supporting Information: Supporting Information may be found in the online version of this article.

Figure 1 .
Figure 1.A probability density distribution of estimated mean entrainment velocity, w ent , for a population of 508 oceanic whitecaps measured from the Air-Sea Interaction Tower (ASIT) at the Martha's Vineyard Coastal Observatory (MVCO) during the SPACE08 campaign.The black line represents the entire data set and the circles are the associated lognormal fit.The gray, red, green and blue lines represent the distributions for Periods I-IV respectively.The associated average wind speeds for these periods were 11.4,8.2, 5.7 and 13.7 m s 1 .Period IV directly followed Period III and is representative of unsteady, rapidly changing conditions.

Figure 2 .
Figure 2. Probability density distributions of estimated values of (a) the total air entrained by individual whitecaps, V tot , (b) the rate of air entrainment per unit whitecap area, Vent , and (c) the rate of air degassing per unit whitecap area, Vdegas .The black lines represent the total data set and the black circles the associated lognormal distribution fits.The gray, red, green and blue lines represent the distributions for Periods I-IV respectively.The associated average wind speeds for these periods were 11.4,8.2, 5.7 and 13.7 m s 1 .Period IV directly followed Period III and is representative of unsteady, rapidly changing conditions.
in Deike et al. (2017) (black dots) and measured laboratory data digitised from Figure12inToba (1961) for wind-driven laboratory breaking waves.While it may be difficult to extrapolate the laboratory data to an open ocean setting, the data have the advantage of being measured (subject to some assumptions) and therefore provide a useful comparison with the present estimates.Indeed, the data fromToba (1961) display a similar trend with increasing wind speed that is seen in both the current work andDeike et al. (2017).The spread of the present Vss estimates at wind speeds below about 10 m s 1 compares well the estimates reported inDeike et al. (2017) albeit with an offset to lower values.At larger wind speeds, the current estimates lie within the scatter of theDeike et al. (2017) estimates.There is more scatter in theDeike et al. (2017) estimates, but these are derived from whitecap measurements from 4 different field campaigns using two measurement approaches, whereas the current estimates are from a single campaign.It should be pointed out thatDeike et al. (2017) stress that wind speed alone cannot fully constrain estimates in air entrainment rates, something I agree with.Indeed, the model for W in Callaghan (2018) is a function of wavefield energy dissipation rate, and estimating W and hence air entrainment rates from a spectral wave model is part of ongoing work, and is expected to show a range of values at any given wind speed.

Figure 3 .
Figure 3. (a) The estimated variation in air entrainment rate per unit sea surface area, Vss , as function of growth phase whitecap coverage, W growth .The solid black line is Equation 9 with α eff = 0.1.The shaded region indicates the associated spear in Vss estimates when α eff is varied from 0.045 to 0.37.The black dotted and dashed lines represent the upper and lower estimates from Deike et al. (2017).(b) The gray circles are estimated values of Vss using the W data from Callaghan et al. (2008) with Equation 10.Also shown are the estimates values from Deike et al. (2017) digitised from their Figure 3d as solid black circles.The black triangles are measured laboratory data from Toba (1961) which have been digitised from their Figure 12.