#### Abstract

This study presents 1D analytical solutions for the ensemble variance of reactive scalars in one-dimensional turbulent flows, in case of stationary conditions, homogeneous mean scalar gradient and turbulence, Dirichlet boundary conditions, and first order kinetics reactions. Simplified solutions and sensitivity analysis are also discussed. These solutions represent both analytical tools for preliminary estimations of the concentration variance and upwind spatial reconstruction schemes for CFD (Computational Fluid Dynamics)—RANS (Reynolds-Averaged Navier-Stokes) codes, which estimate the turbulent fluctuations of reactive scalars.

#### 1. Introduction

Modelling the turbulent fluctuations of a transported scalar interests several industrial processes and environmental phenomena, both in atmosphere ad water bodies, especially where the scalar fly time () is smaller than the Lagrangian integral time scale (microscale dispersion) or if the target parameter, such as damage due to high concentration levels, is nonlinear with respect to the transported scalar. In particular, scalar fluctuations are crucial in modelling pollutant reactions, as they normally depend on the instantaneous concentrations rather than their mean values.

Concentration fluctuations are relevant in several dispersion phenomena: accidental releases, dispersion of reactive pollutants, impact of odours, microscale air quality and water quality, and several industrial processes, such as combustion and process fluid treatments.

In this context, Reynolds’ average or mean concentration () is generally insufficient to represent the time and spatial evolutions of the instantaneous concentration field. Thus, a series of experimental (e.g., [1–5]) and numerical (e.g., [6–8]) studies have investigated concentration fluctuations, mainly focusing on the ensemble variance of concentration (), which is also involved in the definition of the intensity of fluctuations ().

Several numerical methods have been developed in order to evaluate the concentration variance. They are based on Direct Numerical Simulations (DNS; [9]), Large Eddy Simulations (LES) coupled with Lagrangian subgrid schemes [10], probability density function (pdf) models [11, 12], RANS models [13, 14], and Lagrangian micromixing numerical models [15–19].

These methods are somewhat constrained to the respect of the balance equation of the concentration variance, which was first derived by [20], for passive scalars. Reference [21] further discussed its terms and provided a similarity solution. Reference [22] proposed a formulation for the balance equation of , depending on Lagrangian parameters. The influence of the reactive terms is discussed in [23]. So far, the reference analytical solutions for concentration fluctuations were derived in decaying grid turbulence, under homogeneous nonstationary conditions: [15] treated the concentration variance, whereas [24] represented the probability density function of a transported scalar.

In this paper, 1D analytical solutions for the ensemble variance of a reactive scalar are derived, assuming stationary conditions, homogeneous turbulence, and mean scalar gradient, first order kinetics reactions, and Dirichlet boundary conditions. These solutions aim at providing an analytical solver for both preliminary estimations of the concentration variance and upwind schemes for CFD codes, potentially involving CFD models for pollutant dispersion based on the Finite Volume Method [13] or the Finite Difference Method [25, 26].

The paper is organized as follows. Section 2 briefly revises the balance equation of the concentration variance, according to [21], and the uniform solution provided by [15], both simply adapted to represent reactive scalars. In Section 3, this study proposes several analytical solutions for the concentration variance, under stationary conditions. Section 4 discusses a sensitivity analysis on the main nondimensional parameters of the analytical models of Section 3. Finally, Section 5 resumes the main conclusions of the study, whereas Appendix A reports a couple of analytical solutions for the mean scalar gradient, which is one of the key inputs of the analytical models of Section 3.

#### 2. Balance Equation for the Ensemble Variance of a Reactive Scalar in a Turbulent Flow

##### 2.1. Balance Equation of the Ensemble Variance

The balance equation for the concentration variance of a passive scalar dispersed in a turbulent flow was first derived by [20] and then discussed in detail in [21]. Their formulation is here briefly reported and adapted to consider reactive scalars.

The balance equation of the pollutant mass reads where is the instantaneous concentration, the molecular diffusion coefficient, and the velocity vector. Einstein notation applies to the subscript “” hereafter, if not otherwise stated. From left to right, the terms of (1) represent the local rate of change of the instantaneous concentration, the advective term of , the divergence of the molecular diffusion flux of , and the instantaneous reactive term ().

Concentration and velocity can be expressed according to Reynolds decomposition. Reynolds average (over-bar symbol) of (1) provides Hereafter, the apex “” denotes a turbulent fluctuation. Introducing the continuity equations for an incompressible turbulent flow into (2), one can write the balance equation for the mean concentration of a passive scalar: From left to right, the terms of (4) represent the local rate of change of the mean concentration, the advective term of , the divergence of the turbulent and the molecular diffusion fluxes of and the mean reactive term.

