Modeling the diurnal cycle of conserved and reactive species in the convective boundary layer

We have developed a one-dimensional secondorder closure numerical model to study the vertical turbulent transport of trace reactive species in the convective (daytime) planetary boundary layer (CBL), which we call the SecondOrder Model for Conserved and Reactive Unsteady Scalars (SOMCRUS). The temporal variation of the CBL depth is calculated using a simple mixed-layer model with a constant entrainment coefficient and zero-order discontinuity at the CBL top. We then calculate time-varying continuous profiles of mean concentrations and vertical turbulent fluxes, variances, and covariances of both conserved and chemically reactive scalars in a diurnally varying CBL. The set of reactive species is the O3–NO–NO2 triad. The results for both conserved and reactive species are compared with large-eddy simulations (LES) for the same free-convection case using the same boundary and initial conditions. For the conserved species, we compare three cases with different combinations of surface fluxes, and CBL and free-troposphere concentrations. We find good agreement of SOMCRUS with LES for the mean concentrations and fluxes of both conserved and reactive species except near the CBL top, where SOMCRUS predicts a somewhat shallower depth, and has sharp transitions in both the mean and turbulence variables, in contrast to more smeared-out variations in the LES due to horizontal averaging. Furthermore, SOMCRUS generally underestimates the variances and species–species covariances. SOMCRUS predicts temperature–species covariances similar to LES near the surface, but much smaller magnitude peak values near the CBL top, and a change in sign of the covariances very near the CBL top, while the LES predicts a change in sign of the covariances in the lower half of the CBL. SOMCRUS is also able to estimate the intensity of segregation (the ratio of the species–species covariance to the product of their means), which can alter the rates of second-order chemical reactions; however, for the case considered here, this effect is small. The simplicity and extensibility of SOMCRUS means that it can be utilized for a broad range of turbulencemixing scenarios and sets of chemical reactions in the planetary boundary layer; it therefore holds great promise as a tool to incorporate these processes within air quality and climate models.


Introduction
The behavior of trace reactive species in the convective boundary layer (CBL) is of considerable interest for determining the fate of substances emitted by biogenic and anthropogenic sources or entrained into the CBL from the overlying free troposphere (FT).These species may react photochemically or with other species and may be aerosol precursors.If their reaction time constants are between about 0.1 and 10 times the mixing time of the CBL, which we estimate as τ (t) = h/w * , where h(t) is the CBL depth and w * (t) is the convective velocity scale, the species mean and flux profiles may be significantly modified from conserved species profiles.In Eq. (1), g is gravity, T is the mean CBL temperature, and wθ 0 is the surface virtual potential temperature flux.Typical mid-day CBL values are h ≈ 1 km and w * ≈ 1 m s −1 ; thus τ ≈ 1000 s.

D. H. Lenschow et al.: Conserved and reactive species
In order to model the behavior of reactive species correctly, it is important to model both their vertical transport and effective reaction rates since the coupling between the turbulence and the chemistry can have significant impacts on the effective reaction rates and thus on the profiles of these trace species and their products, many of which are important for air quality and climate considerations.One example is the fate of O 3 in the CBL in the presence of other reactive species such as NO and NO 2 .Another is volatile organic compounds emitted by vegetation that react with OH and other oxidants.The interactions among these species are affected by the turbulence in the CBL so that, for example, their flux-gradient relationships are different than for conserved species.Yet, regional air quality and global climate models currently do not take into account these effects even though they may affect the predicted species concentrations.
The effects of chemical reactivity on mean and turbulence statistics of species in the CBL have been investigated previously both with models and observations.An early effort by Lenschow (1982) showed the potential importance of chemical reactivity for the O 3 -NO-NO 2 triad in the surface layer of the CBL.This was followed by a more quantitative analysis of this triad in the surface layer by Fitzjarrald and Lenschow (1983) and an analytical study by Lenschow and Delany (1986).More detailed numerical studies in the surface layer were carried out by Gao et al. (1991), Vilà-Guerau de Arellano and Duynkerke (1992), Vilà-Guerau de Arellano et al. (1995), and Galmarini et al. (1997b).Galmarini et al. (1997a) also used a one-dimensional second-order closure model to study nitrogen oxide chemistry in the nocturnal boundary layer.Donaldson and Hilst (1972) pointed out that locally inhomogeneous mixing of species involved in second-order reactions, as measured by the intensity of segregation (the ratio of the species-species covariance to the product of their means), can change (generally decrease) their reaction rates.Schumann (1989) extended consideration of chemical reactivity effects for two reacting species -one emitted at the surface and the other entrained across the CBL top -to the entire CBL using large-eddy simulation (LES) as a tool and quantified the relationship between the effective reaction rate and intensity of segregation.Sykes et al. (1994) used LES to further study the effects of turbulent mixing on the effective reaction rate between two species, and also compared LES results with a second-order turbulence model using several closures for the triple correlation terms.Krol et al. (2000) used LES with a more detailed chemical scheme that included OH, HO 2 , and a generic hydrocarbon RH in addition to the O 3 -NO-NO 2 triad and obtained a significant reduction in the RH reaction rate in the CBL due to segregation effects, and also showed that nonuniform surface fluxes of RH further slowed its reaction rate.Kim et al. (2012) showed, via LES, that both fair-weather cumulus and the concentration of NO + NO 2 can further modify the reaction rate of isoprene and the O 3 concentration.Vinuesa and Vilà-Guerau de Arel-lano (2003) used LES to elicit more details of terms in the covariance budgets of chemically reactive species and proposed a parameterization for the intensity of segregation of reactive species.
Here we report on continued development of a secondorder closure model of the CBL.The immediate origins of the model -which we call the Second-Order Model for Conserved and Reactive Unsteady Scalars (SOMCRUS) -go back to Verver et al. (1997Verver et al. ( , 2000)), who developed a secondorder closure model to investigate reactive species in the CBL.This work by Verver et al. (1997Verver et al. ( , 2000) ) was subsequently used by Kristensen et al. (2010) as a basis for a simple, one-dimensional second-order closure model to obtain continuous equilibrium profiles of turbulent fluxes and mean concentrations of non-conserved scalars (the O 3 -NO-NO 2 triad) in a steady-state convective boundary layer without shear.The development here combines a simple mixed-layer model (Tennekes, 1973) of the diurnally varying CBL from which we obtain the depth h(t), the mean virtual potential temperature , and the virtual potential temperature difference across the assumed infinitesimally thin CBL top with a second-order model of the turbulence and mean CBL structure for both conserved and reactive species with surface sources and sinks, and turbulent entrainment of FT air across the top of the CBL.SOMCRUS differs from Verver et al. (1997Verver et al. ( , 2000) ) in that it: (1) explicitly calculates h(t) rather than using a prescribed h(t), and (2) does not include parameterized diagnostic equations for the third moments that appear in the second-moment equations.We found that not including the third-moment equations significantly simplified setting up and running the model while not greatly impacting the results.
Here we model a shear-free CBL and use free-convection surface-layer scaling, but our scheme can easily be modified to run other parameterized boundary layers (e.g., incorporating shear and canopy structure).We then apply SOMCRUS first to a conserved species with differing surface and entrainment fluxes, and second to the O 3 -NO-NO 2 triad, and compare the results with LES.

