the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

# Estimating scalar turbulent fluxes with slow-response sensors in the stable atmospheric boundary layer

### Mohammad Allouche

### Vladislav I. Sevostianov

### Einara Zahn

### Mark A. Zondlo

### Nelson Luís Dias

### Gabriel G. Katul

### Jose D. Fuentes

Conventional and recently developed approaches for estimating turbulent scalar fluxes under stable atmospheric conditions are evaluated, with a focus on gases for which fast sensors are not readily available. First, the relaxed eddy accumulation (REA) classical approach and a recently proposed mixing length parameterization, labeled A22, are tested against eddy-covariance computations. Using high-frequency measurements collected from two contrasting sites (the frozen tundra near Utqiaġvik, Alaska, and a sparsely vegetated grassland in Wendell, Idaho, during winter), it is shown that the REA and A22 models outperform the conventional Monin–Obukhov similarity theory (MOST) utilized widely to infer fluxes from mean gradients. Second, scenarios where slow trace gas sensors are the only viable option in field measurements are investigated using digital filtering applied to fast-response sensors to simulate their slow-response counterparts. With a filtered scalar signal, the observed filtered eddy-covariance fluxes are referred to here as large-eddy-covariance (LEC) fluxes. A virtual eddy accumulation (VEA) approach, akin to the REA model but not requiring a mechanical apparatus to separate the gas flows, is also formulated and tested. A22 outperforms VEA and LEC in predicting the observed unfiltered (total) eddy-covariance (EC) fluxes; however, VEA can still capture the LEC fluxes well. This finding motivates the introduction of a sensor response time correction into the VEA formulation to offset the effect of sensor filtering on the underestimated net averaged fluxes. The only needed parameter for this correction is the mean velocity at the instrument height, a surrogate of the advective timescale. The VEA approach is very suitable and simple to use with gas sensors of intermediate speed (∼ 0.5 to 1 Hz) and with conventional open- or closed-path setups.

- Article
(6002 KB) - Full-text XML
- BibTeX
- EndNote

The significance of surface–atmosphere exchanges of trace gases, volatile organic compounds (VOCs), and aerosol species to atmospheric composition and thermodynamics is not in dispute. Increasing concentrations of gases and particles due to natural and anthropogenic sources are modulating the Earth's climate and having deleterious consequences for human health and the environment (Qian et al., 2010; Kolb et al., 2010; Voulgarakis et al., 2015). However, estimating these surface–atmosphere exchanges is particularly challenging in the stable atmospheric boundary layer (ABL), which is characterized by weak mixing and highly anisotropic turbulence (Stull, 1988; Mahrt, 1998; Shah and Bou-Zeid, 2014). Stable ABLs occur at nighttime, in the downdraft region of deep mesoscale convective systems (that transport dry air from the middle troposphere to the surface where it is compressed to higher temperatures) and in polar regions; they persist as one of the least understood regimes in boundary layer meteorology, owing to the inherently complex dynamics and the departure from continuous turbulence towards intermittency (Ansorge and Mellado, 2014, 2016; Shah and Bou-Zeid, 2019; Mahrt and Bou-Zeid, 2020; Allouche et al., 2022). On the sensing side, the so-called flux–gradient or flux–variance relations based on Monin–Obukhov similarity theory (MOST) (Monin and Obukhov, 1954) are challenging to apply due to core assumptions that are tenuous to satisfy in practice for stable ABLs. Specifically, a constant flux surface layer that requires stationarity, planar homogeneity, absence of subsidence, and a high Reynolds number state may not be well established for surface flux measurements under stable conditions. The challenges are exacerbated by surface heterogeneity, such as over surfaces with mixed water and sea ice in polar regions that can accelerate the exchanges of gases, aerosols, and energy between the ocean surface and the atmosphere (Sharma et al., 2012; Fogarty and Bou-Zeid, 2023) and semi-infinite heterogeneity patches, e.g., land–sea interfaces (Allouche et al., 2023a, b). These observational challenges then propagate into theoretical and modeling considerations, escalating the need for improved estimates of scalar fluxes under stable conditions.

To begin addressing these challenges and scientific gaps, turbulence and flux observations using the eddy-covariance (EC) technique for fluxes of heat (an active scalar), momentum, and trace gases (representing passive scalars) are employed here. These EC observations are used to evaluate a series of models that can be applied to parameterize turbulent fluxes of gases for which no high-frequency sensors exist, or such sensors are not available commercially at a moderate price. The formulations considered are (i) the relaxed eddy accumulation (REA) technique (Businger and Oncley, 1990) or (ii) aerodynamic approaches that use mean scalar concentrations that are available in coarse weather or climate models (mixing length–mean gradient models). Here, high-frequency velocity and scalar concentration measurements from two contrasting land-cover types are analyzed (i) over an ice sheet in Utqiaġvik (formerly Barrow), Alaska, and (ii) over a sparsely vegetated grassland downwind of heavy agriculture in Wendell, Idaho. The current work seeks to answer the following research questions: (Q1) What flux/closure models can best reproduce the observed EC fluxes in the stable ABL? Models that best describe the observed fluxes are then tested under scenarios that mimic coarse geophysical variables with mean fields measured using slow-response sensors because fast-response instruments remain largely unavailable for reactive chemical species (mainly those characterized by short atmospheric lifetimes). This motivates the second question: (Q2) Can the models correct for the “unresolved” turbulence scales inherently missed when data are collected using slow-response sensors? For this second question, and since the scalar data are all available and collected using fast sensors, the behavior of a slow-response sensor is simulated by digitally filtering the high-frequency data.

## 2.1 Background and definitions

Any instantaneous flow variable (e.g., *s*) is decomposed as $s=\stackrel{\mathrm{\u203e}}{s}+{s}^{\prime}$, where $\stackrel{\mathrm{\u203e}}{s}$ is an “ensemble mean” quantity, and *s*^{′} is a turbulent perturbation defined as a departure from $\stackrel{\mathrm{\u203e}}{s}$. Operationally, primed variables are determined as excursions from the time-averaged state (hereafter indicated by the overbar). The atmospheric stability is quantified using the dimensionless stability parameter $\mathit{\zeta}=z/L$, where *z* is the wall-normal distance from the surface, and *L* is the Obukhov length (Obukhov, 1971). Under stable conditions, which are the focus here, *ζ* > 0.

The strength of the variability in any flow variable is quantified by ${\mathit{\sigma}}_{\mathrm{s}}=(\stackrel{\mathrm{\u203e}}{{s}^{\prime}{s}^{\prime}}{)}^{\mathrm{1}/\mathrm{2}}$, the root mean squared value of *s*^{′}, while the covariance $\stackrel{\mathrm{\u203e}}{{w}^{\prime}{s}^{\prime}}$ is the average net vertical kinematic scalar flux, and *w*^{′}y is the vertical velocity fluctuation. From definitions, *σ*_{s} is related to $\stackrel{\mathrm{\u203e}}{{w}^{\prime}{s}^{\prime}}$ using the correlation coefficient *R*_{ws} defined as

## 2.2 VCC: variable correlation coefficient flux model

Using Eq. (1), a simplified flux model can be defined based on an empirical parameterization of *R*_{ws} as a function of *ζ*:

The empirical relation (*R*_{ws}(*ζ*), a linear fit here) is, in general, non-generalizable as it may be site-specific and dependent on some other meteorological variables or surface conditions.

## 2.3 ACC: averaged correlation coefficient flux model

Again using Eq. (1), one could also test another simplified model with an averaged correlation coefficient 〈*R*_{ws}〉, taken as the mean over all the available observational periods, yielding

In addition to assuming that the correlation coefficient is stability-independent, the same potential drawbacks of the variable correlation coefficient formulation also apply to this model, and the results could not be extrapolated to other sites, where other factors may be present, such as heterogeneity, seasonality, and the influence of synoptic variability, to name a few.

