Sea surface warming suppresses primary sea spray aerosol number production

Sea spray aerosol (SSA) influences climate through direct and indirect interactions with radiation. However, the magnitudes of these interactions remain poorly constrained, in part due to a lack of understanding of the influences of sea surface temperature (SST) on SSA production. There is no agreed-upon dependence of SSA production on SST despite numerous field, laboratory, and modelling investigations. In this study, we describe a simple theoretical framework relating the interfacial processes and contextualizing previous work. Next, we characterize the connection between SST, bubble concentrations, SSA number concentrations, and SSA emission fluxes using measurements in the Scripps Ocean-Atmosphere Research Simulator (SOARS). This isolated ocean-atmosphere interaction system incorporates wind, waves, and SST controls to produce wave breaking under realistic and controlled conditions. Increasing SST from 2 to 23 °C suppressed total subsurface bubble concentrations (between 6.17 and 830 µm) by a factor of 1.5, SSA number concentrations (between 0.008 and 20 µm) by a factor of 3, and SSA number flux by a factor of 3. Importantly, we identify size-resolved differences including a non-monotonic SST dependence of supermicron number and mass emissions. Using these trends, we derive SST-dependent number and mass emission flux correction factors for SSA source functions in climate models. These controlled wind-wave-SST experiments demonstrate that increasing SST overall suppresses SSA production. Resolving this SST dependence is critical, as it directly alters marine aerosol burdens, cloud condensation nuclei, and radiative forcing, and provides a needed constraint missing from current parameterizations.

1 Introduction

Since the 1980s, global surface temperatures have increased at an estimated 0.2 °C per decade, accelerating to 0.35 °C per decade since about 2010 (Foster and Rahmstorf, 2026; Hansen et al., 2006). More than 90 % of this warming is estimated to be absorbed by the oceans, with warming concentrated at the surface (Llovel and Terray, 2016; Willis et al., 2004), altering the ocean-atmosphere interface (OAI) (Resplandy et al., 2019). Aerosols produced at the OAI play an important role in mitigating global warming directly by scattering solar radiation and indirectly by influencing the micro- and macro-physical properties of clouds (Albrecht, 1989; Haywood and Boucher, 2000; IPCC, 2023; Lohmann and Feichter, 2005; Twomey, 1980). With 71 % of Earth's surface covered by oceans, contributions from both natural and anthropogenic marine aerosols are active areas of study (Lewis and Schwartz, 2004; O'Dowd and De Leeuw, 2007; Quinn et al., 2015; Russell et al., 2023).

Sea spray aerosols (SSA) are the most ubiquitous natural marine aerosol, responsible for an estimated mass flux between 3 and 30 Pg yr⁻¹ (Gliß et al., 2021; Grythe et al., 2014; de Leeuw et al., 2011; Lewis and Schwartz, 2004). SSA are produced during wave breaking, where air is entrained beneath the ocean surface to form bubbles. These subsurface bubbles rise and produce surface whitecaps, which emit SSA primarily through film and jet drop production (Bird et al., 2010; Brasz et al., 2018; de Leeuw et al., 2011; Lewis and Schwartz, 2004; Quinn et al., 2015). Additionally, SSA are produced when surface wind speeds exceed between 7 and 11 m s⁻¹, wherein spume drops are ripped directly from the top of wave crests (Andreas et al., 1995; Mehta et al., 2019; Veron et al., 2012), however, this process is not the focus of the current study.

Quantifying the influence of sea surface temperature (SST) on SSA emissions has been the focus of numerous field measurements, laboratory experiments, and first-principles theories, yet no consensus has emerged. The review by Song et al. (2023) and the introductions of Markuszewski et al. (2024) and Salter et al. (2014) provide recent summaries of this literature, and with no significant advances since. A thorough comparison of previous studies (e.g., methodology, size range, and measured quantity) is beyond the scope of this paper; we therefore summarize their main findings in the Supplement (Sect. S1, Table S1). In brief, some studies report a negative correlation between SST and SSA production (Lehahn et al., 2014; Markuszewski et al., 2024; Nielsen and Bilde, 2020; Sellegri et al., 2023; Sofiev et al., 2011), others a positive correlation (Forestieri et al., 2018; Grythe et al., 2014; Hu et al., 2024; Jaeglé et al., 2011; Liu et al., 2021; Ovadnevaite et al., 2014), and still others non-monotonic relationships (Christiansen et al., 2019; Hu et al., 2024; Mårtensson et al., 2003; Saliba et al., 2019; Salter et al., 2014, 2015; Zábori et al., 2012; Zinke et al., 2022). This disagreement reflects more than a divergence in sign of a single response. Changes in SST simultaneously modifies several properties governing bubble bursting, including seawater viscosity, surface tension, and dissolved gas solubility. These properties in turn influence film and jet drop production in size-dependent ways.

The use of different breaking wave proxies in laboratory studies complicates the interpretation of results in previous studies. The importance of accurate SSA generation mechanisms cannot be overstated (Collins et al., 2014; Fuentes et al., 2010; Lewis and Schwartz, 2004; Prather et al., 2013). A fundamental limitation of previous laboratory studies, and a plausible explanation for conflicting results, is the absence of controlled winds. This study uses the Scripps Ocean-Atmosphere Research Simulator (SOARS) to control SST, wave state, and wind speed to characterize SSA emissions at varied SSTs. Previous measurements of particle number size distributions in SOARS as functions of wind speed and wave breaking intensity show good agreement with field and laboratory estimates (Leibensperger III et al., 2026b). While other previous laboratory proxies offer valuable insights into the dependence of SSA production on SST, SSA production from representative wind-wave interactions will alter bubble plumes, and thus SSA emissions. Measurements in SOARS can be used to correct simulated SSA emissions in Earth System Models (ESMs), which typically incorporate a positive correlation with SST, if there is an included SST dependence at all (Grythe et al., 2014; Jaeglé et al., 2011; Song et al., 2023).

This study presents an integrated investigation of the SST dependencies of marine aerosol production through the conceptualization of a common framework (Sect. 3.1), followed by SST-dependent measurements of whitecap fractions and subsurface bubble populations (Sect. 3.2), SSA number concentrations (Sect. 3.3), a comparison of SSA and subsurface bubble concentrations (Sect. 3.4), and SSA number emission fluxes (Sect. 3.5). These findings culminate with temperature-dependent correction factors derived from the first wind-wave controlled SST experiments for direct integration into climate models (Sect. 3.6) and a discussion of the role of temperature hysteresis on SSA production (Sect. 3.7).

Figure 1 Experimental overview. (a) Sea surface temperature for each experimental day and (b) total aerosol number concentrations during a representative experimental day (11 February 2024, cooling experiment day 1), showing the size measurement periods of a daily experiment: Room Air, Aerosol Mode, Wind Alone, Wind & Waves, Decay. Magenta dashed lines indicate the Production Flux (concentration rise following wave onset) and Steady State Concentration (subsequent plateau), from which SSA number fluxes are derived.

2 Materials and Methods

2.1 Sea Spray Aerosol Generation in SOARS