Basic equations
SOMCRUS is a further development of the model of Kristensen et al. (2010) who carried out similar studies using a second-order closure model to calculate profiles of mean and turbulence statistics, but they considered only steady-state solutions (dh / dt = 0), with the entrainment rate of FT air into the CBL balanced by a mean subsidence velocity.
Here we extend the model of Kristensen et al. (2010)  greatly throughout the day, starting near the surface early in the morning and increasing to a typical depth of a kilometer or more by mid-afternoon.We first solve for h(t), the mean mixed-layer virtual potential temperature (t), and the virtual potential temperature across the inversion at the top of the CBL (t) = h (t) − (t) simultaneously using the mixed-layer approach developed by Tennekes (1973), where γ = ∂ /∂z is the FT lapse rate, θ denotes fluctuations in virtual potential temperature, ∂W/∂z is the largescale CBL subsidence, and is the negative ratio of the virtual potential temperature flux at h to the surface temperature flux.We use the computed h(t) as an input into SOMCRUS.SOMCRUS is a coupled second-order moment system for mean concentrations S i (z, t), fluxes ws i (z, t), speciestemperature covariances θ s i (z, t), and species-species covariances s i s j (z, t) where angle brackets • • • indicate ensemble averaging, which here can be interpreted as averaging over a large enough horizontal domain to obtain stable statistics.The moment equations have the general form of time change + vertical transport + mixing = chemical reaction moments.The relevant equations for this analysis follow Kristensen et al. (2010) and Verver et al. (1997).
The first equation is the mass conservation equation for the concentration of scalars s i (z, t), where s i (z, t) is decomposed into a mean and fluctuation, s i (z, t) = S i (z, t) + s i (z, t), where for simplicity for single variables we use the notation S i = s i .The mean profiles S i (z, t) obey a system of differential equations, Similarly, R i (z, t), which is the rate of concentration change due to reactions with all other species and to photochemistry, is decomposed as where The first-and second-order chemical reaction rates are given by b i j and k i j m , respectively, where the left side contains the reactants and the right side the products: This notation can be extended to higher-order chemical reactions if needed.The reaction rates for a species i are then given by As described in detail by Kristensen et al. (2010), Eq. ( 12) is combined with the three second-moment equations for the flux, temperature-scalar covariance, and scalar-scalar covariance, and to obtain a set of equations that can be solved for the mean and second-order moments.Here we have neglected moments higher than two since Kristensen et al. (2010) found them to be relatively unimportant.Comparing the two systems with and without parameterized third-order moment terms, mathematically the latter is first-order in time and space variables while the former contains second-order derivative terms and requires an additional set of boundary conditions and empirically determined constants.We found, however, that adding the third-moment diagnostic expressions given by Verver et al. (1997) to the second-moment equations reduces the gradients in the mean concentration profiles and improves somewhat the comparison with LES.