## 2.4 REA: relaxed eddy accumulation flux model

Businger and Oncley (1990) proposed the REA method to compute turbulent scalar fluxes. The REA method is ideally appropriate to use when fast-response sensors are available for *w*^{′} (typically from sonic anemometers) but only slow-response measurements are available for the scalar concentration (slow trace gases sensors or even trace gas samples that need to be collected and analyzed subsequently). In such cases, the REA approach offers an enhanced representation of these scalar fluxes (Nie et al., 1995).

The basic idea is inspired by the work of Desjardins (1977), who used conditional sampling techniques to collect scalar information (along with vertical wind speed) in two electronic counters, one for upflow and another for downflow. From linear correlation analysis, the regression slope of ${w}^{\prime}/{\mathit{\sigma}}_{\mathrm{w}}$ against ${s}^{\prime}/{\mathit{\sigma}}_{\mathrm{s}}$ may be estimated from the correlation coefficient following (Baker et al., 1992; Katul et al., 1994, 2018)

where $\stackrel{\mathrm{\u203e}}{{s}^{+}}$ is the conditional average of scalar *s* instantaneously attributed to updraft events (${w}^{\prime}>\mathrm{0}$), and $\stackrel{\mathrm{\u203e}}{{s}^{-}}$ is the conditional average of scalar *s* instantaneously attributed to downdraft events (${w}^{\prime}<\mathrm{0}$), $\mathrm{\Delta}s=\stackrel{\mathrm{\u203e}}{{s}^{+}}-\stackrel{\mathrm{\u203e}}{{s}^{-}}$ reflects this difference in collecting scalar information from the two samples and likewise for the vertical velocity statistics. When this estimate for *R*_{ws} is inserted into Eq. (1), the REA expression emerges as

For a Gaussian-distributed *w*^{′}, it can be shown that ${\mathit{\sigma}}_{\mathrm{w}}/(\stackrel{\mathrm{\u203e}}{{w}^{+}}-\stackrel{\mathrm{\u203e}}{{w}^{-}})=\sqrt{\mathrm{2}\mathit{\pi}}/\mathrm{4}$ (Katul et al., 2018), a constant whose numerical value is 0.63.

Since the linear regression analysis to estimate *R*_{ws} as featured in Eq. (4) is imperfect, an operational REA model is typically used and is expressed as

where *β*_{s} is now treated as an empirical coefficient that corrects for the abovementioned shortcomings of such a slope estimation of *R*_{ws}. Many studies investigated the choice of optimal *β*_{s} over a wide range of atmospheric stabilities, surfaces, and meteorological conditions (Businger and Oncley, 1990; Katul et al., 1996; Milne et al., 1999; Zahn et al., 2016; Vogl et al., 2021). The choice of *β*_{s} is still debatable, yet various studies reported a *β*_{s} ≈ 0.59 (Bowling et al., 1998; Katul et al., 1996), which is not far from a Gaussian prediction derived from *w*^{′} statistics (i.e., 0.63). Hence, a *β*_{s}=0.59 is selected in the current study as a reference baseline in assessing the REA method.

What is less debatable is the expectation of theoretical invariance of *β*_{s} with stability changes: it was recently shown that the required independence of the REA formulation in the limit of free convection from the friction velocity (*u*_{∗}) is not compatible with a stability-dependent *β*_{s} (Zahn et al., 2023), and this stability invariance was in fact reported in many field observational studies. The arguments of Zahn et al. (2023) for a stability-invariant *β*_{s} can be deduced from the dimensionless form of the REA expression

where ${s}_{\ast}=\stackrel{\mathrm{\u203e}}{{w}^{\prime}{s}^{\prime}}/{u}_{\ast}$, with $\stackrel{\mathrm{\u203e}}{{s}^{+}}-\stackrel{\mathrm{\u203e}}{{s}^{-}}\sim {\mathit{\sigma}}_{\mathrm{s}}$. Noting that scalar flux–variance expressions of ${\mathit{\sigma}}_{\mathrm{s}}/{s}_{\ast}$ exhibit opposite scaling exponents with *ζ* compared to ${\mathit{\sigma}}_{\mathrm{w}}/{u}_{\ast}$ across all stability regimes, the dependence of *β*_{s} on *ζ* is likely to be small as the two terms on the right-hand side of Eq. (7) cancel each others' stability dependence. In convective conditions, ${\mathit{\sigma}}_{\mathrm{w}}/{u}_{\ast}\sim |\mathit{\zeta}{|}^{+\mathrm{1}/\mathrm{3}}$, whereas ${\mathit{\sigma}}_{\mathrm{s}}/{s}_{\ast}\sim |\mathit{\zeta}{|}^{-\mathrm{1}/\mathrm{3}}$. For near-neutral conditions, MOST predictions suggest ${\mathit{\sigma}}_{\mathrm{w}}/{u}_{\ast}\sim |\mathit{\zeta}{|}^{\mathrm{0}}$ and ${\mathit{\sigma}}_{\mathrm{s}}/{s}_{\ast}\sim |\mathit{\zeta}{|}^{\mathrm{0}}$ as well, making *β*_{s} also independent of stability in that limit.

Under stable conditions, similar plausibility arguments for a stability-independent *β*_{s} can be made based on the observations of Weaver (1990) that ${\mathit{\sigma}}_{\mathrm{s}}/{s}_{\ast}\sim |\mathit{\zeta}{|}^{\mathrm{0}}$ under very stable conditions. This result was explained by the author based on arguments first presented by Wyngaard (2010) that under very stable conditions the active eddy size scales with *L* rather than *z*, and thus turbulence statistics should become independent of *ζ*. This would then also apply to ${\mathit{\sigma}}_{\mathrm{w}}/{u}_{\ast}$ and by extension to *β*. However, as later shown in the present paper, (i) the practical application of REA uses devices with finite mechanical response time to physically separate the accumulation of the trace gas in updrafts and downdrafts; (ii) a “dead band” is introduced at small *w*^{′}, where the concentrations are counted neither towards $\stackrel{\mathrm{\u203e}}{{s}^{+}}$ nor towards $\stackrel{\mathrm{\u203e}}{{s}^{-}}$; and (iii) the slow response of the scalar sensors may all induce an indirect stability dependence.

## 2.5 VEA: virtual eddy accumulation flux model

Despite its inherent assumptions and some remaining challenges, REA is at present the main tool for estimating fluxes with slow-response gas sensors. However, the device it requires is not simple to design and construct and is not yet available off the shelf. Therefore, we propose and test a version where the separation of air streams into downdrafts and updrafts is done virtually, without a mechanical device. The resulting virtual eddy accumulation (VEA) uses the same REA expression in Eq. (6); however, it assumes that the scalar concentrations are measured using a slow-response gas sensor. In this paper, we simulate this slow signal based on the actual high-frequency scalar concentration measurements using a digital filter; the details will be provided in Sect. 4.2. This method offers the alternative of using cheaper and simpler slow-response scalar sensors instead of setting up an REA apparatus. Sensors operating at measurement speeds slower than the 10 Hz required for eddy covariance but still somewhat rapid (e.g., 0.5 to 1 Hz) are now more widely available and less expensive and cover more gas species (Chen et al., 2021; Burba, 2022). Furthermore, slower gas sampling speeds reduce the needs for massive pumps in closed-path sensors, decreasing overall power usage and bulkiness, thus easing field deployability, especially in remote areas (Burba, 2022). Slower-response sensors also tend to have higher precision, other factors being held equal, as white noise decreases with increasing averaging times (Hodgkinson and Tatam, 2012). In laser-based gas analyzers operating at rates slower than 10 Hz, spectroscopic and temperature corrections can be done in real time, simplifying post-processing and reducing frequency-correction errors (Burba et al., 2019). All these advantages motivate the investigation of such a virtual counterpart to REA that can work with moderate speed sensors.

