Modeling the Diurnal Cycle of Conserved and Reactive Species in the Convective Boundary Layer using SOMCRUS

We have developed a one-dimensional second-order 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 Second-Order Model for Conserved and Reactive Unsteady Scalars (SOMCRUS). The temporal variation of the CBL depth is calculated using a simple mixedlayer model with a constant entrainment coefficient and zero-order discontinuity at the CBL top. We 5 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 10 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 under15 estimates 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), 20 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 turbulence mixing 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. 25


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 30 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 35 temperature flux. Typical mid-day CBL values are h ≈ 1 km and w * ≈ 1 m s −1 ; thus τ ≈ 1000 s.
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 40 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 45 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 50 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. 55 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-60 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 Arellano (2003) used LES to elicit more details 70 on 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 second-order 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 second-75 order 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,80 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 85 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 90 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. Here we extend the model of Kristensen et al. (2010) by considering a diurnally-varying h(t), which typically varies 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 105 approach developed by Tennekes (1973), where γ = ∂Θ/∂z is the FT lapse rate, θ denotes fluctuations in virtual potential temperature, ∂W/∂z is the large-scale CBL subsidence, and is the negative ratio of the virtual potential temperature flux at h to the surface temperature flux. We 115 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), species-temperature 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 120 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 125 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 The first-and second-order chemical reaction rates are given by b i j and k i jm , 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 140 As described in detail by Kristensen et al. (2010), Eq. (12) is combined with the three secondmoment equations for the flux, temperature-scalar covariance, and scalar-scalar covariance, 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 150 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 empiricallydetermined constants. We did find, 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.

155
The chemical moments on the right side of Eqs. (13) -(15) are

160
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 and Wyngaard (1989), and Verver et al. (1997) and use the parameterization 165 where B 0.4 is a dimensionless constant and τ 1 = τ 1 (z) the "return to isotropy" time scale. This parameterization is based on large-eddy simulation 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" time scales τ 4 (z) and τ 3 (z), respectively: We also use the following parameterized second-order moments: 1) the empirical formulation of Lenschow et al. (1980) for w 2 where z * = z/h, and 2) the commonly-accepted empirical formulation (e.g., Tennekes, 1973) for These expressions result from a combination of both observations and laboratory experiments.

180
where a i are dimensionless constants, κ = 0.4 is the von Kármán constant, and τ T KE is the turbulent kinetic energy time scale 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 τ T KE ≈ 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 τ T KE ≈ h/w * . Another difference from Verver et al. (1997) is that, as pointed 185 out by Kristensen et al. (2010), the predicted free-convection surface-layer relationship (Holtslag and Moeng, 1991) for the normalized eddy diffusivity given by

Description of LES model
Due to the enormous complexities associated with real-world observations, we turn to turbulenceresolving 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 subse-200 quently modified by Sullivan et al. (1994Sullivan et al. ( , 1996 (2011); 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;205 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 three-dimensional, wave-cutoff-filtered Boussinesq equations, where a Poisson equation solves for the pressure. In the work described here, a thermodynamic 210 energy equation as well as a conservation equation for each of three passive scalars and three reactive scalars are solved. Unresolved, or subfilter-scale (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 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).

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 ; temperaturespecies covariances, θs i ; and species-species 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 235 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.
We need to impose 3n+n(n+1)/2 boundary conditions (BCs), where n is the number of species.
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, 245 so we would expect a unique well-defined solution throughout the domain {0 < z < h(t)} for the species concentrations and second-order 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 time-dependent CBL due to the singular nature of the parameterized functions; namely, at the lower boundary (z * = 0) the parameterized moment w 2 (z * ), the time scales τ i (z * ), 250 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 surface-layer structure, namely: (1) proper choice and setup of BCs, (2) structure of the solutions, and (3) mathematical and numerical techniques for solving such systems.
255 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.

260
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: 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. 270 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 al-275 lows us to generate the entire SOMCRUS 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) parameterized 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 280 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 285 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.

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 1000 LT, 1200 LT, and 1400 LT (see Table 1 for the meteorological initial 305 and boundary conditions of the variables). Each scalar case (labeled "case A", "case B" and "case C") has different initial conditions (IC) and boundary conditions (BC) 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.

310
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 1000 LT the concentration distribution around the 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 1000 LT at z * ≈ 0.06. Later, at 1200 LT and 1400 LT these differences, although still present, are less 325 pronounced and thus the agreement between SOMCRUS and LES is improved.
Comparing the vertical flux profiles in Fig. 2  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. 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:

350
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 ad- Overall we see from this comparison that the SOMCRUS and LES are in generally good agree-380 ment for concentrations and fluxes, especially at the later times when the differences in the entrainment process, which are most apparent at 1000 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 385 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 mean concentrations for all three species at 1000, 1200 and 1400 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 400 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 stability-corrected 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 405 daytime whose reactive time scale is on 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 1200 and 1400 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. in NO and NO 2 across the top. The SOMCRUS peaks behave similarly, but with much smaller peak 1 In order to maintain the convention of using capital letters for chemical species, we change the notation for mean/fluctuation of chemical species so that roman type represents a mean value and italic type represents a fluctuation.
magnitudes. 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 430 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 θ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 435 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)  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, through-445 out 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.