D. H. Lenschow et al.: Conserved and reactive species
The chemical moments on the right-hand side of Eqs. ( 13)-( 15) are Equations ( 11)-( 18) are formulated for first-and secondorder chemical kinetics, but the moment chemistry scheme could be easily extended to other (higher-order) reactions.Following Kristensen et al. (2010), we assume that the mean virtual potential temperature gradient term in Eq. ( 14) is negligible in the CBL.The constants in Eqs. ( 13)-( 15) are obtained as follows.For the pressure-scalar covariance term in Eq. ( 13) we follow André et al. (1976), Moeng and Wyngaard (1986), Moeng andWyngaard (1989), andVerver et al. (1997) and use the parameterization where B 0.4 is a dimensionless constant and τ 1 = τ 1 (z) the "return to isotropy" timescale.This parameterization is based on LES of the CBL, and is widely used in second-order models of the CBL.Likewise, the viscous terms in Eqs. ( 14) and (15) have been parameterized by "return to isotropy" timescales τ 4 (z) and τ 3 (z), respectively: We also use the following parameterized second-order moments: (1) the empirical formulation of Lenschow et al. (1980) where z * = z/ h, and (2) the commonly accepted empirical formulation (e.g., Tennekes, 1973) for wθ (z), These expressions result from a combination of both observations and laboratory experiments.The time constants in Eqs. ( 13)-( 15) and Eqs. ( 19)-( 21) are parameterized as where a i are dimensionless constants, κ = 0.4 is the von Kármán constant, and τ TKE is the turbulent kinetic energy timescale in mid-CBL.This is similar to Verver et al. (1997), except that we use 18 instead of 10 as the constant in Eq. ( 24).We do this so that τ TKE ≈ 2.8h/w * in mid-CBL, as suggested by the LES results of Moeng and Wyngaard (1989).This differs from Verver et al. (1997), who assumed that τ TKE ≈ h/w * .Another difference from Verver et al. (1997) is that, as pointed out by Kristensen et al. (2010), the predicted free-convection surface-layer relationship (Holtslag and Moeng, 1991) for the normalized eddy diffusivity given by leads to the relation In order to fulfill this condition, we modify the values of {a 1 , a 4 } = {4.85,2.5} given by Verver et al. (1997) to {7.67, 3.96} so as to both maintain the same ratio a 1 /a 4 as Verver et al. (1997) and fulfill Eq. ( 26).The other two constants used here, {a 3 , B} = {2.5, 0.4}, are the same as in Verver et al. (1997).

Description of LES model
Due to the enormous complexities associated with real-world observations, we turn to turbulence-resolving atmospheric LES as a tool to evaluate the ability of SOMCRUS to simulate the time evolution of passive and reactive scalars in the CBL.The National Center for Atmospheric Research's (NCAR) LES was first described in Moeng (1984) and Moeng and Wyngaard (1988), and was subsequently modified by Sullivan et al. (1994), Sullivan et al. (1996), Patton et al. (2005), Vilà-Guerau De Arellano et al. (2005), Sullivan and Patton (2011) and Kim et al. (2012).Over the years, the NCAR LES has proven its ability to simulate observed atmospheric statistics across a wide variety of atmospheric situations and surface characteristics (e.g., Moeng, 1984;Moeng and Wyngaard, 1988;Sullivan et al., 1996;Patton et al., 2003;Vilà-Guerau De Arellano et al., 2005;Beare et al., 2006;Finnigan et al., 2009;Sullivan and Patton, 2011;Lenschow et al., 2012) and has therefore become a close counterpart to field campaigns.Since most of the LES code has been previously described, we present here only a limited discussion of the current code.
The NCAR LES code integrates a set of threedimensional, wave-cutoff-filtered Boussinesq equations, where a Poisson equation solves for the pressure.In the work described here, a thermodynamic energy equation as well as a conservation equation for each of three passive scalars and three reactive scalars are solved.Unresolved, or subfilterscale (SFS) processes, are accounted for by using Deardorff's (1980) 1.5-order TKE model.Reactive scalars are presumed to mix like passive scalars at scales smaller than the filter width.
Horizontal derivatives are estimated using pseudospectral methods (Fox and Orzag, 1973), and vertical derivatives use a second-order centered-in-space finite difference scheme for velocity fields and Koren's (1993) method for all scalar fields.A third-order Runge-Kutta scheme advances the solutions in time (Spalart et al., 1991;Sullivan et al., 1996).
The simulations use 256 × 256 × 256 grid points to resolve a 5.12 × 5.12 × 2.56 km 3 domain.Therefore, the grid resolution is (20, 20, 10) m in the (x, y, z) directions, respectively.Periodic boundary conditions are imposed in the horizontal directions.Klemp and Durran's (1983) radiation boundary condition handles the upper boundary conditions.No-slip conditions are enforced at the ground surface, where the surface stress is calculated following Monin-Obukhov Similarity Theory (MOST) from a prescribed surface roughness length and the velocity or scalar mixing ratio at one-half grid point above the surface, where no modification to MOST is imposed for reactive scalars.
Turbulent fluctuations from the LES are calculated as deviations from the horizontal mean.Turbulence moments are then determined as horizontally averaged fluctuation products which are then time-averaged using a time-evolving vertical coordinate system according to the time-evolving CBL depth.The CBL depth h is estimated using the LES fields as the height of the minimum buoyancy flux.