Subtracting (4) from (1), the balance equation for the concentration fluctuation is obtained: After multiplying (5) times and considering the following equality: one can write Averaging (7), one obtains the balance equation of the concentration variance of a reactive scalar dispersed in a turbulent incompressible flow: where the terms on the left hand side represent the local rate of change of , the advective term and the divergence of the turbulent flux of the concentration variance, respectively. On the right-hand side, (8) shows the production term of (always nonnegative), the dissipation rate of the concentration variance due to molecular diffusion ( always nonpositive), the divergence of the molecular diffusion flux of (negligible) and the reactive term (), which quantifies the direct effects of chemical and physical reactions.

The production term in (8) represents the increase in concentration variance due to nonhomogeneous conditions of the mean concentration field.

The dissipation rate of the concentration variance is ruled by molecular diffusion. Let us consider a fixed point and time: close fluid particles exchange pollutant mass due to molecular diffusion and thus homogenize the instantaneous concentration field. The concentration variance then decreases. This term is not negligible even at very high Reynolds numbers (Re) as the gradient of the instantaneous concentration would tend to infinity.

##### 2.2. Parameterizations of the Turbulent Fluxes

The turbulent fluxes in the balance equation of the concentration variance (8) assume approximated and simpler formulations, according to the “-theory,” which is the theory of the turbulent dispersion coefficients, as reported and further developed by several authors [27]: Here, and are the vectors of the turbulent dispersion coefficients of the scalar mean and variance, respectively. They are usually assumed equal and denoted by .

These parameterizations are commonly used in Eulerian numerical models, although they introduced several limitations. The most important shortcomings emerge in case of strongly nonlinear relationships between the turbulent flux and the mean/variance gradient. Further, (9) loses the information about the concentration-velocity covariance and the velocity autocorrelation. The resulting errors are not negligible at the microscale, even if we just compute the mean concentration.

##### 2.3. Parameterizations of the Dissipation Rate of the Concentration Variance

Several parameterizations are available for the dissipation rate of the concentration variance (). In particular, this study refers to the formula of [15], here reported and discussed. First, one may notice the following equalities: Considering (4) and (5), assuming a turbulent regime and that reaction rates do not affect the dissipation term, (10) can be expressed as follows: IECM (Interaction by the Exchange with the Conditional Mean; [15, 28]) represents a major micromixing formulation, which is alternative to other simpler schemes used in Lagrangian stochastic modelling for pollutant dispersion (e.g., [29]). IECM relies on the following expression for the Lagrangian derivative of the instantaneous concentration: where represents the mean concentration conditioned on the Lagrangian velocity vector () and is the mixing time (i.e., the time scale of the dissipation rate of the concentration variance), which is defined by [15, 30] or [31].

Combining (11) and (12), [15] obtains the following expression: Further, in case of uniform mean scalar gradient and 1D dispersion, [15] reports the following expression for the conditional mean: where stands for the standard deviation of velocity. This parameter can be provided by diagnostic tools, RANS codes (e.g., [32–35]), or LES models, which allow a detailed characterization of the turbulent structure of the atmospheric boundary layer (e.g., [36–38]).

Considering both (14) and (9), the dissipation rate of becomes As we consider 1D dispersion phenomena, the plume spread has the same dimension as the boundary layer height. Under these conditions, the mixing time can be related to the turbulent kinetic energy () and its dissipation rate () via Richardson constant (; [15]): where the Lagrangian integral time scale () is introduced. Its relationship with the dissipation rate of the turbulent kinetic—last equation in (16)—is reported by several authors (e.g., [15]). This relationship can be derived applying Taylor analysis and comparing the formulas of the plume spread alternatively depending on the Lagrangian structure function or the Lagrangian integral time scale.

Taking into account (16), the dissipation term becomes In conclusion, we can refer to a couple of simplified formulations for the dissipation rate of the concentration variance: a 3D generic formulation (13) and a simplified 1D expression in case of uniform mean scalar gradient (17).

##### 2.4. Representation of the Reactive Terms

The role of the reactive term in the balance equation of the concentration variance is here discussed, by alternatively considering 2nd-order kinetics, 1st-order kinetics, and 0th-order kinetics.

In case of 2nd order kinetics, the reactive term is represented by where species and are reactants and is the reaction rate.

The ensemble mean of (18) is equal to and the ratio between the last two terms in (19) defines the segregation coefficient After considering the fluctuation of (18): the reactive term in the balance equation of the variance becomes The first term in the final formulation of (22) depends on the reactant covariance, the second one on the variance of the control reactant () and the last one is a triple correlation term, whose eventual parameterization can refer to [39] The advantage of (23) simply relies on the fact that, in the limit of reactants never coexisting (), (23) guarantees that is exactly zero.