The SOARS instrument is located within the Hydraulics Laboratory at the Scripps Institution of Oceanography at the University of California San Diego in La Jolla, CA, USA. SOARS is a combined wave channel (36 m × 2.4 m × 2.4 m) and recirculating wind tunnel with an air-backed paddle at one end generating user-defined wave fields and a partially submerged ramp at the other end to absorb and dissipate wave energy. Two experiments were conducted from 11 February to 1 March 2024 using SOARS to quantify the effects of increasing SSTs on SSA production. In each experiment, water temperatures in SOARS, hereafter referred to as SSTs for simplicity, were controlled across a ∼ 20 °C range (Fig. 1a). The first experiment (“Cooling”), monotonically decreased SST each day, while the second experiment (“Warming”) did the opposite. Both experiments utilized the same seawater, collected from the Ellen Browning Scripps Memorial Pier (La Jolla, CA; 32°52^′1^′′ N; 117°15^′26^′′ W). The seawater was filtered through 10 µm mesh filters prior to experimental runs to reduce the influence of biology, pollution, and surface-active chemistry on SSA emissions (Kimble et al., 2026). Once filled, the SOARS water was continuously filtered through 1 µm mesh filters and UV sterilizers (λ=264 nm). SSA was produced by combining a wave field, where every fifth wave crest was a gently spilling breaker along the length of the channel, with a 10 m extrapolated wind speed (U₁₀) of 11 m s⁻¹. This wind speed and wave field pairing was chosen to replicate equilibrium seas not under high winds, with more discussion included in Leibensperger III et al. (2026b).

A representative experimental day of operating conditions is shown in Fig. 1b, broken down into six measurement periods: (1) Room Air, (2) Aerosol Mode, (3) Wind Alone, (4) Wind & Waves, (5) Decay, and (6) Room Air. Representative particle number size distributions (PNSDs) from each period on three days with equivalent SSTs are displayed in Fig. S1 in the Supplement. Room air measurements from the beginning and end of each day characterize background aerosol concentrations in the Hydraulics Laboratory. Aerosol mode, wherein the SOARS headspace is cycled through high efficiency air particulate (HEPA), clean carbon MERV 16, and potassium permanganate filter banks, is operated for ∼ 45 min, removing over 85 % of background particles. Once aerosol concentrations are sufficiently low, U₁₀ winds are set to 11 m s⁻¹ for ∼ 45 min (Wind Alone) to remove particles adhered to channel walls and create a stable headspace before SSA generation. During Wind & Waves, a wave field is generated with a peak wave height of 0.42 m, enabling measurements of SSA concentrations and flux estimates (Sect. 2.3–2.4). The features of the aerosol concentration time series explored in the following sections are annotated by the magenta dashed lines in Fig. 1b. SSTs and salinities (Table S2) are measured during Wind & Waves using a YSI Pro 30 (YSI Inc., Yellow Springs, OH). Following Wind & Waves, waves are stopped to allow for estimates of aerosol loss rates during the Decay period.

2.2 Sea Spray Aerosol Sampling

Primary SSA was sampled through two stainless steel inlets (2.27 m length; 10 mm ID) positioned 25 m downstream of wave generation and angled in the direction of wave breaking at a height of 0.6 m above the water surface (Fig. S2). The two parallel sampling systems had total flow rates of 6.2 and 2.6 L min⁻¹. The sampled SSA was dried (relative humidity < 20 %) by passing sample flows through custom-built silica diffusion dryers to minimize the influence of hygroscopic growth on particle size. The sampling system produces size-dependent particle losses that must be quantified and corrected. The Particle Loss Calculator (PLC), described by von der Weiden et al. (2009), was used to estimate the sampling efficiency for each instrument during Wind & Waves conditions (Fig. S3). Loss mechanisms incorporated include diffusion, sedimentation, and inertial deposition in bends and contractions. Sampling efficiency ranged from ∼ 55 % to 100 % over the size range of 8 to 200 nm. Slight enhancements around 200 nm to 4 µm were observed, consistent with previous SOARS studies (Leibensperger III et al., 2026b; Moore et al., 2025). The wind speeds in SOARS produce sub-isokinetic conditions, leading to enhancement in this size range. Beyond 4 µm, sampling efficiency decreases significantly to 20 % at 10 µm. These corrections were applied to all measurements prior to subsequent analysis.

2.3 Sea Spray Aerosol Measurements

Submicron SSA (D<1 µm) were measured using a Scanning Mobility Particle Sizer model 3938 (SMPS; TSI Inc., Shoreview, MN) and Spider-MAGIC model 800 (SM; Aerosol Dynamics Inc., Berkeley, CA). The SMPS consisted of a 0.071 cm impactor, model 3088 advanced aerosol neutralizer, model 3081A differential mobility analyser (DMA) operated in negative polarity, and a model 3789 condensation particle counter (CPC). Aerosol and sheath flow rates of 0.6 and 3 L min⁻¹ enabled measurements of particles with mobility diameters (D_m) from 0.013 to 0.835 µm every minute, with a scan time of 45 s and a purge time of 15 s. The SM, introduced by Amanatidis et al. (2020, 2021) consisted of a Po-210 radioactive charge conditioner (500 µCi; NRD, Grand Island, NY), Spider DMA operated with alternating up and down dual polarity scans, and a Moderated Aerosol Growth with Internal water Cycling (MAGIC) CPC (Hering et al., 2019). Aerosol and sheath flow rates of 0.3 and 0.9 L min⁻¹ enabled measurements of particles with D_m from 0.008 to 0.42 µm. The SM scan time was approximately 30 s for a unidirectional, single polarity voltage ramp, with a complete scan over both polarities and each voltage ramp direction in about two minutes.

Supermicron SSA (D>1 µm) were measured using aerodynamic and optical measurement techniques. An Aerodynamic Particle Sizer model 3321 (APS; TSI Inc., Shoreview, MN) measured particles with aerodynamic diameters (D_a) from 0.52 to 20 µm with a total flow rate of 5 L min⁻¹. An Optical Particle Sizer model 3330 (OPS; TSI Inc., Shoreview, MN) measured particles with optical diameters (D_o) from 0.3 to 10 µm with a sample flow rate of 1 L min⁻¹. The APS and OPS were both operated at one-minute resolution.

Particle number size distributions (PNSDs) from the three measurement techniques were converted to dry, physical diameters (D_p,dry) assuming spherical particles. Equations (1) and (2) have been used previously to combine mobility and aerodynamic measurements (Collins et al., 2014; Lewis and Schwartz, 2004; Stokes et al., 2013). Following Leibensperger III et al. (2026b), D_m is assumed to be equal to D_p,dry (Eq. 1) and D_a is converted to D_p,dry by estimating the effective density (ρ_eff) based on relative humidity, assuming a mixture of pure NaCl salt and water following Fig. 7 from Lewis and Schwartz (2004), and a reference density (ρ₀) of 1 g cm⁻³ (Eq. 2). Finally, D_o is assumed to be equal to D_p,dry (Eq. 3), as previous measurements of the refractive index of dried SSA (m_SSA=1.533) by Liu et al. (2023) were found to be reasonably close to that of the calibration aerosol, polystyrene latex spheres (m_PSL=1.588).

$$\begin{array}{}\text{(1)}& {\displaystyle }{D}_{\mathrm{p},\mathrm{dry}}=D_{\mathrm{m}},\text{(2)}& {\displaystyle }{D}_{\mathrm{p},\mathrm{dry}}={\displaystyle \frac{D_{\mathrm{a}}}{\sqrt{\frac{{\mathit{\rho }}_{\mathrm{eff}}}{{\mathit{\rho }}_{\mathrm{0}}}}}},\text{(3)}& {\displaystyle }{D}_{\mathrm{p},\mathrm{dry}}=D_{\mathrm{o}},\end{array}$$