Another advantage of the VEA is that, in theory, it can be applied over averaging periods comparable to eddy covariance (e.g., as short as 15 min), overcoming some of the challenges of conventional REA or of disjunct eddy covariance (DEC) (another technique for measuring fluxes with a slow gas sensor) (Rinne et al., 2001; Karl et al., 2002; Ruppert et al., 2002; Rinne et al., 2008) that usually require longer averaging time periods of 30 to 60 min to compute fluxes. Such long sampling durations, while improving statistical convergence, might be problematic in the stable ABL, especially when turbulence is intermittent, and thus non-stationarity becomes a concern.

## 2.6 LEC: large-eddy-covariance model

A close analogue of the classic eddy-covariance technique can be obtained when the EC calculations are applied as usual, but here the scalar sensor has a slow physical response time. We refer to this calculation as the large-eddy-covariance (LEC) approach. The sensor may still be sampled at a high rate, equal to that of *w*^{′}, to compute the LEC flux as $\stackrel{\mathrm{\u203e}}{{w}^{\prime}{s}^{\prime}}$, but the method should be cognizant of the inherent physical filtering of the fluxes carried by fast eddies that the slow scalar sensor cannot resolve.

## 2.7 A22: mixing length flux model

Recently proposed models for momentum and heat fluxes based on mixing length analogies (Allouche et al., 2022) that outperformed MOST during stable periods marked with intermittent turbulence dynamics are also tested here. These models were initially formulated using an eddy diffusion representation of fluxes:

The eddy diffusivity (*K*_{s}) was then defined as the product of a characteristic velocity scale (*U*_{char}) and a mixing length scale (*L*_{mix}): *K*_{s}=*U*_{char}*L*_{mix}. Here, *L*_{mix} will be defined differently for momentum (*K*_{u}) and heat (${K}_{{T}_{\mathrm{v}}}$, where *T*_{v} is the virtual temperature), but both use the standard deviation of the vertical velocity (*σ*_{w}) as the characteristic velocity scale (similar to REA); i.e., *K*_{s}=*σ*_{w}*L*_{mix}.

For momentum, *L*_{mix}=*L*_{u}, and *L*_{u} is defined as a harmonic average between two competing shear length scales (*L*_{u1}, a local turbulent shear scale, and *L*_{u2}, the classic bulk shear scale) as follows

In this model, *α*_{u} is an empirical constant; its value is determined as ${\mathit{\alpha}}_{\mathrm{u}}={\mathit{\alpha}}_{{T}_{\mathrm{v}}}$ = 0.35 (same value found to be also adequate for the heat flux model described next). The mean wind speed at the measurement height is given by $\stackrel{\mathrm{\u203e}}{u}$.

Similarly for heat, ${L}_{\mathrm{mix}}={L}_{{T}_{\mathrm{v}}}$, and ${L}_{{T}_{\mathrm{v}}}$ is defined as a harmonic average between two competing length scales. The first is ${L}_{{T}_{\mathrm{v}\mathrm{1}}}$, the Ellison length scale (Ellison, 1957) and the second is ${L}_{{T}_{\mathrm{v}\mathrm{2}}}$, the buoyancy length scale (Stull, 1973; Zeman and Tennekes, 1977). These scales are formulated as

Here,

is the Brunt–Väisälä frequency; *g* is the gravitational acceleration; and TPE is the turbulent potential energy, which is related to *N*_{BV} as shown in ${L}_{{T}_{\mathrm{v}\mathrm{1}}}$ (Zilitinkevich et al., 2013; Katul et al., 2014).

## 2.8 MOST: Monin–Obukhov similarity theory flux model

Based on dimensional analysis, Monin and Obukhov (1954) formulated flux–gradient relations that are still used widely in weather prediction and climate models. MOST has inherent limitations as it applies to planar homogeneous conditions and stationary flows at very high Reynolds number in the absence of subsidence and requires turbulent kinetic energy (TKE) production to be balanced by the TKE dissipation rate. Nevertheless, MOST still serves as a reference for idealized conditions (Foken, 2006). MOST fluxes could still capture the observed fluxes under weakly stable conditions where turbulence is continuously sustained and not intermittently suppressed (Shah and Bou-Zeid, 2014). MOST fluxes here are computed using the Businger–Dyer relations (Businger et al., 1971) as those relations remain pervasively in use today. Such relations are expressed by non-dimensional gradient (diabatic) functions, Ψ_{s}(*ζ*), relating the scalar concentration surface scale ${s}_{\ast}=\stackrel{\mathrm{\u203e}}{{w}^{\prime}{s}^{\prime}}/{u}_{\ast}$ to the gradient following

where *κ* is the von Kármán constant (i.e., 0.4).