1st-order kinetics can be simply represented by imposing in (22): This kind of reactions always decreases both the instantaneous and the mean concentration. Further, a first order kinetics formula is linear with respect to the reactant concentration; thus, it locally decreases the concentration variance (provided the same mean scalar gradient) according to (24).

Finally, a 0th order kinetics (e.g., ) is equivalent to considering both and in (22). In this case, the concentration variance does not depend on the reaction rate: In conclusion, the expressions (22), (24), and (25) alternatively represent the reactive term in the balance equation of the concentration variance, in case of 2nd-order, 1st-order, and 0th-order kinetics reactions, respectively.

##### 2.5. 1D Analytical Solution of Sawford (2004) under Uniform and Nonstationary Conditions

Sawford [15] reported a 1D analytical solution of the concentration variance under homogeneous and nonstationary conditions. Here, this solution is reported and adapted to reactive scalars.

First, consider the 1D balance equation of the concentration variance, as resulting from the combination of (8), (9), (15), and (24), provided a uniform concentration variance: Defining the constant “,” one can write After integration from the initial time () to the generic time , one obtains Finally, considering (38), as explained in the following section, the uniform time-dependent solution for the concentration variance becomes The solution tends to an equilibrium value, when the production and the dissipation terms equalize: This formula highlights the importance of modelling the dissipation term in (26). When this term is absent (this is the case of Lagrangian models without any micromixing scheme or Eulerian models without any dissipation term for ), the concentration variance linearly grows with time, indefinitely: In case of a null mean scalar gradient there is no production term and the variance tends to zero, as follows:

#### 3. 1D Solutions for the Balance Equation of the Concentration Variance under Stationary Conditions

##### 3.1. Main Solution under Stationary Conditions

Provided 1D stationary conditions and homogeneous turbulence, the balance equation of the concentration variance (8) assumes the following form: Introducing the parameterizations for the turbulent fluxes (9) and (17), as well as the expression of the reactive term in case of 1st order kinetics reaction (24), one obtains After assuming the definition of the turbulent kinetic energy in 1D: considering (16) and (36) and dividing by , (35) becomes It is convenient to introduce the relationship between the turbulent dispersion coefficient and the Lagrangian integral time scale, as derived from a simple analysis of these turbulent scale parameters: The system (16), (36), and (38) provides another expression for the turbulent dispersion coefficient and allows writing (37) with no explicit dependency on : Equation (40) represents a 2nd-order inhomogeneous ODE (Ordinary Differential Equation) with the following constant coefficients: Its general solution assumes the form which is completed by the following definition: One can verify that the system (43) and (44) is a general solution of (41), as briefly described in the following. According to (43), the left hand side of (41) becomes The value of can be derived from (44): Replacing and in (45), according to (44) and (46), one obtains which finally verifies that the system (43) and (44) is a general solution of (41): Dirichlet boundary conditions can now be imposed on the left boundary () and the right boundary () of the domain Combining (49) with (50), one can obtain the value of : Thus, the constant from (49) becomes Considering the values of and from (52) and (51), respectively, (43) assumes the following form: Other few and simple algebraic passages finally provide a complete 1D solution of the balance equation of the concentration variance for a reactive scalar in a turbulent flow, under stationary conditions: It is immediate to verify that (54) satisfies the boundary conditions imposed.

The solution (54) is widely discussed and analysed in Section 4, where several examples are available (e.g., Figure 1), in terms of both nondimensional and dimensional physical quantities. Hereafter, we just provide some clarification about the role of the reactive term and the mean scalar gradient.

The reference solution is derived considering 1st-order kinetics reactions. In this case, however, one cannot obtain an exact uniform mean scalar gradient (stationary regime), in the absence of a source term, as we can appreciate from the balance equation of the mean concentration: This means that (54) can only consider a uniform mean gradient as an approximated condition for both 1st and 2nd order kinetics reactions.

On the other hand, an exact uniform scalar mean gradient can relate to a 0th-order kinetics: However, in this case there is no need for a reactive term in the balance equation of the variance, according to (25).

The mean scalar gradient is the key parameter in the production term of the concentration variance and affects the constant in (54). should be approximately uniform and can be calculated by means of simplified analytical solutions Appendix A or using the measured or simulated values of the mean concentration, where they are available in the domain of interest.

##### 3.2. Solution under Stationary Conditions, in the Absence of Mean Velocity

