Analytical Solutions of the Balance Equation for the Scalar Variance in One-Dimensional Turbulent Flows under Stationary Conditions

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.


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.
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 2 , depending on Lagrangian parameters.

Advances in Mathematical Physics
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.

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 2 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 2 , the advective term and the divergence Advances in Mathematical Physics 3 of the turbulent flux of the concentration variance, respectively. On the right-hand side, (8) shows the production term of 2 (always nonnegative), the dissipation rate of the concentration variance due to molecular diffusion ( = −2 ( / ) 2 always nonpositive), the divergence of the molecular diffusion flux of 2 (negligible) and the reactive term ( 2 ≡ 2 ), 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.

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.

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].
Considering both (14) and (9), the dissipation rate of 2 becomes As we consider 1D dispersion phenomena, the plume spread has the same dimension as the boundary layer height. Under 4 Advances in Mathematical Physics 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).

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 2ndorder 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 ( = −1), (23) guarantees that 2 is exactly zero. 1st-order kinetics can be simply represented by imposing = 1 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 = 1 and = 1 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. 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:

1D Analytical Solution of
Defining the constant " , " one can write Advances in Mathematical Physics 5 After integration from the initial time ( = 0) 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 2 ), 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:

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) 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 2 can be derived from (44): Replacing and 2 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 ( = 0) and the right boundary ( = ) of the domain Combining (49) with (50), one can obtain the value of 2 : Thus, the constant 1 from (49) becomes Considering the values of 1 and 2 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.

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 ≡ 1 = − 2 , the first line of (54) becomes Introducing the following hyperbolic functions (and their derivatives) into (58), one can obtain the following expression: 8

Advances in Mathematical Physics
This change in the reference system 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 ( 2 ,0 = 2 , ), as plotted in Figure 1 (blue line). This assumption provides a maximum or a minimum at the domain centre ( = 0), in case ( 2 ,0 < 2 ,eq ) or ( 2 ,0 > 2 ,eq ), respectively, according to the definition of the equilibrium variance (68), as provided in Section 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 2 ("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.

Sensitivity Analysis
The solutions (54) The reference solution (red plot) shows a central region of very low gradients (very weak turbulent and mean transports of 2 ), 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 2 .
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 2 production (low variance) and towards the regions in excess of 2 dissipation (high variance), in order to establish a stationary regime.  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 ( 2 , , = 0), 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 ( 2 , , = 2 ,eq, ), 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 ( 2 , , = 2 2 ,eq, ), the solution is almost symmetric to the first case ( 2 , , = 0) with respect to the equilibrium line, for > /2. 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.  Figure 3 investigates the role of the nondimensional mean velocity. In the absence of a mean flow ( = 0), 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 2 , 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 2 , , in the very inner domain. 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 2 ) 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 2 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.

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. ,eq is 8.694 g 2 /m 2 , 18.90 g 2 /m 2 , 32.507 g 2 /m 2 , 68.66 g 2 /m 2 , and 115 g 2 /m 2 , respectively. (b) Reference dimensional solution (but = 0.001/s) with sensitivity analysis on the scalar mean gradient, which is alternatively equal to 0 g/m 2 (+), 0.1765 g/m 2 (×), 0.353 g/m 2 ( * ), 0.44125 g/m 2 (◻), or 0.5295 g/m 2 (◼)-2 ,eq is 0 g 2 /m 2 , 8.127 g 2 /m 2 , 32.507 g 2 /m 2 , 50.79 g 2 /m 2 , and 73.14 g 2 /m 2 , respectively. where is the time rate of mass injection per unit area. The derivative of (A.1) with respect to is ( ) = − 2 exp (− ) . (A. 2) The curves in Figure 6 show the solutions of (A.1) and (A.2), normalized by the value at the origin, for > 0. 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]: 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.