Implementation of SOMCRUS
The SOMCRUS Eq. ( 6) and Eqs. ( 13)-( 15) contains 3n + n(n + 1)/2 partial differential equations for the following variables: mean concentrations, S i (z, t); vertical eddy fluxes, ws i ; temperature-species covariances, θ s i ; and speciesspecies variances and covariances, s i s j , where 1 ≤ i ≤ j ≤ n, and n is the total number of species.The combined PDE system is configured so that it can be solved in a space-time region consisting of a full or partial diurnal cycle, t 0 < t < t 1 , where t 0 is the initial time (e.g., sunrise, or earlier), and t 1 is the final time (e.g., sunset) with time-dependent spatial boundaries given by the CBL height, 0 < z < h(t), using the mixed-layer Eqs. ( 2)-( 4).
We need to impose 3n + n(n + 1)/2 boundary conditions (BCs).We impose an entrainment relationship for species fluxes across the CBL top, where w e is the entrainment velocity, S i (h + ) is the concentration just above the CBL top, and S i (h − ) the concentration just below the top.We also specify surface values for the temperature and species fluxes as well as for the species variances and temperature-species covariances.
In general, systems like SOMCRUS with top and bottom BCs are well-posed mathematically, so we would ex-pect a unique well-defined solution throughout the domain {0 < z < h(t)} for the species concentrations and secondorder moments.There are, however, some serious mathematical and numerical problems that can have significant impact on the CBL structure and need to be addressed in the timedependent CBL due to the singular nature of the parameterized functions: namely, at the lower boundary (z * = 0) the parameterized moment w 2 (z * ), the timescales τ i (z * ), and many coefficients (e.g., the eddy diffusivity) vanish.This is a well-established feature of surface-layer dynamics (e.g., Stull, 1988) and has important implications for analysis and solutions of CBL systems that attempt to simulate surfacelayer structure, namely: (1) proper choice and setup of BCs, (2) structure of the solutions, and (3) mathematical and numerical techniques for solving such systems.Verver et al. (1997) did not attempt to deal with this problem and thus did not resolve surface-layer structure in a time-varying (diurnal) model as we do here, which may have significant impact on the overlying CBL structure.In the Appendix we lay out our technique for solving the set of Eqs. ( 13) to (15) in a way that allows us to resolve the surface-layer structure and gives an efficient way to solve the moment equations throughout the CBL.
Our boundary conditions (BCs) are similar to those used by Verver et al. (2000).We specify the surface species fluxes ws i 0 (t); the surface variances and covariances are specified based on relations obtained by Wyngaard et al. (1971) from observations in the free-convection regime: s i s j 0 = 1.66 At the lower boundary z = z 0 (z 0 / h is set equal to 10 −3 for numerical calculations; note that z 0 is not the roughness length but a lower boundary condition for solving the differential equation set Eqs. (A9)-(A12), as we assume a free convection boundary layer).Similarly, because of the discontinuity at the top boundary (z = h), which causes numerical difficulties, we actually use z = 0.993h in SOMCRUS; henceforth for simplicity, we redefine h as the height used in SOMCRUS.
We use Mathematica (Wolfram Research, Inc., 2015) at all stages of the model development, implementation, and simulations.The mixed-layer Eqs. ( 2)-(4), are first solved using the Mathematica differential equation solver, and the calculated values for h(t) and (t) are used in SOMCRUS, Eqs.(A9)-(A12).SOMCRUS is designed to cleanly separate the turbulent mixing terms in the moment equations from the chemical reaction terms in the system of Eqs.(A9)-(A12).Mathematica allows us to generate the entire SOM-CRUS system in two steps: (1) using symbolic algebra tools we generate from the basic chemical suite of species and reactions the complete moment chemistry; (2) parameter-ized CBL mixing along with the mixed-layer solution for {h(t), (t)} allows us to generate the turbulent mixing part of the system in regularized form, Eqs.(A9)-(A12).
The next step is to solve Eqs.(A9)-(A12) with the given boundary conditions.The Mathematica solver does this by a proper spatial discretization scheme whose inputs (resolution, difference order, etc.) can be controlled.Thereby a system of partial differential equations is converted into a large (coupled) set of ordinary differential equations solved by time-adaptive numeric codes.The output of the Mathematica solver is a set of interpolating functions over a prescribed space-time range.A single run for a conserved species with a spatial resolution of 100 points in x takes about 30 s of desktop computing time.A system of three reactive species -i.e., the O 3 -NO-NO 2 triad (15 equations) -at the same resolution takes 100-200 s of desktop computing time, depending on the spatial and temporal resolution used in solving the equations.The system size increases with the number of reactive species; e.g., for 10 reactive species, 85 equations must be solved.

Case description
In order to demonstrate the performance of SOMCRUS, we compare SOMCRUS results with those from LES using the same meteorological case as Vilà-Guerau de Arellano et al. ( 2011); namely, 15-day averaged observations from the Tropical Forest and Fire Emission Experiment (TROFFEE, Karl et al., 2007).The initial and boundary conditions in the numerical experiments are presented in Tables 1 and 2. The geostrophic wind is 0 m s −1 (i.e., local free convective conditions).No large-scale forcings (i.e., no horizontal heat and moisture advection, subsidence, nor radiative tendencies) are prescribed.Turbulence is initiated in the LES by imposing a divergence-free random perturbation field on the velocity and temperature fields in the lowest 200 m.The LES results presented in Figs.2-11 represent 1-hour averages centered at the depicted times.The simulation begins at 05:00 local time (LT) and lasts 13 h (sunrise is at 06:00 LT and sunset at 18:00 LT).The depth of the CBL calculated by SOMCRUS and the surface temperature flux are shown in Fig. 1.