In the absence of mean velocity, (44) becomes Under these conditions, we can derive a simpler and alternative formulation of the general solution (54). Considering that , the first line of (54) becomes Introducing the following hyperbolic functions (and their derivatives) into (58), one can obtain the following expression: This change in the reference system implies that Combining (60), (61), and (62), the simplified solution of the 1D stationary balance equation of the concentration variance, in the absence of mean velocity, assumes the following expression: This solution, represented in Figure 1 (e.g., the green plot), is verified, as described in the following. One may first notice that where and are constants, directly defined by comparison with (63).

One can consider the second derivative of the variance in the new reference system: and then write (41), according to (63) and (65): which finally verifies that (63) is the solution of (41), in the absence of mean velocity.

Under these conditions, we can simply write the derivative of the concentration variance: Although it is not immediate to get a general expression to locate the maxima/minima, a simple solution is obtained in case , as plotted in Figure 1 (blue line). This assumption provides a maximum or a minimum at the domain centre (), in case or , respectively, according to the definition of the equilibrium variance (68), as provided in Section 3.3.

##### 3.3. Solution under Uniform and Stationary Conditions

The balance equation of the concentration variance shows a simple formulation under uniform and stationary conditions. In this case, the dissipation term and the reactive term (24) exactly equal the production term of (“equilibrium solution”): This expression also represents the value of the horizontal inflection points in the general solution (54). Under these conditions, the divergence of both the turbulent and the mean transport scalar fluxes is null.

#### 4. Sensitivity Analysis

The solutions (54), (63), and (68) of the 1D balance equation of the concentration variance are here discussed and analysed by means of a sensitivity analysis.

The results are presented in terms of nondimensional parameters: all the physical quantities are normalized (“”) over the reference time, length, and mass scales, respectively: , , and . Thus, the main nondimensional parameters, , and , are always equal to unity, whereas the other nondimensional quantities vary according to the following expressions: , , , and .

The “reference solution” of this analysis is defined by the following arbitrary choice, which simply provides a meaningful and general reference configuration: , , , , , and . Under these conditions, it follows that and .

Discussion first considers some generic features of the solution (54), as plotted in Figure 1 under three different configurations.

The reference solution (red plot) shows a central region of very low gradients (very weak turbulent and mean transports of ), where the production and the dissipation terms almost balance each other.

In this context, the reactive term plays an analogous role to the dissipation term and one may refer to them as “the dissipation terms” of the concentration variance, just for simplicity of notation.

At the edges of the domain, the production (left boundary) or the dissipation (right boundary) terms alternatively dominate. In these lateral regions, the turbulent scalar flux is always negative (directed upstream) and provides a positive (right boundary) or a negative (left boundary) contribution to the rate of change of the concentration variance, thus balancing the excess in the variance dissipation/production. At the same time, a weak nonnegative mean transport is responsible for a local decrease of .

In other words, the production term is uniform because turbulence is homogeneous as well as the mean scalar gradient. At the same time, the dissipation terms linearly grow with the concentration variance. Thus, the turbulent and the mean transport are responsible for driving the concentration variance away from the zone in excess of production (low variance) and towards the regions in excess of dissipation (high variance), in order to establish a stationary regime.

Figure 1 also reports a couple of examples of the solution for passive scalars (63) in the absence of mean velocity (green plot). In particular, the blue line represents the symmetric solution cited at the end of Section 3.2.

Figure 2 reports a sensitivity analysis of the main solution (54) on the choice of the nondimensional boundary values. Let us first examine Figure 2(a), where the left boundary value is null and the right boundary value is alternatively equal to 0, 1, or 2 times the equilibrium value of (68). In the first case (), the production term equals the dissipation terms in the central region of the domain. Elsewhere, the (uniform) production term is higher than the dissipation terms and the turbulent flux transports the excess of the concentration variance away from the domain. The mean velocity makes the solution be slightly asymmetric, with a peak located on the right half of the domain. In the second case (), the solution is almost uniform in the right half of the domain as the equilibrium value is already reached around the domain centre. Finally, when (), the solution is almost symmetric to the first case () with respect to the equilibrium line, for . In this region, the turbulent fluxes balance the deficiency in the variance production by transporting the concentration variance from the right boundary towards the domain centre.

**(a)**

**(b)**

Figure 2(b) considers the uniform solution of (68), established with the choice : the production term is everywhere balanced by the dissipation terms, without any effects of neither turbulent fluxes nor advection. The other solutions refer to () and are almost identical in the right half of the domain. Let us notice that the solution with () is symmetrical to the configuration with (; Figure 2(a)), with respect to the line of the equilibrium value.