Total particle concentrations from 5 nm to 2.5 µm were measured using a MAGIC CPC (Aerosol Dynamic Inc., Berkeley, CA). The MAGIC CPC served as the reference instrument for room air and in-channel number concentrations shown in Fig. 1b. Loss-corrected data from the SMPS and APS were merged using a 2:1 weighting factor favouring the SMPS in the overlapping size range due to low APS counting efficiency below 0.9 µm (Peters and Leith, 2003; Pfeifer et al., 2016). The SM and OPS PNSDs were similarly merged, but with a 1:1 weighting due to a smaller overlapping size range and insufficient data on instrument counting efficiency in this region. Total SSA numbers (N) were estimated by integrating the PNSD. By assuming spherical particles, we also transform the PNSDs into distributions of surface area and volume, similarly estimating total surface area (S) and volume (V) by integrating the respective distributions.

For most of the analysis, measurements were condensed into 5 °C bins covering SSTs from 0 to 25 °C. This amounts to 2 (0–5 °C), 5 (5–10 °C), 2 (10–15 °C), 2 (15–20 °C), and 2 (20–25 °C) experimental runs per temperature bin. Since this process reduces the number of points as a function of SST to n=5, we performed an uncertainty-weighting analysis (see Sect. S4). Applying this weighting factor (Eq. S1) did not dramatically adjust our regressions but increased confidence in our findings (Fig. S4). More information on this analysis and a comparison between the weighted and non-weighted number concentrations can be found in the Supplement. All SST-binned analysis in this study uses the uncertainty-weighting unless otherwise specified, and all following analysis in the main text focuses on measurements from the SMPS and APS covering the size range of 0.01 to 20 µm, while the accompanying analysis of the SM and OPS is included in the Supplement.

2.4 Production Flux Measurements

Estimates of particle number flux distributions (PNFDs) utilized three methods to account for size-dependent influences on SSA production. Certain assumptions are made to simplify the uncharacterized loss mechanisms in SOARS (Leibensperger III et al., 2026b), with all methods assuming a well-mixed reactor and first-order loss rates. Production flux is estimated using Eq. (4) where N is the number concentration at a given aerosol size and time, P is the production flux (# cm⁻³ min⁻¹), and k is the first-order loss rate (min⁻¹).

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\mathrm{d}N}{\mathrm{d}t}}=P-k\cdot N\left(t\right)\end{array}$$

First, submicron PNFDs were estimated by applying a linear regression to the SSA number concentration time series during the initial 5 min of Wind & Waves, referred to as the “initial rise” method. This method assumes a negligible loss rate (k=0) based on the sufficiently long residence times of submicron SSA and sufficiently low background aerosol concentrations, maintaining a large signal-to-noise ratio. The efficacy of this method is evaluated in a separate manuscript, but it showed that reducing the 5 min window to 12 s (the travel time of a particle throughout SOARS) yielded changes to submicron emission fluxes within typical daily variability, indicating that the simplified assumptions do not significantly bias the PNFDs.

Second, PNFDs for supermicron SSA between 1 and 4 µm were estimated for each size bin by applying a power law regression fit to the entire time series of SSA number concentrations during Wind & Waves, referred to as the “power-fit” method. This power law fitting allows P and k to each be freely fit according to each size bin. This method relies on the entire Wind & Waves period to estimate production and loss rates but requires a reduction of the number of size bins to maintain high counting statistics during production. Supermicron SSA reach steady-state concentrations quickly, limiting the validity of this approach for larger particles. Estimates from the power-fit method align well with those from the initial rise method but tend to be lower by about 5 %–10 %.

Third, for SSA above 4 µm, steady state was reached within the first minute of sampling due to short residence times. The PNFDs of these particles were estimated assuming that the production flux is equivalent to the product of the steady-state concentration and the loss rate estimated during the Decay period, referred to as the “steady state” method. This is to say the left-hand side of Eq. (4) becomes 0 and N(t) becomes N_SteadyState, assumed invariable with time. When loss rates could not be directly quantified due to poor counting statistics, a constant loss rate from the last measured size bin was assumed (k is constant). Size- and SST-dependent loss rates exhibit an asymptotic behaviour above 2 µm (approaching a value of 0.55 min⁻¹) and differences within error between experiments (Fig. S12).

The steady state-derived supermicron fluxes are not inherently equal to the initial rise-derived submicron fluxes, but each method is the best approach for the specified size range given the complexities of the SOARS instrument. The loss rate assumptions are likely underestimates, generally leading to underestimates by each method compared to the true interfacial emission flux. However, when the initial rise and steady state methods are compared to one another, the chosen method in each size range represents an overestimate compared to the other method. That is, the initial rise-derived fluxes are larger than the steady state-derived fluxes for submicron SSA (by up to a factor of about 6) and the steady state-derived fluxes are larger than the initial rise-derived fluxes for supermicron SSA (by up to a factor of about 5) (Fig. S5). This variability is of similar magnitude to daily variability of SSA number concentrations in SOARS (Leibensperger III et al., 2026b), reinforcing the validity of these estimates. Critically, this variability is systematic across SST values, such that relative SST-dependent trends are unaffected. A thorough comparison of the flux approximation methods is the focus of a forthcoming manuscript. In the remainder of the text, we combine the three size-differentiated approaches into a single PNFD covering the entire size range (0.01–20 µm) following Eqs. (5)–(7) where P is converted to units of # m⁻² s⁻¹using the headspace volume and water surface area of SOARS:

$$\begin{array}{}\text{(5)}& {\displaystyle }{P}_{D<\mathrm{1}\mathrm{\mu }\mathrm{m}}={\displaystyle \frac{\mathrm{d}N\left(D\right)}{\mathrm{d}t}}\text{(6)}& {\displaystyle }{P}_{\mathrm{1}<D<\mathrm{4}\mathrm{\mu }\mathrm{m}}=k\left(D\right)\cdot N_{\mathrm{Steady}\mathrm{State}}\left(D\right)\text{(7)}& {\displaystyle }{P}_{\mathrm{4}\le D<\mathrm{20}\mathrm{\mu }\mathrm{m}}=k\cdot N_{\mathrm{Steady}\mathrm{State}}\left(D\right)\end{array}$$

Additionally, the PNFDs were used to estimate the particle mass flux distributions (PMFDs). The PNFDs were multiplied by the spherical volume of SSA at each size and assumed a constant density of 2.2 g cm⁻³ (Brock et al., 2019; Wu et al., 2020). This is a reasonable assumption given the water filtration used here to reduce biologic and anthropogenic influences on SSA chemical mixing state. Total number fluxes (F_N) and total mass fluxes (F_M) are the result of integrating the PNFDs and PMFDs. The PMFDs (along with estimates of S and V) are sensitive to supermicron SSA, which are scarcer than submicron SSA. There are known challenges with sampling supermicron SSA in SOARS, so our estimates likely represent lower limits of these quantities (Leibensperger III et al., 2026b).

2.5 Bubble Measurements

Whitecap fractions (WCF) in SOARS were measured via image analysis of wave breaking. During Wind & Waves, a Nikon D3200 SLR camera (shutter speed 1/500 s, f-stop 5.6, ISO 1600) photographed whitecaps from atop the SOARS lid. Each sequence occurred over 12.5 s, captured 27 images, and was repeated for at least ten breaking events. The field of view, located 11 m downwind of wave generation, was illuminated through the side wall of the channel. Background images of still water were subtracted from whitecap images to eliminate the influence of lighting inconsistencies and image deformation. Whitecaps in the resulting images were identified using a thresholding technique. An average WCF value was computed across breaking events for each temperature. Measurements for SSTs below 6 °C were not possible due to significant fogging of the SOARS lid.