Conserved species means and moments
We first compare the mean and moment profiles for three cases of a conserved scalar using both SOMCRUS and LES at 10:00 LT, 12:00 LT, and 14:00 LT (see Table 1 for the meteorological initial and boundary conditions of the variables).Each scalar case (labeled "case A", "case B", and "case C") has different initial conditions and BCs as specified in Table 2.We present these three conserved scalar cases to demonstrate the ability of SOMCRUS to reproduce vertical mixing in the CBL and the influence of surface or entrainment fluxes in the absence of reactivity.
Profiles for case A, which has a surface flux and an initial CBL concentration, but zero concentration in the FT are compared in Fig. 2.This case illustrates the effects of both a surface source and entrainment on the evolving CBL, but since the FT concentration is zero, the total mass of species within the CBL (i.e., the area under the curve) is not affected by entrainment and is the same for both SOMCRUS and LES.We see that particularly at 10:00 LT the concentration distribution around the CBL top is more spread out vertically in the LES than for SOMCRUS, which has a step change in concentration at the CBL top.This smearing out is because the LES resolves horizontal variations in the CBL structurein particular, horizontal variations in the CBL top.The LES also predicts a CBL depth about 150 m higher than SOM-CRUS, which is consistent with the results of Vilà-Guerau de Arellano et al. ( 2011), who used a similar mixed-layer model and made similar comparisons of h with LES for the same case as here.These two features result in a SOMCRUS CBL concentration that is larger than the LES concentration.Furthermore, the LES predicts a smaller gradient throughout the CBL, which increases the difference between the two concentration profiles near the surface as compared to the upper part of the CBL.The maximum difference of about 12 % occurs at 10:00 LT at z * ≈ 0.06.Later, at 12:00 and 14:00 LT these differences, although still present, are less pronounced and thus the agreement between SOMCRUS and LES is improved.
Comparing the vertical flux profiles in Fig. 2 for case A at the same three times, we see that the 10:00 LT LES flux is more spread out vertically, analogous to the concentration, and extends to a higher level than the SOMCRUS flux, with the difference increasing with height up to h.This results in about a 12 % larger flux maximum for SOMCRUS than for the LES.At later times, the LES and SOMCRUS fluxes are in very good agreement, except near the top where the LES flux is again more spread out.The right column of Fig. 2 shows a Table 1.Initial and prescribed values used for SOMCRUS and the LES numerical experiments.The temperature and humidity surface fluxes, and mean profiles are obtained from a simple curve fit to observations from the Tropical Forest and Fire Emission Experiment (TROFFEE), which is the same meteorological case used by Vilà-Guerau de Arellano et al. (2011); see also Karl et al. (2007).All initial conditions are imposed at 05:00 LT, and t is time in seconds.The subscripts ( ) 0 and ( ) h refer to the surface and CBL top, respectively.

Property Value
Initial CBL height, h (m) 200 Surface virtual potential temperature flux (K m s −1 ) wθ 0 = 0.19 sin π(t−8100) comparison of SOMCRUS variances with LES variances for case A. We see that the LES predicts the height of the variance maximum near the CBL top to be about 150 m higher than SOMCRUS, consistent with the predicted higher LES mixed-layer depth.The LES maximum variance is slightly larger than SOMCRUS at 10:00 LT and subsequently decreases more slowly than SOMCRUS so that by 14:00 LT the SOMCRUS variance is only about 17 % of the LES variance.This is likely occurring because the SOMCRUS variance depends explicitly on the CBL growth rate and the jump in concentration across the CBL top, while the LES variance, being a horizontal average, also incorporates contributions from the horizontal variations in CBL height, which are not included in the SOMCRUS results.The SOMCRUS variance is also strongly dependent on the value of a 3 , but adjusting a 3 does not address the more rapid decrease in SOMCRUS variance with time compared with LES; furthermore, decreasing a 3 to obtain a better match to the LES variance near the CBL top also increases the SOMCRUS variance near the surface, which then worsens the comparison of SOMCRUS variance with the LES variance.
Figure 3 shows the variance of the same case A of Fig. 2 at 10:00 LT for the lowest 100 m of the CBL.Here we compare the variance with both the LES and with the local free-convection prediction originally presented by Wyngaard et al. (1971) using dimensional analysis and observational results for temperature variance; later Lenschow et al. (1980) found that this relation, given below, also worked well for humidity variance observations: where s * = ws 0 /w * .Note that the dependency on h cancels out, and we have We see that the SOMCRUS variance agrees well with the LES prediction to within about 40 m of the surface, while the  LES does not capture the z −2/3 dependency close to the surface.We note that Sullivan and Patton (2011) have pointed out that it may be possible for the LES to reproduce this additional near-surface scalar variance if an additional equation for subfilter-scale scalar variance were incorporated akin to that used by Schmidt and Schumann (1989) -a feature not yet implemented in the NCAR LES.The SOMCRUS variance profile has a shape similar to that of the free-convection prediction, but is systematically larger by about 0.2 units 2 .Figure 4 shows the same set of profiles for case B, which has no initial CBL concentration, 6 units FT concentration, and 1 unit m s −1 surface flux.The results are very similar to case A; the combination of surface flux and entrainment results in a CBL concentration remarkably close to case A. Again at 10:00 LT the SOMCRUS concentration is larger than the LES concentration throughout the CBL, with the difference decreasing towards the CBL top, and the LES concentration exceeding the SOMCRUS concentration in the entrainment region near the CBL top.At 12:00 and 14:00 LT, the concentrations are in very good agreement, with the SOMCRUS concentrations slightly exceeding the LES con- centrations near the surface because of a smaller vertical gradient in the LES concentrations.
Comparisons for nonreactive scalar case C at 10:00, 12:00, and 14:00 LT are presented in Fig. 5.This case has no surface flux nor CBL concentration, but an initial FT concentration of 10 units, so it illustrates the effects solely of entrainment on the CBL vertical structure.Here we see almost perfect agreement between the LES and SOMCRUS concentrations, except near the top where the LES variables are again more spread out.The comparison of SOMCRUS variances with LES variances shows that the variance near the CBL top is similar to case A in that the SOMCRUS variance decreases more rapidly with time than the LES variance.In the lowest 200 m of the CBL the SOMCRUS variance becomes negligible since it depends on the surface flux, while the LES variance, particularly at 10:00 LT, is still about 10 % of the maximum variance near the CBL top.Thus, for the LES, variance generated by the entrainment flux is transported all the way down to the surface.
Overall we see from this comparison that the SOMCRUS and LES are in generally good agreement for concentrations and fluxes, especially at the later times when the differences in the entrainment process, which are most apparent at 10:00 LT, have less effect on the overall vertical structure because of the increased CBL depth.However, SOMCRUS significantly underestimates the variances near the CBL top -especially at later times.We also note that SOMCRUS can reproduce the Wyngaard et al. (1971) free-convection prediction for the z −2/3 dependency of scalar variance down to very near the surface.

