Marine aerosols produced by ocean internal waves

. The objective of this work was to investigate the hypothesis that oceanic internal waves could modulate the production of atmospheric aerosols by modulating sea surface roughness. Airborne profiling lidar was used to measure the displacement of the subsurface scattering layer by internal waves, the distribution of surface roughness, and the aerosol concentration over the internal waves. The direct correlation between internal wave displacement and aerosol density was not statistically significant, primarily due to other processes causing mixing in the atmosphere. However, we found that the magnitude of the aerosol spatial power spectral density at the internal wave wavelength was significantly correlated with the internal wave magnitude. We developed a simple model to describe the interactions, which provided excellent agreement with the measurements for two out of three flights. The disagreement for the third flight is thought to be a result of the deeper thermocline on that day, which put the deeper parts of the internal waves out of the depth range covered by the lidar. The conclusion is that internal waves affect aerosol production over the ocean under certain environmental conditions. © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 International License. Distribution or repro-duction of this work in whole or in part requires full attribution of the original publication, including its DOI. [DOI: 10.1117/1.JRS.16.024501]


Introduction
Aerosols produced by the ocean affect the composition of the atmosphere and of the ocean. 1 They influence the climate directly by reflecting sunlight. Aerosols also serve as cloud condensation nuclei that amplify the climate effect. By reducing the sunlight reaching the ocean, they also affect photosynthesis in the upper ocean, thereby affecting its composition. For these reasons, processes in the ocean that affect aerosol production are of interest.
Remote sensing of aerosols is generally done optically. This includes passive measurements of aerosol optical depth from the ground 2 or from the satellite. 3 It includes lidar measurements of aerosol profiles and properties, also made from the ground 4 and from the satellite. 5 Internal waves in the ocean have also been studied for some time. 6 Internal waves are a major source of mixing in the ocean, especially from the pycnocline into the deep ocean. 7 This has climate implications on time scales longer than those of aerosols. We note also that internal waves can propagate large distances across the ocean from their point of origin to where they shoal and break, which affects the horizontal distribution of mixing. 8 The most common technique used for remote sensing of internal waves is microwave radar. The surface currents induced by the internal wave modulate the small-scale surface roughness from which the radar is scattered. 9 While internal waves have been detected by land-based 10 and airborne 11 radars, most have used satellite radar data. [12][13][14] From the right vantage point, the surface-roughness modulation is visible to the naked eye, and optical imagery can also be used as a remote sensing technique. Satellite images of internal waves have been used to detect internal waves [15][16][17] and to infer the heat content of the upper ocean. 17 Airborne detection has been done with photographs, 18,19 hyperspectral imagery, 20 and infrared imagery. 21,22 A related optical technique uses lidar that images the sea surface. These can provide an image based on the modulation of the magnitude of the surface return based on surface roughness 23 or an image of sea-surface height on scales that are modulated by the internal wave currents. 19 Finally, a profiling blue/green lidar can detect internal waves by directly measuring the displacement of a scattering layer. This technique has been used to measure the power spectrum of an internal wave field 24 and an optical estimation theory has been developed. 25 It has been used to investigate the nonlinear characteristics of an internal-wave packet, 26 and the total energy contained within another. 18 The theory behind lidar detection of internal waves has also been developed. 27,28 In this paper, we investigated the hypothesis that internal waves in the ocean can modulate the production of atmospheric aerosols. We did this using an airborne lidar to measure both the internal waves and the aerosols. Section 2 outlines the methods. Section 3 presents the results, including a detailed analysis of one example and a statistical analysis of data from three flights. Section 4 is a theoretical investigation of the processes involved, and a comparison with the measurements. Section 5 summarizes the conclusions.