Figure 3 investigates the role of the nondimensional mean velocity. In the absence of a mean flow (), the equilibrium value is reached at the centre of the domain, where a horizontal inflection point is established (red plot). The solution is symmetric with respect to its central value. Increasing the nondimensional mean velocity, the relative importance of the mean transport grows and the upstream boundary value of influences the solution more and more noticeably; further, the solution is not symmetric anymore and the equilibrium value is located in the right half of the domain.

The effects of the nondimensional standard deviation of velocity, whose square is directly proportional to the nondimensional turbulent kinetic energy, are shown in Figure 4(a). An increase in determines a highest variance production all over the domain. Considering the highest and the lowest values of , the nondimensional variance lies out of the range provided by the Dirichlet boundary values, in the most of the domain. As is relatively low, one can always detect at least a maximum, a minimum, or a horizontal inflection point in the profile of , in the very inner domain.

**(a)**

**(b)**

Contrarily, the dependency of the nondimensional variance on the nondimensional reaction rate has an opposite behaviour. In fact, the reactive term always acts as a dissipation term (linear in ) and its increase tends to lower down the scalar variance all over the domain (Figure 4(b)). In case of extreme values of , one may notice that the most of the domain is characterized by values lying out of the interval provided by the Dirichlet boundary conditions.

The sensitivities on the Lagrangian integral time scale and the mean scalar gradient are finally shown in Figure 5 in terms of dimensional parameters, just to provide a more physical representation of the role of these quantities. An increase in or improves the “production term” with similar effects on the profile of the concentration variance.

**(a)**

**(b)**

#### 5. Conclusions

The study presents 1D analytical solutions for the ensemble variance of reactive scalars in turbulent flows, in case of stationary conditions, homogeneous mean scalar gradient and turbulence, Dirichlet boundary conditions and 1st-order kinetics reactions.

A sensitivity analysis has discussed the role of the different terms in the balance equation of the concentration variance, with focus on the production term, the dissipation rate of the scalar variance, the reactive term, the turbulent transport and advection. The analysis has detected three main scale parameters: the domain length, the Lagrangian integral time scale and the product of the mean scalar gradient times the square of the domain length. Thus, the analysis has investigated the dependency of the solution on the nondimensional boundary values, reaction rate, mean, and standard deviation of velocity.

These analytical solutions represent an immediate tool for preliminary estimates of the concentration variance in several application fields, where concentration fluctuations play a relevant role: accidental releases, pollutant reactions, odour assessment, microscale air quality and water quality, and several industrial processes, such as combustion. Although the solutions referred to 1D configurations (e.g., unstable turbulent boundary layers of plug-flow reactors or pollutant treatment devices), they can still provide approximate estimations for simplified 3D configurations.

Finally, the proposed solutions represent upwind spatial reconstruction schemes for the concentration variance, as they can be implemented in CFD-RANS codes, as well as in prognostic or diagnostic meteorological and ocean models, which simulate the turbulent fluctuations of reactive pollutants.

#### Appendix

#### A. 1D Analytical Solutions for the Mean Scalar Gradient

Several authors (e.g., [40, 41]) presented analytical solutions on the scalar transport in fluid flows. Here, some of them are reported to compute the mean scalar gradient, which is a key input parameter for the analytical solutions of the concentration variance.

##### A.1. Mean Transport, 0th-Order Kinetics, and Continuous Source

In case of a continuous and infinite pollutant source over an horizontal plane and neglecting the along-flow dispersion , the steady-state solution of the balance equation for the mean concentration reads ([40]): where is the time rate of mass injection per unit area. The derivative of (A.1) with respect to is The curves in Figure 6 show the solutions of (A.1) and (A.2), normalized by the value at the origin, for . Given the exponential nature of (A.1) and (A.2), the two normalized solutions coincide.

##### A.2. Turbulent Transport, Instantaneous Source, and Infinite Domain

The solution to the advective-diffusion equation is here reported, in case of instantaneous point source and infinite domain.

Considering a point source in stagnant conditions, the mean concentration can be described as follows [40]: where represents the source strength and is the injection point of the tracer.

For , one can find the maximum value of the mean concentration: If we consider , the concentration tends to an infinite value: .

We can derive (A.4) with respect to time as Instead, the mean scalar gradient can be expressed as The analytical solution of (A.6) can be possibly used as an approximate estimation of the mean scalar gradient, which relevantly influences the analytical solutions for the concentration variance (Section 3). Although (A.6) referred to nonstationary conditions, it can approximate the mean scalar gradient during a late stage of a dispersion process.

#### Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgment

The first author is deeply grateful to Professor Brian Sawford, who provided comments and advices at the beginning of this study.