O 3 -NO-NO 2 means and moments
We now consider the effects of chemical reactivity on the mean and moment profiles for the O 3 -NO-NO 2 triad.The reaction rates are given in Table 3 and the initial conditions in Table 2.These reactions are fast enough (on the order of a hundred seconds around mid-day, increasing at low sun angles) that the reaction time is comparable to the turbulence timescale, h/w * early in the day.The LES surface O 3 flux is specified as a deposition velocity (0.0025 m s −1 ) times the resolved O 3 concentration at the lowest grid level, which for scalars is 5 m above the surface.It is not straightforward to apply this boundary condition directly in SOMCRUS, although it can be done by extrapolating the 5 m O 3 SOM-CRUS concentration down to the lowest level used in the SOMCRUS formulation (z 0 / h = 10 −3 ).Therefore, to ensure as direct a comparison as possible with the LES, we impose a boundary condition for O 3 flux in SOMCRUS that arises via a 30th-order polynomial fit to the time evolution of the horizontally averaged O 3 surface flux predicted by the LES, as shown in Fig. 6.The mean concentrations for all three species at 10:00, 12:00, and 14:00 LT are shown in Fig. 7.We see that the agreement between SOMCRUS and LES is very good for O 3 , again subject to the effects of a smaller CBL depth h for SOMCRUS compared to that predicted by LES, but for NO + NO 2 -i.e., for the total odd nitrogen which is conserved -the LES predicts a higher concentration than SOMCRUS.This is because the LES imposes a rough-wall stabilitycorrected boundary condition that treats reactive scalars as passive; that is, no reactivity is permitted between the surface and the first grid point in the domain.As a result, for reactive species such as NO, NO 2 , and O 3 during daytime whose reactive timescale is of the order of a minute or two, the LES domain produces a surface flux, in this case an NO surface flux, that appears slightly larger than that imposed.
The LES also predicts a larger vertical gradient for NO than SOMCRUS for 12:00 and 14:00 LT.This is somewhat puzzling since NO should be in approximate chemical equilibrium throughout most of the mixed layer, but with positive surface and entrainment fluxes.
Figure 8 shows a comparison of SOMCRUS species flux profiles in the CBL (blue lines) with LES predictions (red lines) for the O 3 -NO-NO 2 triad.The SOMCRUS produces the non-linearity in the vertical flux profiles resulting from the chemical reactions, similar to the LES.We also note the effects of the greater vertical spread over which the entrainment processes occur in the LES similar to what was observed for the conserved scalar cases.Both models produce about the same curvature in the lower half of the CBL, and because NO + NO 2 is conserved, the sum of the NO and NO 2 fluxes is a straight line.
A comparison of the θ s i covariance profiles at 12:00 LT in Fig. 9 shows that near the surface, the LES and SOMCRUS profiles are very similar.Since the surface flux of ozone is negative and the temperature flux positive, θ s i is negative; the NO flux is positive at the surface and the NO 2 flux is positive just above the surface (due to chemical reaction), thus θ NO1 and θ NO 2 are both positive near the surface.The SOMCRUS covariances decrease in magnitude throughout the mixed layer and change sign near the CBL top, while the LES covariances change sign about midway up, with a large positive θ O 3 peak at the CBL top because of the positive jumps in both and O 3 across the top, and large negative peaks in both θ NO and θ NO 2 because of the negative jumps in NO and NO 2 across the top.The SOMCRUS peaks behave similarly, but with much smaller peak magnitudes.2. Top panel is at 10:00, the middle panel at 12:00, and the bottom panel at 14:00 LT.
We note that in the θ s i covariance equations, the generation term is a sink for θO 3 and a source for θ NO and θ NO 2 throughout most of the CBL.On the other hand, the result of the SOMCRUS assumption of a zero gradient in virtual potential temperature means that the term is neglected in SOMCRUS, while in the LES, for ∂ /∂z > 0, this is a source for θO 3 , and a sink for θ NO and    2.
Top panel is at 10:00 LT, the middle panel at 12:00 LT, and the bottom panel at 14:00 LT.
θ NO 2 .Thus we conclude that SOMCRUS may have some shortcomings in realistically modeling this process compared to the LES; one possibility to address this may be to incorporate a modeled virtual potential temperature gradient in SOMCRUS.
The species variances are compared in Fig. 10, and we see that the LES variances are consistently larger than the SOMCRUS variances throughout the CBL.Near the surface, the SOMCRUS species variances are negligible, as in the conserved case C (Fig. 5) with no surface flux, because the surface flux for NO 2 is zero, and the O 3 and NO surface fluxes are not large enough to generate variances comparable to those generated by entrainment near the CBL top.On the other hand, the LES is able to transport this entrainment-generated variance down to the surface, particularly at 10:00 LT.
A comparison of the s i s j covariances in Fig. 11 shows that SOMCRUS generates generally smaller species peak covariances in the entrainment region than the LES, and a more rapid decrease with time as the entrainment rate decreases.As with the variance and the θ s i covariances, throughout most of the CBL the SOMCRUS s i s j covariances are considerably smaller than the LES.In the entrainment region, SOMCRUS second moments are generated by the entrainment flux and do not include contributions from the undulating capping inversion that are present in the LES because of horizontal averaging.Covariances of two species involved in a second-order chemical reaction can alter the effective reaction rate since the rate is proportional to the concentration of both species.For O 3 NO , however, the covariance may be significant near the surface, but is not large enough to significantly impact the chemical reaction rate throughout the bulk of the mixed layer.This is because the chemical reaction timescale (of order 100 s) is much less than the mixing timescale h/w * ; but for second-order reactions that may occur on timescales comparable to h/w * , the covariances can significantly affect the reaction rates throughout the CBL (e.g., Schumann, 1989).