Methods
The data were collected with the National Oceanic and Atmospheric Administration (NOAA) oceanographic lidar. It has been described previously, 29-31 but the main characteristics will be briefly reviewed. The laser was a linearly polarized, flashlamp-pumped, Q-switched, frequencydoubled Nd:YAG that generated 100 mJ, 12-ns pulses of 532-nm light at a repetition rate of 30 Hz. A diverging lens was used to expand the beam to 17-mrad full angle and a pair of steering mirrors was used to overlap the beam with the receiver telescopes at a distance of 300 m.
Two Kepler telescopes collected light at two orthogonal polarizations, one of which was copolarized with the transmitted beam. The telescope diameters were 6 cm for the co-polarized receiver and 15 cm for the cross-polarized receiver. The polarization for each telescope was selected by a sheet of polaroid material on the front. An aperture at the primary focus of each telescope set the field of view of each receiver to be slightly larger than the transmitter divergence. A 3-nm wide interference filter reduced the amount of background light reaching the photomultiplier-tube photodetectors. The detected signals were input to a logarithmic amplifier that provided 80 dB of dynamic range when coupled with the eight-bit, giga-sample-per-second digitizer. A computer recorded 2000 samples of each pulse, global-positioning system information, and the photomultiplier gains.
We collected lidar data from a Basler BT-67 aircraft over the New York Bight. This is a region where internal waves have been studied previously, 12,32,33 so we were confident we would find internal waves there. To reduce specular reflections from the surface, the lidar was installed so the incidence angle was 15 deg from nadir. Flights were made on 3 days, June 25, 2020, July 30, 2020, and August 11, 2020 at altitudes between 150 and 1700 m and a speed of about 40 m s −1 . The location was in the New York Bight (Fig. 1) just south of the Hudson Canyon. The pattern was a series of parallel transects across the shelf break with a spacing of about 7 km.
Conversion of digitization levels to attenuated backscatter values was done using laboratory measurements of the various system component characteristics. These include the transmitter pulse energy and pulse length; the collecting area, field of view, and transmission of each telescope; the photocathode response and gain curve for each photomultiplier tube; and the response of each log-amp and digitizer. The recorded photomultiplier gain settings and the recorded aircraft altitude were also used in these calculations. The maximum value of the calibrated signal in each shot was used to find the sea surface. This value is the volume backscattering coefficient at the surface.
After locating the position of the surface in each lidar return, the data were corrected for the effects of signal attenuation. For the atmospheric data, the effects of scattering and absorption are negligible, except in the presence of low clouds or fog. Any data with low clouds or fog were removed from consideration, and the atmospheric data were corrected for geometric loss.
Because the penetration depth of the lidar was much less than the aircraft altitude, geometric losses were neglected in the oceanic data. The oceanic data were also investigated for possible effects of after-pulsing from the surface return. None were observed, probably because the attenuated signal at the relevant depth was greater than the after-pulse signal. The effects of the exponential attenuation in the water were removed using a perturbation approach described previously. 34 A recent improvement to this approach uses a Klett algorithm to obtain profiles of attenuation and the perturbation approach to obtain profiles of backscattering. 35 Since we are only interested in profiles of backscattering for this work, the perturbation approach is sufficient. The result of this processing was a data set that was free of the effects of fog and low clouds, referenced to the surface of the ocean, and corrected for attenuation in the atmosphere and the ocean.
The processed data were then inspected visually for obvious internal wave signatures. Several examples were further processed to look at correlations between the co-polarized atmospheric return at a height of 5 m above the surface, co-and cross-polarized surface returns, and the depth of the scattering layer inferred from the cross-polarized oceanic return. To allow for possible effects of currents and winds, the correlations were calculated for various delays between one parameter and the other. The spatial power-spectral density of each of the four parameters was also calculated. The results for one of these example internal wave packets are presented in this paper. The next step was to automate the spectral analysis. The surface was taken to be the point of maximum signal in each shot. The layer depth was taken to be at the maximum signal below the surface. Files were removed from consideration by visual inspection if clouds were misidentified as sea surface or if there were double layers. In the latter case, the algorithm would identify one layer some of the time and the other layer the rest of the time producing depth variations that could be mistaken for waves. The depths were filtered with a nine-shot median filter. Spectral analysis was performed on time series of layer-depth data that were 200-shots long. Within each data segment, the data were detrended, a Hann window was applied, and the power spectral density was calculated. Successive data segments were separated by 100 shots. If an internal wave train ran off the end of one segment, this overlap should put it near the center of the next one, so the full extent can be captured.
Within each window, the spatial frequency with the maximum spectral density was assumed to be that of the dominant internal wave in that window. The power spectral density was calculated in the same way for three parameters-the co-polarized surface return, the crosspolarized surface return, and the co-polarized return from the atmosphere at a height of 5 m above the surface. This height was chosen to be high enough that no contamination from surface returns was possible. From each of the three power spectra, the value at the dominant internalwave frequency was extracted. The reasoning behind this was that the energy in the spatial frequency was not sensitive to the effects of wind and currents, and the correlations between parameters at the same spectral component would not depend as much on other sources of variability.
We also considered the effects of stratification, as measured by Airborne eXpendable BathyThermographs (AXBT) deployed by the aircraft. Two profiles were measured on June 25, 2020, three on July 30, 2020, and three on August 11, 2020. As a measure of the strength of stratification, we used the density difference, Δρ, between the surface and a depth of 100 m. Density was calculated from the measured temperature using 36 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 5 9 2 (1) where ρ is the density, T is the temperature in°C, and z is the depth in meters. A salinity of 35 practical salinity unit has been assumed in this equation. As a measure of the depth of the mixed layer, we used the depth where the gradient of temperature was greatest. This corresponds to the depth of maximum Brunt-Väisälä frequency. Error bars on the AXBT data were calculated as the standard deviation of the two or three casts on each day. Meteorological data were downloaded from NOAA station 44017 at 40.693°N 72.049°W (location shown in Fig. 1, Ref. 37). We considered wind speed and direction and air-sea temperature difference, which is a measure of atmospheric stability, for the three flight days.
For comparison, we also obtained wind speed and air temperature from an aircraft-integrated meteorological measurement system (AIMMS) and sea-surface temperature from a thermal radiometer on the aircraft. To represent values near the surface, AIMMS values were taken from the lowest flight altitude of 46 m.

