(8.34)
In this section the expected (mean) value of a random quantity X (v, Ψ) is denoted by brackets:
(8.35)
and the fluctuating component is denoted by a prime: X′ =X – ⟨X⟩.
The evolution of the turbulent field at a fixed point of the flow is simulated by solving a transport equation for the joint PDF. Fox (2003) provides details on the derivation of this equation. In short, the equation can be established by equating two independent expressions for the expected value of the total derivative of an arbitrary scalar function F(v, Ψ):
(8.36)
and
(8.37)
The total derivatives of the fluid velocity components (i = 1, 3) and the concentration components (n = 1, N) that appear in (8.37) are substituted by expressions derived from the Navier–Stokes and continuity equations, respectively
(8.38)
(8.39)
where ν [m2 s–1] is the kinematic viscosity coefficient and Dn [m2 s–1] is a molecular diffusion coefficient. Quantity gi denotes the gravitational acceleration and ρa the air density. The expected values of the different terms appearing in expressions (8.38) and (8.39) are expressed as a function of the joint PDF, applying (8.35) for the different variables. Noting that the equality between (8.36) and (8.37) must hold for any arbitrary choice of scalar function F, one derives the following transport equation (Fox, 2003; Cassiani et al., 2010)
(8.40)
where the conditional fluxes are
(8.41)
and
(8.42)
The first terms in (8.40) refer successively to the local rate of change of the joint PDF and to its transport in the geometric space (xi). The first term on the right-hand side of the equation accounts for the effects of viscous stress, pressure fluctuations and gravity, and the second term for the transport by molecular fluxes in the composition space and for chemical sources. Terms that require multiple-point information (gradient and Laplacian terms) are not closed, and closure expressions based on empirical information or DNS data must be added to the system.
The numerical solution of the transport equation (8.40) provides the value of the joint PDF at a specified point r and time t, from which mean physical quantities can be derived by integration over velocities and concentrations following (8.35). An advantage of the stochastic method over the classic Reynolds decomposition method is that it provides not only the mean values of the velocity and concentrations, but also the full PDF from which to compute higher moments of these quantities, including the covariances between fluctuating quantities. In particular, no closure assumption is needed to treat the chemical source term involving products of concentrations.
Solution of (8.40) is computationally burdensome and more economical methods are usually adopted. A common approach is to consider a simpler composition PDF obtained by integrating the joint velocity-composition PDF over the entire velocity phase space:
(8.43)
In this case, one assumes that the stochastic wind fields are provided by an external turbulence model. The transport equation for the composition PDF is derived by integrating (8.40) over the entire velocity phase space (see Fox, 2003 for a complete derivation). One finds:
(8.44)
The terms on the left-hand side of (8.44) represent again the time evolution of the PDF and its advection by the mean wind field 〈vi〉 in the spatial space xi (sometimes referred to as macro-mixing). The terms on the right-hand side include the effect of unresolved concentration-conditioned velocity fluctuations (meso-mixing) and the transport in the composition space due to molecular mixing (micro-mixing). Again, the chemical term involves only single-point information, and is described in a closed form. One usually assumes that the term related to the concentration-conditioned velocity can be parameterized by
(8.45)
where K is an eddy diffusion coefficient, supposed here to be isotropic. Similarly, the conditional term that accounts for the dissipation of the concentration fluctuations by molecular diffusion is also unclosed, and can be parameterized by a linear dissipation toward the mean concentration state 〈Ψn〉.
It has been shown (Valiño, 1998) that the solution of the transport equation for the composition PDF (8.44) is equivalent to the solution provided by an ensemble of M stochastic partial differential equations (PDEs)
(8.46)
where m denotes one of the M possible realizations and Tmix a mixing timescale. The two first terms on the right-hand side of this equation represent the effects of the large-scale advection and of mixing by eddy diffusion. The third term introduces the effects of random motions with dξi denoting the random increment of an uncorrelated Wiener process with a zero average and a variance equal to dt (see Section 4.11.2). The fourth term evaluates the chemical sources and sinks; it implicitly accounts for the effects of fluctuations (chemical eddies) on the reaction rates. The last term introduces a relaxation toward the mean concentration. The average of the concentration field as well as the covariances are obtained from the ensemble of solutions in the M different realizations.
Stochastic model modules, if inserted into a coarser chemical transport modeling system, provide an approach to simulate the influence of chemical and dynamical subgrid processes on large-scale dynamics and chemistry, including, for example, the dispersion of reactive plumes in the PBL and the effects of subgrid heterogeneities in surface emissions. The formulation discussed above is Eulerian because the PDFs are calculated at fixed points in space. The evolution of fluid particle properties in turbulent reacting flows can also be simulated by Lagrangian stochastic models (see Section 4.11.2) in which the flow is represented by a large number of particles, each being characterized by its own position, velocity, and chemical composition. These different properties evolve according to stochastic Lagrangian model equations, and each sampled particle is assumed to be representative of a different realization of the flow.
8.6 Numerical Solution of the Diffusion Equation
Diffusion equations are often used in atmospheric models to parameterize the effects of small-scale turbulence. They must be solved numerically except in idealized cases (see Box 8.2). Here, we first consider numerical methods adopted to solve the 1-D diffusion equation
(8.47)
which is frequently encountered in chemical transport models.
Box 8.2 Analytic Solution for the 1-D Diffusion Equation
Consider the 1-D diffusion equation for an inert tracer concentration Ψ(x, t)
∂Ψ/∂t = K∂2Ψ/∂x2
where K is a constant diffusion coefficient. We assume at the initial time (t = 0) that the entire mass of the tracer is concentrated at the spatial origin (x = 0). Thus, the initial value of the concentration is Ψ(x, 0) = S δ(x) where δ(x) is the Dirac function (this function equals zero everywhere except at x = 0 and its integral value is 1). For an infinite space with Ψ = 0 for x → ± ∞, the analytic solution is the Gaussian expression
Ψ(x, t) = {S/(4πK t)1/2} exp {–[x2/(4 K t)] }
This expression provides the basis for Gaussian plume models. The plume dilutes with a variance σ2 = 2 K t that increases linearly with time t. See Box 8.2 Figure 1 and Section 4.12.1.
Box 8.2 Figure 1 Dilution of a plume.
8.6.1 Explicit Schemes for the 1-D Diffusion Equation
If the diffusion coefficient K is constant, the partial differential parabolic equation written now as
(8.48)
can be approximated by the forward-in-time, centered-in-space (FTCS) scheme:
(8.49)
where, as in previous chapters, Δx and Δt represent the grid interval and the time step, respectively. This expression is easy to solve due to its explicit nature. It is, however, only first-order accurate in time. Following von Neumann’s stability analysis considered in Chapter 7, in which the evolution of a wave with wavenumber k is examined, the amplification factor |g(k)| associated with the FTCS method is
|g(k)| = 1 − 2β(1 − cos (kΔx))
> (8.50)
For the method to be stable, this quantity must be less than 1 for any value of k. This is achieved if the dimensionless parameter
(8.51)
is less than ½. Thus, to ensure stability of the solution, the integration time step Δt must be smaller than Δx2/2K. This timescale characterizes the diffusion of air parcels across the grid cell width Δx. However, even though this criterion guarantees stability, it may not be sufficient to ensure a correct simulation of the shortest resolved waves, and specifically the 2Δx wave mode. Durran (2010) shows that a more appropriate condition is
(8.52)
In the explicit Richardson scheme, the time derivative in the FTCS scheme is replaced by a centered difference:
(8.53)
This scheme is second-order accurate in time but is unconditionally unstable. Unconditional stability, however, can be achieved by splitting the term into . The resulting DuFort–Frankel algorithm is:
(8.54)
Although the scheme appears implicit, the solution is easily derived if the dependent variables are known at the two previous time steps tn and tn–1:
(8.55)
The amplification factor is:
(8.56)
In spite of its unconditional stability, the DuFort–Frankel scheme does not correctly damp the short-wavelength components of the solution when implemented with an excessively large time step. It is therefore recommended to adopt a value of β smaller than ½.
Accuracy of the solution can be improved by considering a weighted time derivative over three time levels tn–1, tn, and tn+1, and a more elaborate formulation of the space derivative:
(8.57)
where γ and θ are adjustable parameters. For θ = 0, this expression corresponds to the explicit FTCS scheme if γ = 0, to the Richardson scheme if γ = ½, and to the DuFort–Frankel scheme if γ = –½ + β. Stability conditions vary with the value of parameter γ: the stability limit for β varies from about 0.35 to 0.5 when γ increases from 0 to 6. The stencils corresponding to different explicit algorithms are shown in Figure 8.6.
Figure 8.6 Stencils for different explicit schemes to approximate the 1-D diffusion equation. Indices j and n refer to the spatial and temporal dimensions, respectively.
If, as is often the case in atmospheric applications, the diffusion coefficient K is not a constant parameter, but depends on the spatial variable x, a more general form of the diffusion equation (8.47) must be considered. In this case, the explicit FTCS algorithm, for example, takes the form
(8.58)
or
(8.59)
Index j + ½ corresponds to the point located at equal distance between grid-points xj and xj+1. If the value of K is not known at such intermediate locations but only at the grid points, one can use a simple interpolation as Kj+1/2 = (Kj + Kj+1)/2. This scheme is stable only if
(8.60)
for all values of j. However, as indicated above, a step twice as small should be preferred to ensure a proper treatment of the shortest resolved waves. If both the diffusion coefficient and the grid spacing vary in space, we approximate the diffusion term by a second-order finite difference expression:
(8.61)
or
(8.62)
where
with
Δxj = xj − xj − 1, Δxj + 1 = xj + 1 − xj and Δxj + Δxj + 1 = xj + 1 − xj − 1.
8.6.2 Implicit Schemes for the 1-D Diffusion Equation
The implicit form of the FTCS scheme
(8.63)
provides solutions that are unconditionally stable because the amplification factor
(8.64)
is less than 1 for all values of β. The errors on the solution can, however, be substantial for large values of the time step because this algorithm is only first-order accurate in time.
Expression (8.63) represents a tridiagonal system of J – 1 linear equations:
(8.65)
The solution requires that boundary conditions be specified for j = 0 and j = J. The Thomas algorithm, often adopted to solve tridiagonal systems, is given in Box 4.4.
The accuracy of the fully implicit FTCS algorithm can be improved by using a combination of spatial derivatives at time tn and tn+1:
(8.66)
This algorithm, called the Crank–Nicholson scheme, is unconditionally stable since the associated amplification factor
(8.67)
is always less than or equal to 1. It is second-order accurate in time and space. As in the fully implicit method, a system of linear equations involving a tridiagonal matrix must be solved.
A general approach is to consider a weighted time differencing over three time levels together with a combination of implicit and explicit approximations for the space derivatives:
(8.68)
where γ and θ are again adjustable parameters. The Crank–Nicholson scheme corresponds to the case where γ = 0 and θ = 1/2. For γ = 1/2 and θ = 1, one obtains the following three-level scheme:
(8.69)
which is second-order accurate in space and time, and is unconditionally stable. Figure 8.7 gives schematic representations of the different implicit schemes.
Figure 8.7 Stencils for different implicit schemes to approximate the 1-D diffusion equation. Indices j and n refer to the spatial and temporal dimensions, respectively.
If the diffusion coefficient K varies in space, expressions similar to (8.58) can easily be established for the fully implicit and for the Crank–Nicholson schemes; in both cases, they lead to a tridiagonal system of linear equations. These algorithms are usually stable for any value of time step Δt. More elaborate and accurate forms can be adopted for the discretization of the second derivative in space (see Table 4.2). In this case, the matrix associated with the system of linear equations may not be tridiagonal anymore. Techniques to solve this system, such as the Gauss elimination algorithm, are often prohibitively expensive.
8.6.3 Numerical Algorithms for the Multidimensional Diffusion Equation
The explicit FTCS method described earlier can easily be generalized to more than one dimension. However, because of the constraints on the adopted time step, implicit methods are often preferred. Again, a system of coupled algebraic equations can easily be established, but the matrix of this system, although sparse, is no longer tridiagonal.
An alternative to the fully implicit approach is to apply operator splitting and treat the problem as a succession of 1-D problems. In the popular alternating-direction implicit (ADI) method, the time step is divided into sub-steps at which a different dimension is treated implicitly (Figure 8.8). For example, in the 2-D case, the diffusion equation expressed as
(8.70)
is integrated from time tn to tn+1 by considering the following sequence
(8.71)
(8.72)
where Ψ* is an intermediate value for function Ψ. Assuming constant diffusion coefficients and defining
(8.73)
and
(8.74)
the resulting algebraic equations are:
(8.75)
(8.76)
where j and k are the indices referring to the spatial discretization in the x- and y-direction, respectively. Each step requires solving an implicit equation in one dimension; the first step is solved for parameter k fixed and the second step for parameter j fixed. This approach requires that at each time step two tridiagonal systems be solved; it is thus computationally efficient. The von Neumann stability analysis can be applied for each step and for the steps combined to yield:
(8.77)
Figure 8.8 Stencil representing the two successive steps in the 2-D alternating direction method.
This method is unconditionally stable because the amplification factor is always less than or equal to 1. It is second-order accurate in space and time. We close by noting that diffusive schemes, as those discussed above, are not necessarily used to represent physical processes, but rather to stabilize algorithms adopted for solving the advection equation. Furt
her, unwanted numerical diffusion can be produced by low-order advection schemes (see Box 8.3).
Box 8.3 Numerical Diffusion in Atmospheric Models
In addition to the diffusion parameterization used to describe atmospheric turbulence, spurious numerical diffusion may be produced in Eulerian models by the adopted advection schemes (see Chapter 7) and by the use of numerical filters, which are often applied to avoid numerical instability when solving the transport or dynamical equations. This additional (unphysical) diffusion leads to excessive mixing and can be reduced by using a higher-order advection scheme and a finer grid to resolve spatial gradients. Numerical diffusion is a particular problem at dynamical barriers through which physical transport is restricted (see Section 8.12). In Lagrangian models, on the other hand, air parcels moving with the flow are considered to be isolated from other air parcels, and no mixing occurs. Without an appropriate parameterization to account for mixing processes between neighboring air parcels, Lagrangian models may overestimate spatial gradients and produce unrealistic small-scale structures (Collins et al., 1997; McKenna et al., 2002). To account for small-scale mixing processes in a Lagrangian model, an expression can be added that brings the mixing ratio μi of a given air parcel closer to the background mixing ratio . Collins et al. (1997) use a relaxation term where the degree of exchange d is considerably faster in the troposphere than in the stratosphere. The background mixing ratio in a given area is defined as the average for all air parcels located within that area. In an alternate approach, McKenna et al. (2002) entirely mix two air parcels when their spatial separation is smaller than a given threshold value.
Modeling of Atmospheric Chemistry Page 44