Total subsurface bubble concentrations (n) with radii from 6.17 to 830 µm were measured acoustically. Measurements were based on Vagle and Farmer (1998), employing two submerged transducers. The first emits sound between 4 and 500 kHz during a 3 s interval, while the second receives the attenuated signal. Signal attenuations at several frequencies are inverted into a bubble size distribution using the resonant bubble approximation (Czerski, 2012). The bubble size distribution was then integrated across the size range to estimate total concentrations.

3 Results and Discussion

3.1 A Simple Framework of SSA Production

To begin, we present a simple framework to contextualize our measurements and summarize the current understanding of SSA production, using notation that is conserved throughout the remainder of this work. This framework functions as a conceptual interpretation tool, where the equations will not be solved empirically as more measurements or assumptions are necessary for full solutions. This framework also serves to highlight areas for future research where the role of SST in controlling bubble-bursting dynamics is understudied and difficult to constrain across bubble sizes and small-scale processes. This framework is one interpretation to describe the relevant processes and motivate how we approach the problem. Several terms included in the equations below are functions of variables in addition to SST (e.g., salinity, air-water temperature differential, diversity and concentrations of surfactants) but will not be thoroughly described as such. We aim to keep the framework simplified to ensure a consistent view of the SST dependencies of the OAI.

The change in net amount of SSA at a given D_p,dry that enters the atmosphere over a given ocean area (N₀, [#_SSA, m⁻² µm⁻¹]) produced from surface bubbles (ñ, [#_bubble m⁻² µm⁻¹]) over time is characterized by Eq. (8):

$$\begin{array}{}\text{(8)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}{N}_{\mathrm{0}}(D_{\mathrm{p},\mathrm{dry}},t)}{\mathrm{d}t}}& =\underset{r_{\min }}{\overset{r_{\max }}{\int }}E(D_{\mathrm{p},\mathrm{dry}},r)\times B\left(r\right)\times \mathit{ñ}(r,t)\times \mathrm{d}r\\ & -D_{\mathrm{1}}\left(D_{\mathrm{p},\mathrm{dry}}\right)\times N_{\mathrm{0}}(D_{\mathrm{p},\mathrm{dry}},t),\end{array}\end{array}$$