Example
An example is presented of data collected at 1816 UTC near the continental shelf break 40-km southwest of the Hudson Canyon. Figure 2 shows both atmospheric and oceanic returns from this region. The subsurface return in the figure shows the cross-polarized channel. The plankton layer that the internal wave is perturbing is at a depth of about z ¼ −15.1 m. The internal wave itself has an average wavelength of 120 m. After detrending, the peak-to-peak excursion was Fig. 2 Example of lidar return as a function of distance, z, from the sea surface and distance, d , along the flight track. The color scales for the atmospheric return (β atm , z > 0) and oceanic return (β sea , z < 0) are on the right. The smallest bin in the oceanic color scale (β sea < 3 × 10 −8 m −1 sr −1 ) has been set to white, but this is not evident in the color bar.
12.9 m, and the root-mean-square (RMS) variation was 2.3 m. The strong asymmetry of the displacement is characteristic of a strongly nonlinear wave. The cross-polarized surface return (at z ≤ 0) can be seen to align with the negative displacement of the internal wave. The atmospheric return in the figure shows the co-polarized channel. The strong surface return contaminates the atmospheric return right near the surface, so the bottom 3 m of the atmospheric return have been blanked out.
The relationship between the internal wave displacement and the surface return is clear from Fig. 2. The correlation was R ¼ 0.17 (P < 0.001) between the depth and the co-polarized surface return and R ¼ 0.13 (P < 0.01) between the depth and the cross-polarized return. However, if the surface return is delayed relative to the displacement, the correlation can be larger. For the copolarized surface return, the largest correlation was R ¼ 0.32 (P < 10 −11 ) at a shift of nine shots, which corresponds to about 12 m or about 10% of the wavelength of the observed internal wave. For the cross-polarized surface return, the largest correlation was R ¼ 0.38 (P < 10 −16 ) at a shift of 11 shots (about 15 m).
The relationships with the atmospheric return are not as clear. There is no significant correlation between the atmospheric return at 5 m and the internal-wave displacement or either of the surface returns without shifting any of the data. Delaying the atmospheric return relative to the surface return by 79 shots (105 m) produces a correlation of R ¼ 0.12 (P ¼ 0.01) for both the co-and cross-polarized surface returns. This is only marginally significant. Combining a 12 to 15 m delay from the internal-wave displacement to the surface return with a 105 m delay from the surface return to the atmospheric return suggests that there should be a significant correlation between the internal-wave displacement and the atmospheric return if the latter were delayed by 117 to 120 m, but the correlation values in the region were not significant.
The low correlations between the ocean parameters and the atmosphere are not surprising because of the turbulent mixing of aerosols in the atmosphere. However, we might still expect to see an increase in the spatial power spectral density, S, in the atmospheric return that corresponds with the peak in the spatial power spectral density in the internal-wave displacement. Figure 3 does show such a peak, although the correspondence between the spectra of displacement and surface returns is greater. The three spectral bins between 11.1 and 13.3 km −1 contain 28% of the variance of the internal-wave displacement. While 23% of the variance of both surface returns is within this same spectral range, only 2% of the atmospheric return is within this range. However, only 6% of the atmospheric spectral peaks are larger than the one that corresponds to the internalwave peak. Fig. 3 Spatial power spectral density, S, normalized by its maximum value, as a function of spatial wavelength, K , for the internal-wave displacement (black line with circles), co-polarized surface return (red line with squares), cross-polarized surface return (blue line with diamonds), and copolarized atmospheric return (green line with triangles).