Intensity of Segregation
Intensity of segregation, defined as quantifies the change in effective reaction rate resulting from the covariance of two species involved in a second-order chemical reaction.Therefore, for the triad, the covariance O 3 NO can change the effective reaction rate for these two species, according to the relationship given by, e.g., Sykes et al. (1994), Reaction (R2) in Table 3 is first order, and therefore the other two species-species covariances do not affect the reaction rates.
For the triad case modeled here, O 3 NO is relatively small near the surface (Fig. 12) because the surface fluxes of both O 3 and NO are relatively small.Therefore, the turbulence makes little change to the reaction rate near the surface in both the SOMCRUS and LES results, although for SOM-CRUS the O 3 NO intensity of segregation increases negatively very near the surface, as it should for species with surface fluxes of opposite sign.Similarly, the O 3 NO 2 intensity 991 of segregation also shows a negative increase approaching the surface.This results from the negative O 3 flux producing negative fluctuations in NO 2 via chemical reactivity.Similarly, the positive NO flux produces positive NO 2 flux, which produces positive NO NO 2 intensity of segregation near the surface.
The entrainment flux also generates species-species covariances that are transported down to the surface, and here the covariances are relatively large in magnitude so the intensity of segregation also becomes large in magnitude.The Fig. 12 plots are cut off at the top of the SOMCRUSpredicted h -i.e., about 150 m below the LES top -since above about this level, the LES intensities of segregation become ill defined because the mean concentrations of NO and NO 2 are zero in the FT.For this case, at 10:00 LT O 3 NO reduces the reaction rate in both the SOMCRUS and the LES results by as much as 5 % near the entrainment zone.
The effects of the intensity of segregation on the effective chemical reaction rates are not included in, e.g., the boundary-layer parameterizations of the Weather Research and Forecasting model coupled with Chemistry (WRF-Chem, Grell et al., 2005), which is used to simulate the emission, transport, mixing, and chemical transformation of trace gases and aerosols simultaneously with meteorology for investigation of regional-scale air quality, field program analyses, and cloud-scale interactions between clouds and chemistry; nor in the mixed-layer model described by Vilà-Guerau de Arellano et al. ( 2009) which examines the evolution of isoprene in the CBL.We also note that if we were to use a more complete chemical mechanism such as Model for Ozone and Related chemical Tracers, version 4 (MOZART-4, Emmons et al., 2010), the influence of the intensities of segregation might be enhanced/reduced as a result of in situ species production via alternate chemical production.

Eddy diffusivity
The concept of an eddy diffusivity is often used in simplified models involving diffusion in the CBL to parameterize turbulent mixing.We therefore examine one obvious approach to this by applying the equations implemented in SOMCRUS to derive an explicit formula for the eddy-diffusivity function K(z, t) = − ws /(∂S/∂z). (37) For a conserved scalar, using Eqs.( 13) and ( 14) we have For steady-state conditions, ∂ ∂t ws = ∂ ∂t θ s = 0, and Eqs. ( 38) and ( 39) can be solved for ws and θ s : Then the eddy diffusivity is Kristensen et al. ( 2010) considered the stationary case where the CBL depth did not change with time because the buoyancy-driven entrainment rate was balanced by the mean subsidence.In that case, Eqs. ( 40) and ( 41) are exact.Here, however, the time changes are not zero, so there is no reason to expect a priori that the stationary relation Eq. ( 42) correctly describes the dynamic case under consideration.Interestingly, the "quasi-stationary" flux-gradient relation Eq. ( 37) holds consistently at all times t.To demonstrate this, we use as an example a case with the same meteorological conditions as the previous case, but with the following differences in the scalar variable: no initial concentration and a surface flux of ws 0 = 0.05 units m s −1 .We still use the same Mathematica implementation scheme, including the changes in variables.Figure 13 shows that there is little difference between two sets of profiles.
We might expect, therefore, that we could use Eq. ( 42) to calculate the S(z, t) profiles for the dynamic case considered here by solving the eddy-diffusion equation However, unlike SOMCRUS, whose solutions are almost completely independent of z 0 , the eddy-diffusion approach is very sensitive to z 0 because of the singular surface boundary condition, with K(z, t) ∼ O(z 4/3 ).In Fig. 14 we see that the eddy diffusion approximation can capture the behavior of the concentration and flux profiles for this test case, but it requires a high-resolution calculation in Mathematica because this singular surface boundary condition creates a large gradient in the concentration near the surface resulting in a sensitive dependence of computed profiles on surface flux and system discretization.Figure 14 shows that 100-point numerical resolution significantly underestimates both the surface flux and concentration, but that both can be adequately resolved with 1000 point resolution.SOMCRUS, however, is very stable to boundary conditions at the surface because the flux and concentration equations are separate and the flux equation is regular at z = 0, while in the explicit diffusivity formulation, the two equations are linked.Another advantage of SOMCRUS, of course, is that it generates second-order moments and intensity of segregation.Although it may seem more straightforward to use an eddy diffusivity, we point out that this does not save computational time compared to SOMCRUS.