450
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 time scale (of order 100 s) is much less than the mixing time scale h/w * ; but for second-order reactions that may occur on time scales 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 460 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.

465
For the triad case modeled here, O 3 NO is relatively small near the surface (Fig. 12)  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 485 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).

500
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 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 515 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 eddydiffusion 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 cal-525 culation in Mathematica because this singular surface boundary condition creates a large gradient in the concentration near the surface. 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 530 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.

535
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 SOMCRUS with a 540 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) SOMCRUS mostly under pre-545 dicts 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-toisotropy 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 550 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 species-species 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 time scale and the chemical reaction time scale 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, 560 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 largeeddy 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 565 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 570 accurately simulate 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-scale numerical model. SOMCRUS can be obtained in the currently-reported O 3 -NO-NO 2 Mathematica notebook configuration by requesting 575 a copy from lenschow@ucar.edu.

Appendix A: Numeric Implementation and SOMCRUS Solutions in Mathematica
The standard technique for solving singular boundary-value problems known as matched asymptotic expansions (Nayfeh, 2008) calls for approximate "inner" (surface layer) and "outer" solutions, as series expansions whose coefficients are matched in the intervening transitional layer. Our approach 580 here is simpler and more efficient than the matched asymptotic expansion. In the context of free convection in the CBL we use the known asymptotic behavior of the following variables as z → 0 to write them as products of scaling factors and regular functions of z: where S i (z, t), θs i , s i s j , are all now regular functions of z at 0, and fluxes ws i are already regular functions. This singular behavior makes it difficult to implement and run SOMCRUS, even when the singular boundary condition at z = 0 is replaced with a positive value that is regular at Here we propose a regularization scheme for SOMCRUS that allows us to compute solutions more efficiently than was the case for Kristensen et al. (2010), using the standard built-in numeric differential equation solvers of Mathematica. The idea is to change variables (independent z and dependent S i , ws i , θs i , s i s j ) to make the system "regular" (or less singular) using a technique similar to the Method of Strained Coordinates (e.g., Nayfeh, 2008, ch. 3), as an alternative 595 to matched asymptotic expansion. Indeed, the asymptotic form in Eqs. (A1) to (A3) suggests a proper change of variables, as well as the choice of surface boundary conditions for ( ws i , θs i , s i s j ); specifically, we replace z by the dimensionless variable x = [z/h(t)] 2/3 (0 < x < 1), and {S, ws i , θs i , s i s j } by the regularized variableŝ Having a fixed range 0 < x < 1 is also an important feature in the standard Mathematica solvers.
The regularized system of variables Eqs. (A4) -(A7) requires replacement of the standard partial 605 derivatives (∂ t , ∂ z ) in SOMCRUS with differential operators For a conserved scalar, the resulting system of equations takes the form Here w 2 (x, t) = w 2 (z, t) , wθ (x, t) = wθ (z, t), and ws i (x, t) = ws i (z, t); furthermore 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).