Articles | Volume 15, issue 23
Atmos. Meas. Tech., 15, 7171–7194, 2022

Special issue: Cloud, Aerosol and Monsoon Processes Philippines Experiment...

Atmos. Meas. Tech., 15, 7171–7194, 2022
Research article
 | Highlight paper
13 Dec 2022
Research article  | Highlight paper | 13 Dec 2022

In situ particle sampling relationships to surface and turbulent fluxes using large eddy simulations with Lagrangian particles

In situ particle sampling relationships to surface and turbulent fluxes using large eddy simulations with Lagrangian particles
Hyungwon John Park1, Jeffrey S. Reid2, Livia S. Freire3, Christopher Jackson4, and David H. Richter5 Hyungwon John Park et al.
  • 1National Research Council Postdoctoral Fellow, US Naval Research Laboratory, Monterey, CA, USA
  • 2Marine Meteorology Division, US Naval Research Laboratory, Monterey, CA, USA
  • 3Instituto de Ciências Matemáticas e de Computação, University of São Paulo, São Carlos, Brazil
  • 4Consultant, Global Science and Technology Inc., College Park, MD, USA
  • 5Department of Civil and Environmental Engineering & Earth Sciences, University of Notre Dame, Notre Dame, IN, USA

Correspondence: Hyungwon John Park (


Source functions for mechanically driven coarse-mode sea spray and dust aerosol particles span orders of magnitude owing to a combination of physical sensitivity in the system and large measurement uncertainty. Outside special idealized settings (such as wind tunnels), aerosol particle fluxes are largely inferred from a host of methods, including local eddy correlation, gradient methods, and dry deposition methods. In all of these methods, it is difficult to relate point measurements from towers, ships, or aircraft to a general representative flux of aerosol particles. This difficulty is from the particles' inhomogeneous distribution due to multiple spatiotemporal scales of an evolving marine environment. We hypothesize that the current representation of a point in situ measurement of sea spray or dust particles is a likely contributor to the unrealistic range of flux and concentration outcomes in the literature. This paper aims to help the interpretation of field data: we conduct a series of high-resolution, cloud-free large eddy simulations (LESs) with Lagrangian particles to better understand the temporal evolution and volumetric variability of coarse- to giant-mode marine aerosol particles and their relationship to turbulent transport. The study begins by describing the Lagrangian LES model framework and simulates flux measurements that were made using numerical analogs to field practices such as the eddy covariance method. Using these methods, turbulent flux sampling is quantified based on key features such as coherent structures within the marine atmospheric boundary layer (MABL) and aerosol particle size. We show that for an unstable atmospheric stability, the MABL exhibits large coherent eddy structures, and as a consequence, the flux measurement outcome becomes strongly tied to spatial length scales and relative sampling of crosswise and streamwise sampling. For example, through the use of ogive curves, a given sampling duration of a fixed numerical sampling instrument is found to capture 80 % of the aerosol flux given a sampling rate of zf/w 0.2, whereas a spanwise moving instrument results in a 95 % capture. These coherent structures and other canonical features contribute to the lack of convergence to the true aerosol vertical flux at any height. As expected, sampling all of the flow features results in a statistically robust flux signal. Analysis of a neutral boundary layer configuration results in a lower predictive range due to weak or no vertical roll structures compared to the unstable boundary layer setting. Finally, we take the results of each approach and compare their surface flux variability: a baseline metric used in regional and global aerosol models.

1 Introduction

Aerosol particles in the atmosphere contribute to environmental processes related to cloud physics, radiative forcing, and geochemical cycles. Coarse- and giant-mode sea salt and sea spray particles affect the mass, momentum, and energy exchange between the atmosphere and ocean (Lewis and Schwartz2004; Veron2015) based on their diameter, total mass, and composition. Sea spray can be equally as important as fine-mode aerosol particles for indirect effects (Jensen and Lee2008), can perturb radiation budgets in the visible and infrared (Hulst and van de Hulst1981), and can even impact latent heat transport within the marine boundary layer (Fairall et al.1994; Andreas et al.2015). These aerosol particles may also contribute to air pollution through new particle formation (Lee et al.2019; Wu et al.2021) and inhibit lightning frequency (Pan et al.2022). Despite their importance and extensive study, sea salt and spray physics remain poorly quantified, as a large range of uncertainty in their spatiotemporal distribution continues to be observed in models and field observations (Bian et al.2019; Watson-Parris et al.2019); similar large uncertainty extends to any mechanically generated coarse- or giant-mode aerosol species, such as dust. These large uncertainties are due to the strong physical sensitivities of mechanically generated aerosol particles to their surface and surface layer turbulent properties as well as large length scale and timescale histories of atmospheric turbulence that local measurements do not necessarily capture.

The beginning of the primary marine aerosol (PMA) life cycle process is particle generation by the entrainment of air into water from breaking waves, the formation of bubble rafts, and subsequent bubble bursting with the formation of film and jet drops (Monahan et al.1986; de Leeuw et al.2011; Lewis and Schwartz2004; Veron2015). Larger droplets, known as spume, can also be torn from breaking waves. Known as the sea spray generation function (SSGF), the rate of aerosol loading into the atmosphere has persisted with high uncertainty, especially for large droplets and aerosols. Thus, the SSGF is an active area of research, as breaking waves (Toba and Chaen1973; Wu1992) and surface–bubble dynamics (Keene et al.2017; Deike et al.2018) are examples of the large number of air–sea processes that challenge an accurate SSGF quantification (Sutherland and Melville2015). Indeed, large differences persist between measurement datasets (e.g., Andreas1998; Lewis and Schwartz2004; Reid et al.2006; Grythe et al.2014). Fundamental measurements have also been questioned with the use of different instrumentation (stationary flux towers, airborne measurements, lidars, ships, differing in situ probe and/or measurement types), which have led to persistently uncertain concentrations and/or size-dependent profiles within field observational studies (Blanchard et al.1984; Porter and Clarke1997; Reid et al.2006). Even the methods of flux calculations themselves are highly varied, ranging from local correlations and vertical profiles to box model approaches and eddy correlation – each with their advantages and disadvantages (see, e.g., the review by Lewis and Schwartz2004). When these flux deductions are then applied in models, the models might then propagate incorrect forecasting and prediction of their spatiotemporal distribution. Thus, a biased SSGF (or any other aerosol source function) may improve or diminish model skill independent of the flux parameter's true skill (see, e.g., the discussion in Sessions et al.2015). There are numerous factors that are hypothesized and still debated as to SSGF dependencies that may account for some of the differences in the literature and with ongoing field research, including wave state, sea surface temperature, atmospheric stability, and biological activity including diurnal effects (e.g., Irwin and Binkowski1981; Kapustin et al.2012; Grythe et al.2014; Keene et al.2017; Srivastava and Sharan2017).

As a case in point, the NASA Cloud, Aerosol, Monsoon Processes Philippines Experiment (CAMP2Ex) P3 aircraft and Office of Naval Research Propagation of Intraseasonal Oscillations (PISTON) R/V Sally Ride tried to quantify both turbulence and aerosol properties over the Sulu Sea and subtropical western Pacific (WESTPAC) within the Maritime Continent's boreal summertime southwest monsoon. Like many other field-based studies (Lauros et al.2011; Li et al.2017), these platforms have obtained the vertical distribution of aerosol particles in the marine atmospheric boundary layer (MABL) (Blanchard et al.1984; de Leeuw1986; Gong et al.1997; Reid et al.2001). Both the P3 and the R/V Sally Ride made point or curtain measurements of atmospheric aerosol and turbulence properties from which aspects of the aerosol life cycle are hoped to be inferred, including source functions and particle lifetimes. However, such measurements represent a small set of parcels that have both local and long-range aerosol influences: from individual gust fronts to as far away as the Indian Ocean (Reid et al.2013, 2015). Sea salt data in particular are difficult to interpret due to both measurement shortcomings and the nature of MABL flows (e.g., Reid et al.2006). Therefore, when the P3 or Sally Ride makes a measurement of sea salt it is often difficult to know what that measurement actually represents.

CAMP2Ex found that the monsoonal MABL is quite heterogeneous, with features such as convergence lines, land breezes, cold pools, and detainment features within a large littoral environment. The observations led this study to hypothesize that MABL heterogeneity is a further complication of SSGF quantification that may also help explain the differences between flux measurements that are likewise overwhelmingly littoral. That is, high volumetric variability and strong gradients within the MABL, ranging from changes in wind along a sampled trajectory over hundreds to thousands of kilometers to scales of the order of 1 km, can lead to further diversity between individual flux measurements. Flux measurement methods assume some form of atmospheric steady state and spatial homogeneity. For example, eddy correlation or vertical profile methods measure the local net flux, which is often assumed to be the total production flux; however, when there are high winds upstream of the measurement point, that may result in a net downward flux even if active but less production is being taken place. Box model methods, whereby upstream and downstream states are measured, likewise measure a net flux, but if regions and conditions are chosen properly, uncertainty in the downward and dry deposition flux can be sequestered (e.g., Reid et al.2001; Kapustin et al.2012). Even so, in these, knowledge of the air mass history and an invariant atmospheric structure are assumed. In all cases one must be cognizant of uncertainties associated with high-resolution atmospheric features and their association with the sea spray aerosol distribution, including the state of the surface layer, the presence of MABL roll features and diurnal (i.e., nonstationary) turbulence (e.g., LeMone1973; Wurman and Winslow1998; Kapustin et al.2012; Monaldo et al.2014; de Szoeke et al.2021; Prajapati et al.2021), or the prevalence of cold pools (Reid et al.2015, 2016; de Szoeke et al.2017).

Examples of a simple and common MABL feature observed as surface winds speeds increase to whitecapping levels are MABL roll structures. During the CAMP2Ex field campaign, the Sentinel-1 synthetic aperture radar (SAR) obtained 20 m resolution imagery over the Sulu Sea on 1 September 2019, as shown in Fig. 1. The Sentinel-1 normalized cross section values were converted to 250 m sampled ocean surface wind speeds using the CMOD4 geophysical model function and wind directions obtained from the Global Forecast System (Monaldo et al.2013). On this day were typical southwest monsoon conditions with surface wind features ranging from  5 to 9 m s−1 over an approximately 30 km× 30 km area (see black box in Fig. 1 for scale reference). As mean winds reach 7 m s−1 within a small sector, the surface wind speed minimum and maximum organize along lines  2 km apart, parallel to the direction of the wind as provided by the ERA5 reanalysis (model summary in Hersbach et al.2020). The configuration represents the formation of coherent roll structures and vortices in the MABL. In comparison, the ERA5 reanalysis wind speeds (not shown) provide a smooth wind analysis ranging from 5 m s−1 in the lower left corner to 6 m s−1 in the upper right, below the threshold for whitecaps (Monahan et al.1986; Wu1992). This SAR wind imagery demonstrates some of the challenges facing field campaigns: the MABL is not uniform and exhibits horizontal and vertical structure of ranging wind speeds. Instantaneous snapshots by ship or aircraft lead to questions regarding the local representativeness of in situ concentration and fluxes made by ships or aircraft. How does a regional SSGF relate to in situ ship or aircraft measurements? How much sampling (temporally or spatially) is required to achieve convergence between the measurement and regional flux? What are the potential measurement sensitivities to unresolved upstream sea spray production? Are there measurement practices that can aid in the measurement of turbulent aerosol fluxes?

Figure 1A 250 m resolution image taken from the Sentinel-1 synthetic aperture radar (SAR), post-processed to show surface winds at the Sulu Sea on 1 September 2019. For scale reference, the black box represents a square distance of 11 km in each cardinal direction. The black arrow represents the orientation of the wind direction as given by the ERA5 reanalysis.

To assist in hypothesis evaluation and establishing a predictive range of retrieved SSGFs in field settings, numerical models can be used to help understand volumetric evolution of the MABL from which flux measurements can in turn be simulated (Chen et al.2018; Eaton and Fessler1994; Li et al.2016; Wei et al.2018). In this study we examine the effect of turbulence aerosol particle transport through large eddy simulations (LESs) (Moeng1984; Sullivan and Patton2011) combined with Lagrangian point particles, representing a simple cloud-free marine atmospheric boundary layer with prescribed aerosol particle sources. From this numerical model, uncertainties in simulated aerosol particle flux measurements can be quantified in full 3-D spatiotemporal detail from which flux disaggregation can be achieved in an idealized manner. We study the relation between sampling methodologies, such as stationary and airborne data acquisition; this approach has been previously used to record turbulence statistics in LES-generated settings (Scipión et al.2008; Wainwright et al.2015). Here we expand upon previous work with additional features of boundary layer dynamics – turbulent coherent structures and aerosol particle size, as well as their effect on the surface flux of aerosol particles. Finally, we use one-dimensional column models (Kind1992; Freire et al.2016; Nissanka et al.2018) to estimate net surface emissions and measure their performance against the stationary and airborne data acquisition. Combined, we show the effect of factors which affect the surface flux, including directionality with respect to roll features, aerosol particle size, and height in the boundary layer, and compare how sampling methodologies determine the rate of convergence of the net aerosol flux. These insights aim to assist in both interpretation of field observations of aerosol flux and experimental design.

2 Methodology

2.1 Large eddy simulations with Lagrangian particles

This study uses the National Center for Atmospheric Research (NCAR) large eddy simulation (LES) model (Moeng1984) that is combined with Lagrangian particles (NTLP model, see Park et al.2020; Richter et al.2021). The Eulerian fields of mass, momentum, and energy are solved using the filtered Navier–Stokes equations with the Boussinesq approximation:


where ũi is the resolved velocity, θ̃ is the resolved potential temperature, τij is the subgrid stress, f is the Coriolis parameter, τθi is the subgrid turbulent flux of potential temperature, and p̃ is the resolved dynamic pressure containing both hydrostatic and normal forces, which are used to satisfy the divergence-free condition. The Eulerian subgrid-scale turbulent fluxes are parameterized with a prognostic equation that solves for the subgrid-scale turbulent kinetic energy, which is then used to define a mixing length (Deardorff1980).

We drive the flow by imposing a constant, unidirectional geostrophic wind Ug. The Eulerian field is assumed periodic in the horizontal (x and y) directions and resolved using a uniform grid in all three directions. At the start of the simulation, an inversion is imposed at the midpoint of the vertical extent (zinv 510 m), with a radiation condition at the top of the domain (Klemp and Durran1983). A pseudo-spectral discretization is used for spatial gradients in the horizontal directions, while a second-order finite-difference scheme is used in the vertical direction. Time integration is done using a third-order Runge–Kutta method, and a divergence-free filtered velocity field is achieved via a fractional step method. The lower boundary condition is prescribed by a rough-wall Monin–Obukhov similarity relation, and the surface is assumed flat with a constant aerodynamic roughness that is representative of the open ocean (z0= 0.001 m). Further details of the LES model can be found in Moeng (1984) or Sullivan et al. (1996).

In addition to the integration of the filtered Navier–Stokes equations for boundary layer turbulence, airborne aerosol particles are represented as Lagrangian point particles – a computational method which assumes that the particles are smaller than the smallest resolved scales of the flow (Balachandar and Eaton2010) and which allows for straightforward treatment of gravitational settling, inertia, and thermodynamic effects (the latter two are not considered here). As opposed to an Eulerian framework, the Lagrangian approach allows for a natural approach of sampling discrete particles as opposed to a concentration field, given that average concentration and flux profiles remain the same. The Lagrangian point particles are governed by the following equations for their position and velocity.


Here xp is the position of particle p, vp is the particle's velocity, uf is the resolved fluid velocity interpolated to the particle position using sixth-order Lagrange interpolation, and ws=τpg is the terminal settling velocity of the particle, where τp is the Stokes timescale (Wang and Maxey1993). In Eq. (4), the final two terms on the right-hand side account for particle transport done by unresolved, subgrid scales. The parameter η is an independent and identically distributed random value from a normal distribution, while K(xp) is the average subgrid momentum diffusivity obtained from the LES model interpolated using trilinear interpolating from the surrounding grid points. Overbars refer to averaging in the horizontal directions. The fourth and final term, zK(xp)Δt, takes into account the vertical transport that is caused by spatial gradients of the mean subgrid diffusivity, which is required to conserve the well-mixed condition that would otherwise be artificially violated (Delay et al.2005).

Lagrangian–LESs are computationally expensive for sufficient aerosol flux statistics, and some simplifying assumptions need to be made. The simulated Lagrangian sea spray aerosol particles maintain a constant size, meaning that condensation–evaporation or aerosol swell are not considered (Winkler1988). Likewise we are not allowing condensation or cloud formation – a topic of ongoing research. We also ignore the momentum and energy exchange between the aerosol particles and the air; that is, we neglect the two-way coupling of droplets within the surface or cloud layers (Peng and Richter2019; Mellado2017). Lastly, as mentioned above, the LES assumes a flat surface with a prescribed aerodynamic roughness length, although in the open ocean the moving surface waves may play a substantial role in the transport of sea spray aerosol particles in the wave boundary layer (Richter et al.2019). These simplifications do not impact our primary goal of this study, which is to investigate sampling strategies and statistical convergence of the aerosol source function, and not necessarily the thermodynamic evolution and subsequent turbulence of the fully coupled system. This study is also dependent on the realism of the LES, including the limited domain extent and periodic horizontal boundary conditions; however, we can represent features such as those shown in Fig. 1 and subsequently discuss the importance of fine-scale structure to interpreting field data.

2.2 Simulation configuration

The domain size and number of grid points are held fixed at 10 km× 5 km× 1 km (x×y×z) and 320 × 320 × 320, respectively. This corresponds to 30 m streamwise × 15 m spanwise horizontal resolution, and roughly 3 m vertical is similar or finer in resolution compared to other convective atmospheric boundary layer studies (Moeng and Sullivan1994; Sullivan and Patton2011). We define our inversion height zinv as the location where the planar average of the potential temperature gradient is at its maximum. The inversion height allows us to define a largest eddy timescale as Teddy=zinv/w, where w=[gT0Qzinv]1/3 is the convective velocity scale which is a function of the surface heat flux Q, the inversion height zinv, and a reference temperature T0 (Deardorff1972). The time step is set to 0.5 s with an initial temperature inversion of 0.15 K m−1 beginning at approximately 510 m. The use of such a strong inversion is to limit entrainment and boundary layer growth in order to provide nearly statistically stationary conditions. The minimal growth also maintains a relatively constant (< 5 % difference) value of Teddy 700 s. The total time of the simulations is 15 000 s, or t/Teddy 21.5.

We emphasize that our goal is not to literally reproduce real meteorological phenomena, as in Fig. 1, but to obtain a setup that displays a heterogeneous distribution of turbulence and investigate how this affects the transport and subsequent sampling of aerosol particles. This said, general environmental parameters used as boundary conditions of the model do represent the basic features of the wind speed enhancement in Fig. 1. A unidirectional geostrophic wind (Ug=10ms-1) and a surface heat flux Q of 0.02 K m s−1 are applied to force the MABL dynamics (resulting in  6 m s−1 mean surface wind), representing a typical unstable convective boundary layer corresponding to an air–sea temperature difference of roughly 1.5 K, with a sea surface temperature θs of 283 K. We note here that the simulation setup is similar to that of Moeng and Sullivan (1994). The friction velocity u is around 0.35 m s−1. After running the simulation for 1 h for turbulence spin-up, particles sampled from a monodisperse distribution are generated randomly along an xy plane at the surface (z= 0 m). Simulations are conducted for two particle sizes with diameters of 10 µm and 50 µm. The 10 µm size is chosen as it largely accounts for the sea salt coarse mode (Reid et al.2006) and has an appreciably small settling velocity relative to the mixing velocity of turbulence (ws< 0.3 cm s−1). Likewise, 50 µm is also chosen as the settling velocities become appreciable (ws 7 cm s−1) and is a common ending size for wing-mounted probes, such as the forward-scattering spectrometer probe (FSSP). Bubble-bursted film or jet droplets can be in this range (Lewis and Schwartz2004), as can other sources of “giant” aerosols where production is the most uncertain (Ryder et al.2019; Jensen and Lee2008). The source flux is set to Φss= 0.0002 m-2s-1, or 10 000 Lagrangian particles over the whole domain per time step for the 10 µm, and it is against this known source flux that the estimates are compared. Sensitivity tests are conducted for 50 000 and 100 0000 Lagrangian particles per second, and although turbulent statistics and concentration profiles remain relatively insensitive to injection rate, the distributions of particle concentration for the simulated instrumentation (explained in Sect. 3.3) are unaltered for the 50 µm sizes at Φss= 50 000 particles per second. If the Lagrangian particles fall below the surface (xp,3<0), the particle is removed from the simulation, representing dry deposition.

3 Results

3.1 Particle-laden flow characteristics

We first examine the evolution of the aerosol particle concentration in time with the goal of understanding the accuracy of surface flux estimates using common measurement strategies. The normalized vertical concentration profiles are shown in Fig. 2 for 10 and 50 µm diameter particles at multiple times tgen/Teddy, where Cr is the reference concentration, defined to be the concentration at the first grid point (z= 1.6 m); tgen represents the time since the first generation of particles; and Teddy is the largest eddy timescale, being around 12 min for the current setup. As expected, the maximum aerosol concentration is located within the surface layer (z≲0.1zinv) at all times, consistent with results observed in the field (Blanchard et al.1984; Bian et al.2019; Schlosser et al.2020). Above the surface layer, the concentration becomes more uniform, meaning that there is little average variation in height over much of the central portion of the simulation domain. Then, for times tgen/Teddy2, the concentration rapidly drops at the inversion, where detrainment is inhibited by the strong inversion and weakened turbulence. In the context of in situ measurements, for z/zinv>0.2 (or 100 m for a zinv of 500 m), C/Cr is nearly constant. Yet, 100 m is the minimum altitude at which the NASA P3 and many other research aircraft fly. Under some circumstances, the minimum altitude can be 30 to 50 m, but not under high wind conditions. Therefore, more common in situ ship, tower, or aircraft measurements are at the maximum and minimum values of the concentration profile; observations between the two require special technology such as towed instruments (Yamaguchi et al.2022) or unmanned aerial vehicle measurements (Kemppinen et al.2020; Mehta et al.2021).

Figure 2Normalized concentration profiles for (a) 10 µm and (b) 50 µm diameter aerosols at intervals of tgen/Teddy. The term tgen is the time since the initiation of surface aerosol generation; Teddy is the largest eddy timescale, a function of the inversion height and convective velocity scale: it is approximately 12 min.


Although the mass concentration profile is steadily increasing in time, the fluxes (both resolved and subgrid) are expected to approach statistical stationarity soon after particles are introduced. This rapid approach to stationarity of the fluxes and other turbulence statistics, despite the continued growth of aerosol concentration, has been shown in other LES studies (Freire et al.2016; Nissanka et al.2018). For simplicity, we remove the effects of initial transients by performing all analysis after tgen/Teddy=2. To obtain a baseline quantity to compare against all future turbulent flux measurements, we calculate the net aerosol flux Φ(z) by recording the net number of aerosol particles that have traversed through a plane z= constant at each time step, then divide by the time step and horizontal extent of the domain; these fluxes are calculated at the same Eulerian grid heights (in z) where the flow information is stored. The net aerosol flux is the combination between the turbulent flux against the gravitational settling flux.

Based on the bottom-up diffusion concept for a conservative scalar in a convective planetary boundary layer (Wyngaard and Brost1984), we expect the 10 µm aerosol particles (and likewise any smaller particles) to behave largely like passive scalars and therefore exhibit a linear flux with height. This is verified in Fig. 3a, which shows the time-averaged net turbulent flux in blue. The translucent profiles in the background are instantaneous snapshots of the net turbulent flux at the start, middle, and end times of the averaging (only for 10 µm), illustrating that the flux statistics are stationary despite the concentration increasing in time. Although the 50 µm diameter particles (in red) exhibit the same endpoint values, a slight curvature in the flux profile is observed. This is due to the additional presence of a non-negligible settling flux wsC that modifies the balance for a typical conserved scalar (Nissanka et al.2018).

Figure 3(a) Vertical flux profile for 10 and 50 µm particles normalized by the respective flux at the first grid point (Φ1). The three faded profiles are the instantaneous net fluxes at the start, middle, and end of the averaging time for 10 µm particles. (b) Comparison of mean particle velocity for both sizes and the standard deviation of the flow vertical velocity at different heights (the left vertical axis is scaled by the inversion height, and the right side is the dimensional height). The dashed line is the gravitational settling velocity for 50 µm particles.


To further investigate the influence of particle size on the flux profiles, we also consider the mean aerosol particle vertical velocity; this can be used as an indicator of gravitational transport over turbulent advection on the aerosol particles. Figure 3b shows the average vertical velocity of the particles at each height (wp) and the fluid vertical velocity variance (w2). For the 10 µm particles (and smaller), the negligible settling velocity throughout the boundary layer suggests a small deposition rate (thus a long lifetime), for which negative average velocities are reported only up to 1.6 m (the first grid point). On the other hand, heavy 50 µm particles have a negative mean velocity up to 20 m from the surface. This negative mean velocity will lead to an overall deposition of particles this size. The presence of the lower surface and the strong inversion at the top of the boundary layer cause w2 (in yellow) to be at its weakest in these regions. As a result, the mean vertical velocities of both 10 and 50 µm aerosol particle sizes approach their respective gravitational settling velocities ws, since this would be the mean velocity in the absence of the surrounding flow (the vertical dashed line indicates the value of ws for the 50 µm diameter particles).

To provide a more qualitative view of the boundary layer turbulence, Fig. 4 shows instantaneous contours of turbulent vertical velocity in multiple planes cut through the domain. The combined influence of shear- and buoyancy-generated turbulence creates coherent streaks at low levels that are a part of convective roll structures at higher levels; this interaction and the resulting roll structures are well documented in both observations (LeMone1973; Grossman1982; Weckwerth et al.1997) and numerical simulations (Moeng and Sullivan1994; Khanna and Brasseur1998; Salesky and Anderson2018).

Figure 4Image of instantaneous vertical velocity contours for an xy top-down view (a, b), xz streamwise side view (c), and yz spanwise side view (d). The xy top-down views are in panel (a) z/zinv=0.083 (or 50 m) and in panel (b) z/zinv=0.525 (or 300 m). x is the streamwise direction of the flow. The solid black lines across (c, d) are the heights of panel (a) and panel (b) plane visualizations.


Considering Fig. 4a, an xy plane near the upper portion of the surface layer at 50 m (z/zinv=0.083) exhibits elongated streaks of positive vertical velocities in the streamwise direction. As gravitational settling becomes stronger (e.g., 50 µm), particles require persistent upward velocities to counteract their settling tendency; hence, these streaks are critical for the vertical flux of heavy particles. In Fig. 4b, the assumed “well-mixed” region at z= 300 m (z/zinv=0.525) exhibits streaks that have grown in spanwise width, with updrafts and downdrafts at their peak values; this height is associated with the highest vertical velocity variance in Fig. 3. From an observational perspective, this horizontal heterogeneity means that if sampling were restricted to only within one of these streaks, then the measured turbulent flux of aerosol would be biased as an overprediction (and vice versa). These boundary layer structures necessitate strategies for achieving full statistical representation. Figure 4c and d show the same velocity contours but for spanwise and streamwise planes (the solid black lines across the xz and yz planes are the heights of the xy plane visualizations in Fig. 4a and b). The surface layer is characterized by a short spatial periodicity of vertical velocity and vice versa in the well-mixed layer.

To further quantify the turbulent coherent structures, we compute their two-point spatial correlations of the vertical velocity in both streamwise and spanwise directions. As noted above, we perform this calculation only after waiting for statistical stationarity, which develops beyond Teddy 2. These are plotted in Fig. 5 at the same two heights in Fig. 4. Figure 5a is the two-point spatial correlation in the streamwise direction. Near the center of the convective boundary layer (red), the streamwise vertical velocity correlation is strong at multiple inversion heights, while it is significantly weaker in the upper regions of the surface layer (blue). Figure 4b visualizes this effect, whereby instantaneous turbulent structures at 300 m exhibit prominent streamwise coherence compared to 50 m.

Figure 5Two-point correlation of vertical velocity variance Rww in the (a) streamwise and (b) spanwise direction for z/zinv=0.088 (red) and z/zinv=0.525 (blue). Δ is the spatial distance dependent on the coordinate direction, nondimensionalized by the inversion height, zinv.


Figure 5b shows the two-point correlation in the spanwise direction. Near the top of the surface layer (blue), the correlation is again quickly reduced to zero, while the mixed layer exhibits oscillations which are characteristic of alternating streaks of high and low vertical velocity; this can be seen in Fig. 4d where the vertical velocities alternate with distance y. Later, we will show that sampling along these structures causes substantial underprediction or overprediction when estimating a representative aerosol flux.

3.2 Idealized direct flux measurements

Now that the overall LES has been described, we introduce idealized sampling studies to establish a baseline uncertainty of sampled concentration fluxes. To start, we calculate idealized, “true” direct flux measurements from the simulation data as a function of sampling area. Once these direct flux measurements are obtained in height, we will use them to estimate an emission flux and compare it with the known surface flux prescribed in the simulations. That is, to provide an estimate of how representative a limited domain is for a given time span, we systematically decreased the sample domain and compared it to the entire LES to determine the spatial sampling to convergence. While these statistics are not representative of any measurement technology, they are illustrative of the spatial and temporal scales that must be considered when interpreting flux data. Later, these measurements will be used in comparison to the simulated instrumentation and indirect estimate results, tying the analysis to real-world applications and providing a theoretical best-case scenario for computing and evaluating concentration fluxes.

We begin by exploring the effect of limiting the spatial extent of the sampling regions. This is done by directly calculating the net concentration flux in multiple subregions of size 3.1 km2, or 1/16 of the horizontal domain. This 1/16 unit area is represented by dark blue in Fig. 6a. We denote this as Af, where f is the number of these 3.1 km2 equal-area rectangles used to create a larger subregion of the total area. In Fig. 6a, the red area, for example, would be denoted as A4 and uses 1/4 of the total horizontal area to subsample, while the green area would be denoted as A9. A16 would represent the total area of the domain. While there is no measurement analog for area sampling, this sort of flux disaggregation approach has proved useful in previous numerical studies (Hutjes et al.2010; Sühring et al.2019) to conceptualize and compare the subsampled net fluxes ΦAf(z) with the domain-integrated value Φ(z).

Figure 6Illustration of the wall-normal plane decomposed into 16 equal-sized rectangles. Panel (a) shows how the sampling region is varied by size; the dark blue rectangle is an example of A1, the red rectangle is for A4, and green shows A9. Panel (b) illustrates the varying directional coverage with a fixed area A4. The minimum rectangular unit covers 3.1 km2, or 1/16 of the horizontal domain area.


We employ these subregions in two ways: first to examine how quickly the subsampled flux ΦA1(z) approaches Φ(z) as a rectangular subsampled area approaches the full domain and second as the same subsampled area A4 is configured into various directional patterns. The directional-based sampling of the net flux will be susceptible to variability based on the area's alignment with coherent boundary layer structures. With periodic boundary conditions in both horizontal directions and allowing partial overlap of the subregions, there are 16 possible configurations for the rectangular regions of Fig. 6a; the spanwise and streamwise sampling has only four possible configurations, and the diagonal sampling has eight configurations.

We first calculate the flux ΦA1(z), which occurs over each time step of Δt= 0.5 s, after statistical stationarity has been achieved for all 16 regions. These are presented in Fig. 7 as a probability distribution function. Each panel corresponds to a different height, and the sampled particle flux is normalized by the true, domain-wide Φ(z) for the same aerosol particle size. The blue bars represent the probability of 10 µm, and the red represents the 50 µm aerosol particles. In all figures, the expected value is shown by the vertical dashed–dotted line at ΦA1/Φ(z)= 1, meaning that the average of the sampled fractional areas converges to Φ(z).

Figure 7Probability density function of 10 µm (red) and 50 µm (blue) aerosol particle fluxes at z=10m (a), z=50m (b), z=300m (c), and z=500m (d) for A1. The pink probabilities signify the overlap of fluxes between the two sizes.


Due to uniform particle sourcing at the surface, Fig. 7a shows that ΦA1 exhibits the narrowest range of flux for both particle sizes near the bottom of the domain. The heavy 50 µm aerosol particles exhibit a wider variation based on A1 due to their gravitational settling and their stronger reliance on higher fluid velocities for upward transport. The small 10 µm flux distribution is nearly all a net positive at z= 10 m, but then exhibits negative turbulent fluxes in Fig. 7b, as the turbulent flux broadens for both aerosol particle sizes. The subsequent panels show that aerosol particles that reach higher altitudes will have a larger net flux range for ΦA1 all the way to the inversion height (note the x-scale change as one considers higher altitudes).

Given this initial analysis, we now quantify in Fig. 8 the range of fluxes recorded at each height for the various subsampled areas (A1, A4, and A9). We take a sample time relevant to field observations (Geever et al.2005; Norris et al.2008, 2012): 12 min, or ∼1Teddy.

Figure 8Vertical profiles of flux ranges in sampled (a, b) 10 µm and (c, d) 50 µm aerosol particle fluxes for a time-averaged duration of 1Teddy. The directional sampling (b, d) is for a sampling area of A4. Colors from Fig. 6 correspond to those in the current figure. The dashed black line is the zero value reference of net flux, while the solid black line is the reference value for the sampled aerosol particle flux equal to the true, domain-wide LES aerosol particle flux.


For both 10 and 50 µm particles, the ranges between the minimum and maximum net fluxes are shown in Fig. 8. Figure 8a and c are the rectangular configurations, while Fig. 8b and d are for directional-based configurations. The dashed black line in all panels represents the zero flux value (which distinguishes an overall upward or downward net flux), while the blue vertical line represents the mean net flux. The multiple colors in each panel correspond to the subregions of sampling previously described in Fig. 6.

For any rectangular configuration (regardless of size) in Fig. 8a and c, the range in the subsampled net flux increases with height, centered about the mean value of 1. This result can be attributed to the increasing fluid length scales of the turbulent coherent structures as one approaches the center of the mixed layer, as quantitatively shown in Fig. 5 and qualitatively in Fig. 4. The coloring past the zero value line (from a reference point of ΦAf(z)/Φ(z)=1) shows the heights at which a sampled aerosol flux can erroneously report a net downward flux. As expected, sampling from a larger subregion (green being the largest) results in a reduction of the variability in the measured value ΦAf(z).

In Fig. 8b, we quantify the effect of flux sampling based on the directional alignment described in Fig. 6. Again, sampled net fluxes become less representative of Φ(z) with height; either underprediction or overprediction can be sampled from the spread. The streamwise sampling has the largest net flux variability and is attributed to the persistent turbulent structures of high vertical velocity seen in Fig. 4. The diagonal subsampling and streamwise subsampling exhibit similar levels in aerosol net flux variability. Finally, the spanwise sampling (green) in the crosswind direction provides the least variability for the 1Teddy sampling time considered. These variations in net fluxes, although considering the same subregion A4, are due to the coherent turbulent structures present in the boundary layer.

Shifting the focus to large 50 µm aerosol particles using the same methodology, Fig. 8c exhibits many of the same features as the 10 µm particles. However, the significant gravitational settling provides a few key differences in the results. The predictive range for both subsampling strategies is nearly doubled compared to the 10 µm counterparts: for example, there is roughly a 100 % increase in predictive range given ΦA1 at z/zinv=0.80 and a 300 % increase given the streamwise sampling at z/zinv=0.60. Since Φ(z) decreases to zero at the top of the boundary layer, this means fewer particles reach that height, resulting in large underpredictions or overpredictions of turbulent flux: this effect is shown in all panels. A net downward flux is more probable when considering the larger range of prediction to the left of the zero reference line. When considering the directional sampling combinations for the 50 µm sizes in Fig. 8d, the same conclusions are made with respect to the 10 µm aerosol particles, namely that the spanwise subsampling provides the fastest convergence towards the true flux, although the flux range increases substantially at the top of the boundary layer.

The disaggregation technique employed here (similar to Sühring et al.2019, and Hutjes et al.2010) demonstrates the importance of areal coverage and directional sampling when calculating aerosol mass flux; even in an idealized numerical setting, there is a large range of possible flux values when only sampling a subregion of the turbulent boundary layer. In all cases, this highlights the fact that even if a perfect measurement of flux were to be made at various limited-area heights in the MABL, there is a high likelihood of this measurement not capturing the true surface flux below – in many cases it may even get the sign incorrect. Due to the coherent turbulent structures, the directional sampling causes substantial differences in the net flux, and perhaps intuitively, sampling the domain across the spanwise direction exhibits the fastest convergence towards the true horizontal mean.

While the analysis above discusses the accuracy of sampled fluxes at specific heights under various subsampling strategies, it is often surface aerosol fluxes which are of interest, however, and not necessarily those at some elevated point within the boundary layer. We now use the results above to investigate how this sampling range propagates into inferred surface fluxes.

In the limit of negligible gravitational settling, the aerosol concentration behaves as a conserved scalar and its flux exhibits a linear shape with height; gravity causes this to deviate (see Fig. 3). Regardless of size, however, we fit a linear profile to each profile Φ(z), then extrapolate to the surface to retrieve an inferred surface flux Φs,z:


where we assume the zero flux at the inversion, Φ(zinv)=0. We perform this extrapolation from each height by considering the largest overprediction or underprediction of fluxes, similar to the analysis of Fig. 8. Figure 9 shows the inferred effective surface flux Φs,z, which is normalized by the time-averaged true surface flux Φs at the first grid point. Figure 9a and b represent the 10 µm aerosol particles, and Fig. 9c and d represent the 50 µm aerosol particles. For 10 µm aerosol particles, much of the profile is very similar to Fig. 8a and b, which is to be expected owing to the linearity of the flux Φ(z) of these particles in Fig. 3a. For heavy aerosol particles, however, not only do they have a larger range of ΦAf based on subsampling, but the downward bend of the overall turbulent flux (Fig. 3a) also violates the assumption of a linearly decreasing flux profile. This violation causes the curvature in the surface flux profile, with maximum uncertainty at the upper boundary layer heights at z/zinv 0.7. Thus, for 50 µm aerosol particles, Φs,z shows an even wider range of possible values due to the compounding effect of significant settling, leading to a departure of an assumed linear flux profile with height. Note that it is possible to incorporate gravitational settling into the particle flux balance (Nissanka et al.2018), and this is considered below in Sect. 3.4.

Figure 9Inferred effective surface fluxes Φs,z normalized by the true effective surface flux Φs as a function of sampling height for 10 µm (a, b) and 50 µm (c, d) aerosol particles. Elevated fluxes are used to infer surface fluxes via linear extrapolation. The left panels (a, c) illustrate variability due to area coverage, and the right panels (b, d) illustrate directional sampling for the same total area A4.


Overall, the inferred surface fluxes demonstrate the same level of predictive range as the subsampled net fluxes of Fig. 8. The resulting flux range is amplified when extrapolating to obtain a surface flux – a quantity that is required in regional and global aerosol models. In the next section, we use a similar framework using constraints analogous to field observations.

3.3 Simulated instrumentation

Thus far, we have conducted geometrical subsampling to quantitatively determine the range of estimated net fluxes (both local in height ΦAf(z) and extrapolated to the surface, Φs,z) as it is influenced by sampling area. In the interest of better interpreting the observational constraints of field observations, we now turn our attention to simulated instrumentation. In particular, the previous section analyzed the predictive range when the true flux, as computed by the mass of particles crossing a horizontal plane, was calculated exactly. This kind of information, however, is unavailable in practice and only represents a best-case scenario of observations.

Therefore, we now introduce a virtual “probe” to obtain a time series of concentration C(t) and vertical velocity w(t), from which we compute an eddy covariance wc within the local LES flow field – a technique more representative of field observations of aerosol fluxes (Norris et al.2012). Sampling LES-generated fields through simulated instrumentation has been previously performed to evaluate in measurements and identify its potential biases and limitations (Wainwright et al.2015; Sühring et al.2019). The results of the eddy covariance for predicting the surface flux Φs,z are shown in Fig. 14 in comparison to the other turbulent flux methods.

As we have seen in the previous subsection, sampling across subregions of the flow results in a potentially large spread in estimates of ΦAf(z) and subsequently Φs,z due to the coherent turbulent structures and their spatial and temporal correlations. Therefore, it stands to reason that sampling in a manner that captures a representative distribution of the turbulence will result in less uncertainty of aerosol particle flux. To test this hypothesis, an ogive curve can be used to quantify the coherent turbulent structures' effect on the aerosol particle flux. A common technique in field measurements (Desjardins et al.1989; Friehe1991; Norris et al.2012), an ogive curve Ogwc(f0) represents a running integral from the highest frequencies to a reference frequency of a co-spectrum:

(8) Og wc ( f 0 ) = f 0 Co wc ( f ) d f ,

where Cowc is the co-spectrum of the net aerosol flux wc. In other words, ogive curves show the convergence rate of a net flux given a defined sampling time. If the reference frequency is set to the lowest possible value based on sampling time, the cumulative integral of the entire spectrum results in the covariance (i.e., the average net flux over that time interval).

To compute an ogive curve, we obtain C(t) and w(t) using two different sampling strategies. First, the virtual probe can be moving (mimicking an aircraft or moving ship) or stationary (say a tower or a stationary ship). Further, each probe has its own sampling volumes and constraints – far too many to simulate here and in fact unnecessary for the point that is being made in this study. For now, consider a probe that has perfect fidelity over a given sampling volume. For each time step (0.5 s), moving probes are translated at a user-defined speed of 50 m s−1 (typical air speeds of a Twin Otter research aircraft) and record both vertical velocity and aerosol concentration over a spherical volume of approximately 1500 m3, or a fraction of 3 × 10−8 of the domain volume. At the same time, 10 stationary probes are uniformly spaced along the cross-flow direction within the domain, with rows being positioned at z= 50 m and z= 300 m (so 20 stationary probes total). These stationary probes have 9 times the sampling volume compared to the moving probes for statistical sampling convergence. A schematic of the moving and stationary probes in the field of turbulence is presented in Fig. 10b. We perform analysis for sampling durations of 1 Teddy as done above. A single moving probe is recorded throughout the whole simulation, resulting in six segments of length 1 Teddy, emulating a total distance of 35 km of aerosol particle sampling.

Figure 10Calculated ogive curves of the net aerosol flux for simulated instrumentation. Moving instrumentation is shown in dashed and dash-dotted lines, while stationary instrumentation is shown in dotted lines at (a) z/zinv=0.083 and (b) z/zinv=0.525. Within the moving probe, streamwise (dash-dotted) and spanwise (dash–dash) are presented. Colors are subsequently used to distinguish the aerosol particle size difference. The sampling duration is 1Teddy. A small schematic figure containing detail of the simulated instrumentation is provided: the left dots represent the spanwise-placed stationary probes, and the upper right dot represents the aircraft sampling in the specified direction.


The net flux convergence through the ogive curves is shown in Fig. 10 for all virtual probes. Consider reading any line in the top panel beginning from the right end of the abscissa to its left end. For any given frequency, the ogive curves accumulate the co-spectrum of the net flux beginning from the highest frequency until the lowest frequency point represents the covariance (normalized to 1). The top panel is from just above the surface layer (z/zinv=0.088 or z=50 m), and the bottom panel is from within the mixed layer (z/zinv=0.525 or z=300 m). Legends identify the particle size, moving or stationary, and the direction in which the moving sampling is retrieved. Averages are taken over all 1 Teddy samples for each probe type and height. Each ogive curve is normalized by its calculated respective covariance value at 2 Hz (Ogwc,f=2Hz), while at the highest frequencies the value of Ogwc approaches zero by definition.

For Fig. 10a, there are negligible differences between the net fluxes of different aerosol sizes, as shown by the overlap in each respective virtual instrument. When comparing different instrumentation, however, the convergence rates of the net flux are not the same. For stationary probes, an 80 % net flux convergence requires a lower normalized sampling rate of zf/w0.2, where higher frequencies (zf/w>20) provide a negligible contribution to the total net flux. Since stationary probes (like buoys and flux towers) by definition sample in the streamwise flow direction, they are constrained to the periodicity of the coherent turbulent structures that drift in the spanwise direction over time. Thus, stationary probes require longer sampling times for the traversal of the full range of turbulent structures. Even so, the nonzero slope of the ogive curve at the lowest sampling rate suggests that stationary probes require additional lower frequencies to adequately capture all the scales of the net flux (Norris et al.2008) (i.e., a sampling duration of 1 Teddy is not sufficient for a stationary probe). That is, for a tower or stationary ship (often pointed into the wind), turbulent structures align with the flow and are preferentially sampled, thus requiring much longer to achieve convergence with area-average fluxes.

For completeness, we inspect the convergence rate of the net flux with respect to the moving speed by running additional simulations. By varying the moving speed of the probes 100 and 200 m s−1, we consider only the spanwise direction and 10 µm particles at z=300 m for 1 Teddy. The computed ogive curves (not shown) demonstrate a change in the convergence rate of net flux to shift toward higher normalized frequencies. In other words, probes that traverse faster throughout the boundary layer (given a fixed sampling frequency) require less evaluation time to adequately capture the total net flux of aerosol particles.

As we would expect in the presence of turbulent structures, the ogive behavior is different for streamwise and spanwise moving probes; the leveling-off at low frequencies is indicative of a fully resolved net flux. The same sampling rate of zf/w20 results in a 95 % net flux convergence for spanwise sampling probes. Sampling the coherent turbulent structures across the spanwise direction reduces the required total sampling time, since it captures the overall statistical distribution of the net flux as seen in the previous section. Streamwise sampling achieves around a 60 % net flux convergence for zf/w20. Finally, the nonzero contribution of net flux begins at different sampling frequencies for all simulated instrumentation, with the spanwise sampling beginning at the highest frequencies (zf/w60).

Figure 10b shows the net flux convergence within the mixed layer, where there are a few differences compared to eddy covariance sampling just above the surface layer. First, there are visible, albeit small, discrepancies between particle sizes. Larger particle sizes shift the net flux convergence to lower frequencies (i.e., longer sampling times are required to capture net flux for larger sizes). However, this effect is negligible for stationary probes, for which size again plays a negligible role in the convergence of the net flux compared to the long time required to achieve adequate statistical sampling of the background flow. Second, the start of a nonzero net flux covariance is roughly at the same normalized frequency for both directions of the moving probe, regardless of size.

Overall, the frequencies determining the net flux convergence are unique for all moving and stationary measurements. The streamwise and spanwise sampling of the moving probe requires a shorter sampling time to capture all the necessary frequencies for a net flux convergence. This consideration is critical during the field observation planning of airborne and stationary instrumentation. How these results impact the estimation of Φs,z will be discussed below.

3.4 Theoretical flux profile method

Before fast-response sensors were available for eddy covariance measurements, indirect surface flux estimations were made through so-called flux–profile relationships that relate vertical profiles of mean particle concentration to its respective surface flux (Gillette et al.1972; Gillies and Berkofsky2004). These equations are based on the Reynolds-averaged conservation of mass of monodisperse particles under horizontally homogeneous conditions, which (by neglecting molecular diffusivity and applying an eddy diffusivity closure) results in

(9) - K C d C d z - w s C = Φ s - β z .

In this equation, C is the horizontal mean particle concentration, KC is the eddy diffusivity, ws=dp2ρpg/18µ is the settling velocity (from Stokes' law for a spherical particle, where μ is the dynamic viscosity of air), and β=Φs/zinv, the surface flux divided by the computed inversion height. Note that for now, in order to introduce and discuss the flux profile method in its ideal form, we use the full horizontally averaged concentration C(z) known from the LES; however, as above, this true average is not a known quantity in the field. Below, the flux profile method will be subjected to the same types of sampling uncertainties as described in previous sections.

Equation (9) represents the vertical flux budget, combining the turbulent transport (first term) and the gravitational settling (second term on the left-hand side) into a total vertical flux which is assumed to decrease linearly from the effective surface flux Φs to zero at zinv (right-hand side). While this assumption is clearly valid for small particles (see fluxes in Fig. 3a), it may lead to non-negligible error for the larger 50 µm particles.

When representing the surface layer, Eq. (9) can be further simplified by setting β=0, which corresponds to the constant-flux layer assumption (Wyngaard2010). In addition, Eq. (9) requires a closure for the eddy diffusivity parameter, which can be obtained by assuming KC to be proportional to the eddy viscosity parameter (Km) from the log law of the wall (Davidson2004), i.e., KC=κzu/Sct (where κ 0.4 is the von Kármán constant and Sct=Km/KC is the turbulent Schmidt number). Using these assumptions, Chamberlain (1967) and Kind (1992) obtained a general relationship between the effective surface flux and the mean concentration profile of particles, namely

(10) C C r = Φ s C r w s + 1 z z r - w s Sc t / κ u - Φ s C r w s ,

where Cr is the mean concentration at a reference height zr (note that Sct=1 in their original work). Freire et al. (2016) extended this result to include effects of atmospheric stability by using the Monin–Obukhov similarity theory's result KC=κzu/ϕ(ζ)Sct, obtaining

(11) C C r = Φ s C r w s + 1 z z r - w s Sc t / κ u × exp w s Sc t κ u ψ ( ζ ) - Φ s C r w s ,

where ζ=z/L is the stability parameter (L=-u3θs/(κgwθs) is the Obukhov length), ϕ(ζ) represents stability functions, and the terms in ψ(ζ)zr/Lz/L1-ϕ(x)xdx are defined as

(12) ψ ( ζ ) = 2 ln 1 + ( 1 - 16 ζ ) 1 / 2 1 + ( 1 - 16 ζ r ) 1 / 2 , if ζ < 0 , - 5 ζ + 5 ζ r , if ζ > 0 , 0 , if ζ = 0 .

Note that due to the assumptions adopted for the eddy diffusivity and β, these equations are expected to be valid in the surface layer only. In addition, Eq. (11) tends to Eq. (10) when stability goes to neutral (ζ→0). From both models, an analytical value of Φs can be obtained from two measurements of C at different heights – hence, surface fluxes can be estimated from elevated mean concentrations.

In order to extend these results to the entire ABL, Nissanka et al. (2018) proposed the use of Eq. (9) with β=Φs/zinv and

(13) K C ( z ) = κ u z ϕ ( ζ ) Sc t , if z < z b , a κ u z ϕ ( ζ ) Sc t 1 - z z inv 2 , if z z b ,

where zb is the surface layer height (taken as 0.1 zinv) and a=1/(1-zb/zinv)2, which creates a continuous transition between the two parts of Eq. (13). Although this model does not provide an analytic solution for Φs, it can be used to estimate Φs by numerically fitting the model to the data. Note that, for z<zb, the difference between the Nissanka et al. (2018) and Freire et al. (2016) models comes from the β term only.

As in previous works, a value of Sct=1.3 was adopted to account for numerical differences in the diffusivity of particles and momentum in the simulation (Freire et al.2016; Nissanka et al.2018). Figure 11a compares the vertical concentration profiles from different theoretical approaches with the simulation result. It is clear that, in the case of non-neutral conditions, the effect of atmospheric stability on the theoretical profile is not negligible, as the Kind (1992) model does not represent the simulation results well. As expected, the Freire et al. (2016) model is similar to the simulation up to z∼0.15 zinv, whereas the Nissanka et al. (2018) model accounts for the decrease in concentration at zzinv due to the assumption of a linearly decreasing flux profile.

Figure 11Vertical profiles of the (a) mean particle concentration normalized by its surface value and the (b) estimated surface flux from concentrations at zr and z (Φs,z) normalized by the imposed surface flux in the LES for 10 µm diameter particles. Theoretical profiles from the Kind (1992) model (yellow lines), Freire et al. (2016) model (red lines), and Nissanka et al. (2018) model (blue lines), with LES results in circles. Concentration data were obtained from the time and horizontal average of t/Teddy1.

By taking the reference concentration Cr and a second value at a different height C(z), it is possible to estimate the effective surface flux Φs,z from the theoretical profiles at any height z, as it would be performed in the field if mean concentration values at two different heights were available. The result is presented in Fig. 11b for each height z where a value of C(z) is available from the LES. From the Kind (1992) model, an underestimation between  20 % and 70 % is obtained depending on the height of the second concentration value, confirming that the Kind (1992) model is not appropriate for unstable conditions. On the other hand, the Freire et al. (2016) and Nissanka et al. (2018) models provide an error of less than 20 % for heights up to 0.9 zinv, being similar when compared to each other. While the advantage of using the Freire et al. (2016) model is the analytic solution for Φs, the Nissanka et al. (2018) model provides an estimate of Φs with less than 10 % error for zzinv. Note that, under neutral conditions, the Kind (1992) and Freire et al. (2016) models provide the exact same result.

We again emphasize that in this flux profile analysis, the calculation of C was taken as a horizontal average over the entire domain. As in the previous section, the convergence of this mean to the true horizontal mean (the theory of Freire et al.2016, and related flux–profile relationships assume horizontal homogeneity) will be heavily dependent on limited spatial or temporal sampling in much the same way that the fluxes were subject to the influence of the large coherent boundary layer structures. This uncertainty is discussed further in Sect. 4.1, where we use the concentration of the simulated instrumentation for the flux profile method rather than the domain horizontal average.

3.5 Neutral boundary layer

Up to this point we have analyzed the effect of turbulent coherent structures in a typical unstable boundary layer. As atmospheric conditions can widely vary, for completeness we perform additional analysis for a shear-driven, neutrally stratified boundary layer to identify the key differences and similarities. Consideration of a second atmospheric state serves as a helpful foundation for inferring the effects of different stabilities and the resulting sampled net flux. In the case of a neutral boundary layer, the wind shear is solely responsible for the mechanical generation of turbulence (Stull1988; Jacob and Anderson2017). We provide a briefer synopsis of this configuration, since much of the analysis follows the same patterns as the unstable boundary layer.

The neutral simulations use an identical grid configuration, boundary conditions (both Eulerian and Lagrangian), and geostrophic wind speed (details in Sect. 2.2), but now with Q=0. As a result we define an alternative large eddy timescale which is denoted as Tneut=zinv/u, where u is the friction velocity recorded at tgen/Tneut>2. This different large eddy timescale is found to be around 2.5 times longer than Teddy: around 33 min. Due to the lack of nonlocal, buoyancy-driven transport, the time for turbulence to reach a well-mixed, quasi-steady state takes longer than its unstable counterpart ( 1.5 tgen in Sect. 3.1). As a result, particle generation is initiated at tgen,neut= 6 Tneut (Φss= 10 000 s−1), and time averaging is done over tgen,neut+ 2 Tneut.

In the neutral boundary layer, the lack of surface heat flux (and thus buoyant convection) results in a different and less organized turbulence structure. Figure 12a shows the contours of vertical velocity at z/zinv=0.525. The range of the vertical velocities is lower than the unstable boundary layer, leading to less vertical transport of particles throughout the mixed layer and hence the longer time to reach quasi-steady state. Although there are coherent structures visible in the horizontal u and v contours (not shown) qualitatively consistent with previous neutral boundary layer simulations (Jacob and Anderson2017), the vertical velocity exhibits a less structured distribution. The vertical velocity correlation lengths in Fig. 12b and c also demonstrate the short spatial coherence of the fluid turbulence relative to the unstable case. Both the streamwise and spanwise correlations drop to zero at the normalized distance Δ/zinv 0.5 or Δ= 300 m, independent of boundary layer height. These differences in the spatial coherence lead to substantial differences in the particles fluxes shown in Fig. 12d, where the temporal and planar-averaged vertical profile of 10 µm (cyan) and 50 µm (green) aerosol particle flux is displayed. The same quantities shown in Fig. 3a are provided as dashed lines for reference. Contrary to the nearly linear profiles of Fig. 3, Φ(z) decreases more quickly from the surface; this effect is again a function of particle diameter. The lack of strongly coherent structures in the neutral boundary layer leads to a more disorganized and thus less efficient upwards transport of particles.

Figure 12Characteristics of the simulated neutral boundary layer. The top panel (a) is the snapshot vertical velocity field at z/zinv=0.525. The bottom left panels are the spatial correlation at heights z/zinv=0.083 and z/zinv=0.525 for both streamwise (b) and spanwise (c) directions. The bottom right panel (d) is the planar-averaged vertical flux profile for 10 µm and 50 µm diameter aerosol particles. “Neut” stands for the neutral boundary layer configuration, while the dashed references represent the previous unstable boundary layer configuration (abbreviated as Uns).


We perform the same analysis as Sect. 3.2, namely by taking fractional areas Af as well as directional sampling given A4 and comparing them to the true Φ(z). Given an averaging time of t/Tneut=1, Fig. 13 shows the results for 10 µm (Fig. 13a and b) and 50 µm (Fig. 13c and d) particles. Again, the left panels show the different Af for both particle sizes, and the right panels use various directional subsampling.

Figure 13Vertical profiles of predictive flux ranges for 10 µm (a, b) and 50 µm (c, d) aerosol particle sizes. Panels (a, c) are different subsampling sizes, and panels (b, d) represent directional subsampling given A4. The panel layout and configurations are identical to Fig. 8.


As before, larger coverage areas provide less variation to Φ(z). For small 10 µm particles (Fig. 13a and b), any subsampling technique applied in this study (with the exception at the top of the inversion height) results in a non-negative net flux result, given the dashed black zero reference. The large oscillation of flux range at high boundary layer levels, especially for larger sizes, is attributed to low particle count due to weak vertical mixing: low counts lead to lacking statistical convergence. In the case of the neutral boundary layer, and in sharp contrast to the unstable boundary layer, Fig. 13b confirms that directional subsampling is inconsequential as a predictor of the net flux, as all three predictions are nearly identical to each other. This is because the lack of coherent boundary layer structure no longer affects the spatial distribution of flux sampling. For the 50 µm particles in Fig. 13c and d, the sampled net flux range is small near the center of the boundary layer when comparing with Fig. 8c and d. The direction again has a minimal difference in accuracy when comparing streamwise, spanwise, or diagonal flow.

Overall, subsampling areas, height in the boundary layer, and aerosol particle size still contribute to the accuracy of representative net fluxes, although the neutral boundary layer estimates of flux experience much less variation. The most significant difference is that directional-based sampling is negligible in predicting Φ(z), as the vertical velocity (associated with its correlation lengths) does not exhibit the same spatial coherence as found in the unstable boundary layer. This difference causes higher aerosol particle concentrations toward the surface layer and lower toward the inversion height. Practically speaking, the results here imply that flux measurements for neutral conditions are representative given a local sampling region: a stark difference to that of the coherent turbulent structures exhibited by an unstable boundary layer.

4 Discussion

Our results conclude that stability, sampling direction, height in the boundary layer, and aerosol particle size are all critical factors to consider when retrieving aerosol fluxes. For unstable conditions, the consideration of convective roll structures is paramount in achieving statistical convergence of net flux measurements. In the context of field observations, impacts of these large coherent structures under commonplace unstable conditions are shown through the unique convergence rates of single-point and moving point measurements (Sect. 3.3). The theoretical profiles in Sect. 3.4 were agnostic to these convective roll structures, but only because they were applied using the true horizontal average over the domain. Finally, a neutral boundary layer state removes the impact of directional sampling, as the vertical velocity roll structures are significantly less pronounced and prove to be inconsequential to aerosol flux measurements. For neutral conditions, measurement challenges and overall biases are greatly reduced. This simple observation from model data under different parameter spaces demonstrates how changes in sampling method can result in integer factor differences in the value of any given flux data point and may thus explain some of the extreme variance in sea salt flux reports in the literature. Nevertheless, in the interest of bridging the gap between all approaches, we compile and summarize the predictive range of Φs,z associated with each of the methods.

Summary of prediction ranges

The main purpose of this study is to create a baseline summary of sampled surface flux estimates from different numerical approaches. Doing so bounds the sampled surface flux in predictive range: a critical tool in constricting the order of magnitude variation in surface flux inferences shown in the sea salt literature, thereby informing field observation practices and improving field study design. Due to the large parameter space, we continue to limit our analysis at two heights and two aerosol particle sizes along with both stability cases. We compile the estimated surface flux ranges, including uncertainties from spatial sampling, (1) direct flux measurements, (2) simulated instrumentation which computes an eddy covariance, and (3) the theoretical model.

Before summarizing the prediction ranges of surface fluxes, we briefly mention a couple technical details: one regarding the eddy covariance method and one regarding the theoretical flux profile method. From the simulated instrumentation, we obtain a time series of concentration C(t) and vertical velocity w(t) to compute the net aerosol flux wc for 1 Teddy (recall the discussion of the ogive curves in Sect. 3.3). From this we repeat the procedure described in Sect. 3.2: fit a linear profile to wc using a zero flux at the inversion height and extrapolate to the surface to retrieve an inferred surface flux Φs,z. We use the largest overprediction or underprediction of fluxes as the basis for the error bars in the following discussion. For the flux profile method, we use the same 1 Teddy segments of concentration for the two specified heights along with the full horizontal surface concentration (at the first grid point) and apply the theoretical profile of Nissanka et al. (2018) (henceforth denoted as N18) to obtain a surface flux. Again, we use the largest variation of surface flux as the error bars.

The summarized information is presented in Fig. 14, where surface flux estimates are calculated from methods discussed in the previous sections: elevated, idealized flux measurements extrapolated to the surface; eddy covariance using either stationary or moving probes extrapolated to the surface; and the theoretical flux profile method based on the simulated instrumentation aloft combined with the reference surface concentration. Each panel shows both aerosol particle sizes, separated by the slightly transparent center black bar. The first legend item is the idealized direct flux measurement subject to the smallest sampling area. The legend items with EC are based on the numerical instrumentation using eddy covariance. The next legend items are the application of the theoretical profile from N18. Finally, for completeness, the last legend item represents the idealized direct flux measurements under a neutral stability condition for t/Tneut=1. The dashed reference line represents Φs,z/Φs=1, where the estimated surface concentration fluxes match the true LES surface flux.

Figure 14Summary of the various methods for inferring effective surface fluxes from (a) z/zinv=0.088 and (b) z/zinv=0.525. The left side of the abscissa in each plot is for 10 µm diameter aerosol particles, and the right side represents 50 µm sizes. The first legend entry corresponds to the subregion sampling analysis in Sect. 3.2. The second set of legend entries is eddy covariance (EC) and streamwise (as stream), spanwise (as span), and stationary (as stat), corresponding to the simulated instrumentation. The next set of legends corresponds to the Nissanka et al. (2018) model as N18, in which the elevated concentration is constrained by the simulated instrumentation. The last legend item is the surface flux variation from a neutral stability condition for t/Tneut=1. Aside from the last legends, the data are based on an averaging time of 1Teddy. The dashed reference line represents Φs,z/Φs=1, meaning that the inferred surface flux matches the LES surface flux.

When considering Fig. 14a for the 10 µm particles, the idealized flux measurements at z= 50 m confirm the surface flux variation as we have seen in Sect. 3.2. The eddy covariance results show that the spanwise moving probe provides the closest and most reliable retrieval of Φs,z out of any sampling method. In comparison, the large uncertainty shown by the stationary probe is in agreement with the ogive curve discussion in Sect. 3.3: an insufficient sampling time results in a large surface flux range. The N18 model fed with mean concentrations retrieved from simulated instrumentation exhibits a similar uncertainty range, regardless of whether the concentration is measured from a moving or stationary probe. The neutral stratification, subregional sampling shows near-zero variation: this low variation results from sufficiently sampling the turbulent transport of aerosol particles given the relative absence of large, coherent structures spanning the entire depth of the boundary layer.

For the 50 µm particles (right side of the center black bar), all measurements increase in uncertainty, with the stationary probe resulting in the largest overprediction of surface flux, while the N18 model exhibits the greatest negative prediction (or underprediction). From inspection of the concentration time series of the stationary probe for 1 Teddy (figure not shown), the updrafts carrying aerosol particles cause large positive concentration fluctuations, thereby causing significant ranges of positive surface flux. On the other hand, the N18 results in large negative surface flux prediction (underprediction) due to the violation of the linearly decreasing flux at these large sizes (Fig. 3a). Lastly, the neutral stability of the boundary layer widens in surface flux range for 50 µm particles compared to their 10 µm counterparts, although it remains the most robust of the various estimates.

Near the center of the boundary layer (Fig. 14b), elongated turbulent coherent structures that exhibit the largest spatial correlation can lead to subsampling that is considerably worse in representing the overall surface flux for the unstable condition. Thus, the predictive range for any sampling method is greater when compared to the surface layer flux subsampling. Again, heavier 50 µm particles worsen the predictive range, as the range of estimated surface fluxes increases for all numerical sampling strategies. Given an averaging time of 1 Teddy, heavier particles and sampling with this height may lead to predictions of a negative surface flux. For the neutral boundary layer, the lack of these turbulent coherent structures causes the variation of inferred surface fluxes to be inconsequential in height (with exception near the inversion), as seen in Fig. 13c and d.

5 Conclusions

This study investigates the fundamental problem of how well a point and/or moving in situ measurement of sea spray aerosol or dust particle can estimate or represent the true surface flux in a turbulent atmospheric boundary layer. The recent data acquisition of the CAMP2Ex field campaign (like Fig. 1) exemplifies the real-world challenges of sampling with regards to the local representativeness of concentration and fluxes obtained. Through the use of large eddy simulations combined with Lagrangian point particles, we simulate a configuration of coherent turbulent structures that analogously demonstrates the heterogeneous nature of aerosol particle transport.

These numerical simulations, which provide full 3-D spatiotemporal detail, are then used to measure aerosol particle surface fluxes through three approaches. The first approach was through a direct method whereby particles are counted crossing a horizontal boundary using multiple configurations of subregions to compare their net aerosol flux with the domain-integrated value. This approach results in a net flux range that is found to be sensitive to multiple parameters: aerosol particle size, spatial sampling strategy, and height in the boundary layer. Indeed, the spatial correlation of the roll structures is dependent on height (Fig. 4), causing differences in the predictive range of net aerosol flux. Gravitational settling based on size competes against turbulent transport that would otherwise lead to tracer-like behavior of aerosol particles (Fig. 3b). These aerosol fluxes were then used to estimate a surface flux which could be compared to the known emission flux, since surface fluxes are required in regional and global aerosol models for forecasting. The surface flux variability was found to be nearly monotonic in height (Fig. 9), with a nearly 50 % increase when considering 50 µm sizes from the small tracer particles. The directional subsampling (given a sampling subregion A4) also supports this and offers the lowest predictive range in the spanwise subsampling direction.

Using simulated instrumentation similar to field practices as the second approach, moving airborne and stationary probes were used to sample the simulated LES fields in computing an eddy covariance. These were also presented in the form of ogive curves (Fig. 10), which suggest that the manner of observation (whether moving or stationary), directionality (streamwise or spanwise sampling for the moving probe), and aerosol particle size determine the important frequencies of the aerosol flux: our results suggest that, given a sampling duration of 1Teddy, stationary probes require additional lower frequencies (meaning a longer sampling time) to adequately capture all of the scales of the turbulent flux. As expected, sampling in the spanwise direction provides the least amount of time required in sampling the true flux. The eddy covariance technique was then also used to estimate a surface flux. The results suggest that the method of aerosol particle sampling results in unique surface flux variability due to coherent roll structures with the marine atmospheric boundary layer. The sampling method for surface flux variability can result in integer factor differences in the value of any given flux data point and may thus explain some of the extreme variance in sea salt flux reports in the literature, especially for large particle sizes.

The third approach uses the flux profile method for surface flux extraction. This approach is used in two ways: one procedure, meant to provide a theoretical best-case prediction, uses horizontally averaged elevated concentrations to compare the accuracy of the theoretical flux profiles in calculating surface flux between Kind (1992), Freire et al. (2016), and Nissanka et al. (2018). The N18 model provides a small underestimation of less than 10 % for 10 µm aerosol particle sizes. The second procedure applies the N18 model using concentration segments provided by the simulated instrumentation for streamwise or spanwise moving and stationary probes. While the N18 model records small surface flux variability for 10 µm particles, substantial variability is found for 50 µm particles.

Although the direct flux approach is purely numerical, the approach suggests that potentially large uncertainty in measurements extrapolating the surface flux is inevitable, even with perfect sampling. Again, this result is attributed to the convective roll features. This attribution is further validated through the converse, in which the neutral boundary layer has a significantly lower range of surface flux variability based on lack of vertical coherence: the flow field demonstrates a more local distribution of turbulent transport. In the context of more realistic point or moving sampling strategies in the atmospheric boundary layer, these results suggest that in an unstable boundary layer, the recommended frequency and duration of sampling measurements are constrained to the manner of the observation.

We conclude that, given all the assumptions and constraints initially provided in the study (such as mesoscale and synoptic-scale negligence and a realistically based sampling duration), there is surface flux variability due to coherent turbulent structures that affect sampling across boundary layer height and aerosol particle size. To better constrain this flux variability, the heterogeneous patterns of turbulent flow due to the convective roll features must be taken into consideration. Although the focus of this study emphasizes sea spray aerosol particles over the open ocean, this analysis can in principle be applied to dust or any other kind of particle that is transported from the surface aloft. As this is the first of several studies to better interpret field data, this study acts as a foundation in bridging the gap between numerical simulations and in situ acquisition of aerosol particle transport information.

Code and data availability

The NTLP model specifically designed for this study is publicly available at the GitHub repository: (Park and Richter2022). The data generated from this model are available from the authors upon request.

Author contributions

HJP, JSR, and DHR developed or contributed to all the work in this paper and wrote most of it. LSF analyzed, wrote, and provided informatics regarding the theoretical profile section. CJ provided the SAR imagery and analysis surrounding it.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Special issue statement

This article is part of the special issue “Cloud, Aerosol and Monsoon Processes Philippines Experiment (CAMP2Ex) (ACP/AMT inter-journal SI)”. It is not associated with a conference.


The authors would like to thank Coda Phillips and Charlotte Wainwright for providing insightful discussions on the development of this study. Computational resources were partially provided by the Computing Research Center at the University of Notre Dame and through the High Performance Computing Modernization Program (HPCMP) by ONR, specifically the Topaz, Koehr, and Onyx machines.

Financial support

Hyungwon John Park and David H. Richter were supported by the Office of Naval Research (ONR; grant no. N00014-16-1-2472). Hyungwon John Park was also supported by the National Research Council through the Research Associateship Programs. Jeffrey S. Reid was supported by the NRL base program and ONR (grant no. 322). Livia S. Freire was funded by the São Paulo Research Foundation (FAPESP, Brazil; grant no. 2018/24284-1).

Review statement

This paper was edited by Jing Wei and reviewed by two anonymous referees.


Andreas, E. L.: A new spray generation function for wind speeds up to 32 m s−1, J. Phys. Oceanogr., 28, 10,<2175:ANSSGF>2.0.CO;2, 1998. a

Andreas, E. L., Mahrt, L., and Vickers, D.: An improved bulk air–sea surface flux algorithm, including spray-mediated transfer, Q. J. Roy. Meteor. Soc., 141, 642–654,, 2015. a

Balachandar, S. and Eaton, J. K.: Turbulent dispersed multiphase flow, Annu. Rev. Fluid Mech., 42, 111–133,, 2010. a

Bian, H., Froyd, K., Murphy, D. M., Dibb, J., Darmenov, A., Chin, M., Colarco, P. R., da Silva, A., Kucsera, T. L., Schill, G., Yu, H., Bui, P., Dollner, M., Weinzierl, B., and Smirnov, A.: Observationally constrained analysis of sea salt aerosol in the marine atmosphere, Atmos. Chem. Phys., 19, 10773–10785,, 2019. a, b

Blanchard, D. C., Woodcock, A. H., and Cipriano, R. J.: The vertical distribution of the concentration of sea salt in the marine atmosphere near Hawaii, Tellus B, 36 B, 118–125,, 1984. a, b, c

Chamberlain, A. C.: Transport of lycopodium spores and other small particles to rough surfaces, P. Roy. Soc. Lond. A-Math., 296, 45–70,, 1967. a

Chen, S., Yau, M.-K., Bartello, P., and Xue, L.: Bridging the condensation–collision size gap: a direct numerical simulation of continuous droplet growth in turbulent clouds, Atmos. Chem. Phys., 18, 7251–7262,, 2018. a

Davidson, P.: Turbulence, an introduction for scientists and engineers, Oxford University Press, Oxford,, 2004. a

de Leeuw, G.: Vertical profiles of giant particles close above the sea surface, Tellus B, 38, 51–61,, 1986. a

de Leeuw, G., Andreas, E., Anguelova, M., Fairall, C., Ernie, R., O'Dowd, C., Schulz, M., and Schwartz, S.: Production Flux of Sea Spray Aerosol, Rev. Geophys., 49, RG2001,, 2011. a

de Szoeke, S. P., Skyllingstad, E. D., Zuidema, P., and Chandra, A. S.: Cold pools and their influence on the tropical marine boundary layer, J. Atmos. Sci., 74, 1149–1168,, 2017. a

Deardorff, J. W.: Numerical Investigation of Neutral and Unstable Planetary Boundary Layers, J. Atmos. Sci., 29, 91–115, 1972. a

Deardorff, J. W.: Stratocumulus-capped mixed layers derived from a three-dimensional model, Bound.-Lay. Meteorol., 18, 495–527,, 1980. a

Deike, L., Ghabache, E., Liger-Belair, G., Das, A. K., Zaleski, S., Popinet, S., and Séon, T.: Dynamics of jets produced by bursting bubbles, Physical Review Fluids, 3, 1–20,, 2018. a

Delay, F., Ackerer, P., and Danquigny, C.: Simulating solute transport in porous or fractured formations using random walk particle tracking: A review, Vadose Zone J., 4, 360–379,, 2005. a

Desjardins, R. L., MacPherson, J. I., Schuepp, P. H., and Karanja, F.: An evaluation of aircraft flux measurements of CO2, water vapor and sensible heat, Bound.-Lay. Meteorol., 47, 55–69,, 1989. a

de Szoeke, S. P., Marke, T., and Brewer, W. A.: Diurnal Ocean Surface Warming Drives Convective Turbulence and Clouds in the Atmosphere, Geophys. Res. Lett., 48, e2020GL091299,, 2021. a

Eaton, J. K. and Fessler, J. R.: Preferential Concentration of Particles by Turbulence, Int. J. Multiphas. Flow, 20, 169–209, 1994. a

Fairall, C., Kepert, J., and Holland, G.: The effect of sea spray on surface energy transports over the ocean, Global Atmos. Ocean Syst., 2, 121–142, 1994. a

Freire, L. S., Chamecki, M., and Gillies, J. A.: Flux-profile relationship for dust concentration in the stratified atmospheric surface layer, Bound.-Lay. Meteorol., 160, 249–267,, 2016. a, b, c, d, e, f, g, h, i, j, k, l

Friehe, C. A.: Air–sea fluxes and surface layer turbulence around a sea surface temperature front, J. Geophys. Res., 96, 8593–8609,, 1991. a

Geever, M., O'Dowd, C. D., van Ekeren, S., Flanagan, R., Nilsson, E. D., de Leeuw, G., and Rannik, Ü.: Submicron sea spray fluxes, Geophys. Res. Lett., 32, 2–5,, 2005. a

Gillette, D. A., Blifford, I. H., and Fenster, C. R.: Measurements of aerosol size distributions and vertical fluxes of aerosols on land subject to wind erosion, J. Appl. Meteorol., 11, 977–987,<0977:MOASDA>2.0.CO;2, 1972. a

Gillies, J. and Berkofsky, L.: Eolian suspension above the saltation layer, the concentration profile, J. Sediment. Res., 74, 176–183,, 2004. a

Gong, S. L., Barrie, L. A., Prospero, J. M., Savoie, D. L., Ayers, G. P., Blanchet, J.-P., and Spacek, L.: Modeling sea-salt aerosol in the atmosphere 2. Atmospheric concentrations and fluxes, J. Geophys. Res., 102, 3819–3830,, 1997. a

Grossman, R. L.: An analysis of vertical velocity spectra obtained in the bomex fair-weather, trade-wind boundary layer, Bound.-Lay. Meteorol., 23, 323–357,, 1982. a

Grythe, H., Ström, J., Krejci, R., Quinn, P., and Stohl, A.: A review of sea-spray aerosol source functions using a large global set of sea salt aerosol concentration measurements, Atmos. Chem. Phys., 14, 1277–1297,, 2014. a, b

Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz-Sabater, J., Nicolas, J., Peubey, C., Radu, R., Schepers, D., Simmons, A., Soci, C., Abdalla, S., Abellan, X., Balsamo, G., Bechtold, P., Biavati, G., Bidlot, J., Bonavita, M., De Chiara, G., Dahlgren, P., Dee, D., Diamantakis, M., Dragani, R., Flemming, J., Forbes, R., Fuentes, M., Geer, A., Haimberger, L., Healy, S., Hogan, R. J., Hólm, E., Janisková, M., Keeley, S., Laloyaux, P., Lopez, P., Lupu, C., Radnoti, G., de Rosnay, P., Rozum, I., Vamborg, F., Villaume, S., and Thépaut, J. N.: The ERA5 global reanalysis, Q. J. Roy. Meteor. Soc., 146, 1999–2049,, 2020. a

Hulst, H. C. and van de Hulst, H. C.: Light scattering by small particles, Courier Corporation, New York,, 1981. a

Hutjes, R. W., Vellinga, O. S., Gioli, B., and Miglietta, F.: Dis-aggregation of airborne flux measurements using footprint analysis, Agr. Forest Meteorol., 150, 966–983,, 2010. a, b

Irwin, J. S. and Binkowski, F. S.: Estimation of the Monin–Obukhov scaling length using on-site instrumentation, Atmos. Environ. (1967), 15, 1091–1094,, 1981. a

Jacob, C. and Anderson, W.: Conditionally Averaged Large-Scale Motions in the Neutral Atmospheric Boundary Layer: Insights for Aeolian Processes, Bound.-Lay. Meteorol., 162, 21–41,, 2017. a, b

Jensen, J. B. and Lee, S.: Giant sea-salt aerosols and warm rain formation in marine stratocumulus, J. Atmos. Sci., 65, 3678–3694,, 2008. a, b

Kapustin, V. N., Clarke, A. D., Howell, S. G., Mcnaughton, C. S., Brekhovskikh, V. L., and Zhou, J.: Evaluating primary marine aerosol production and atmospheric roll structures in Hawaii's natural oceanic wind tunnel, J. Atmos. Ocean. Tech., 29, 668–682,, 2012. a, b, c

Keene, W. C., Long, M. S., Reid, J. S., Frossard, A. A., Kieber, D. J., Maben, J. R., Russell, L. M., Kinsey, J. D., Quinn, P. K., and Bates, T. S.: Factors that modulate properties of primary marine aerosol generated from ambient seawater on ships at sea, J. Geophys. Res.-Atmos., 122, 11,961–11,990,, 2017. a, b

Kemppinen, O., Laning, J. C., Mersmann, R. D., Videen, G., and Berg, M. J.: Imaging atmospheric aerosol particles from a UAV with digital holography, Sci. Rep.-UK, 10, 1–12,, 2020. a

Khanna, S. and Brasseur, J. G.: Three-dimensional buoyancy- and shear-induced local structure of the atmospheric boundary layer, J. Atmos. Sci., 55, 710–743,<0710:TDBASI>2.0.CO;2, 1998. a

Kind, R. J.: One-dimensional aeolian suspension above beds of loose particles-a new concentration-profile equation, Atmos. Environ. A-Gen., 26, 927–931,, 1992. a, b, c, d, e, f, g, h

Klemp, J. B. and Durran, D. R.: An Upper Boundary Condition Permitting Internal Gravity Wave Radiation in Numerical Mesoscale Models, Mon. Weather Rev., 111, 430–444,<0430:aubcpi>;2, 1983. a

Lauros, J., Sogachev, A., Smolander, S., Vuollekoski, H., Sihto, S.-L., Mammarella, I., Laakso, L., Rannik, Ü., and Boy, M.: Particle concentration and flux dynamics in the atmospheric boundary layer as the indicator of formation mechanism, Atmos. Chem. Phys., 11, 5591–5601,, 2011. a

Lee, S. H., Gordon, H., Yu, H., Lehtipalo, K., Haley, R., Li, Y., and Zhang, R.: New Particle Formation in the Atmosphere: From Molecular Clusters to Global Climate, J. Geophys. Res.-Atmos., 124, 7098–7146,, 2019. a

LeMone, M.: The Structure and Dynamics of Horizontal Roll Vortices in the Planetary Boundary Layer, J. Atmos. Sci., 30, 1077–1091,<1077:TSADOH>2.0.CO;2, 1973. a, b

Lewis, E. R. and Schwartz, S. E.: Sea salt aerosol production: mechanisms, methods, measurements, and models–A critical review, American Geophysical Union, Washington, DC,, 2004. a, b, c, d, e

Li, D., Wei, A., Luo, K., and Fan, J.: Direct numerical simulation of a particle-laden flow in a flat plate boundary layer, Int. J. Multiphas. Flow, 79, 124–143,, 2016. a

Li, Z., Guo, J., Ding, A., Liao, H., Liu, J., Sun, Y., Wang, T., Xue, H., Zhang, H., and Zhu, B.: Aerosol and boundary-layer interactions and impact on air quality, Natl. Sci. Rev., 4, 810–833,, 2017. a

Mehta, M., Khushboo, R., Raj, R., and Singh, N.: Spaceborne observations of aerosol vertical distribution over Indian mainland (2009–2018), Atmos. Environ., 244, 117902,, 2021. a

Mellado, J. P.: Cloud-top entrainment in stratocumulus clouds, Annu. Rev. Fluid Mech., 49, 145–169,, 2017. a

Moeng, C.-H.: A large-eddy-simulation model for the study of planetary boundary-layer Turbulence, J. Atmos. Sci., 41, 2052–2062,<2052:ALESMF>2.0.CO;2, 1984. a, b, c

Moeng, C.-H. and Sullivan, P. P.: A comparison of shear- and buoyancy-driven planetary boundary layer flows, J. Atmos. Sci., pp. 999–1022,<0999:ACOSAB>2.0.CO;2, 1994. a, b, c

Monahan, E. C., Spiel, D. E., and Davidson, K. L.: A model of marine aerosol generation via whitecaps and wave disruption, Springer, Dordrecht, 167–174,, 1986. a, b

Monaldo, F. M., Jackson, C. R., and Pichel, W. G.: Seasat to Radarsat-2: Research to operations, Oceanography, 26, 34–45, 2013. a

Monaldo, F. M., Li, X., Pichel, W. G., and Jackson, C. R.: Ocean wind speed climatology from spaceborne SAR imagery, B. Am. Meteorol. Soc., 95, 565–569,, 2014. a

Nissanka, I. D., Park, H. J., Freire, L. S., Chamecki, M., Reid, J. S., and Richter, D. H.: Parameterized vertical concentration profiles for aerosols in the marine atmospheric boundary layer, J. Geophys. Res.-Atmos., 123, 9688–9702,, 2018. a, b, c, d, e, f, g, h, i, j, k, l, m, n

Norris, S. J., Brooks, I. M., de Leeuw, G., Smith, M. H., Moerman, M., and Lingard, J. J. N.: Eddy covariance measurements of sea spray particles over the Atlantic Ocean, Atmos. Chem. Phys., 8, 555–563,, 2008. a, b

Norris, S. J., Brooks, I. M., Hill, M. K., Brooks, B. J., Smith, M. H., and Sproson, D. A.: Eddy covariance measurements of the sea spray aerosol flux over the open ocean, J. Geophys. Res.-Atmos., 117, 1–15,, 2012. a, b, c

Pan, Z., Mao, F., Rosenfeld, D., Zhu, Y., Zang, L., Lu, X., Thornton, J. A., Holzworth, R. H., Yin, J., Efraim, A., and Gong, W.: Coarse sea spray inhibits lightning, Nat. Commun., 13, 1–7,, 2022. a

Park, H. and Richter, D.: In situ particle sampling relationships to surface and turbulent fluxes using large eddy simulations with Lagrangian particles software code (1.0.0), Zenodo [code],, 2022. a

Park, H. J., Sherman, T., Freire, L. S., Wang, G., Bolster, D., Xian, P., Sorooshian, A., Reid, J. S., and Richter, D. H.: Predicting vertical concentration profiles in the marine atmospheric boundary layer with a Markov chain random walk model, J. Geophys. Res.-Atmos., 125, e2020JD032731,, 2020. a

Peng, T. and Richter, D.: Sea spray and its feedback effects: Assessing bulk algorithms of air–sea heat fluxes via direct numerical simulations, J. Phys. Oceanogr., 49, 1403–1421,, 2019. a

Porter, J. N. and Clarke, A. D.: Aerosol size distribution models based on in situ measurements, J. Geophys. Res., 102, 6035–6045, 1997. a

Prajapati, J., Shukla, B. P., and Kumar, R.: Orientation of marine atmospheric rolls in a SAR imagery using wavelet transform: A case study over Bay of Bengal, J. Earth Syst. Sci., 130, 34,, 2021. a

Reid, J. S., Jonsson, H. H., Smith, M. H., and Smirnov, A.: Evolution of the vertical profile and flux of large sea-salt particles in a coastal zone, J. Geophys. Res.-Atmos., 106, 12039–12053,, 2001. a, b

Reid, J. S., Brooks, B., Crahan, K. K., Hegg, D. A., Eck, T. F., O'Neill, N., de Leeuw, G., Reid, E. A., and Anderson, K. D.: Reconciliation of coarse mode sea-salt aerosol particle size measurements and parameterizations at a subtropical ocean receptor site, J. Geophys. Res.-Atmos., 111, 1–26,, 2006. a, b, c, d

Reid, J. S., Hyer, E. J., Johnson, R. S., Holben, B. N., Yokelson, R. J., Zhang, J., Campbell, J. R., Christopher, S. A., Di Girolamo, L., Giglio, L., Holz, R. E., Kearney, C., Miettinen, J., Reid, E. A., Turk, F. J., Wang, J., Xian, P., Zhao, G., Balasubramanian, R., Chew, B. N., Janjai, S., Lagrosas, N., Lestari, P., Lin, N. H., Mahmud, M., Nguyen, A. X., Norris, B., Oanh, N. T., Oo, M., Salinas, S. V., Welton, E. J., and Liew, S. C.: Observing and understanding the Southeast Asian aerosol system by remote sensing: An initial review and analysis for the Seven Southeast Asian Studies (7SEAS) program, Atmos. Res., 122, 403–468,, 2013. a

Reid, J. S., Lagrosas, N. D., Jonsson, H. H., Reid, E. A., Sessions, W. R., Simpas, J. B., Uy, S. N., Boyd, T. J., Atwood, S. A., Blake, D. R., Campbell, J. R., Cliff, S. S., Holben, B. N., Holz, R. E., Hyer, E. J., Lynch, P., Meinardi, S., Posselt, D. J., Richardson, K. A., Salinas, S. V., Smirnov, A., Wang, Q., Yu, L., and Zhang, J.: Observations of the temporal variability in aerosol properties and their relationships to meteorology in the summer monsoonal South China Sea/East Sea: the scale-dependent role of monsoonal flows, the Madden–Julian Oscillation, tropical cyclones, squall lines and cold pools, Atmos. Chem. Phys., 15, 1745–1768,, 2015. a, b

Reid, J. S., Lagrosas, N. D., Jonsson, H. H., Reid, E. A., Atwood, S. A., Boyd, T. J., Ghate, V. P., Xian, P., Posselt, D. J., Simpas, J. B., Uy, S. N., Zaiger, K., Blake, D. R., Bucholtz, A., Campbell, J. R., Chew, B. N., Cliff, S. S., Holben, B. N., Holz, R. E., Hyer, E. J., Kreidenweis, S. M., Kuciauskas, A. P., Lolli, S., Oo, M., Perry, K. D., Salinas, S. V., Sessions, W. R., Smirnov, A., Walker, A. L., Wang, Q., Yu, L., Zhang, J., and Zhao, Y.: Aerosol meteorology of Maritime Continent for the 2012 7SEAS southwest monsoon intensive study – Part 2: Philippine receptor observations of fine-scale aerosol behavior, Atmos. Chem. Phys., 16, 14057–14078,, 2016. a

Richter, D. H., Dempsey, A. E., and Sullivan, P. P.: Turbulent Transport of Spray Droplets in the Vicinity of Moving Surface Waves, J. Phys. Oceanogr., 49, 1789–1807,, 2019. a

Richter, D. H., MacMillan, T., and Wainwright, C.: A Lagrangian Cloud Model for the Study of Marine Fog, Bound.-Lay. Meteorol., 181, 523–542,, 2021. a

Ryder, C. L., Highwood, E. J., Walser, A., Seibert, P., Philipp, A., and Weinzierl, B.: Coarse and giant particles are ubiquitous in Saharan dust export regions and are radiatively significant over the Sahara, Atmos. Chem. Phys., 19, 15353–15376,, 2019. a

Salesky, S. T. and Anderson, W.: Buoyancy effects on large-scale motions in convective atmospheric boundary layers: Implications for modulation of near-wall processes, J. Fluid Mech., 856, 135–168,, 2018. a

Schlosser, J. S., Dadashazar, H., Edwards, E., Hossein Mardi, A., Prabhakar, G., Stahl, C., Jonsson, H. H., and Sorooshian, A.: Relationships Between Supermicrometer Sea Salt Aerosol and Marine Boundary Layer Conditions: Insights From Repeated Identical Flight Patterns, J. Geophys. Res.-Atmos., 125, 1–17,, 2020. a

Scipión, D. E., Chilson, P. B., Fedorovich, E., and Palmer, R. D.: Evaluation of an LES-based wind profiler simulator for observations of a daytime atmospheric convective boundary layer, J. Atmos. Ocean. Tech., 25, 1423–1436,, 2008. a

Sessions, W. R., Reid, J. S., Benedetti, A., Colarco, P. R., da Silva, A., Lu, S., Sekiyama, T., Tanaka, T. Y., Baldasano, J. M., Basart, S., Brooks, M. E., Eck, T. F., Iredell, M., Hansen, J. A., Jorba, O. C., Juang, H.-M. H., Lynch, P., Morcrette, J.-J., Moorthi, S., Mulcahy, J., Pradhan, Y., Razinger, M., Sampson, C. B., Wang, J., and Westphal, D. L.: Development towards a global operational aerosol consensus: basic climatological characteristics of the International Cooperative for Aerosol Prediction Multi-Model Ensemble (ICAP-MME), Atmos. Chem. Phys., 15, 335–362,, 2015. a

Srivastava, P. and Sharan, M.: An Analytical Formulation of the Monin–Obukhov Stability Parameter in the Atmospheric Surface Layer Under Unstable Conditions, Bound.-Lay. Meteorol., 165, 371–384,, 2017. a

Stull, R.: An introduction to boundary layer meteorology, Atmospheric and Oceanographic Sciences Library, Springer Netherlands, Dordrecht,, 1988. a

Sühring, M., Metzger, S., Xu, K., Durden, D., and Desai, A.: Trade-Offs in Flux Disaggregation: A Large-Eddy Simulation Study, Bound.-Lay. Meteorol., 170, 69–93,, 2019. a, b, c

Sullivan, P. P. and Patton, E. G.: The effect of mesh resolution on convective boundary layer statistics and structures generated by large-eddy simulation, J. Atmos. Sci., 68, 2395–2415,, 2011. a, b

Sullivan, P. P., McWilliams, J., and Moeng, C.-H.: A grid nesting method for large-eddy simulation of planetary boundary-layer flows, Bound.-Lay. Meteorol., 80, 167–202,, 1996. a

Sutherland, P. and Melville, W. K.: Field measurements of surface and near-surface turbulence in the presence of breaking waves, J. Phys. Oceanogr., 45, 943–965,, 2015. a

Toba, Y. and Chaen, M.: Quantitative expression of the breaking of wind waves on the sea surface, Rec. Oceanogr. Works Jpn., 12, 11 pp., 1973. a

Veron, F.: Ocean spray, Annu. Rev. Fluid Mech., 47, 507–538,, 2015. a, b

Wainwright, C. E., Bonin, T. A., Chilson, P. B., Gibbs, J. A., Fedorovich, E., and Palmer, R. D.: Methods for Evaluating the Temperature Structure-Function Parameter Using Unmanned Aerial Systems and Large-Eddy Simulation, Bound.-Lay. Meteorol., 155, 189–208,, 2015. a, b

Wang, L.-P. and Maxey, M. R.: Settling velocity and concentration distribution of heavy particles in homogeneous isotropic turbulence, J. Fluid Mech., 256, 27–68,, 1993. a

Watson-Parris, D., Schutgens, N., Reddington, C., Pringle, K. J., Liu, D., Allan, J. D., Coe, H., Carslaw, K. S., and Stier, P.: In situ constraints on the vertical distribution of global aerosol, Atmos. Chem. Phys., 19, 11765–11790,, 2019. a

Weckwerth, T. M., Wilson, J. W., Wakimoto, R. M., and Crook, N. A.: Horizontal convective rolls: Determining the environmental conditions supporting their existence and characteristics, Mon. Weather Rev., 125, 505–526,<0505:HCRDTE>2.0.CO;2, 1997.  a

Wei, W., Zhang, H., Wu, B., Huang, Y., Cai, X., Song, Y., and Li, J.: Intermittent turbulence contributes to vertical dispersion of PM2.5 in the North China Plain: cases from Tianjin, Atmos. Chem. Phys., 18, 12953–12967,, 2018. a

Winkler, P.: The growth of atmospheric aerosol particles with relative humidity, Phys. Scripta, 37, 223–230,, 1988. a

Wu, H., Li, Z., Jiang, M., Liang, C., Zhang, D., Wu, T., Wang, Y., and Cribb, M.: Contributions of traffic emissions and new particle formation to the ultrafine particle size distribution in the megacity of Beijing, Atmos. Environ., 262, 118652,, 2021. a

Wu, J.: Individual characteristics of whitecaps and volumetric description of bubbles, IEEE J. Oceanic Eng., 17, 150–158,, 1992. a, b

Wurman, J. and Winslow, J.: Intense sub-kilometer-scale boundary layer rolls observed in hurricane Fran, Science, 280, 555–557,, 1998. a

Wyngaard, J. C.: Turbulence in the atmosphere, Cambridge University Press, Cambridge,, 2010. a

Wyngaard, J. C. and Brost, R. A.: Top-Down and Bottom-Up Diffusion of a Scalar in the Convective Boundary Layer, J. Atmos. Sci., 41, 102–112,, 1984. a

Yamaguchi, R. T., Wang Q., Kalogiros, J., Ruiz-Plancarte, J., Bucholtz, A., Woods, R., Ortiz-Suslow, D., and Barge, J.: Spatial Variability of Optical Turbulence in the Coastal Region from In-situ Measurements Onboard a Novel Low-Flying Aerial Platform, in: Directed Energy Student Workshop, Mobile, AL, 25–29 April 2022, Directed Energy Professional Society, 22-EduWksp-014,, last access: 6 December 2022. a

Executive editor
Next to tackling the fundamental problem in the representativeness of point measurements of particles and retrieving the surface flux in a turbulent atmosphere, this manuscript provides a general framework to plan sampling strategies for aerosol field campaigns and provides tools to quantify related uncertainties.
Short summary
We use numerical models to study field measurements of sea spray aerosol particles and conclude that both the atmospheric state and the methods of instrument sampling are causes for the variation in the production rate of aerosol particles: a critical metric to learn the aerosol's effect on processes like cloud physics and radiation. This work helps field observers improve their experimental design and interpretation of measurements because of turbulence in the atmosphere.