In this study, data from two field experiments are analyzed. One data set is collected over the frozen tundra near Utqiaġvik, Alaska (U09), as part of the OASIS-2009 (Ocean–Atmosphere–Sea Ice–Snowpack) field campaign (Staebler et al., 2009; Perrie et al., 2012). The second data set is collected from November 2022 to January 2023 at a sparsely vegetated grassland in Wendell, Idaho (W22). At Utqiaġvik, four sonic anemometers were mounted on a 10 m tall tower at 0.58, 1.8, 3.2, and 6.2 m above the snowpack, and the data analyzed herein correspond to *z*_{m} = 1.8 m. The lowest anemometer (model TR90-AH, Kaijo Denki, Japan) had a 5 cm pathlength and provided data at a frequency of 20 Hz, while the other three sonic anemometers (model CSAT3, Campbell Scientific Inc., Logan, Utah) each had a 10 cm pathlength and provided data at a frequency of 10 Hz. Three-dimensional velocity (*u*, *v*, *w*: longitudinal, lateral, and vertical components) and sonic virtual temperature (*T*_{s}≈*T*_{v}, where *T*_{v} is the true virtual temperature) measurements were recorded. At the Wendell site, data were acquired only at one height (*z*_{m} = 2.4 m) above the ground surface. Chemical scalar concentration (carbon dioxide, CO_{2}; ammonia, NH_{3}; and water vapor, H_{2}O), in addition to the three-dimensional velocity and temperature measurements, were recorded using a commercial open-path analyzer (${\mathrm{CO}}_{\mathrm{2}}/{\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ 7500A, LiCor Inc., Lincoln, NE), a custom-made open-path sensor with a quantum cascade laser (NH_{3}), and an R.M. Young 81000 sonic anemometer (Miller et al., 2014; Sun et al., 2015; Pan et al., 2021). Network Time Protocol (NTP) was used to ensure all the gas, environmental, and meteorological sensors were synchronized by GPS (Global Positioning System). For both sites, instantaneous molar density measurements of the chemical species were converted to mass concentrations using the pressure sensor on the LiCor 7500A, which has an accuracy of ± 0.4 kPa from 50 to 110 kPa and a resolution of 0.006 kPa (Edson et al., 2011). The gas fluxes were then calculated based on their mass concentrations, in lieu of applying the so-called WPL (Webb, Pearman, and Leuning) density corrections to the fluxes after processing (Webb et al., 1980; Detto and Katul, 2007).

The sampling frequency at both sites was set to *f*_{s} = 10 Hz, and the post-processing involved (i) de-spiking, (ii) linear detrending, and (iii) double rotation of wind components (Wilczak et al., 2001). Fluxes and other required statistics were then computed for various scalar quantities (i.e., *T*_{v}, CO_{2}, NH_{3}, and H_{2}O; *u* and its associated momentum flux were also tested here for comparison). Analysis periods were set to 15 min in U09 (the 15 min Reynolds average choice here is selected because U09 periods reveal strong intermittent behavior) and 30 min in W22. These were then the periods used for double rotation and time averaging throughout. Details of the data quality control for these data sets can be found elsewhere (Allouche et al., 2022).

An earlier study (Zahn et al., 2023) investigated REA under non-ideal unstable conditions and concluded that the REA method outperforms MOST flux models. One of the main aims of the present study is to examine whether REA outperforms MOST under stably stratified conditions as well. MOST is used as a reference for comparison as it reflects the “state of the science” in climate models. Since A22 established the limitations of MOST under stable conditions for the Utqiaġvik data set and further proposed the closure models detailed previously that also outperform MOST, REA will then be compared to A22. Other model details in Sect. 2 will serve as additional benchmarks to understand model performance, but the analyses focus on REA and A22 and the VEA approach that we will introduce later.

## 4.1 Model inter-comparison using high-frequency measurements

All introduced models pertinent to high-frequency measurements are now tested at the Utqiaġvik site because it has the multiple levels that are required for testing the A22 and MOST models. The middle panel subplots of Fig. 1, corresponding to the REA and A22 models, depict the strongest correlation between modeled and observed EC fluxes for both momentum and heat. In addition, the ACC model, which incorporates stability information (Eq. 2), performs slightly better than the constant ACC model (Eq. 3). Although these models at first may appear to be better approximations of the eddy-covariance fluxes than REA, the REA model performance is in practice superior, benefiting from the cancellation of the effect of stability in the model coefficient *β*_{s}, as detailed previously. This agrees with prior findings (Zahn et al., 2023) for unstable conditions. All proposed models outperform MOST (Fig. 1c–f). With their superior performance established, the REA and A22 models will be the only ones retained in the subsequent analyses. When filtering is applied, the VEA model replaces the REA model.

## 4.2 Simulating a slow scalar sensor for model testing

To address the limited bandwidth of many trace gas sensors, we simulate the output of a real slow sensor $\stackrel{\mathrm{\u0303}}{s}$ measuring a variable *s* as the numerical solution to the first-order ODE in Eq. (13). Here, *s* would be the “fast” turbulent sensor, and $\stackrel{\mathrm{\u203e}}{\mathrm{\Delta}}$ is the timescale of the filter width (the response timescale of the slow sensor), which we selected to vary in the range [0–5 s] with increments of 1 s. Thus, solving the following equation numerically (using explicit forward Euler time advancement scheme) converts the *s* time series from 10 Hz to a lower frequency down to 0.2 Hz when $\stackrel{\mathrm{\u203e}}{\mathrm{\Delta}}$ = 5 s:

The observed filtered EC fluxes computed with the filtered scalar signal $\stackrel{\mathrm{\u0303}}{s}$ are referred to as the LEC fluxes. We also tested another filter type, a low-pass Gaussian filter, and the obtained results were not sensitive to the filter type, but the ODE solution is a more accurate model for a first-order slow-sensing system. The filtering is not applied to the vertical velocity (*w*) as high-frequency anemometers are readily available. Hence, a tilde denotes the filtered virtual temperature ($\stackrel{\mathrm{\u0303}}{{T}_{\mathrm{v}}}$), scalar concentration ($\stackrel{\mathrm{\u0303}}{s}$), and streamwise velocity ($\stackrel{\mathrm{\u0303}}{u}$). All three-dimensional velocity components (*u*, *v*, and *w*) are available in high-frequency output of the sonic anemometer, but only *u* is filtered ($\stackrel{\mathrm{\u0303}}{u}$) to compare its kinematic momentum fluxes to those of scalars.

In theory, an REA system should not suffer from the slow response of the trace gas sensor since it only requires the mean measurements of $\stackrel{\mathrm{\u203e}}{{s}^{+}}$ and $\stackrel{\mathrm{\u203e}}{{s}^{-}}$. However, this requires a mechanism to separate the gas streams that has a 10 Hz response time as well. While such systems have been constructed and used, they are typically operated using a certain threshold value that defines a dead-band velocity *w*_{0} below which the system does not switch intake or keeps both intakes closed. This is designed to avoid an excessive number of movements and to guarantee larger individual air samples. In addition, since these systems are typically custom made, there is still a possibility of some latency in the response of the mechanics of some models. In the present implementation, a dynamic dead band that is linked to the turbulence conditions in each period is adopted such that similar amounts of air are sampled in updrafts and downdrafts. We used the empirical finding of Grelle and Keck (2021), as depicted in Eq. (14), which yielded roughly similar amounts of air in the updraft and downdraft reservoirs. Sampling is only activated if the vertical wind exceeds this threshold value in *w*_{0}, i.e., |*w*|>*w*_{0}.

We should underline the fact that the computations of $\stackrel{\mathrm{\u203e}}{{s}^{+}}$ and $\stackrel{\mathrm{\u203e}}{{s}^{-}}$ for the REA are first done with the unfiltered *s* signal (mimicking a flawless and very fast mechanical system), and the dead band of vertical velocity defined above is then the only effective filter that applies to the REA computations. However, when exploring the VEA approach that we will detail later, we will compute the VEA fluxes with $\stackrel{\mathrm{\u203e}}{{s}^{+}}$ and $\stackrel{\mathrm{\u203e}}{{s}^{-}}$ from the filtered signal produced by Eq. (13).

## 4.3 LEC, VEA, and A22 model evaluation using simulated slow-sensor data

Again focusing on Utqiaġvik data set with its multilevel measurements, the models are now tested using inputs from a slow sensor. The top and middle panel subplots of Fig. 2 show that the LEC and VEA methods significantly, and in similar fashion, underestimate the observed (unfiltered) heat and momentum EC fluxes based on filtered quantities $\stackrel{\mathrm{\u0303}}{{T}_{\mathrm{v}}}$ and $\stackrel{\mathrm{\u0303}}{u}$, respectively, as $\stackrel{\mathrm{\u203e}}{\mathrm{\Delta}}$ increases. This underestimation, however, is not surprising because under stable conditions small eddies carry a significant proportion of the fluxes, especially when the background flow is laminarizing and intermittent, as shown elsewhere (Ansorge and Mellado, 2014, 2016; Allouche et al., 2022; Issaev et al., 2023). These small eddies are filtered appreciably by the sensor's slow response. We should underline here that an REA system with fast-response mechanical valves will give results equivalent to the unfiltered VEA (VEA recovers an REA system for $\stackrel{\mathrm{\u203e}}{\mathrm{\Delta}}$ = 0) and will be in good agreement with the actual EC fluxes. However, any latency in the mechanical response of the REA device will introduce some type of filtering to the signal that depends on the device design; this underlines the importance of the mechanical design of these systems. Remarkably, the bottom-panel subplots of Fig. 2c and f show that the A22 model's performance is not sensitive to the signal filtering and provides good estimates of the observed (unfiltered) heat and momentum EC fluxes, even when filtered quantities $\stackrel{\mathrm{\u0303}}{{T}_{\mathrm{v}}}$ and $\stackrel{\mathrm{\u0303}}{u}$ are used, and up to the highest filter width $\stackrel{\mathrm{\u203e}}{\mathrm{\Delta}}$ = 5 s. The A22 model relies on multilevel means and variances in computing the fluxes, which tend to be carried by larger scales than the actual fluxes; this may explain the independence of the model performance from sensor response.

Given the results above, a followup question is whether the VEA model estimates with a filtered signal actually correspond to the fluxes that would be computed using eddy covariances of the filtered scalar signal, i.e., the LEC fluxes. Figure 3a and b indeed show that the VEA model is still a reliable method to capture these filtered observed heat and momentum LEC fluxes at Utqiaġvik as $\stackrel{\mathrm{\u203e}}{\mathrm{\Delta}}$ increases. All quantities here, including the fluxes, are computed based on filtered series $\stackrel{\mathrm{\u0303}}{{T}_{\mathrm{v}}}$ and $\stackrel{\mathrm{\u0303}}{u}$. This implies that VEA flux estimates are broadly comparable to the filtering of the LEC fluxes by slow-response scalar sensors.

This match between VEA and LEC fluxes is also observed at the Wendell site for heat and momentum, as well as for the other scalars available at that site (CO_{2}, NH_{3}, and H_{2}O). As in Utqiaġvik's site, Figs. 4b, d and 5b, d, and f show that VEA performs better when evaluated against the LEC fluxes, compared to Figs. 4a, c and 5a, c, and e that compare it to total EC fluxes. Note that since only one level of measurements at Wendell is available, gradients of first-order moments could not be computed, which precludes testing of MOST or A22 models.

These similar VEA findings among the two different sites indicate that, under stable conditions, the VEA captures the fluxes of the “resolved” eddies. Any slow-response filtering will cause the method to significantly underestimate the needed high-frequency, full observed fluxes. An important question (addressed in the next subsection) that follows is whether the VEA fluxes can be “corrected” under stable conditions to recover the missed fluxes. The proposed modifications, outlined below for the VEA methodological framework, offer a correction for VEA flux estimates to recover their corresponding total EC fluxes. It is to be noted that under unstable conditions at the Wendell site, where flux carrying eddies are of much larger scales than under stable conditions, the VEA method is found to be almost insensitive to the considered filter widths ($\stackrel{\mathrm{\u203e}}{\mathrm{\Delta}}$'s) (refer to Appendix A). Therefore, VEA performs well and captures the observed (high- and low-frequency) fluxes under unstable convective regimes, and hence biases in predicting scalar fluxes are expected to be minimal (Figs. A1 and A2). This also implies that for conventional REA, the mechanical response speed is not as critical under unstable conditions. In the next subsection, we aim to develop an effective correction for the under-resolved fluxes under stable conditions that will provide a method to obtain continuous accurate fluxes using VEA over the whole diurnal cycle.

## 4.4 A sensor-response correction for the optimal VEA coefficient *β*_{v}

Figure 3 reveals that the scatter between VEA and EC or LEC results is larger when the fluxes are small, which would translate into larger errors and scatter in the values of *β*_{s}. Analyses not shown here also reveal that *β*_{s} values calculated for each period converge well towards the 0.59 value when the correlation *R*_{ws} increases; i.e., *R*_{ws} > 0.2, indicating larger fluxes, with more scatter for lower correlation values. However this scatter is randomly distributed around the 1:1 line, indicating some error cancellation when the fluxes are integrated over longer periods of time. Discussions on these random variations in *β*_{s} have linked them to the effect of height above canopies (Gao, 1995) and to the energy content influence of the associated eddy motions (Katul et al., 1996), among others, and are not a focus of the present paper.

Figures 2a and d and 4a and c, on the other hand, reveal that the chosen value of *β*_{s} = 0.59 in the modeled fluxes is a good estimate in recovering the observed EC fluxes when the signal is not filtered ($\stackrel{\mathrm{\u203e}}{\mathrm{\Delta}}$ = 0 s). Missing smaller eddies (when the signal is filtered) that contribute significantly to fluxes under stable, but not unstable, conditions thus requires a larger *β*_{s} to predict the correct fluxes using VEA. Such underestimation was attributed to the filtering operation, dictated here by the choice of $\stackrel{\mathrm{\u203e}}{\mathrm{\Delta}}$. This motivates a model development for *β*_{s} that incorporates the effect of filtering of the fast eddies, which we will express as *β*_{v}.

Furthermore, the agreement between the LEC and VEA predictions when the scalar signal is filtered also opens the possibility of applying the VEA method as a surrogate for REA, without an actual device that separates air streams from downdrafts and updrafts. If a slow-response sensor (open- or closed-path) is available, the VEA method that uses the REA equation (Eq. 6) can be applied, with $\stackrel{\mathrm{\u203e}}{{s}^{+}}$ and $\stackrel{\mathrm{\u203e}}{{s}^{-}}$ computed using conditional averaging of the scalar time series based on the sign of *w* – similar to what is done to the simulated REA measurements – but without the dead band given in Eq. (14).

For this purpose, we first compute each period's optimal *β*_{v} (i.e., by inverting the expression in Eq. 6) that causes the REA-predicted fluxes, equivalent to VEA computed with the raw signal without filtering, to match the exact observed fluxes when $\stackrel{\mathrm{\u203e}}{\mathrm{\Delta}}$ = 0 s (no filtering) across the two contrasting sites. The top-panel subplots, Fig. 6a for heat, Fig. 6c for momentum, and Fig. 6e for all three passive trace gases (CO_{2}, NH_{3}, and H_{2}O), show a scatter plot of these exact *β*_{v}'s relative to the integral timescale *τ*_{int}(*s*) for each period's co-spectrum, where the data points are colored with the MOST stability parameter. The integral timescale (*τ*_{int}) is determined by integrating the autocorrelation function (${\mathit{\rho}}_{{w}^{\prime}{s}^{\prime}}\left(\mathit{\tau}\right)$) up to the first zero-crossing from the measured ${w}^{\prime}{s}^{\prime}$ instantaneous time series. We observe *β*_{v} values that depart from *β*_{s} = 0.59 under all stabilities, but in general there is no clear *β*_{v} dependence on *τ*_{int} as depicted here.

For $\stackrel{\mathrm{\u203e}}{\mathrm{\Delta}}\ne \mathrm{0}$ s, the corresponding bottom-panel subplots, Fig. 6b, d, and f, show these same *β*_{v}'s relative to $\stackrel{\mathrm{\u203e}}{\mathrm{\Delta}}/{\mathit{\tau}}_{\mathrm{a}}$ for different $\stackrel{\mathrm{\u203e}}{\mathrm{\Delta}}$'s, where *τ*_{a} is an advective timescale. After experimenting with various choices of timescales for normalizing the filter scale $\stackrel{\mathrm{\u203e}}{\mathrm{\Delta}}$, an advective timescale formed by *z* and $\stackrel{\mathrm{\u203e}}{u}$, hereafter labeled as ${\mathit{\tau}}_{\mathrm{a}}=(\mathit{\kappa}z/\stackrel{\mathrm{\u203e}}{u})$, appears to provide the best scaling for the variations of *β*_{v} with the filter size (*κ* is the von Kármán constant). This converges with the work of Horst (1997), who formulated corrections to estimate the attenuation of scalar flux measurements by the slow response of sensors (akin to LEC). They used the sensor response frequency and a normalized frequency formulated based on *τ*_{a}, which is the frequency of the peak of the logarithmic cospectrum ${f}_{\mathrm{m}}={n}_{\mathrm{m}}\mathit{\kappa}/{\mathit{\tau}}_{\mathrm{a}}$, to correct for the missed fluxes. In that work, *n*_{m} is the dimensionless frequency at the cospectral maximum, where it is estimated from observations of its behavior as a function of atmospheric stability *ζ*. The advective timescale was similarly found to be a plausible choice in describing the drift and non-linear diffusion terms of a proposed non-linear Langevin equation to model the turbulent kinetic energy in stably stratified ABL (Allouche et al., 2021). This characteristic timescale, *τ*_{a}, measures the advection time of the attached eddies, of size *κ**z*, past a fixed sensor.

As depicted in Fig. 6, an empirical fit that relates *β*_{v} to $\stackrel{\mathrm{\u203e}}{\mathrm{\Delta}}/{\mathit{\tau}}_{\mathrm{a}}$ i.e., ${\mathit{\beta}}_{\mathrm{v}}=f(\stackrel{\mathrm{\u203e}}{\mathrm{\Delta}}/{\mathit{\tau}}_{\mathrm{a}})$, is proposed here to recover the real observed (unfiltered) EC fluxes using the VEA method with ($\stackrel{\mathrm{\u0303}}{s}$) measurements. This relation is best described using a power-law model, ${\mathit{\beta}}_{\mathrm{v}}=a(\stackrel{\mathrm{\u203e}}{\mathrm{\Delta}}/{\mathit{\tau}}_{\mathrm{a}}{)}^{b}+c$, and Table 1 summarizes these empirical coefficients (*a*,*b*) for momentum, heat (active scalar), and passive scalar (CO_{2}, NH_{3}, and H_{2}O) fluxes at both sites, where $c={\mathit{\beta}}_{\mathrm{v}}(\stackrel{\mathrm{\u203e}}{\mathrm{\Delta}}=\mathrm{0})={\mathit{\beta}}_{\mathrm{s}}$. The reported *b* exponent for all fluxes, which describes how *β*_{v} scales with $\stackrel{\mathrm{\u203e}}{\mathrm{\Delta}}/{\mathit{\tau}}_{\mathrm{a}}$, varies between active and passive scalars. Note that if the dead-band criterion is removed, *b* ≈ 0.7 and does not vary much among all scalars as shown in Table 2, hinting at the possible universality of such dependence for a VEA model that removes the need for a complex mechanical REA system. Nevertheless, further exploration at disparate sites and analyses of observational data for different scalars across wider stability ranges and different non-ideal surfaces are needed to have increased statistical confidence in the reported values of (*a*,*b*) and their generalizability, especially for the Utqiaġvik site, which does not have measurements of the passive scalars.

A common feature for all fluxes, as depicted in the bottom panel of Fig. 6, is that the proposed model becomes less certain as $\stackrel{\mathrm{\u203e}}{\mathrm{\Delta}}$ increases (as expected). Therefore, such a model becomes less reliable with very slow sensors that cannot resolve most of the scales of turbulence and is ideally suited for moderately slow sensors that can still resolve the larger eddies. Otherwise, to compute *β*_{v}, the only needed inputs are (i) $\stackrel{\mathrm{\u203e}}{\mathrm{\Delta}}$ (ideally provided by the slow-sensor manufacturer) and (ii) *τ*_{a} (computed from mean wind measurements). The reported fitting parameters in Table 1 are obtained using the least absolute residual (LAR) method, so that extreme values, which occur less frequently and may be related to unusual conditions or measurement errors, have a lesser influence on the fit. Given the variability in the optimal values of *β*_{v}, the proposed models would be mostly suited for quantification of long-term aggregates (net averaged fluxes) of the scalars cycle.

Conventional and novel formulations are assessed to predict scalar fluxes under stable conditions. These conditions pose measurement challenges because the turbulence timescales are shorter, thereby amplifying the effects of any spatial or temporal averaging by sensors – or even instrument separation between *w*^{′} and *c*^{′}. The models tested include the MOST-based flux–gradient model; the A22 gradient- and variance-based flux model; the conventional REA approach; and VEA, a virtual version of REA requiring no flow separation and collection device that can be used with moderately slow gas sensors. The test used measurements collected at two different sites (Utqiaġvik and Wendell).

It was found that the REA and A22 models outperform the conventional models (e.g., MOST) and are thus less sensitive to departures from ideal flow conditions of homogeneity, steadiness, negligible vertical transport of covariances, and TKE production–dissipation balance. The A22 model was also found to perform well and provide good estimates of the observed (unfiltered) heat and momentum EC fluxes, even when filtered quantities $\stackrel{\mathrm{\u0303}}{{T}_{\mathrm{v}}}$ and $\stackrel{\mathrm{\u0303}}{u}$ are used and up to the highest filter width $\stackrel{\mathrm{\u203e}}{\mathrm{\Delta}}$ = 5 s. This is because the A22 model is insensitive to filtering operations of the turbulent scales as it relies on multilevel means and variances in computing the fluxes, which tend to be carried by larger scales than the actual fluxes; this may explain the independence of the model performance from sensor response.

The REA model, on the other hand, would be sensitive to any latency in the response of the mechanical air flow separation system but is not affected by the *w*_{0} dead band. This underlines the importance of the design and construction of REA devices, which remain inaccessible commercially and are rather custom-made by researchers who use them. This prompted us to examine the possibility of using a virtual eddy accumulation (VEA) version where separating the updrafts and downdrafts to obtain the conditional mean concentrations needed for REA is done digitally, only using the sign of *w* and without a need for a mechanical device.

With numerically simulated slow sensors, however, it was noted that the A22 model outperforms VEA in predicting the observed (unfiltered) EC fluxes. VEA was instead reproducing the EC fluxes computed with a filtered scalar signal, which we called the large-eddy-covariance (LEC) fluxes here. This suggests that a VEA approach can plausibly be utilized without a physical device to separate the updraft and downdraft air streams but requires a correction for the missed (filtered) scales. To correct the underestimated VEA fluxes, relative to the observed (unfiltered) EC fluxes, a model for the *β*_{v} coefficients of the VEA model that incorporates the effect of filtering is proposed. The results indicate that the corrected VEA model remains robust in terms of reproducing long-term averages of the scalar fluxes across their ecosystem lifetime cycle but becomes less certain over individual periods as the sensors response time $\stackrel{\mathrm{\u203e}}{\mathrm{\Delta}}$ increases. Use of this model, along with the observations reported in Appendix A that reveal an insensitivity of *β*_{v} to sensor response under unstable regimes (*β*_{v}=*β*_{s}), suggests that VEA can be a robust framework for estimating turbulent fluxes when only single-level measurements are available but with sensors that can still resolve the larger scales of turbulence (i.e., with response frequencies ≈ 0.5 to 1 Hz). The A22 can be alternatively used, under stable conditions, when multilevel measurements are available to compute needed mean gradients, even with slow-response sensors.

The data set of all the observational data for the two field experiments (Utqiaġvik (formerly Barrow) and Wendell) is publicly available at https://doi.org/10.5281/zenodo.10073726 (Allouche et al., 2023c).

MA and EBZ designed the framework of this study; MA, EBZ, ND, EZ, and GGK developed the analysis methodologies; MA performed the formal analysis and visualization; and VIS and MAZ collected and curated the data from Wendell, working with April Leytem from the U.S. Department of Agriculture and John Walker and Ryan Fulgham from the U.S. Environmental Protection Agency. JDF collected and curated the data from Utqiaġvik, working with collaborating Environment Canada scientist Ralf Staebler. MA and EBZ developed the first draft of the paper; all of the authors assisted in interpreting the results and contributed to writing the paper.

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 made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

Mohammad Allouche and Elie Bou-Zeid were supported by the Cooperative Institute for Modeling the Earth System at Princeton University under award no. NA18OAR4320123 from the National Oceanic and Atmospheric Administration and by the US National Science Foundation under award no. AGS 2128345. Vladislav I. Sevostianov was supported by the National Defense Science and Engineering Graduate Fellowship from the U.S. Department of Defense and Army Research Office. Mark A. Zondlo was supported by the U.S. Environmental Protection Agency (grant no. 68HERH21D0006). Gabriel G. Katul was supported by the U.S. National Science Foundation (grant no. NSF-AGS-2028633) and the Department of Energy (grant no. DE-SC0022072). Jose D. Fuentes was supported by the National Science Foundation to complete the PHOXMELT field studies (grant no. PLR-1417914) to collect the data. Mohammad Allouche was also funded, in part, by the Lawrence Livermore National Laboratory. This work was performed under the auspices of the U.S. Department of Energy by the Lawrence Livermore National Laboratory under contract DE-AC52-07NA27344 and LLNL-JRNL-866478.

This paper was edited by Thijs Heus and reviewed by two anonymous referees.

Allouche, M., Katul, G. G., Fuentes, J. D., and Bou-Zeid, E.: Probability law of turbulent kinetic energy in the atmospheric surface layer, Physical Review Fluids, 6, 074601, https://doi.org/10.1103/PhysRevFluids.6.074601, 2021. a

Allouche, M., Bou-Zeid, E., Ansorge, C., Katul, G. G., Chamecki, M., Acevedo, O., Thanekar, S., and Fuentes, J. D.: The Detection, Genesis, and Modeling of Turbulence Intermittency in the Stable Atmospheric Surface Layer, J. Atmos. Sci., 79, 1171–1190, 2022. a, b, c, d

Allouche, M., Bou-Zeid, E., and Iipponen, J.: The influence of synoptic wind on land–sea breezes, Q. J. Roy. Meteor. Soc., 149, 3198–3219, https://doi.org/10.1002/qj.4552, 2023a. a

Allouche, M., Bou-Zeid, E., and Iipponen, J.: Unsteady Land-Sea Breeze Circulations in the Presence of a Synoptic Pressure Forcing, arXiv [preprint], https://doi.org/10.48550/arXiv.2401.00863, 2023b. a

Allouche, M., Sevostianov, V. I., Zahn, E., Zondlo, M., Dias, N. L., Katul, G. G., Fuentes, J. D., and Bou-Zeid, E.: Data Sets: Estimating scalar turbulent fluxes with slow-response sensors in the stable atmospheric boundary layer, Version v1, Zenodo [data set], https://doi.org/10.5281/zenodo.10073726, 2023c. a

Ansorge, C. and Mellado, J. P.: Global intermittency and collapsing turbulence in the stratified planetary boundary layer, Bound.-Lay. Meteorol., 153, 89–116, 2014. a, b

Ansorge, C. and Mellado, J. P.: Analyses of external and global intermittency in the logarithmic layer of Ekman flow, J. Fluid Mech., 805, 611–635, 2016. a, b

Baker, J., Norman, J., and Bland, W.: Field-scale application of flux measurement by conditional sampling, Agr. Forest Meteorol., 62, 31–52, 1992. a

Bowling, D., Turnipseed, A., Delany, A., Baldocchi, D., Greenberg, J., and Monson, R.: The use of relaxed eddy accumulation to measure biosphere-atmosphere exchange of isoprene and other biological trace gases, Oecologia, 116, 306–315, 1998. a

Burba, G.: Eddy Covariance Method for Scientific, Regulatory, and Commercial Applications, LI-COR Biosciences, Lincoln, Nebraska, USA, 2022. a, b

Burba, G., Anderson, T., and Komissarov, A.: Accounting for spectroscopic effects in laser-based open-path eddy covariance flux measurements, Glob. Change Biol., 25, 2189–2202, 2019. a

Businger, J. A. and Oncley, S. P.: Flux measurement with conditional sampling, J. Atmos. Ocean. Tech., 7, 349–352, 1990. a, b, c

Businger, J. A., Wyngaard, J. C., Izumi, Y., and Bradley, E. F.: Flux-profile relationships in the atmospheric surface layer, J. Atmos. Sci., 28, 181–189, 1971. a

Chen, W., Venables, D. S., and Sigrist, M. W.: Advances in Spectroscopic Monitoring of the Atmosphere, Elsevier, Amsterdam, Netherlands, 632 pp., ISBN: 9780128150146, ISBN: 9780128156896, 2021. a

Desjardins, R.: Description and evaluation of a sensible heat flux detector, Bound.-Lay. Meteorol., 11, 147–154, 1977. a

Detto, M. and Katul, G.: Simplified expressions for adjusting higher-order turbulent statistics obtained from open path gas analyzers, Bound.-Lay. Meteorol., 122, 205–216, 2007. a

Edson, J., Fairall, C., Bariteau, L., Zappa, C. J., Cifuentes-Lorenzen, A., McGillis, W. R., Pezoa, S., Hare, J., and Helmig, D.: Direct covariance measurement of CO_{2} gas transfer velocity during the 2008 Southern Ocean Gas Exchange Experiment: Wind speed dependency, J. Geophys. Res.-Oceans, 116, C00F10, https://doi.org/10.1029/2011JC007022, 2011. a

Ellison, T.: Turbulent transport of heat and momentum from an infinite rough plane, J. Fluid Mech., 2, 456–466, 1957. a

Fogarty, J. and Bou-Zeid, E.: The Atmospheric Boundary Layer Above the Marginal Ice Zone: Scaling, Surface Fluxes, and Secondary Circulations, Bound.-Lay. Meteorol., 189, 53–76, 2023. a

Foken, T.: 50 years of the Monin–Obukhov similarity theory, Bound.-Lay. Meteorol., 119, 431–447, 2006. a

Gao, W.: The vertical change of coefficient *b*, used in the relaxed eddy accumulation method for flux measurement above and within a forest canopy, Atmos. Environ., 29, 2339–2347, 1995. a

Grelle, A. and Keck, H.: Affordable relaxed eddy accumulation system to measure fluxes of H_{2}O, CO_{2}, CH_{4} and N_{2}O from ecosystems, Agr. Forest Meteorol., 307, 108514, 2021. a

Hodgkinson, J. and Tatam, R. P.: Optical gas sensing: a review, Meas. Sci. Technol., 24, 012004, https://doi.org/10.1088/0957-0233/24/1/012004, 2012. a

Horst, T.: A simple formula for attenuation of eddy fluxes measured with first-order-response scalar sensors, Bound.-Lay. Meteorol., 82, 219–233, 1997. a

Issaev, V., Williamson, N., and Armfield, S.: Intermittency and critical mixing in internally heated stratified channel flow, J. Fluid Mech., 963, A5, https://doi.org/10.1017/jfm.2023.303, 2023. a

Karl, T. G., Spirig, C., Rinne, J., Stroud, C., Prevost, P., Greenberg, J., Fall, R., and Guenther, A.: Virtual disjunct eddy covariance measurements of organic compound fluxes from a subalpine forest using proton transfer reaction mass spectrometry, Atmos. Chem. Phys., 2, 279–291, https://doi.org/10.5194/acp-2-279-2002, 2002. a

Katul, G., Albertson, J., Chu, C.-R., Parlange, M., Stricker, H., and Tyler, S.: Sensible and latent heat flux predictions using conditional sampling methods, Water Resour. Res., 30, 3053–3059, 1994. a

Katul, G., Peltola, O., Grönholm, T., Launiainen, S., Mammarella, I., and Vesala, T.: Ejective and sweeping motions above a peatland and their role in relaxed-eddy-accumulation measurements and turbulent transport modelling, Bound.-Lay. Meteorol., 169, 163–184, 2018. a, b

Katul, G. G., Finkelstein, P. L., Clarke, J. F., and Ellestad, T. G.: An investigation of the conditional sampling method used to estimate fluxes of active, reactive, and passive scalars, J. Appl. Meteorol. Clim., 35, 1835–1845, 1996. a, b, c

Katul, G. G., Porporato, A., Shah, S., and Bou-Zeid, E.: Two phenomenological constants explain similarity laws in stably stratified turbulence, Phys. Rev. E, 89, 023007, https://doi.org/10.1103/PhysRevE.89.023007, 2014. a

Kolb, C. E., Cox, R. A., Abbatt, J. P. D., Ammann, M., Davis, E. J., Donaldson, D. J., Garrett, B. C., George, C., Griffiths, P. T., Hanson, D. R., Kulmala, M., McFiggans, G., Pöschl, U., Riipinen, I., Rossi, M. J., Rudich, Y., Wagner, P. E., Winkler, P. M., Worsnop, D. R., and O' Dowd, C. D.: An overview of current issues in the uptake of atmospheric trace gases by aerosols and clouds, Atmos. Chem. Phys., 10, 10561–10605, https://doi.org/10.5194/acp-10-10561-2010, 2010. a

Mahrt, L.: Nocturnal boundary-layer regimes, Bound.-Lay. Meteorol., 88, 255–278, 1998. a

Mahrt, L. and Bou-Zeid, E.: Non-stationary boundary layers, Bound.-Lay. Meteorol., 177, 189–204, 2020. a

Miller, D. J., Sun, K., Tao, L., Khan, M. A., and Zondlo, M. A.: Open-path, quantum cascade-laser-based sensor for high-resolution atmospheric ammonia measurements, Atmos. Meas. Tech., 7, 81–93, https://doi.org/10.5194/amt-7-81-2014, 2014. a

Milne, R., Beverland, I., Hargreaves, K., and Moncrieff, J.: Variation of the *β* coefficient in the relaxed eddy accumulation method, Bound.-Lay. Meteorol., 93, 211–225, 1999. a

Monin, A. and Obukhov, A.: Basic laws of turbulent mixing in the surface layer of the atmosphere, Dokl. Akad. Nauk SSSR+, 24, 163–187, 1954. a, b

Nie, D., Kleindienst, T., Arnts, R., and Sickles, J.: The design and testing of a relaxed eddy accumulation system, J. Geophys. Res.-Atmos., 100, 11415–11423, 1995. a

Obukhov, A.: Turbulence in an atmosphere with a non-uniform temperature, Bound.-Lay. Meteorol., 2, 7–29, 1971. a

Pan, D., Benedict, K. B., Golston, L. M., Wang, R., Collett Jr., J. L., Tao, L., Sun, K., Guo, X., Ham, J., Prenni, A. J., Schichtel, B. A., Mikoviny, T., Müller, M., Wisthaler, A., and Zondlo, M. A.: Ammonia dry deposition in an alpine ecosystem traced to agricultural emission hotpots, Environ. Sci. Technol., 7776–7785, 2021. a

Perrie, W., Long, Z., Hung, H., Cole, A., Steffen, A., Dastoor, A., Durnford, D., Ma, J., Bottenheim, J. W., Netcheva, S., Staebler, R., Drummond, J. R., and O’Neill, N. T.: Selected topics in arctic atmosphere and climate, Climatic Change, 115, 35–58, 2012. a

Qian, Y., Gustafson Jr., W. I., and Fast, J. D.: An investigation of the sub-grid variability of trace gases and aerosols for global climate modeling, Atmos. Chem. Phys., 10, 6917–6946, https://doi.org/10.5194/acp-10-6917-2010, 2010. a

Rinne, H., Guenther, A., Warneke, C., De Gouw, J., and Luxembourg, S.: Disjunct eddy covariance technique for trace gas flux measurements, Geophys. Res. Lett., 28, 3139–3142, 2001. a

Rinne, J., Douffet, T., Prigent, Y., and Durand, P.: Field comparison of disjunct and conventional eddy covariance techniques for trace gas flux measurements, Environ. Pollut., 152, 630–635, 2008. a

Ruppert, J., Wichura, B., Delany, A., and Foken, T.: 2.8 Eddy sampling methods, A comparison using simulation results, in: 15th Symposium on Boundary Layers and Turbulence, Wageningen, The Netherlands, 15–19 July 2002, American Meteorological Society, 15, p. 27, 2002. a

Shah, S. and Bou-Zeid, E.: Direct numerical simulations of turbulent Ekman layers with increasing static stability: modifications to the bulk structure and second-order statistics, J. Fluid Mech., 760, 494–539, https://doi.org/10.1017/jfm.2014.597, 2014. a, b

Shah, S. and Bou-Zeid, E.: Rate of decay of turbulent kinetic energy in abruptly stabilized Ekman boundary layers, Physical Review Fluids, 4, 074602, https://doi.org/10.1103/PhysRevFluids.4.074602, 2019. a

Sharma, S., Chan, E., Ishizawa, M., Toom-Sauntry, D., Gong, S. L., Li, S., Tarasick, D., Leaitch, W., Norman, A., Quinn, P., Bates, T. S., Levasseur, M., Barrie, L. A., and Maenhaut, W.: Influence of transport and ocean ice extent on biogenic aerosol sulfur in the Arctic atmosphere, J. Geophys. Res.-Atmos., 117, D12209, https://doi.org/10.1029/2011JD017074, 2012. a

Staebler, R., Fuentes, J., and Bottenheim, J.: The role of surface and boundary layer dynamics in Arctic ozone depletion episodes, in: 2009 AGU Fall Meeting, San Francisco, CA, 14–18 December 2009, 2009, A24B–05, 2009. a

Stull, R. B.: Inversion rise model based on penetrative convection, J. Atmos. Sci., 30, 1092–1099, 1973. a

Stull, R. B.: An introduction to boundary layer meteorology, vol. 13, Springer Science & Business Media, https://doi.org/10.1007/978-94-009-3027-8_9, 1988. a

Sun, K., Tao, L., Miller, D. J., Zondlo, M. A., Shonkwiler, K. B., Nash, C., and Ham, J. M.: Open-path eddy covariance measurements of ammonia fluxes from a beef cattle feedlot, Agr. Forest Meteorol., 213, 193–202, 2015. a

Vogl, T., Hrdina, A., and Thomas, C. K.: Choosing an optimal *β* factor for relaxed eddy accumulation applications across vegetated and non-vegetated surfaces, Biogeosciences, 18, 5097–5115, https://doi.org/10.5194/bg-18-5097-2021, 2021. a

Voulgarakis, A., Marlier, M. E., Faluvegi, G., Shindell, D. T., Tsigaridis, K., and Mangeon, S.: Interannual variability of tropospheric trace gases and aerosols: The role of biomass burning emissions, J. Geophys. Res.-Atmos., 120, 7157–7173, 2015. a

Weaver, H. L.: Temperature and humidity flux-variance relations determined by one-dimensional eddy correlation, Bound.-Lay. Meteorol., 53, 77–91, https://doi.org/10.1007/BF00122464, 1990. a

Webb, E. K., Pearman, G. I., and Leuning, R.: Correction of flux measurements for density effects due to heat and water vapour transfer, Q. J. Roy. Meteor. Soc., 106, 85–100, 1980. a

Wilczak, J. M., Oncley, S. P., and Stage, S. A.: Sonic anemometer tilt correction algorithms, Bound.-Lay. Meteorol., 99, 127–150, 2001. a

Wyngaard, J.: Turbulence in the Atmosphere, Cambridge University Press, 215–240, https://doi.org/10.1017/CBO9780511840524.011, 2010. a

Zahn, E., Dias, N. L., Araújo, A., Sá, L. D. A., Sörgel, M., Trebs, I., Wolff, S., and Manzi, A.: Scalar turbulent behavior in the roughness sublayer of an Amazonian forest, Atmos. Chem. Phys., 16, 11349–11366, https://doi.org/10.5194/acp-16-11349-2016, 2016. a

Zahn, E., Bou-Zeid, E., and Dias, N. L.: Relaxed Eddy Accumulation Outperforms Monin–Obukhov Flux Models Under Non-Ideal Conditions, Geophys. Res. Lett., 50, e2023GL103099, https://doi.org/10.1029/2023GL103099, 2023. a, b, c, d

Zeman, O. and Tennekes, H.: Parameterization of the turbulent energy budget at the top of the daytime atmospheric boundary layer, J. Atmos. Sci., 34, 111–123, 1977. a

Zilitinkevich, S., Elperin, T., Kleeorin, N., Rogachevskii, I., and Esau, I.: A hierarchy of energy-and flux-budget (EFB) turbulence closure models for stably-stratified geophysical flows, Bound.-Lay. Meteorol., 146, 341–373, 2013. a