Statistics
The various correlations for July 30, 2020, above have been summarized in Table 1. This table includes the correlations between the peak of the depth spectra of the individual segments and the values of the spectra of the co-polarized surface return, the cross-polarized surface return, and the co-polarized atmospheric return at a height of 5 m. Values of all parameters were taken at the same spatial frequency, which was defined by the peak of the depth spectrum. These correlations are all highly statistically significant, with P < 10 −10 except for the correlation between the depth parameter and the atmospheric parameter, which was P ¼ 10 −4 .
The correlation between depth and atmospheric return was calculated for different heights in the atmosphere (Fig. 4). Each data series was separated into ten segments, and the standard deviation of the mean value of these provided the error bars in the figure. The correlation coefficients are significant (P < 0.005) to a height of 90 m, but the value at 100 m is not statistically significant (P ¼ 0.16). A linear regression of correlation versus height shows a significant decrease (R 2 ¼ 0.79, P ¼ 4.7 × 10 −5 ). The decrease in correlation with height is consistent with turbulent mixing of aerosols; aerosols at the internal wave frequency are spread over frequencies by turbulence as they move away from the source.
We should point out that the value at 5 m in the plot is different from the value in Table 1. This is because the correlation is not a linear operator, so the correlation of the full data set is not Table 1 Correlation matrix for the various parameters of the spatial spectra of the 200-shot data segments. These parameters are the peak value of the spectrum of the scattering layer depth (depth), the value of the co-polarized surface return at that same spatial frequency (co-surf), the value of the cross-polarized surface return at the same frequency (X -surf), and the value of the atmospheric return at the height of 5 m at the same frequency [atm (5 m)].  Fig. 4 Correlation, R, between the value of the maximum spectral bin in the power spectral density of the depth of the layer and the normalized atmospheric return in that same spectral bin at height, H. Error bars are the standard deviation of the mean of values from ten subsets of the full data set. The solid line is the linear regression.
necessarily equal to the average over 10 subsets. The value in the table is well within the error bars in the figure, however.
On the other two flight days, the relationships were weaker (Fig. 5). On June 25 (day of year 177), the correlations between the scattering layer depth variations and the surface returns at the two different polarizations were statistically significant, but that between the surface layer depth and the atmospheric return was not. On August 11 (day of year 224), none of the correlations with scattering layer depth were significant.

In Situ Data
The AXBT temperature profiles from the 3 days [ Fig. 6(a)] show that the depth of the mixed layer and the strength of the stratification have increased between the first flight and the latter flights. This is clearer in Fig. 6(b), which shows the depth and strength directly as a function of  the day of the year. The first day has a shallow thermocline, but a relatively weak stratification. The second day has a thermocline that is slightly deeper, with much stronger stratification. This condition is optimal for lidar detection, and many layers were detected on this day. The third day still has strong stratification, but the thermocline is much deeper, making lidar detection of layers more difficult.
The meteorological data ( Table 2) show that the conditions on the 3 days are slightly stable, with air temperatures warmer than the sea surface temperatures. For the NOAA station data, the air/sea temperature differences were between 0.2 deg and 0.6 deg. The aircraft data have a bigger variation, from 0 deg to 1.1 deg. Generally, the air temperatures at the aircraft were higher than those at the surface. The water temperatures measured by the radiometer on the aircraft were also higher than those measured at the NOAA station. Wind speeds differed by <0.5 m s −1 , except for the first flight.