where E is the number of SSA emitted per bubble [#_SSA #${}_{\mathrm{bubble}}^{-\mathrm{1}}$], B is the burst rate of surface bubbles [s⁻¹], r is the bubble radius [µm], D₁ is the deposition rate of SSA [s⁻¹], and t is time [s]. Liu and Yang (2022) include a useful schematic relating changes in SST to potential changes in whitecaps, which is related to the ñ term in Eq. (8). Our framework extends that schematic to include the entire OAI system through a simple mathematical lens. Importantly, the N₀ term in Eq. (8) is related to, but distinct from, the aerosol concentration (N) and aerosol flux (F_N) measurements discussed in the following sections.

Based on previous work, we expect that E, B, and ñ will vary as a function of SST due to changes in seawater and bubble properties (Nielsen and Bilde, 2020; Sellegri et al., 2023; Song et al., 2023). These properties include: the bubble film drainage (φ), bubble film thickness (η), bubble film brittleness (β), air entrainment (ϵ), kinematic viscosity (ν), surface tension (γ), and the solubility of gases (κ), each of which are functions of temperature (Deike et al., 2022).

The number of SSA ejected per burst bubble (E) is known to be proportional to β while inversely proportional to η, φ, and ν, and can vary by a factor of 58 over an SST range of 0 to 19 °C (Dubitsky et al., 2024; Nielsen and Bilde, 2020; Walls and Bird, 2017):

$$\begin{array}{}\text{(9)}& E(D_{\mathrm{p},\mathrm{dry}},r)\propto {\displaystyle \frac{\mathit{\beta }\left(r\right)}{\mathit{\eta }\left(r\right)\mathit{\phi }\left(r\right)\mathit{\nu }}},\end{array}$$

The number of SSA ejected per bubble is different for film and jet drop production, and thus dependent on D_p,dry of SSA and r of bubbles. With increasing SSTs, jet drop production is likely to decrease while film drop production is likely to increase, leading to aerosol- and bubble-size-dependent shifts to E (Deike et al., 2022). Film drop enhancement is due to increasing η with increasing SSTs, although previous laboratory studies focus on very warm water temperatures (10–77 °C) (Poulain et al., 2018). Future studies should therefore focus on thorough characterization of film and jet drop production across realistic SSTs (0–40 °C).

Surface bubble burst rate (B) is proportional to surface tension and inversely proportional to viscosity (Nielsen and Bilde, 2020; Poulain et al., 2018):

$$\begin{array}{}\text{(10)}& B\left(r\right)\propto {\displaystyle \frac{\mathit{\gamma }}{\mathit{\nu }}},\end{array}$$

In addition to bursting, surface bubbles can coalesce into rafts wherein neighbouring bubbles may merge, altering the surface bubble size distribution, thereby altering the resulting SSA produced (Néel and Deike, 2021). The surface bubble concentrations (ñ) will vary non-linearly with SST, depending on the indirect effects of water properties, following the relations given in Eqs. (S2)–(S8). Subsurface bubbles also interact with one another, altering their size-dependent concentrations and likely influencing SSA production (Jiang et al., 2024).

The lack of consensus on the SST dependence of SSA production stems in part from the competing physical processes described in Eqs. (8)–(10) and in the supplement (Sect. S5). For example, increasing SST simultaneously enhances bubble burst rates (Eq. 10) and reduces the number of SSA ejected per bubble (Eq. 9), and these opposing effects are functions of bubble size. Without isolating and quantifying each process independently, observations of SSA production will reflect the sum of these competing influences, making it difficult to attribute trends to any single mechanism (Sellegri et al., 2023; Sofieva et al., 2022). Previous mechanistic frameworks of SSA production provide greater detail of the processes summarized here, but lack experimental support for SST effects (Deike et al., 2022). Characterizing the relationship between surface and subsurface bubbles with changing SST is therefore a critical next step toward disentangling these effects.

This framework is a starting point for the following discussion as the processes and typical measurements span different length scales: SSA are microscopic particles (10 nm–50 µm) produced from bubbles (1 µm–10 mm) contained within wave breaking events (0.1–10 m) (Deike, 2022). While we point to the importance of continued research at the smallest scales (e.g., single-bubble bursting experiments), here we focus on the sum of the processes across entire breaking wave events. Breaking waves in SOARS span the width of the channel (2.4 m) and break along most of the length of the channel (approximately 25–30 m). We therefore analyse the aerosols produced over whole breaking events through measurements of surface and subsurface total bubble concentrations, SSA concentrations, and SSA emission flux.

Figure 2 Trends in whitecap fraction and subsurface bubble concentrations. (a) Whitecap fraction and (b) subsurface bubble concentrations as a function of SST.

3.2 Surface and Subsurface Bubbles

Whitecap fractions (WCFs) were measured to estimate the effect of changing SST on the surface expression of subsurface bubble plumes at SSTs above 6 °C. The WCFs displayed here relate to the surface bubble concentrations (ñ) in Eqs. (8) and (S2), although the ñ term is explicitly the concentration of surface bubbles, which is a function of bubble radius. Whitecaps consist of layers of bubbles near the surface contributing to the high albedo, the measured quantity, whereas the ñ term focuses on estimates of the top-most layer of bubbles directly contributing to SSA production. In Fig. 2a, different results were obtained during warming and cooling experiments and no clear dependence is seen between WCF and SST (R²≤0.53, p>0.1), which is consistent with previous findings (Jia and Zhao, 2019; Liu and Yang, 2022). On the other hand, a previous laboratory study using plunging waterfalls found an increase in air entrainment and surface foam with increasing SSTs over a similar range (6 to 30 °C) (Callaghan et al., 2014). While these results approximately align with the warming results in Fig. 2a, our laboratory proxy utilizes wind-wave interactions, likely contributing to the differing results during the cooling experiment.

Subsurface bubble concentrations (n, Eqs. S2–S3) provide additional insight into how SST modifies the OAI. Figure 2b shows that n decreased by a factor of 1.5 over our SST range, regardless of whether the SST was increasing or decreasing, and the trend is captured by a linear fit (R²=0.83, p≪0.01). These results show the same trend as the work of Thorpe et al. (1992), who found that n decreased by 50 % for every 10 °C increase in SST over the bubble-size range of 5 to 145 µm. However, our observed decrease in n is lower and may reflect the larger size range measured here (6.17 to 830 µm). Notably, this decrease in subsurface bubble concentrations does not correspond to a clear change in WCF, different from the consistent findings of Callaghan et al. (2014). This suggests that SSA production under wind-wave interactions may be governed by processes occurring beneath the surface (Chu et al., 2025) or by bubbles smaller than 6 µm.

Figure 3 SST trends of SST-binned, uncertainty-weighted PNSDs. (a) Mean PNSDs for five SST ranges measured by the SMPS-APS. (b) Total Aitken, (c) accumulation, and (d) supermicron mode number concentrations. The shaded regions in a and error bars in (b)–(d) indicate one standard deviation.

3.3 Particle Number Size Distributions

Mean PNSDs were calculated during the last hour of each Wind & Waves run and binned in 5 °C increments from 0 to 25 °C. Increasing SST from 2 to 23 °C produced a clear suppression of submicron SSA concentrations (Fig. 3), consistent with the diffuser experiments of Christiansen et al. (2019). However, diffuser-generated PNSDs exhibit mode diameters approximately 50 nm smaller than ours, a discrepancy attributable to differing SSA generation mechanisms. In contrast, the plunging jet system of Christiansen et al. (2019) produces mode diameters like ours but yields a U-shaped SST dependence for submicron SSA concentrations. These method-dependent inconsistencies reinforce the critical importance of representative SSA production mechanisms (Collins et al., 2014). By generating SSA through realistic wind-wave interactions, SOARS minimizes the artifacts introduced by other laboratory proxies, lending greater confidence to the monotonic suppression observed here for submicron SSA. These trends are further corroborated by independent PNSD measurements from the SM-OPS system (Fig. S6, Table S3).

Beyond the overall suppression in number concentration, we evaluated whether SST also shifts the shape of the PNSD by examining the temperature dependence of geometric mean diameter (D_g), geometric standard deviation (σ_g), and mode diameter (D_mode) (Fig. S7). None of these parameters exhibit a clear trend with SST, in contrast to the strong linear decrease observed in total number concentration (N, R²=0.99, p<0.05). This indicates that increasing SST reduces the amount of SSA produced without substantially altering the size distribution at steady state. The potential for SST to modulate the chemical composition and mixing state of emitted SSA, a subtler effect that would not be captured by these bulk size metrics, is the focus of ongoing single particle analysis to be the focus of a forthcoming manuscript.

Integrated number concentrations over three size ranges clarify size-dependent relationships with SST that are not captured by evaluating the complete PNSD. Suppression of SSA with increasing SST for both Aitken ($\mathrm{0.02}<D_{\mathrm{p},\mathrm{dry}}<\mathrm{0.1}$ µm) and accumulation ($\mathrm{0.1}<D_{\mathrm{p},\mathrm{dry}}<\mathrm{1}$ µm) modes is well represented by a linear regression (R²≥0.97, p<0.05 Fig. 3b–c). On the other hand, supermicron SSA exhibits a non-linear relationship with SST, which is well represented by a quadratic regression (R²=0.98, p<0.05). These trends were consistent across cooling and warming experiments (Fig.  S10) and suggest a mechanistically driven (i.e., film vs. jet drop) size-dependent shift in production. This indicates that increasing SSTs have a complex influence on the resulting SSA and emphasizes the importance of size-resolved measurements. For example, film and jet drop production may have different SST dependencies, causing size-binned N to vary with SST while the shape of the distribution (a consequence of the ratio of total film to jet drop production) remains constant. More work is needed to isolate the SST dependencies of film versus jet drop production. Subtle shifts in supermicron SSA production are not captured when assessments are made on bulk changes to SSA number concentrations. Although supermicron concentrations are limited in SOARS (Leibensperger III et al., 2026b), supermicron N estimates are dominated by SSA between 1–2 µm, where counting statistics and counting efficiency remain high, providing robust results.

Figure 4 Extreme SST-binned mean PNSD modal decomposition. Average PNSD for the coldest and warmest binned SSTs. Each is decomposed into three log-normal modes and reconstructed.

A tri-modal log-normal fit was applied to the coldest and warmest PNSDs to further characterize the changes in the modal distributions as a function of SST (Fig. 4, Table S4). From the lowest to the highest SSTs, the greatest change occurs in the accumulation mode number concentration, which decreases by 75 %. Beyond a decrease in magnitude, each mode experiences increases in D_g and σ_g. Aitken mode concentrations decrease by 35 %, D_g increases by 22 %, and σ_g increases by 4 %. The accumulation mode N decreases by 75 %, D_g increases by 18 %, and σ_g increases by 1 %. The supermicron mode N decreases by 36 %, D_g increases by 2 %, and σ_g increases by 5 %. The less pronounced changes in D_g and σ_g again suggest SST primarily influences the number of SSA at steady state rather than the shape of the distribution.

Figure 5 Higher order moments of the SST-binned, uncertainty-weighted PNSDs. (a) Total number, (c) surface area, and (e) volume of SSA in SOARS for each temperature range, with a dashed line representing a linear least squares best fit line. (b) Number, (d) surface area, and (f) volume size distribution for each temperature range.

Surface area and volume size distributions were derived from PNSDs, assuming spherical particles, to assess the influence of changing SST on higher-order moments of the size distribution. Figure 5 shows that increasing SST drives a decrease in total SSA number, surface area, and volume by more than 50 %, across the SST range of 2 to 23 °C. As with modal number concentrations, suppression of S and V was observed during both experiments (Fig. S11). At the SST-binned scale, all three moments exhibit strong, negative correlations with SST, which are well represented by linear regressions (R²≥0.98, p<0.05).

The higher-order moments of the size distribution reveal a structural shift that the number distribution alone does not capture (Fig. 5). As SST increases from 2 to 7 °C, the surface area distribution transitions from a broad, log-normal mode centred near 0.3 µm to a sharper mode at 0.55 µm, which then persists across all higher SSTs (Fig. 5d). The volume distribution undergoes a parallel transition, with a secondary mode emerging near 1 µm over the same temperature range (Fig. 5f). Although this secondary mode falls near the SMPS-APS merging point, it also appears in the SM-OPS PNSDs (Fig. S6) where merging occurs at 0.4 µm, confirming it is not an artifact of measurement stitching. These modal shifts are also present in the probability density functions of the number, surface area, and volume of the coldest and warmest SSTs (Fig. S8).

These shifts carry direct implications for climate; particle surface area underpins several modelling parameterizations for ice nucleating particle concentrations (DeMott et al., 2016; Moore et al., 2022; Phillips et al., 2008), so a redistribution of SSA surface area into larger sizes could alter the available sites for ice nucleation and, in turn, the abundance of ice and mixed-phase clouds. Simultaneously, the suppression of SSA volume combined with shifts in mode diameter could reduce direct scattering of solar radiation, potentially enhancing poleward heat transport as described by Kay et al. (2016).

The abruptness of these transitions between 2 and 7 °C, followed by relative stability at higher SSTs, points to a possible physicochemical threshold followed by a gradual scaling with temperature. One potential explanation is the influence of marine gels. Despite water treatment to minimize biological contributions (Sect. 2.1), dissolved organic matter and its gel-phase products can persist. Since marine gels are temperature sensitive, with significant hydrodynamic diameter growth over the SST range covered here (Chang et al., 2022; Verdugo, 2012), they are likely to become most prevalent in the seawater (and therefore enriched in SSA) at the coldest SSTs. Nielsen and Bilde (2020) reported that removing gel-stabilizing divalent cations (Ca²⁺, Mg²⁺) from sea salt mixtures produced a similar suppression of particle concentrations as SSTs increased. Another explanation is a shift in physical production mechanisms modulated by SST. Callaghan et al. (2014) reported that temperature had the greatest effect on subsurface bubbles (n) larger than 1 mm, which are expected to produce submicron film drops (de Leeuw et al., 2011). The rapid decay of larger bubbles suggests that bubble rise times at low SSTs are accelerated. Furthermore, Nielsen and Bilde (2020) support the conclusion that the number of SSA produced per bubble is highest at low SSTs. However, no previous laboratory studies have incorporated wind, which significantly influences SSA concentrations (Leibensperger III et al., 2026b). Consequently, these results are expected to be the most representative of open ocean conditions and SOARS represents an opportunity to test these hypotheses in the future.

Figure 6 Evaluation of the bubble-SSA correlation using SST-binned, uncertainty-weighted SSA concentrations. Comparison of SSA total number (N) to subsurface bubble concentrations (n), separated by SSA mode: (a) Aitken, (b) accumulation, and (c) supermicron. Weighted least squares regression and resulting best fit lines are shown for linear (dashed) and quadratic (dotted) fits.

3.4 Connecting SSA and Bubble Concentrations

We compare SSA number concentrations to subsurface bubble concentrations to contextualize which quantities exhibit similar SST-dependencies (Fig. 6). The sub-millimetre subsurface bubble concentrations (n) are likely to be related to jet droplet production, controlling supermicron SSA production and up to 43 % of submicron SSA production (Wang et al., 2017). However, the highest correlation occurs between the accumulation mode SSA and subsurface bubbles (R²=0.99, p<0.05, Fig. 6b). Supermicron SSA exhibits a quadratic relationship with subsurface bubbles (R²=0.93), but the result is not significant at the 95 % confidence level (p=0.07). More work is needed to understand the non-linear processes connecting changing SSTs to SSA production via size-dependent bubble processes. This analysis provides a cursory look into the temperature-dependent relationship between subsurface bubbles and SSA.

Figure 7 Extreme SST-binned mean PNFD modal decomposition and SST-binned, uncertainty-weighted F_N trends. (a) Average PNFD for 2 and 23 °C. Each PNFD is decomposed into three log-normal modes and reconstructed. Weighted least squares linear (dashed lines) and quadratic (dotted lines) regressions of total number flux for the (b) Aitken, (c) accumulation, and (d) supermicron modes.

While SSA production has been investigated with various breaking wave analogues in the laboratory (Bowyer et al., 1990; Forestieri et al., 2018; Salter et al., 2014), experiments on SSA production should be conducted with an accompanying emphasis on understanding the dynamics and evolution of the bubble plume, similar to Callaghan et al. (2014). The structure of the subsurface bubble plume is governed by the fluid dynamics within the entire water column, which is dependent on winds and waves (Crawford and Farmer, 1987). The air-side fluid dynamics are also dependent on the method of air entrainment (plunging jet vs. spilling breaker, e.g., Yang et al., 2018), and the spilling breakers in SOARS likely generate dynamic differences from the plunging sheets used previously (Callaghan et al., 2014; Forestieri et al., 2018). The fluid dynamics of each system will vary owing to the system dimensions and presence or absence of wind, influencing the SSA deposition rate, D₁ (Eq. 8).

Table 1 Least squares regression and coefficients of determination for linear and quadratic fits to total number fluxes. Astricts indicate a statistical significance of p-values exceeding the 95 % (^*) and 99 % confidence levels (${}^{**}$).

3.5 Particle Number Flux Distributions

Emission fluxes of SSA in SOARS also decrease with increasing SST, but with a striking size-dependence. Size-resolved PNFDs at the lowest and highest SSTs (Fig. 7a) reveal a crossover: SSA smaller than 0.6 µm are emitted at lower rates at high SSTs, while SSA larger than 0.75 µm show the opposite trend. Tri-modal fitting quantifies this divergence, where total number flux (F_N) for the Aitken and accumulation modes decreases by 82 % and 58 %, respectively, while supermicron F_N increases by more than a factor of three (Table S4). This size-dependent behaviour suggests that the mechanisms governing film drop production (predominantly submicron) and jet drop production (predominantly supermicron) respond differently to changes in SST. Linear regressions capture the Aitken and accumulation mode trends well (R²≥0.99, p≪0.01), while the supermicron mode requires a quadratic fit, consistent with its non-monotonic SST dependence (Table 1, Eqs. 11–12). These regressions were evaluated separately for each experiment and for all measurements combined (Fig. S13) to account for the limited number of temperature bins (n=5).

$$\begin{array}{}\text{(11)}& {\displaystyle }{F}_{\mathrm{N}}=a_{\mathrm{1}}\times \mathrm{SST}+a_{\mathrm{2}},\text{(12)}& {\displaystyle }{F}_{\mathrm{N}}=b_{\mathrm{1}}\times {\mathrm{SST}}^{\mathrm{2}}+b_{\mathrm{2}}\times \mathrm{SST}+b_{\mathrm{3}},\end{array}$$

Aitken and accumulation mode SSA number flux SST dependencies exhibit trends well captured by both linear and quadratic regressions (R²≥0.99, p≪0.01), while the SST dependence of supermicron SSA is only sufficiently captured by a quadratic regression (R²=0.96, p<0.05). SSA number emissions are dominated by submicron SSA (Fig. 7a), so a first-order SST dependency is recommended for computational efficiency during model integration. Further, a monotonically decreasing SST dependence for total SSA number flux agrees with previous estimates (Grythe et al., 2014; Lehahn et al., 2014; Markuszewski et al., 2024). This first-order relationship forms the basis of the temperature-dependent number flux correction factor (CF_N), developed in Sect. 3.6 for direct integration into ESMs.

To ensure compatibility with SSA modelling schemes that parameterize emissions by mass and compare to previous measurements of mass fluxes, the PNFDs and F_N values are converted to particle mass flux distributions (PMFDs) and total mass flux (F_M), assuming spherical particles and a density of 2.2 g cm⁻³ (Brock et al., 2019; Wu et al., 2020) (Figs. S9, S14). The mass flux results diverge markedly from the number flux results. Because particle mass scales with the cube of diameter, F_M is dominated by the largest SSA in each mode – particles whose number emission rates increase with SST (Fig. 7a). As a result, while the Aitken mode PMFD mirrors its PNFD counterpart, both the accumulation and supermicron mode PMFDs exhibit non-monotonic SST dependencies with maxima near 13 °C, resembling the supermicron PNFD behaviour. Note that these PMFD trends differ from the $\mathrm{d}V/\mathrm{dlog}{D}_{\mathrm{p},\mathrm{dry}}$ curves in Fig. 5f, which reflect steady-state concentrations rather than emission fluxes.

Figure 8 Correction factors applied to the complete sets of F_N and F_M measurements. (a) Total number emission flux (F_N) and (b) total mass emission flux (F_M) as a function of SST, normalized by the average flux at 14 °C. The black diamond depicts flux at 14 °C, and the vertical error bars represent one standard deviation. The horizontal error bars represent the annual SST range of the coastal seawater sourced for this study.

While a direct comparison of SSA flux estimates from SOARS to literature values is evaluated in a separate publication, it's important to contextualize the SST-dependent flux estimates determined in this work. The F_M estimates vary between $\mathrm{0.63}\times {\mathrm{10}}^{-\mathrm{8}}$ kg m⁻² s⁻¹ at 2 °C to $\mathrm{3.4}\times {\mathrm{10}}^{-\mathrm{8}}$ kg m⁻² s⁻¹ at 14 °C, corresponding to approximate annual global SSA fluxes of 72–383 Pg yr⁻¹. These values lie within the bounds of previous studies (Grythe et al., 2014), but are greater than most estimates. Therefore, we normalize the SST-dependent fluxes by the estimate at moderate SSTs to capture the SST trends. The opposing SST dependencies of F_N and F_M – one monotonically decreasing, the other non-monotonic – necessitate separate correction factors for number and mass emissions in models. A temperature-dependent mass flux correction factor (CF_M) is developed alongside CF_N in Sect. 3.6.

3.6 Correction Factors (CFs) for Model Integration of SST Dependence

The uncertainties of simulated SSA in ESMs can be reduced by incorporating SST-dependent emission correction factors derived from the SOARS total number flux (CF_N) and mass flux (CF_M) measurements. Critically, these correction factors arise from regressions fit through normalized production fluxes to remove any dependence on absolute emission rates in SOARS, making them applicable to any SSA emission scheme (Fig. 8). F_N at each SST is normalized by the average F_N at 14±3 °C ($\mathrm{8}\pm \mathrm{2}\times {\mathrm{10}}^{\mathrm{6}}$ #_SSA m⁻² s⁻¹), a reference temperature chosen to represent the annual SST range of the local seawater used in this study (Espinosa-Carreón et al., 2001). CF_M is generated through a parallel normalization of F_M at 14±3 °C ($\mathrm{3.4}\pm \mathrm{0.6}\times {\mathrm{10}}^{-\mathrm{8}}$ kg m⁻² s⁻¹). CF_N and CF_M are represented by Eqs. (13) and (14), respectively:

$$\begin{array}{}\text{(13)}& {\displaystyle }{\mathrm{CF}}_{\mathrm{N}}=-\mathrm{0.064}\times \mathrm{SST}+\mathrm{1.950}\text{(14)}& {\displaystyle }{\mathrm{CF}}_{\mathrm{M}}=-\mathrm{0.006}\times {\mathrm{SST}}^{\mathrm{2}}+\mathrm{0.152}\times \mathrm{SST}-\mathrm{0.035}\end{array}$$

These correction factors are chosen as simple, yet encompassing, representations of the SST dependencies of total F_N and F_M. They rely on all the F_N and F_M measurements made across the two experiments and are not binned by SST. The linear CF_N shows good correlation (R²=0.88, p≪0.01) and is valid across the SST values ranging from freezing (approximately −1.8 °C) to 30 °C. On the other hand, the quadratic CF_M shows adequate correlation (R²=0.61, p<0.01) and is valid for SST values between 0.25 and 25 °C. These validity bounds are determined by when the CFs become negative; however, we advise only adapting them to SSTs between 2 and 25 °C, consistent with the conditions tested here. These CFs mirror the trends of the dominant modes, where F_N is dominated by the accumulation mode, and F_M is dominated by the supermicron mode. While they do not account for the size-dependent differences highlighted in this work, they represent a first step towards integrating SOARS measurements into ESMs to improve the simulation of SSA.

Incorporating these correction factors into ESMs would substantially alter simulated SSA burdens, consistent with sensitivity analyses in previous work (Grythe et al., 2014; Jaeglé et al., 2011; Ovadnevaite et al., 2014; Sofiev et al., 2011). The magnitude of this impact is underscored by the fact that only half of the emission schemes that currently incorporate an SST dependence predict the correct sign as measured here (Grythe et al., 2014; Mårtensson et al., 2003; Salter et al., 2015). The other emission schemes simulate increasing SSA production with SST. The CF_M from SOARS introduces an additional sensitivity at moderate SSTs, where small deviations from approximately 13 °C will lead to large decreases in mass flux, a regime that encompasses much of the global ocean. Because Eqs. (13)–(14) are derived from SSA produced by realistic wind-wave interactions, they represent the most physically grounded, laboratory-based SST correction factors currently available for ESM integration. Once integrated, these CFs would deviate from the currently utilized SST dependencies of Jaeglé et al. (2011) and Sofiev et al. (2011), referred to here by the abbreviations J11 and S11. The CF_M of J11 is similar to our CF_M for SSTs less than 13 °C. For warmer SSTs, our CF_M would drastically decrease mass production, likely decreasing sea salt aerosol optical depth (SSAOD) over equatorial regions. The S11 CFs mimic our findings in the Aitken mode (increasing F_N with decreasing SST) but deviate in the accumulation and supermicron modes showing an increase in F_N with increasing SSTs. Therefore, replacing the S11 CFs with ours would decrease SSA number concentrations over warm water, likely outside most of the regions they identified for comparison to satellite retrievals of SSAOD.

These CFs are designed so that the SOARS-derived SST relationships can be adopted to emission schemes regardless of whether they are estimating number or mass (volume) flux. However, one limitation of the presented CFs is that they are not mutually consistent when applied to the same emission scheme. That is to say that applying the CF_N to an SSA number flux emission scheme and then calculating SSA mass produces a different answer than converting the PNFD to a PMFD and then applying the CF_M. This is due to the drastically different relationships of SSA number and SSA mass. These mismatches can be as large as a factor of 8 for SSTs around 2 °C, decrease to a minimum of about 1.05 for SSTs around 18 °C, and increase to about a factor of 3 for the highest SSTs. These differences highlight the variability arising from number versus mass measurements and the importance of retaining quantity consistency (i.e., number, surface area, volume, or mass) between measurements, emission schemes, and correction factors. However, these differences are within typical variability across SSA emission schemes and SST dependencies. Therefore, the CFs improve emission schemes without an SST dependence or with one of the wrong sign; users should apply the CF that matches their emission scheme quantity (number or mass) and not cross-convert.

3.7 Temperature Hysteresis Effects

We also characterize the effect of SST hysteresis (whether the SST dependence is different when SSTs are increasing or decreasing across a common range) on the relationship to SSA production. Section 3.3 discussed the influence of SST on SSA concentrations for SST-binned data, demonstrating a linear relationship between concentrations of Aitken and accumulation mode SSA and a quadratic relationship for supermicron SSA. Figure S10 shows the same analysis for the cooling and warming experiments. Binning data by SST produced an average increase of R² by 0.1 for linear models and 0.04 for quadratic models and did not alter which regression type produced the best result. Furthermore, all regressions remained statistically significant at the 95 % confidence level. Regressions comparing N, S, and V also showed slight improvement when data was binned compared to cooling and warming experiments independently (Fig. S11).

The loss rates used in the PNFD calculations were similarly separated by each experiment (Fig. S12). To our knowledge, this is the first estimate of the effect of varied SSTs on aerosol loss rates, a process that is best investigated in laboratory studies. The mean loss rate between the experiments and across sizes was found to be 0.30 min⁻¹, with the warming and cooling experiments approaching asymptotes of 0.38 and 0.30 min⁻¹, respectively. The warming experiment exhibits larger loss rates on average (by approximately 20 %–30 %), though the cooling and warming experiments remain within one standard deviation of each other across all sizes. The small variations in loss rates are three to four times smaller than the variations in F_N, and are therefore not expected to drive the larger-scale findings of the PNFD analysis and are not explored further here. Future studies should target SST-dependent loss rate variability to better understand the influence on flux and concentration estimates.

The SSA flux analyses (Sect. 3.5) were repeated for the cooling experiment, the warming experiment, and both experiments combined (Figs. S13–S14). Interestingly, the warming experiment exhibited better regression statistics with linear fits for accumulation (R²=0.96, p≪0.01) and total (R²=0.91, p≪0.01) F_N compared to the cooling experiment (accumulation: R²=0.84, p<0.05; total: R²=0.87, p<0.05). The supermicron mode quadratic regression exhibited similar R² values between the experiments (warming: 0.79, cooling: 0.80) accompanied by significant p-values (warming: 0.02, cooling: 0.04). Conversely, for F_M, the cooling experiment exhibited better fitting regressions than the warming experiment, with a linear regression fitting the Aitken mode (cooling: R²=0.98, p≪0.01; warming: R²=0.70, p<0.01) and quadratic regressions in the accumulation mode (cooling: R²=0.61, p=0.1; warming: R²=0.46, p>0.1), supermicron mode (cooling: R²=0.99, p≪0.01; warming: R²=0.86, p=0.01), and total SSA mass (cooling: R²=0.97, p<0.01; warming: R²=0.80, p<0.05). While this indicates that temperature hysteresis may influence the strength of the fitted SST dependencies of SSA emission flux, the trends are consistent despite hysteresis. Future work should identify if F_N continues to be better correlated under warming SSTs while F_M is better correlated under cooling SSTs.

We further evaluate the impact of temperature hysteresis on the CFs described in the previous section by generating linear and quadratic CFs for the cooling and warming experiments individually (Fig. S15). The linear form of CF_N remains a statistically significant good fit for both the cooling experiment (R²=0.87, p<0.05) and the warming experiment (R²=0.91, p≪0.01). The quadratic form of CF_M similarly remains a good fit and significant for both experiments: R²=0.97, p<0.01 and R²=0.80, p<0.05 for the cooling and warming experiments, respectively. This reaffirms that the functional forms of the CFs in Eqs. (13)–(14) are correct based on all available SOARS data; however, the coefficients of the CFs vary between each experiment and the general CFs in Fig. 8. The similarity of the regression coefficients was evaluated using the Chow test. The p-value of the linear CF_N regression was 0.79, indicating that the coefficients were not significantly different between the experiments. On the other hand, the p-value of the quadratic CF_M regressions was found to be 0.03, indicating that there were statistically significant differences in the coefficients. Temperature hysteresis may therefore lead to varied mass CFs. However, SST hysteresis is difficult to apply in real world fashions due to variability on timescales of days to weeks and substantial noise in the variability. We therefore recommend the generalized CF_M in Eq. (14) as the best approximation to account for the SST dependence of SOARS-derived SSA mass fluxes in models.

4 Conclusions

This is the first study to generate SSA through wind-wave interactions in controlled laboratory conditions under varied SSTs, using a comprehensive set of measurement techniques to probe the uncertain relationship between SST and SSA production. The results presented here incorporate measurements of whitecap fraction, subsurface bubble concentrations, steady-state SSA concentrations, and SSA emission fluxes to build a robust framework for characterizing the SST dependence of SSA. Wind-driven SSA production leads to SST dependencies inconsistent with previous laboratory studies that replicate wave breaking conditions. This reinforces the importance of producing SSA from wind-wave interactions to clarify the size-resolved, mechanistically driven SST dependencies.

Increasing SST in SOARS from 2 to 23 °C suppresses subsurface bubble concentrations (Fig. 2b), total SSA number concentration (Fig. 3), and total SSA number emission flux (Fig. 7) by factors of approximately 1.5, 3, and 3, respectively. Critically, this general suppression is not size-independent; the shift of the SST dependence of F_N between 0.6 and 0.75 µm emphasizes that no single monotonic SST dependence describes SSA production. The size-resolved SST dependencies differed in strength and shape, suggesting the occurrence of a physicochemical or mechanistic shift. Surface area and volume distributions (Fig. 5) showed pronounced shifts while supermicron SSA fluxes and most SSA mass fluxes show non-monotonic SST dependencies sensitive to the SST value. These findings depend on realistic SSA production methods to corroborate previous estimates (Christiansen et al., 2019; Lehahn et al., 2014; Markuszewski et al., 2024; Nielsen and Bilde, 2020; Salter et al., 2014, 2015; Sellegri et al., 2023).

Although there are system-specific considerations for interpreting these measurements, the use of wind-wave interactions to produce SSA are the most representative of conditions outside of SOARS to date (Leibensperger III et al., 2026b). The general suppression of subsurface bubble concentrations, total SSA number concentrations, and total SSA number fluxes with increasing SST is common across experiments, independent of SST hysteresis. Supermicron F_N and F_M for SSA larger than 0.1 µm are suppressed as SSTs diverge from approximately 13 °C. SST hysteresis can lead to meaningful changes in CF_M coefficients, indicating an open area of future research. These results indicate that aerosol number exhibits a stronger dependence on the magnitude of SST than hysteresis, while aerosol mass (and the supermicron mass-containing aerosols) depend on the direction of SST change coupled with the magnitude of the SST.

Subsequent work should pursue three avenues to better understand the temperature dependence of SSA production and the associated impacts. The first avenue should characterize the physicochemical mixing state of SSA at varied SSTs. Changes in mixing state are likely to accompany the shifts in production identified here and will directly depend on water composition. Aerosol composition shifts with temperature may explain the variability from previous studies and partially reconcile differences between laboratory and field estimates of the SST dependence of SSA production.

The second avenue should seek to better characterize the SST-dependent SSA production at the single-bubble level. We hypothesize that changes in the emission rate of SSA per bubble (E) or in subsurface bubble concentrations (n) are the cause of changes in SSA number, with strong dependencies on aerosol and bubble sizes. While the aerosol size-dependent SST trends are well documented here, more focus on the bubble size-dependent SST dependencies of the terms in Eq. (8) is vital for future work. Laboratory experiments focused on the role of SST in determining single-bubble bursting characteristics are currently missing from the literature and represent a knowledge gap that inhibits our understanding of interfacial fluxes. As mentioned in Sect. 3.3., single-bubble experiments can also isolate the SST dependencies of film drop production and jet drop production, which are needed to better understand the mechanistic drivers of the SST dependence of bulk SSA measurements.

The third avenue should focus on improved simulation of the SSA-SST dependence. Simulations should span the scales from direct numerical simulation of single-bubble processes to Earth system modelling using the CF_N and CF_M described here. Direct numerical simulation of single bubbles and bubble plumes can connect SST-dependencies across scales and improve the framework described in Sect. 3.1. Recent direct numerical simulation of breaking waves emphasizes the need to perform fine-resolution simulations (Yang et al., 2018), where SST-dependent systems should be a primary focus. Incorporation of the SOARS-based CFs could dramatically change simulated SSA production, likely decreasing marine-derived SSAOD over the warmest equatorial water. This change in SSAOD is likely to be even stronger under historical and projected scenarios in which SSTs are likely to change significantly. Improving simulated SSA dependence on SST is critical to understand associated changes in marine aerosol burdens, cloud condensation nuclei concentrations, and radiative forcing estimates.