HYNHallis In dip-kalita Model fop tap-Stagnant Ralsttii Happison Has 17, 1987 ------- Report No. EPA-910/9-87-179 WYNDvalley An Air-Quality Model for Near-Stagnant Flows in Constricted Terrain Halstead Harrison WYNDsoft, Inc. 6333 77th Avenue SE Mercer Island, Washington 98040 May 17, 1987 Prepared for U S Environmental Protection Agency Region 10 1200 Sixth Avenue Seattle, Washington 98101-3183 ------- On Air-Quality Modeling in Complex Terrain, with Emphasis on WYNDvalley, A Dispersion Model for Near-Stagnant Flows Halatcad Harrison WYNDeoft Inc. 6333 77th Avenue SE Mercer Island Washington 98040 (206)-232-1810 The mathematics of air-quality modeling is discussed with special attention to problems associated with complex terrain. Practical solutions to some of these problems are suggested and an implementation is described as WYNDvalley, a dispersion model adapted to the study of trace species concentrations in near-stagnant conditions within complex valley drainage systems. The model is illustrated by application to the Methow, Twisp, and Chewack valleys of Okanogan County, Washington. I. Introduction To guide wise management of the air we live in, we wish to develop models and simulations to permit accurate estimates of the concentrations of air pollutants, appropriately averaged over time and space. This is a difficult task. The motions of the air are characterized by erratic vascillations that display temporal scales between seconds and seasons, and spatial scales between millimeters and thousands of kilometers. Limitations of cost affect the density of observations that we may make to specify the initial state of the systems we wish to model, and to follow the evolution of the winds and emissions throughout the period that interests us. Perhaps surprisingly in these times of exponentially improving computer technology, we are also limited by the computer speeds and memory that would be needed for near-rigorous solutions to the equations we pose. For these reasons practical progress requires approximation. In the following section I shall explore alternate choices for these approximations, starting from common assumptions and with moderate rigor. Here I shall emphasize practical implementation of computer models to be useful under condi- tions when the flows are dominated by constraints imposed by local topography. In section III I shall illustrate these choices in more detail through a specific model, WYNDvalley, and discuss some results from a specific simulation, the sen- sitivities of these results to the data and to the model's assumptions, accuracies ------- - 2- and limitations, ways the model might be improved, and the likely costs and benefits of effort to effect these improvements. I conclude with a summary and documentary appendices. II. Some Mathematics of Tracer Dbpersion in the Atmosphere Readers anxious to get on with specific details of WYNDvalley may turn immedi- ately to section ID. My purpose in including the present discussion is to make explicit a rather long series of approximations that must be introduced along the way towards solutions of equations over which ... as it seems to me ... unnecessary arguments sometimes erupt about the merits of alternate and disparate seeming algorithms. The point is that excessive partisanship over 'Gaussian Plume1 or 'Box' models is ... or ought to be ... somewhat vitiated by an appreciation of the many severe approxi- mations that are common to both these useful approaches. In the sub-sections that follow I shall highlight those approximations that are common to both Gaussian and box models, and those separate to each. I shall conclude this section with a discussion of the virtues and limitations of the two approaches to air-quality modeling, with emphasis on matching models to the questions that are to be asked of them. A. The Dispersion Equation: General Approximations The central concept of tracer dispersion is, very simply, that mass must be con- served. The governing equation is: In [l] C(x,y ,z,t) is the concentration of a trace chemical species whose motions and dispersion we wish to model as functions of the spatial coordinates * ,y ,z, and time f; y is a generalized representation of the gradient operator which, in Cartesian coordinates, becomes d/dx + d/dy + d/d«. The wind velocity is u(x ,y ,z ,t) and Q (x ,y ,z ,t) is a source or sink. The bold fonts emphasize that v and u are vectors [while C and Q are scalars] and the central dot between the vectors denotes the scalar product. ------- - 3 - One can imagine solving [l] directly by numerical integration: set up a three dimensional grid, specify C(x ,y ,z u(x ,y ,z ,t) and Q (* ,y ,z ,t), and then use Euler's rule .. or other convenient algorithm .. to solve for C(x,y ,x ,tt +dt). Problems arise, however, because very small increments must be specified in each of the coordinates. We have all seen near-neighbor wind flags blowing in disparate directions, and we have all experienced wind speeds and directions that change abruptly within a few seconds. Suppose then, for illustration, that we specify a grid with 10 meter resolution over a 10 km x 10 km x 1 km box. This would require 108 grid elements just to specify the instantaneous concentration field. I understand that there are parallel machines that can do this, barely, but for our purposes the most apt comparison may be to the CRAY I, that in one common configuration can address 'only' 7 x io6 random-access memory locations. Ignoring this, just for the illustration, with 10 second time resolution for the simulated fields, and 10 Mhz computational clock speed with 10 cycles per 4-byte multiplication, this brute-force integration would require approximately 10 com- puter minutes per simulated hour. That is, the computer would be running only six times faster than the world, for a 10 km x 10 km x 1 km box that often would be too small for our interest, anyway. This illustration is silly for a further important reason: we simply cannot afford the density of measuring instruments necessary to specify the initial concentra- tions and winds in those 10® grid elements, nor to verify the precision of our cal- culations, after the fact. For these reasons, let us proceed to condition equation [1] for approximations that are forced upon us anyway by the intrinsic variability of our noisy atmo- sphere. First let's try to stretch out the time unit for each integration step. Notice that in the brute-force example above, a ten second time step is limited to wind speeds of less than one meter per second, across the grid lengths of ten meters. But most measurements will report winds and tracer concentrations that have been averaged over tens of minutes, or hours, and often at speeds exceeding 1 m/s. ------- - 4 - To spread out the time steps, let's define: +r/2 fa ,v ,*)s — / fa >y >* J W [2] r -r/2 where { may be either the concentration, C ,the source, Q, or any component of the wind velocity, u. The time, r is a somewhat imprecisely defined subscale over which we want to smooth our equation variables to match more closely the real averaging times of our measurements. With [2] we next define primed quantities such that: 0{x,y ,',t) s 0{x,y,z) + e {x,y,z,t) [3c] Q(x,1,i,t)ea Q(x,y,z) + J (x,y,z,t) [3q] + [3u] V (x ,y ,z,t)ss v(x ,y,z) +r/ (z ,y ,z ,t) [3v] u>(x,y,t,t)s3iv{x,ylz)+u! (x ,y ,z ,t) [3w] In [3] the barred quantities are time-averaged over r and the primed quantities express the more slowly varying fluctuations at times greater than r, « , v, and w are the three scalar components of the vectorial wind, respectively parallel to x, y, and z. Note that by definition: Fr=qr=u7=tr=w7=0 [2'] Substituting [3u] back into [l], while temporarily neglecting the coordinates y and z: T? —T, [(5 + t' Xr + * >]+<« + * ' IU1 Expanding and with a second overbar to denote averaging the operators of [la] over r, + ) + Q + 9" [lb] The middle two terms within the [ ] average to zero by definition [2'], as does q7. Then dropping double overbars and noting that HC = uC, —^[i®+r?l + 5 M ------- - 5 - The effect of this manipulation is to lump into the ««• all the temporal fluctua- tions that occur in times that are short with respect to r, the measurement time. Since by definition this term is no*-measured, it has to be modeled. What shall we do with it? We could set it equal to zero: that after all is a model, too. But this vigorous approximation excessively degrades the accuracy of what we're trying to do. One useful approach is to write additional differential equations of the form du -c „ [—5- -j — 1 — = F[n-J but this expansion rapidly increases the complexity of the equation set, and clo- sure problems intrude with higher-order moments. By far the most common approach is to model the ire7 in terms of the quantities that are measured at the time scale, r, namely a and C. This is the natural 'next higher approximation'. Drawing heavily upon analogies to the diffusion of heat in rigid bodies [an imperfect analogy that will causes trouble later, because the air is not rigid], Batchelor, Prandl, Reynolds, Richardson, and many others have explored the utility of postulating: For obscure reasons this is called the 'exchange' approximation (or the 'austasch' approximation by those who wish to make it even more obscure). The Kt defined by [4] is ultimately an empirical coefficient, called a diffusivity. It has dimensions of lengtlAime"1, and it may be thought of as the product of a root- mean-square 'eddy velocity' times a scale length over which these eddies are effective in transporting momentum. The subscript 2 denotes «,hat K, is a scalar component of the vector K. Some authors generalize [4] still further by defining K as a tensor, but this is really just mathematical exhibitionism: remember that the components of K are ultimately empirical and can be measured in ways that force a vectorial interpretation. Now we can proceed by substituting [4] into [lc], hereafter dropping the overbar notation as excessively busy; remember, however, that all the measurements and modeled variables remain implicitly as averaged over t. We get: ------- - 6- *£s=[_l + _L + ±) dt [ dx dy dzJ -UC+K— + K— 4-K — uC + K' 17+ K' 9F + ^~ dz + C) = Q dx dx dy dy dx dz [5a] and dx dy dz „ dC „ iC „ dC] ' dx ' dy ' dz if dC j, dC jr dC jV« iV¦ TJ"" A* 1 dx ' dy ' dz j. dC „ dC j, dC ' dx ' dy ' dz is diagonal (5b] Condition [5a] is just a statement that the iPiF, Pc7, and Ore7 do not affect the orthogonality of the * -y -z coordinate system. [This is not equivalent to stating that the turbulence is isotropic.] Condition [5b] again states that the diffusivities coordinate system, and [5b] again states that the diffusivities need not be treated as tensors. In comparison to other approximations that follow, neither [5a] or [5b] are very severe. With these approximations [id] simplifies to: dC a(nC) d(vC) d(ivC) dt dx dy dz [le] d_ dx dc K*~di dy k¦ dC dx j, dC K'-d7 Q In [le] the K,,Kt, and K, have so far been left within the exterior spatial differential operators to emphasize that they are in fact functions of scale. This property of the atmosphere makes eddy diffusion fundamentally different from that of heat in rigid solids. Measurements over widely differing scales show that atmospheric diffusivities are larger for bigger puffs than smaller [See for example: Gifford, 1982].' At scales larger than about 30 km eddy-diffusivity coefficients scale roughly as X2, where X is a curvature scale of a diffusing puff, that may or may not be isotropic. ------- - 7 - At smaller scales K,,Kt, and K, vary with diminishing powers of X,, Ar> and X,, but they do not approach scale independence until puff sizes of millimeters, or less. To proceed with any elegance beyond [le] therefore requires differentiation at all scales, a need perhaps best assuaged by spectral expansions of the variables in Fourier series. Alternately, with less elegance but perhaps more common sense, one may ignore the scale dependence of the K( at this step, and proceed with the differentiations as if the they were scale independent; then at a later step empiri- cal scale corrections can be imposed a postertori. It is further a common approxi- mation to neglect the spatial gradients of «, v, and w at this step, though in fact this is not strictly required ... nor is it a very good idea. However, for simplicity and because Gaussian models do require this further approximation, I shall make it here to achieve a 'final1 working version of this series of equations [1]: dC dt —1» dC dC - XD dc dx ,dy bz + *'-0 + *-0 + *-0+« M Let us pause for a moment to reflect on what... if anything ... useful that we have done by recasting [1] into [if]. Recall that the central idea was to subsume all the temporal fluctuations that occur over time scales smaller than the meas- urements' averaging times, r, into terms of the form JFc7, and then through the 'exchange' approximation to parameterize these into K(. This may permit us both larger time steps and a wider spatial net, because at constant wind speeds each parcel travels further during the longer times. Unfortunately, however, this latter advantage is not automatic, as a new and serious difficulty intrudes when we attempt to integrate [If] numerically with a relaxed grid. This trouble is perhaps best illustrated with a sketch, as in figure 1. Suppose that an initial puff concentration is uniformly spread between two adjacent grid points, as with the prism at the left of the figure. With a wind blowing to the right, then after a time increment dt the puff might be translated to the position of the center prism. But with a finite grid the mass distribution of this puff would be interpreted as split between the two adjacent grid elements, as with the pair of prisms at the right. Suddenly the puff has 'diffused' to greater width, purely as an artifact of the finite grid. Unless the translation fortuitously places the puff exactly over a neighboring grid, which is hard to achieve generally when ------- uotsnj.^iQ xeaxaauinN U U 1 i i r BB ^ T 3^inE> x j ------- - 8 - the winds vary in direction and speed, the numerical diffusion may easily exceed the physical diffusion. One may minimize this artifact by requiring tight grids, but this is a step backwards in the whole approach from [l] towards [if]. What is to be done about this new problem? Two approaches are common: brute force, or a clever transform. Both have virtues and defects. The WYNDvalley model that will be described in the next section is essentially a brute-force approach: set up a moderately fine grid (100's of meters) with moderately small integration steps (100's of seconds), accept moderate numerical diffusion, but compensate for it by diminishing the external diffusion parameters, K(. The advantages and limitations of this approach will be described in appropriate detail later. Because the transform approach is the more commonly used, how- ever, it is proper here to continue its development for a few more paragraphs, to illuminate the comparison. B. Gaussian Models: Special Approximations In this sub-section I amplify briefly on an additional set of approximations inherent in Gaussian solutions to [If]. If «, v, and w are approximately constant over short times and spaces, then [if] is a linear second order differential equation. Several of the approximations lead- ing to [If] where consciously adapted to achieve linearity, so as to take advantage of a convenient property of linear equations, that two independent solutions of [if] may be added to one another to make a a third. Thus if we can derive an analytical solution to [if] for an infinitesimal puff, then we may superimpose puff solutions to achieve a general solution. This approach leads to the justly famous 'Gaussian7 models. While the derivation of a Gaussian solution to [if] is straight forward, it is also rather lengthy, and it would overbalance the present text. Interested readers should therefore consult other references [see for example: Pasquill and Smith., 1983; Longhetto, 1980]. For present purposes permit me to sketch the derivation, only, for the more common 'plume' version. A more complex 'puff version may be derived similarly, with somewhat fewer restrictions. ------- - 9 - 1. Assume a steady-state for times less than r. Temporal changes at times greater than r will be accounted for by specifying slower variations in the winds and diffusivities. 2. Neglect diffusion parallel to the wind: that is, set K, = 0. 3. Set dC /dt — 0 and seek a solution for the remaining partial derivatives in 2, p, *, of the form X(x)Y(y)Z{z). 4. Substituting XYZ into [if], and rearranging one may separate the partial differential equation into three total differential equations, one for x, one for y, and one for z. These have particular solutions of the form X(s)=Cle-b Y(y) = <7asin(nf y) + Cjcos(n, y) Z (*) = C^sii^n, x) + Cecos(nt x) 5. Set C3 and CA equal to zero, to force transverse symmetry. 6. Set the lateral and upper boundary conditions at infinity. 7. Evaluate the general solution by integrating over n, and n, with Fourier transforms. The expansion coefficients of these transforms are evaluated by assuming a profile for Y(y) and Z(x) at a specified z. 8. Choose Gaussian functions for these profiles. This choice takes advantage of a unique property of Fourier transforms, that the cosine transform of a Gaussian is a Gaussian also. [Other choices are possible, but most require solutions involving sums of infinite series.] Then if the real plume is a Gaus- sian at any x, it will be a Gaussian at all other x. 9. Evaluate k by conservation of mass. The result of these assumptions and manipulations can be cast in several closely related ways. Perhaps most quoted of these [Turner, D3., 1967] is: In [7] Q has a slightly different meaning from that in [if]: here it refers to a total current [kg/s]; there it referred to a current per unit volume. H is a plume center-line height above the surface (at r = 0); x as before is a downwind dis- tance from the source of Q; y is the transverse distance across the plume, and av ------- - 10- and a, are dispersion lengths equal respectively to \/2Kf s /o and y/2K, x /« . The coordinate z is measured upwards from the local surface, and the term in x +H accounts for perfect reflection of the plume from that surface. Usually ar and |