Theoretical Considerations
The modulation depth, M, of the aerosol lidar return at a particular spatial wavelength, K, can be expressed as where S is the power spectral density of the lidar signal. From the lidar equation, 38 we can see that the modulation depth is given as where β is the volume scattering function at the lidar scattering angle of 180 deg, δβ is the RMS variation of β at the internal wave wavelength, and β 0 is the mean value of β. Equation (4) follows from Eq. (3) because the Fourier transform is a linear operator on the lidar signal and the lidar equation shows that the lidar signal is linearly dependent on β. Thus, the square root of the power spectral density at any wavelength is proportional to the RMS variation in β at that wavelength.
The lidar extinction-to-backscatter ratio can be used to relate β to the scattering coefficient, b, for sea-salt aerosols because the absorption can be neglected compared with the scattering. In this case, β ¼ b∕L, where L is the lidar ratio. While values for the lidar ratio of sea-salt aerosols have been published, 39-41 the value is not important here. L cancels out, and we have For a particular situation, the scattering coefficient is proportional to the dry aerosol mass concentration, m, where the constant of proportionality depends on the composition of the aerosols, their size distribution, and the humidity. [42][43][44] Because of this linear relationship, the constant of proportionality cancels, and we have E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 6 ; 6 5 2 A common model of sea salt aerosol production over the ocean is one with an exponential dependence on wind speed, 45 where m 1 is the value at zero wind speed, a m is a coefficient, and U 10 is the wind speed at 10 m above the surface. Assuming that δU 10 is much less than the average wind speed, we can expand the exponential and get Here, the coefficient becomes important. A number of values have been reported for this parameter. 51 These include several measurements of 0.23 to 0.28 s m −1 46,47,49 and others at 0.59 to 0.62 s m −1 . 45,49 Following Lewis and Schwartz, 51 we will use a value of 0.23.
Next, we will assume that the relevant wind speed variations are not the result of actual variations in the wind, but of variations in the internal wave-induced surface currents that change the wind speed relative to the surface. These surface current variations are the result of the orbital motion of water induced by the internal wave. They are generally less than, but of the same order of magnitude as the phase speed of the wave. 12,33,52 For this analysis, we will assume that δU 10 ¼ C p , where C p is the phase speed. Using the formula for the phase speed of linearly propagating internal waves in a two-layer system, we have 18 where g is the acceleration due to gravity, Δρ is the density difference between the lower and upper layers, H is the thickness of the upper layer, and ρ is the density of seawater. Lidar measurements for the three flight days (Fig. 7) show excellent agreement with this model for the first 2 days, but not for the third day. This difference cannot be explained by the meteorological conditions, which were similar for all three flights. Mean wind speeds were 2.8, 6.2, and 4.2 m s −1 . Mean wind directions were 223 deg, 220 deg, and 209 deg. Mean air-sea temperature differences were 0.6 deg, 0.2 deg, and 0.4 deg, so the atmosphere was near neutral to slightly stable. The quantity that did change from flight to flight was the pycnocline depth, which went from 5.7 to 10.3 to 19.4 m. With a pycnocline depth of 19.4 m, the lidar is likely not capable of measuring the full extent of the downward portion of strong internal waves, leading to an underestimation of the modulation depth on August 11, 2020.
The penetration depth of the lidar depends on several factors. To illustrate, we will consider the average values from the August 11, 2020 flight. First, the saturation level of the digitizer corresponds to a photocathode current of 1.63 μA for the photomultiplier gain that was used on this day. This gain has to be low enough that the largest surface return would be less than this value. The mean surface return was 0.14 μA, well below saturation. The mean volumetric return at the surface was less than the total surface return because of specular reflections, foam, and bubbles. It was 97.0 nA. The mean background light signal was 1.79 nA, so the range from background to saturation was 911 (59 dB). This is less than the full range of the log-amp because of the background light signal. The mean attenuation coefficient for this day was 0.054 m −1 and the mean penetration depth was 31.3 m. The signal at the penetration depth is then 97 nA expð−2 × 0.054 × 31.3Þ ¼ 3.30 nA. This is above the background light signal because we used the criterion that the signal should be five standard deviations of the noise above the background light level to obtain the penetration depth.
From this analysis, it would appear that the lidar should be capable of detecting a layer at the pycnocline depth of 19.4 m. However, highly nonlinear internal waves such as those found in the New York Bight can easily extend to twice the depth of the pycnocline. In this case, the lidar would underestimate the amplitude of the internal wave and therefore the aerosol modulation depth.
In this development, we have assumed that variations in aerosol scattering are due to variations in the relative wind speed. We know that there are many other factors, and a study of optical scattering concluded that about 20% of the variations in aerosol optical scattering were related to variations in local wind speed. 53 In light of this study, we would not expect our correlations between surface scattering and aerosol scattering at the internal wave wavelength to be higher than measured. Our results are also consistent with those of Ortiz-Suslow et al., 54 who found that momentum flux in the lower atmosphere was affected by internal waves, but the effect was small.

Conclusions
The primary conclusion is that oceanic internal waves can affect the local production of aerosols. The effect is small but statistically significant. The fact that it is small is not surprising, since there are other processes that affect the distribution of aerosols in the atmosphere. The fact that the effect decreases with height above the ocean provides evidence of these competing effects.