Conclusions
We have extended the model of Kristensen et al. (2010) to treat the behavior of conserved and reactive species in the diurnally varying CBL by using: (1) the Tennekes (1973) mixed-layer model to calculate mixed-layer height, mean virtual potential temperature, and virtual potential temperature jump across the CBL top, and (2) a second-order moment closure model to calculate mean and turbulence statistics of reactive species throughout the daytime.Comparing SOM-CRUS with a turbulence-resolving LES for a free-convection case, we note that SOMCRUS has a discontinuous jump across the CBL top, while horizontal averaging of the LES output smears out the variables across the top.We also found: (1) generally good agreement for concentrations and fluxes of both conserved and reactive species throughout most of the mixed layer, including the curvature in the flux profiles throughout the CBL due to chemical reactions; and (2) SOM-CRUS mostly underpredicts the variances and covariances compared to LES, indicating that the time constants used in the second-moment equations in SOMCRUS for parameterizing the rates of dissipation and return-to-isotropy terms may not be optimal.SOMCRUS is able to model the rapid changes in concentrations, variances, and covariances in the surface layer to within a few meters of the surface, as predicted by free-convection similarity theory.We also show that using an eddy-diffusivity formulation for vertical transport is problematical for a time-varying CBL because of the inherent singularity as the diffusivity goes to zero approaching the surface, which is not an issue for SOMCRUS because the flux and concentration equations are separate and the flux equation is regular at z = 0.Because SOMCRUS includes equations for speciesspecies covariances, it can be used to calculate intensities of segregation which can modify the reaction rates for second-order chemical reactions.Although not very important throughout most of the mixed layer for the case considered here (because of the disparity between the turbulence mixing timescale and the chemical reaction timescale for the O 3 -NO-NO 2 triad), this effect can be significant for other reactive species in the CBL (e.g., Krol et al., 2000).
We have shown that SOMCRUS provides a simple and robust tool for predicting concentration, variance, and flux profiles of trace reactive species in the CBL.SOMCRUS is intermediate in ease of use between simple mixed-layer models (e.g., Vilà-Guerau de Arellano et al., 2009) and large-eddy simulation models.SOMCRUS also provides considerably more detail of the vertical variation of first-and second-order species statistics than a mixed-layer model.Furthermore, it is portable, requires little time to run on a PC or laptop using Mathematica, and it is easy to change and to quickly make runs with different scenarios.
SOMCRUS can easily be extended to include more complicated chemistry, such as schemes involving isoprene and related reactions, and to incorporate parameterizations for different surface boundary conditions and meteorological regimes.Examples of this include a parameterized canopy layer and surface stress.We believe that this tool has possibilities for use in air quality models to more accurately simulate the behavior of reactive species in the CBL.We note that software tools exist to convert Mathematica code to Fortran and C++ (e.g., https://store.wolfram.com/view/app/mathcodef90) and that the SOMCRUS code contains separate turbulent mixing and chemistry modules that could in principle be independently incorporated into a larger-

Table 2 .
Surface moisture flux (g kg −1 m s −1 ) wq 0 = 0.13 sin π(t−3600) Specifications for the conserved tracers and the O 3 -NO-NO 2 triad in the numerical experiments with SOMCRUS and LES.The free-troposphere (FT) concentration is constant in time; the convective boundary layer (CBL) concentration and the height h vary with time.

Figure 3 .
Figure 3.Comparison of SOMCRUS (blue curve) with the local free-convection prediction of Lenschow et al. (1980) (green dashed curve) and with LES (red dots) for conserved scalar case A at 10:00 LT.Each dot denotes a layer-averaged LES value.

Figure 6 .
Figure 6.30th-order least squares polynomial fit to the LES surface flux of O 3 .

Figure 7 .
Figure 7.Comparison of SOMCRUS mean concentrations (blue lines) with LES concentrations (red lines) of O 3 , NO, and NO 2 .Initial and boundary conditions are given in Table2.Top panel is at 10:00, the middle panel at 12:00, and the bottom panel at 14:00 LT.

Figure 8 .
Figure 8.Comparison of SOMCRUS fluxes (blue lines) with LES concentrations (red lines) of O 3 , NO, and NO 2 .Initial and boundary conditions are given in Table2.Top panel is at 10:00 LT, the middle panel at 12:00 LT, and the bottom panel at 14:00 LT.

Figure 9 .
Figure 9.Comparison of SOMCRUS θ -species covariances (blue lines) with LES (red lines) of O 3 , NO, and NO 2 at 12:00 LT.Initial and boundary conditions are given in Table2.Top panel covers the entire CBL, while the bottom panel is up to 1 km to accentuate the region below the CBL top.

Figure 10 .
Figure 10.Comparison of SOMCRUS species variances (blue lines) with LES (red lines) of O 3 , NO, and NO 2 at 10:00, 12:00, and 14:00 LT.Initial and boundary conditions are given in Table2.Top panel is at 10:00 LT, the middle panel at 12:00 LT, and the bottom panel at 14:00 LT.

Figure 13 .
Figure 13.A comparison of the flux-gradient profiles for the dynamic SOMCRUS case considered here (red lines) versus the quasistationary diffusivity K(z, t) derived from the SOMCRUS parameterizations (blue lines).

Table 3 .
The chemical reaction scheme used for the O 3 -NO-NO 2 triad in the numerical experiments with SOMCRUS and LES.χ is the zenith angle.