ANALOG COMPUTING TECHNIQUES APPLIED TO ATMOSPHERIC DIFFUSION: CONTINUOUS POINT SOURCE by Fred. V. k' E. Wendell ~ewson Technical Report No. 3 ORA Project 05652 NATIONAL SCIENCE FOUNDATION GRANT G-114o4 College of Engineering Institute of Science & Technology Department of Engineering Mechan'ics Great Lakes Research Division Meteorological Laboratories Special Report No. 12 TH UNIVERSITY OF MICHIGAN ANN ARBOR, MICHIGAN Jul 1961TI r-

%E2OA (AU4o 5-3K ACKNOWLEDGMENT The authors wish to acknowledge the wholehearted cooperation and assistance of Mrs. Anne. C. Rivette.

TABLE OF CONTENTS Page LIST OF ILLUSTRATIONS v ABSTRACT ix 1. INTRODUCTION 1 )N 2. FINITE DIFFERENCE MODEL 5. ANALOG MODELS iLS 4. COMPUTATION AND RESULTS 10 A. Isotropic Diffusion with an 84-Cell Model 10 B. Isotropic Diffusion with a 56-Cell Model 15 C. Nonisotropic Diffusion 19 D. Eddy Diffusivity Related to Plume Dimensions 22 E. Diffusion in Transitional States 2 F. Wind Direction Shear52 G. Large-Scale Eddies 4 H. A Possible Combination of Two of the Models 50 5. SUMMARY AND CONCLUSIONS.51 APPENDIX'55 REFERENCES 5

LIST OF ILLUSTRATIONS Table Page I Cell Boundaries in the Horizontal and the Vertical for the First Quadrant of the 84-Cell Model 12 A-I Comparison of Plume Concentrations for the Meandering Plume and the Straight Plume Figure 1 The 84-cell, finite difference model of isotropic diffusion from a point source. 6 2 The passive network analog of diffusion from a point source. 9 Computing circuit for the first quadrant of the 84-cell model. 1 4 Typical cell from the circuit of Fig. 3 showing the connections to the integrating amplifier. 14 5 Concentration plots for each cell as a function of distance downwind from a point source using the 84-cell model in isotropic conditions. 16 6 Concentration plots for each cell as a function of distance downwind from a point source using the 84-cell model in isotropic conditions. 17 7 Computing circuit for the first qiuadrant of the 56-cell model shown in (a). An integrating amplifier for a t'ypical cell is shown in (b). 18 8 Comparison of concentration plots for the source cell using the 84- and the 56-cell models with isotropic diffusion conditions. 20 9 Comparison of concentration plots for some cells using the 56-cell model. 21 10 Computing circuit for the first quadrant of the 56-cell model

LIST OF ILLUSTRATIONS (Continued) Figu~re Page 11 Computing circuit for the generation of the special forms of the coefficient of eddy diffusivity. 26 12 Plot of the modifying parameters of the coefficient of eddy diffuvisity. 28 13 Plot of the concentrations in the source cell as a function of the modifying parameter used. 29 14 Plot of the concentration in cell 1, 2 as a function of the modifying parameter used. 0 15 Plot of the concentration in cell 1, 3 as a function of the modifying parameter used. 1 16 The 25-cell, finite difference model of diffusion from a point source. o 17 A schematic of the wind shear component V as a function of height in the 25-cell finite difference model. 4 18 Computing circuits for the 25-cell model for the case with wind shear. 36 19 Cross section of the plume at X = 1 showing lines of constant concentration in arbitrary units. 57 20 Comparison of the concentration plots for cells to the right and the left of the source cell in the presence of wind shear. 58..21 Concentration plots for.some cells in the 25-cell model in the presence of wind she'ar. 59 22 Concentration plots for some cells in the 25-cell model in the presence of wind shear. 40 25 The 15-cell finite difference model of diffusion from a point source for use in the case of diffusion in the presence of large scale eddies. 42 24 Computing circuit for the 15-cell model including provision

LIST OF ILLUSTRATIONS (Concluded) Figure Pg 25 Computing circuit for the generation of the sinusoidal horizontal wind component. 45 26 Typical cell from the circuit of Fig. 23. 46 27 Plot of the concentrations for cells 0,2, 0,3, and 0,4 in the presence of large-scale eddies. 47 28 Plot of the concentration for cells -1,3 and +1,3 in the presence of large-scale eddies. 48 29 Plot of the concentrations for cells -1,2; +1,2; -1,4; and +14 in the presence of large-scale eddies. 49 0 Schematic of the fluctuating plume model showing the sinusoidal nature of the plume axis.

ABPSTRACT The analog computer has been used to obtain solutions to the diffusion;euation for the case of a continuous polnt source aloft. The basic model used is that of a source locate.d between an inversion and the ground where the wind speed and eddy diffusivity are constant over the diffusing region. Solutions were also ob.tained with several special features added to the basic model: an increase of lateral diffusivity over the vertical; representation of the eddy diffusivity as a function of plume dimensions or as a function of distance downstream; presence of wind direction shear in the diffusing layers; and the presence of large-scale eddies which cause the meandering of the plume. The solutions are presented as deviations from the reference solution, that solution derived from the basic model, to show graphically the effects of each parameter change. ix

1. INTRODUCTION Turbulent diffusion from an elevated point source is of primary concern in air-pollution meteorology. Although such diffusion is the central component of a common problem, it is by no means a well-understood phenomenon. This results in part from the lack of a sufficiently general mathematical model with a complete set of solutions. No matter how well diffusion is documented in the fieldY it cannot be completely and compactly expressed and understood except in reference to such a model. But atmospheric diffusion is such an intricate mechanism that no one model is ever likely to satisfy all damands. One model which is perhaps sufficiently general is the parabolic diffusion ecquation dt ~X_X'D C y y ~y 6z z - where dt dx yd, and X=average concentration x,y,z = space coordinates u,v,w = average wind components

This model has two main shortcomings: it yields only average concentration, and a complete set of solutions has not been found. In fact, such analytical solutions as are available apply only to special cases of limited application12 Another deficiency is that the model implies an infinite diffusing velocity, but in cases where this is important, the hyperbolic diffusion equation can be usedo3 In the case of steady-state diffusion, the velocity of diffuson is not important. The validity of the application of the parabolic diffusion equation to atmospheric diffusion has been questioned,5 but the lack of adequate solutions has made verification difficult. The use of electronic computing machines mak~ies solutions available for any type of parameter variation but here another difficulty arises. Since analytical solutions are not readily available, the functional dependence of some of the parameters has not been determined. For example, Richardson4 has suggested that the coefficient of eddy diffusivity be a function of plume size, but the exact form.of this function has not been indicated. The purpose of this paper is to show how the diffusilon equation may be implemented on the electronic analog computer and to show the qualitative nature of the solutions. Means for obtaining improved solution accuracy are also indicated-6 The reasons for using the analog computer in this application are that it'is easy to program and that it presents the solution in graph'ical form. The problems chosen for solution'..are based up.on the case of diffusion

of the ground and the base of an inversion. The special cases considered here were selected from the many problems of interest described elsewhere. 6-8 The methods outlined are perfectly general, however, and with an analog computer of sufficient capacity, a similar analysis may be made for nearly any meteorological conditions Added to the basic model are several functions of the average wind components and the coefficients of eddy diffusivity so as to consider increased lateral over vertical diffusivity, the case of a plume passing from one turbulent regime to another, the theoretical question of making the diffusivity a function of cloud size, the case of wind direction shear with height in the diffusing region, and the case of the meandering plume. The model is arranged so that the origin is at the source, the x-axis oriented downwind, the y-axis crosswind, and the z-axis vertical The basic model includes no provision for variations in the wind speed with height or in the eddy diffusivity in the y- or z-directions. This can be readily provided9 but the emphasis here is on special techniques. 2. FINITE DIFFERENCE MODEL The model presented in Eq. (1) can be simplified by deleting the terms that do not apply, to u ~x Ky Y2 - K _-z2~ (2) where the diffusion term in x was neglected in comparison to the transport

X + X as x,y,z + H H Kz z + 0 as z +for x> + where H is the height above ground of the base of the inversion layer These conditions state that the concentration is infinite at the source and that the ground and the inversion are impervious to the diffusing substance. The continuity condition is +H %-c 2 u x(x,y,z) dz dy -00 2 for x. O where is the rate of emission from the source. It is convenient to state (2) in nondimensional form: X K ~Y2 ~Z2 6 where S = X/x'o Y = y/H Z = Z/H x = H u x/Kz V = Hv7/Kz and the symbol X is reserved for later use. The constant will be de

The mathematical model to be used, Eq. (4), must be put in a form suitble for programming on the analog computer. Since the computer can integrate only with respect to one independent variable, computer time, two of the inependent variables must be replaced with finite differences. The derivaives in Y and Z will be expressed as finite differences as followso a2S 1 Sml Sm Sm - S ml S Sm+1 - Sm-l (5) )y2 5 (AY)m (AY)m (AY)mlj Y 2(AY)m m m 62S Sn+l - 2Sn + Sn-1 6Z2 2 (AZ)n2 n he expressions for the second derivatives in Y and Z are different because, in all the cases studied here, uniform differences are taken in Z but nonuniform differences are taken in Y o If uniform differences were taken both ways, the two expressions would be identical. Equation (4) may now be stated as a set simultaneous, ordinary differential equations: dS Ky 1 SM.n-Mn Sm,n- Smdx - Kz (Ym (AYl,n m-ln m,n 2 2 + Sm,n+l 2Srn,n +SM nl- S m+1 n -Sm..1f (6 (AZ)n2 2(AY)m where the subscripts m and n refer to discrete values in Y and Z,9 respectively. The physical model corresponding to Eq. (6) is shown in Fig.

z H INVERSION.500 J,428 C0 357 to.286.214 to.143 \N I 714 — ~~~X Xo 0 cI XO~. 0714.1934 C GROUND m = I 2 3 Fig. 1. The 84-cell, finite difference model of isotropic diffusion from a point source. 6

array of cells of rectangular cross section with width AY, height AZ and with one end at x = 0 and the other at x = Xo The cells are stacked vertically from the ground to the inversion base and horizontally from - to + The number of cells in Y is finite since the end cells have width AY = co The source cell is in the center of the array and is designated 0,0. All the properties of each cell are lumped in the center of the cross section sothat the model is an array of rods with one end at x = 0 and the other at x = o The source is no longer a point source but a vertical area source. In this model, the continuity Eq. (3) holds for all values of and by evaluating it at x = 0, we can define o = Y2 Z2'-Q = / ~ H2 u S Xo dZ dY Y1 Z1 The integration at x= 0 is carried over the source cell only, and in a given cell all the pa-rameters in the integrand are constant. Since the value of S at x = 0 in the source set at unity, Q - y~0 (AY)O (AZ). and therefore =o H 2 (Ay)o (AZ)0 The solutions that will be obtained with this model willl be the average con

3. ANALOG MODELS In analog computation, it is helpful to realize that there exists a correspondence between the components of the circuit set up to solve Eq. (6) and the physical model shown in Fig. 1. This can be shown by constructing on paper a passive network analog as shown in Fig. 2. In this circuit, there is a capacitor for each cell and a network of resistors connecting them. In this analog, the voltage at each node is a function of time and corresponds to the concentration as a function of distance downstream. The current into each capacitor is given by dV dX i = c- = u Ay Az.dt dx ~~and ~c = u Ay Az A y _Az By K Az'R - K Ay y In this circuit the boundary conditions are implemented directly. Since there must be no flux through the ground or the inversion., the, counterpart of flux, current, is made zero by simply making no connections across the boundaries.'The initial concentration (voltage) in the source cell is provided by disconnecting the capaci-tor for that cell., charging it to arbitrary voltage, and switching it back into the network at x (time) equals zero. Since the end cells'in y are infinite, the concentration (voltage) is always zero there so these nodes are connected to ground

INVERSION S- s s~~~~~~ I I S-.0s l I l S_, I I ~~~~~~~~~~~~~~~~~I I - I I So,o I'ZR I _?I I_ ~ ~ ~~s I 1i I'I _Is I I I IO,I 1~~~~~~ I I — I I I'I H_ — >! I I~~ ~~~~~~~~~~~~~~~~~~ I_ ~~-l,0 I /V - I ~~~~~~~~~~~~~~~~~~~I I GROUND Fig. 2. The passive network analog of diffusion fronm a point source. The center cell is the source cell and its capacitor is initially changed to represent the source concentration. <:? I - - I~~

sistors and capacitors involved. If any of these parameters are a function of x, then they must be made to change in time. This is not readily accomplished and constitutes one of the disadvantages of the passive network. In general, the passive network is much less flexible than the active circuit. An electronic analog circuit will be shown later for each case. They are constructed by assigning to each component a mathematical operation and these are connected so as to satisfy the governing set of differential equations. Identical circuits can be derived by simulating the passive network analog; thus each circuit has a dual nature. It is a series of mathematical operations designed to solve an equation and it is also an analogy of the system being studied where components and voltages represent physical parts of the atmospheric diffusion model. Both the active and the passive analogs are physical systems designed to behave like another., idealized physical system. The electronic analog circuits will be explored as they appear. 4. COMPUTATION AND RESULTS Since,each of the-cases considered is qluite different., the f orm of Eq. (6) appropriate to each will be given along with the computing circuit and the results. A. ISOTROPIC DIFFUSION WITH AN 84-CELL MODEL T'he fPi-rst ncaselI (is diffusQion f'ro-m thepint sorc -in isotropic conI

the coordinate axis so only the cells in the first quadrant need be simulated and of these, the ones in the outer column, m = 3, are always zero so that there are only 14 active cells in the simulation. In this case, Eq (6) becomes dS 1 Sm+l mn S m,,n - Sm-l,n dX m n 196 (AY)m (AY)m +1 (AY)m 1_l m,n 2 -2 Sm,n+l - 2Sm,n + Sm,n-l 7 - (7 196 (AZ)2 where X = 196 x and the boundary conditions are S l,n = S-l,n Sm,1 = Sm, 1 s3,n s The factor 196 was inserted to make the coefficients more convenient. The cell boundaries given in Table I show that the spacing in the vertical is linear whereas the spacing in the horizontal is exponential and the source cell has a square cross section.

TABLE I CELL BOUNDARIES IN THE HORIZONTAL AND THE VERTICAL FOR THE FIRST QUADRANT OF THE 84-CELL MODEL Y = y/H Z = z/H O 0 0.0714 0.0714 0.1934 0.1428 00 0.2142 0.2856 0.3570 O.4284 O. 5000 Since AZ is a constant for each cell, Eq.o (7) can be conveniently expressed as two sets of difference-differential equations for the columns used. ldX0.8 28(52n - Sl,,n) + (Sl.,n+l- 2S51,n + Sl,n-l) dS o2n.4957 S1,n - 0.7410 s2, + (S2,n+1 - 2S2,,n + S2 n-1.) (8) The computer circuit for Eq. (8) which is shown in Fig. 5 has one amplifier for each of the 14 differential equation,s represented by (8). The top row of amplifiers simulates the first column of cells and the top equation in (8). The bottom equatilon is solved by the lower row of amplifiers. The Sl, 1 amplifier in the top row represents the source cell and has an initial condition

+ 0~~~~~~~ o & -o co U) U)~~~~~~~~~~~~~~~4)Q rd) 00 IH rd r —0 0 0b ao U)~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~fr rd rd + 0)~~~~~~~~~~~~~~~~~~~~~~~44-) 4-'W CH~ co 4-P) o 0' -co CO U) U)b0 2 4-. 10 U)CI) E-i C UDO 0 0 10 M~~~~~~~~~~~~~4HU)~ ~ H-

-Sm-l,n IL__4 + Sm, n - Sm,n-I -Sm, n+l * - I /if ImI HIGH-GAIN d-c AMPLIFIER Fi.4 yialcl rmte ici fFi.5soigth oncin

Solution curves for the first row of cells are shown in Figs. 5 and 6. B. ISOTROPIC DIFFUSION WITH A 56-CELL MODEL To facilitate the simulation, the number of cells can be reduced. In this case, only 4 cells in Y are to be used instead of six. The purpose of this reduction is that more involved problems to be treated are more readily handled with a simpler basic model. Not so much equipment is needed in this case but the reduction in the number of cells imposes a severe accuracy penalty. Before introducing the more interesting cases, the magnitude of the error introduced can be shown by solving the same problem employing the same techniques. Cell boundaries in Y for the 56-cell model are Y = 0, 0.0714, while the cell boundaries in Z are as stated in Table I. In the first quadrant, there are still 7 active cells in Z but only 1 in Y; thus there is 1 active column and Eq. (8) simplifies to ldnX -0.7220 S1,n + (Sl,n+l - 2Sl1,n + Sil,n-1) (9) where X= 196 -x and the boundary conditions are S S 5 1,91(0) =1.00 n =n 0 Rl7=~,

1.00 0.90 x1 96K x 0.80 - _ 0.70 0.60 0,) z 0 0.50 z w z 0 00.40 0.300.20 0 DISTANCE DOWNWIND. X

0.20 0.18,2 SZ X/X0 0.16 196Kg uH2 0.14 ) 0.12 z 0 a: 0.10 z w 0 z 0 o 0.08 0 24 14 68121 DISTANCE DOWNWIND X

e~~~~~~~~~~~~~~~~~~~~~~~~~~ I -I UT~~~~~~~~~~~~~~~~~,-v~~~~~~~~~~~~~~~~~E o (I)~~~~ rA ) S~~~~~~~ ~~~ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 51(0') j+ 0 4-1 0 4-A'') - ~ ~ ~ - 0 0 0 CH CH 0 b.O C-p..~4.r-I +~~~~~~~~~~~~~~~~~~~~~ b.D 0 C~~~~~~~~~~~~~~~~~ T CH — -U:'H Or-. 0 co 0~~~~~~~~~~~~~~~~H0~~~~~~~~~~~~~~0~

one integrating amplifier for each active cell. The solutions for the source cell obtained from the 56-cell and the 84-cell models are compared in Fig. 8. Evidently the approximations are not bad for small values of X For all the following work, the absolute error for values of X is large, the objective is to show deviations from the reference solutions in each case. ~~~C. N~ONISOTROPIC DIFFUSION Changing the ratio of the eddy diffusivity in the Y-direction to that in the Z-direction produces a very marked alteration in the solutions. As )an example, set Kz in the 56-cell model. Equation (9) may now be written ~~~dSl~ n+l - 2Sl,n + Sl,n-1l - 0~7220 Sl,n + 10 1,n where X= 1960 -x and the boundary conditions are l,1n -1l,n 1,1 1 1 S 1,1(o) =1.00 1,97 1,8 The-computer circuit for Eq. (10) is identical to the one shown in Fig. 7

1.00 0.80 0.60 0.40 /96KX 0.20 LU 0'n 4 - CELL MODEL z ) 0.10 0.08 0.06 2 56- CELL MOE 0 0.04 z 0 0.0 0.021 000 1 2 34 56 DISTANCE DOWNWIND X Fig. 8. Comparison of concentration plots for the source cell using the 84- and the 56-cell models with isotropic diffusion conditions.

I.00 1, / 1.10 C,)~~~~~~~~~~~~~~~~~~~~Oo o 0.01 z I z 0 0.001 0.00OI0 0.1 0.2 0.3 0.4 0.5 060 7

has a pronounced effect. The magnitude of the effect is enhanced by the limited vertical diffusion due to the presence of the inversion. D. EDDY DIFFUSIVITY RELATED TO PLUME DIMENSIONS Plume dimensions may influence the magnitude of the effective eddy diffusivity acting on the plume since, as the cloud grows, larger eddies are able to operate to diffuse the materials. The eddy diffusivity might be made a function of the plume size and of the power spectrum of eddies present. In the usual parabolic diffusion equation, the plume has no easily efined boundaries except those imposed upon it, such as the surface of the ground. Therefore the diffusivity cannot be related directly to plume size but must be specified in some arbitrary manner. For example, it could be made proportional to the standard deviation of the concentration distribution taken perpendicular to the plume axis. If this distribution is normal, the standard deviation is inversely proportional to —the central maxima or axial concentration, so that the diffusivity can be set inversely proportional to the axial plume concentration. Since the power spectra would never be continuous indefinitely, the diffusivity could not become infinite as the axial concentrat'ion went to zero but would reach some finite value and remain constant. In the finite difference model., the source is finite so that the diffusivity wou-ld increase from one small finite value to a larger one inversely as the concentratilon in the cell containing the source. Represent the diffusivity as O(S1)K where K is the diffusivity for infilnite plume dimen

For example, set 0(S) = 1 0.7 S (1) since 1.00 S,1 1 >- 0 Setting Ky = Kz, Eq. (9) becomes dX lS n = (Sl)(Sl,n+l - 2S1,n + Sl,n-l dX 1., n 1,~ n 1,n-1-0720OS),n12 with the same boundary conditions as before. The computing circuit, shown partially in Fig. 10, is similar to that of Fig.. 7 in principle. This one uses 3 amplifiers, one integrator and two summers, per cell to isolate all the inputs to the integrating amplifier. With this configuration, it is only necessary to introduce the modifying parameter (S1) once per cell. Since O(S1) is a function of the axial concentration, its form is determined by a feedback look defined.by Eqs. (11) and ('15). The function Sl is generated as the solu'tion proceeds. The computing circuits used in this problem Are virtually identical with those of the next section. Therefore the description of the additional components Idequired for the generation of O(l and the presentati on of the solutions will be deferred to that section. E. DIFFUSION IN TRANSITIONAL STATES A problem of special interest arises when the trajectory of the diffusing plume crosses the boundary from one tu-rbulent reg ime to another. The physical

the thermal turbulence decrease markedly as it does so. The associated diffusion is referred to in this research as diffusion in a transitional state; the inversion breakup fumigationio is the result of a corresponding time transition in diffusion. In the mathematical model, a change in the eddy diffusivity accompanies the transition from one turbulent regime to the other. This will be implemented lby applying another modifying parameter 0(X) to the constant K Then 0(X) will be unity in the first regime and fall to some lesser value in the second. For the present preliminary study, this function was simulated directly on the computer to conform to an inltuitive impression of what it should be without using a mathematical expression. It is possible to contrive a mathematical description, which is: O(X) =0,52 + 0.4~8{%[ - h(X -Xl + cos(X- Xl) FhX - Xl) - h(X -l - (15) where h(a) is the unit step function which has the following properties: h(a) =0, a< 0 -1, a~>0, and X is the value of X where the turbulence starts to decrease. The diffusion equation is just Eq. (12) with O(S1) replaced with 0(X) The computing circuit is also that of Fig. 10 and the circuit for generation

-S 1,2 -Si, 1.72-2 O~~~~~ I INTEGRATING AMPLIFIER +S, 2 -S +-S, 22/ | ~~~~~SUMMING AMPLIFIER S 14 1,3 -S,3 +SI, 5 1.4 _~~f~ SI.4 2.722 i +,+S +SI, 5,5S1, 2.7922 +S~~~~~~~~~~~5 2.722 C) ll~~~~~~ ~(,X) ss DsSERVO Si, 1,7S IU MULTIPLIER 1.722 I MECHANICAL LINKAGE TO FOLLOW-UP POTENTIOMETERS Fig. 10. Computing circuit for the first quadrant of the 56-cell model which permits the inclusion of special forms of the coefficient of eddy diffusivity. This configuraton permitted te generation o solutions whee the diffusio coefficients wre a function f plume size ad in the case f diffLision i tran~~~~~~~~~~~~~~~~~~~Sitoa tts h ote,ie niaemehnclcnetos 25~~~~~S,

0 0 ~ ~ ~ ~ ~ 0 o(I /) ~~p LL~ ~ ~ ~~~ix *H r~ 4-) 0 0 =__ -.- "'"0'.0, 0'HO \ < O~~~rd LLJ P-1 ~~LL 0 Ic I o, 0> / o - ~~~~~~~~~~~~~~~~~~-P *H o 00 0 OCH 0co [9i7{7~~~~~~~~~~~~~~~~~~~~~~~4)Q $:ird4 F- ~'qH~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~L- Q 0 000 0 r-l0~~~~~~~ ~. 5' ~0 o rd 0 0~~~~~ LJLJ 0 +-. 26 ~l 0~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~=. "r) i~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ — $:i IL ~~~~~~~~~~~-

are different despite a somewhat deceptive similarity in the governing equation used. The differences between the problems are indicated in the form of the modifying parameter 0 o It is only because of the similar way in which 0 is implemented that the problems are handled together. It is possible to introduce ooth effects in the problems. In this case, a new parameter will be defined $uch that 0(S1,X) = 0(S1) 0(X) o Figure 11 shows how each of three functions is generated. When switch 1 is on -100 and switch 2 is closed, the outpuit of amplifier 4 is 0(S1) as defined by Eq. (11) and it passes through servomultiplier 1 unchanged to control servomultiplier 2. When switch 1 is on 0(X) and switch 2 is open, the function 0(X), generated by amplifiers 1., 2, and 3, controls servomultiplier 2. The function 0(S1,X) is generated when switch 1 is on 0(X) and switch 2 is closed. The functions 0(S1), O(X) and 0(S1,X) are plotted in Fig. 12 as they were generated by the computler. Solutions for S1.,1 s1yl2'and S1,5 are presented in Figs.115, 14., and 15, respectively, for 4 cases: 0 = 1., 0 = 0(51), 0 = 0(X) and 0 = 0(S1,X). The case0=I is the reference solution. Each modifying parameter acts to change the scale in X, that is~, the new solutions can be found by translation in X from the reference solLution. This is shown best in Figs. 14 and 15 where the maxima have been trailslated and expanded but the value of the.maxima remains the same. The firot modifying parameter 0(S1 is significantly different from unity only i~or small values of X so the greatest effect on-the solutions should occ~ur there. Each solution shows large t-rans

1.00_.90-.80 -.70U.60U ar-.50 z LL >~~~~~~~~~~ I 0 I 2 3 4 5 6 7 DISTANCE DOWNWIND X fusiity.~(s1 isa fuctio of lum ( siz wil (x) ersnstas ito rmoetruetrgm o ante.40s,)i iml utpi

1.00 \.90.80.70 _.60 z o.0,.~ - 2 - z~~~~~~~~~~~~~~S X c rw 0 ~ ~ ~ ~ ()__, — z ~ ~ ~ ~ () 64 0~~~~~~~~~~~~ O~~~~~~~~~Is 0 D(SiAC DONWN X Fi.13. P lo of h cnrtOS i h orc ela fnto f the~~~~~ moifin p=araetrwsd ~L ty) 9

.18r.16 -.14..12 (I, cj z 0.10z w 0 _~.08 ).(s),X).06(Sd).04.02 I IIIII 0 I 2 3 4 5 6 7 DISTANCE DOWNWIND X

.10-.09.08.07 u).06z 0 H z.05 z 0.04.03.02-.01 0 I 2 3 4 5 6 7 DISTANCE DOWNWIND X

along the plume axis and to delay the occurrence of the maxima off the axis. The solutions for 0(1) and 0(Sl) do not converge for large X but reach a constant translation in X. But since the gradient is small for large X the absolute value of the differences between the solutions is small, Since the function 0(X) can be determined independently, the solution could be found from the reference solution by plotting it on a new scale defined by X X* = 0(x) dX o Jox The observed effect corresponds to this point of view in that 0(X) has no effect on the solution for X Z 1 but acts to translate the solution for X >1. F. WIND DIRECTION SHEAR The presence of wind direction shear in the diffusing layers causes the plume cross section to be skewed, and the amount of skewness can be determined by inserting terms in the diffusion equation for cross-wind components. Equation (6) contains the nondimensional parameter V = H 7/K for these cross-wind components. Since the diffusion pattern is not symmetrical,. it will not suffice to simulate only the first quadrant. In the interest of simplicity,, a new 25-cell model is used as shown in Fig. 16. The source cell is 0,1 and there are 15 active cells. Figure 17 shows VJ(Z), which is an exponential function of height, in the finite difference model in relation

+ + ~ + + +1 0 0 0.3 4d - - - - - U)~~~~ — + + ~ + + +.rd 0) 0c CH CH Ord - CH0 0 0 0 0 0 0~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 0 (0 N ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~rd Hrd CH~ 0) 3rd Ur. U. 4VbE-' X! *H0 0'C) N -~~~~~~~~~~~~~~~~~~~~~~D U. U. U. U.~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

x~~~~ 0" 0 0 - 0 INC~~~~~~C 4-) 0 rOH a ~rd 0 0 Cd~ ~rd'H -H N IC) r~~~~~~~~~~~~~~~) ~ CjCd LO r1r) ~ ~ ~ ~ ~ 1) co C\) 0~~~~~~~~ld b'H H-P

written dS |_ 1 S m+l,n-sm Sm,in dXm n 25(AY)m (AY)m+ (AY)ml 50 (/\Y)m LSm+1,n - Sm-l,nJ + LSmn+1 2Sm,n + SminJ (14) S2,n ~0 SmwQ =Sml 1 Sm5 = Sm,6 so 3() = 1.00 The computing curcuit is shown in Fig. 18 and the solutions are in Figs. 20, 21, and 22. Each solution is compared to the reference solution, where V(Z) = 0.We can compare the concentrations in the cells for two cases at some distance downstream., say at X = 1.* The results, as shown in Fig. 19, indicate that the wind direction shear produces a marked skewness in the plume. The skewed plume is to be expected but a more interesting question is whether the average concentrations in the plume are reduced because of the skewness. Figures 20, 21, and 22 indicate that the concentrations-are-reduced, *In making this comparison, note that changes in wind direction were provided by adding v components while holding u constant which increased the mag+nite ln~ of the tota vecto^ —-I% r ban a%- Tr= ve ragerz -1.6%. ( - V SinceXI is a funcion- of __

Sn U, /' + 2A2I Sn Sn _ Sn U) c'H + + 0.'~~~~~~~~1 G)4 U, + + + 4- I, —. 0 1* - 0 (V-C U) Cf) Cf) cn d~~~~~U + + I 0 Ul) 6f' t,/'U) 0))t + ie Li\ 01).4Uf) Uf) U) U/) + + Ca) 4-" U)P C) ) o ~'-H N 9~~~~~~~~~~~~~~~4-CM...7 N - N~~~~~~~~~~~~~~~~~~~~~~~~- ci~ o ~~~~~~~~~~~~~~~0 I f U) 6 _ U)+ 0' co I 36

__-I 0 +1 n 5 g 9~~~~~~~~~~~ ii 4 3 2 39 (a) WITHOUT WIND SHEAR -I 0 +I i\ i\i i ii4 Fig. 19 Crs eto ofth9p0m atX=1soiglie fcntn concntrtio inarbirar unts.Thi fiuesosterfrneslto 57~~~~

0.020 0.018 0.016 0.014 0.012 z 0 n~ 0.010 z w z 0.008 \ 0.006 0.004 0.002 0 I 0 2 4 6 8 10 12 14 DISTANCE DOWNWIND X

0.020 0.018 0.016 0.014 0.012 zO0tw o0 w z + L)00D08 0.004 0.002 0 2 4 6 8 10 12 14 DISTANCE DOWNWIND X Fig.oo -i2. Concentration posfr mecl inth 2 m i the prsneofwn har h oidlnsae eeec sltos

0.020 0.018 0.016 0.014 c 0.012 z c 0.010 / ~ L #1,1 Z I N' wN z I*0.' 00.008 I/. -/,5 N 0.006 N~~/00 0.004 1.i Q00 2 1 21 0~~~~~~~~4

by about 4% in the center of the plume at X = 1. The integrated concentration over the plume at X = 1 was reduced about 2%. Thus it seems valid to conclude that wind direction shear does increase the diffusing propetties of the atmosphere. It is planned to analyze the influence of wind directon shear in concentrations downwind from a continuous point source using an analog computer of much greater capacity. G. LARGE-SCALE EDDIES 3S The presence of large-scale eddies or, since we are considering the steady state, a standing wave in the horizontal wind components can be represented as a sinusoid v = a sin bx. A stationary plume might be an example of this kind of process. The small-scale eddies which operate to diffuse the plume are represented by the coefficient of eddy diffusivity, whereas the large-scale eddies act only to transport the plume in the y-direction. Equipment limitations at the time of the research dictated the choice of a more coarse finite difference model-as shown in Fig. 25. The mathematical model., based on Eqj. (6), is dS -= 5m+l.9n + 5m-1l,7n + 5m,Jn+l + SM,n-1 -S dX m,n VXX) [Sm+l,n - Sm-1,](5 Ha where V(X) = -7 sin BX and the boundary conditions.are

0 0) _ __ (D 0(0 ~~~~~CH I~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~C+ 00 I ~~~~~~~~~~~~~~~~~~~~~~~~~~~~rI(12 ro ~0 rd - 1 ________ rE') 00 rf) 1 0~~~~~~~~~~~~~~~~~~~~~~~~~~~1 N ~~~~~~~~~~~~~~I ~~~ — 0- — 0 —~ ------ G..~..0-0-(0 0 P, —I 0 (D I rf) -.' (0 r() Yo'H o 0~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~1 d ~ ~ ~~~~ ~~~~~~~~~~~~~~~~~~~~~ 4-c),,o~~~~~~~~~c ~~~~~~1 q~~~~~~~~~~~~H, 0 ~~~~~~I'HO I H 0 0 0 - ~~~~~~~~I r4' ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0 tOir 0~~~~~~~~~~~~~c I N > ~~~~~~~~~~~~~~~~~~~~~C\i a() E z I b0 rd I'HO~~~~~~~~ — rd ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~pi -1riII w E tp I =~

Sm 1 = Sm,2 SmE5 = Sm,4 and AY = AZ for the active cells. If we set 5A = Ha(AZ)/2K, (15) becomes dS I1 F Smn- 4Sm,n + (1 - 5A sin BX) Sm+l dX m,n 10 nn + (1-+ 5A sin BX) Sm-l,n] (16) where X = 90 x. The computing circuit is shown in Figs. 24 and 25, and a typical cell is shown in Fig. 26. The solutions are plotted in Figs. 27, 28, and 29 along with the reference solution. This system exhibits a phenomena like resonance in a forced second-order system. It is not pure resonance but the sytem does respond markedly to changes in frequency of the cross-wind component. Changes in amplitude do not affect the qualitative nature of the solution. The plume meander lags behind the cross-wind components. The point at which the plume axis crosses the middle cell can be found by noting the point where the concentration of opposite pairs of cells is equal. This occurs for X = 1.506, 2.448, and.531, while the cross-wind goes to zero at X = 1.047, 2.094, and 3.141. The phase lage is not constant for the above three points but is 78.9 deg, 60.9 deg, and 67.1 deg, respectively. It is of interest to compare the concentrations in the plume, for sampling periods short enough to avoid plume shifting, when it meanders with when it d.oes not. At an arbitrary point d.ownwind., at X = 1.6, the total concentraion reported. in the case with no meand.er was 161o higher than in the case

x x CSX rn Go X ~~~~~~~~<~ 0 0 0 0 c c C\J C\j C U) ii I -~~~~~~~~~~~~U)V II I ~~~~~~~~~~~~~0 I CH + I I I 0 _ ~ ~ ~~I Un~ w ~I I'I Ii ~~~~~~~~~~~~~0 l lI rd l l to U) l)+ I + U) 0~~~~~ ~~~~~~~~~ I 0) U A-''-I+ X - U X II If\ 0 ii~~~~~~~~~~~~~ I, I ~~~~~~~~~~~~~~~~~~~~~rd 0H to - Ito ~~~~~~~~~~~~~~ L..... UH~~~~~~~~~tL U) I.\'U)\I + II~~~~~~~~~~~~~~~~~ I I~~~~~~~~~~~~~~~~~~~~~~~~~~~0C II bO~~~~~~~~~C- Ii *Hc~~~~5 II~~~~~~~~~~~~~\ -

cr. LiJ -4 t5 ~x 5 m 0 co~~~~ cr~~~~~~~~~~c >u < 0 LLJ N~ H O o cd.r-I 0 0!X x ~ xN co Q Q 0 0 U) ~~~~~~~~~~~~~~~0 4-, 133 c~~~~~x + Q~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~) to 40) 0 0 4-= H F-i LuL 0.rI d H.H I~~~~~~~ 0 2 2 F __ <: L):J~~~~~~~~~~~L

Zc~ E C -I I4 o~1 t: 0 ZL 0 Al U.H'c -P ~~~~~LL *~~~~~~~~~~~~~~H btD z~~~ 0H!Zl <~~~~~~~~~~~~~~~~~~~F~* C: ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~~~~~CU + 0~~~~~~~~~~~~~~~~ E

1.00 ~~\ V(X) S/N 3X - ~~\ 0,3 0.10' Z _! C) z 0 0, 0.0 0.001 0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 DISTANCE DOWNWIND X

1.00 V(X) =SN 3K 0.10 z 0~~~~~~~~~~~~~~~~~10 II —, Iz o _ \,.i o #1,3'/,3 /~~~~~~~~~~~ 0.01 J 0.001IIIIII 0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 DISTANCE DOWNWIND X Fi.28 lo f'te ocetato frcel -, ad+15inte rsec oflrg-cleedes hesld ieiste eeenesluin

1.00 V(X) = S/N 3X 0.10 -,;-, z 0 l'I. —. z w z 0 0 0.01 0.001IIII I 0 0.5 1.0 1.5 2.0 2.5 I3.0 3.5 DISTANCE DOWNWIND X fi.2.Pltothcocnrtosfrcls-,;+2;l.1.; n. +) intepeec_flre-cl_ddis h oid iei h eeec solution. ~ ~ ~ ~ ~ V() -.)~~~~~~~~~~~~~~4

where the plume meanders. As shown in the appendix, this decrease in concentration in the meandering plume is due entirely to the longer path length followed by the plume. If one knew the total concentration across a straight plume as a function of distance downwind. one could then compute the long-term concentration at a station subjected to a plume in meandering but otherwise identical conditions. The calculations may be so tedious as to be impractical since they would involve computing the length of the trajectories over the given interval, adjusting the concentration and integrating them. It would still be of interest to perform this calculation on the computer using random large-scale eddies, and it is intended that this topic will be the subject of a future report. H. A POSSIBLE COMBINATION OF TWO OF THE MODELS Two of the models presented above are related to some extent. The model in which eddy diffusivity increases with the plume dimensions can be related to that including large-scale eddies. In the first case, the complete spectrum of eddies can be included and the effect on the plume partially determined. The solution is partially complete in that, while the diffusing effects of all the eddies are included, the transport due to eddies much larger than the plume dimensions is neglected. The model simply ignores eddies of a given scale until the plume dimensions grow to that scale. The transport phenomenon is shown in the latter model where the large-scale eddies are

A combination of the two models might be an improvement. Let the wind components represent eddies larger than the plume while all smaller eddies act through the eddy diffusivity. The wind components would become smaller with distance downstream whereas the eddy diffusivity would increase. The wind components incorporated in the model are defined as averages which raises the question of how much they should be allowed to fluctuate or to possess time or distance-variant properties. An averaged wind speed could be defined as the average or mean value obtained by the use of a sampling time just large enough to exclude all the turbulent eddies which contribute to the diffusion process. If we allow the assumption that the scale of eddies involved depends upon the scale of the plume, then the sampling time in the above definition would depend upon the plume size also. The wind-speed components could be represented as a Fourier series whose higher-frequency components are gradually eliminated as the plume grows., and conversely, the coefficient of eddy diffusilvity would increase as more turbulent eddies were added to its domain. This is a rather involved model even for analog computation and no attempt, as yet has been made to set it up. 5.SUMMARY AND CONCLUSIONS As the model has become more complicated, it has been necessary to use fewer cells and to tolerate greater error. The amount of equipment needed is proportional to the number of cells and to the complexity of the model,

is to use a larger machine which would permit simulation of more cells. The errors are proportional to the fourth partial space derivative so that the error can be reduced by solving for a function related to the solution which possesses a smaller fourth derivative. Such a function could be found by perturbation techniques. The parabolic diffusion equation has been shown to be of general utility in the variety of problems covered. It is a model which is well adapted to analog computation. With sufficient equipment the actual programming seldom presents any difficulty even in the most involved problems discussed here. The trapping situation was used throughout this report to limit the scope of the research,and to facilitate handling the problems on a computer of given size. Extension of these methods to nontrapping situations can be readily done using the technique reported earlier.1

APPENDIX In the case where the meandering plume was treated, there were two possible mechanisms which could affect the diffusion of the plume. The one noted in the main body of the report was the increased path length of the trajectory. The other is that large-scale eddies might change the shape of the plume, thereby influencing the diffusion process. An attempt was made to evaluate this mechanism by correcting for the first one. When there is a cross-wind which causes meandering, the magnitude of the total wind vector and the path length along the plume axis, as shown in Fig., has been increased. To compensate for this, the concentration for no cross-wind should be taken at a point further downstream in inverse proportion to the relative decrease in the magnitude of the wind vector. Taking the first zero crossing., at X, = 1.506, the compensating factor is =2 xi BXin where sin2 BX - -$ sin2 BX dX = 7[X1 - sin BX1 cox BXj 0. 4790 Then X2 = 1.506./T.47 = 1.851. The following table gives the concentration in both cases.

t t \s \. \) o Z 0.4/ J /,U0 /1/ / * ~ / +.r-' / d - o ~~~~~~~~~~~~~~~0 /~~ Q. c o 54~~~~~~~~~~P i/i// CH /~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'' Ht /~~~~~~~~~~~~~~~~~~m Ocd i,~~~I ~~~~~~~~~~~~~bA'o 5".

TABLE A-I COMPARISON OF PLUME CONCENTRATIONS FOR THE MEANDERING PLUME AND THE STRAIGHT PLUME Cell Straight Plume Meandering Plume x = 1.831 x = 1.56 O, 2 0. 0595 0.0599 0,3 0.0596 0.0617 0,4 0.0595 0.0599 -1,2 0.0412 0.o406 -1, 3 o. 0422 o. 0420 -1, 4 o 0O412 o.o406 +1,2 0.0412 0.0406 +1, 3 o.0422 o.0420 +1, 4 o 0.0412 o. 0406 Total concentration o. 4278 0. 4279 Evidently the plume shape is a little different since the individual cells do not report identical concentrations, but the sume of the concentrations for all the cells is not significantly different.

REFERENCES G1. Sutton, 0 G., Micrometeorology, New York, McGraw-Hill Book Co., Inc., 19575. 2. Godson, W. L., "The diffusion of particulatematter from an elevated source Arch. Meteorolo, Geophys. und Bioklimatol., Serie A, 10, No. 4 (1958).. Goldstein, S., "On diffusion by discontinuous movements and on the telegraph equation," Quart. JO Mech. and Appl. Mech., 4:129-56 (1951). 4. Richardson, L. F., "Atmospheric diffusion shown on a distance-neighbour graph" Proc. R Soc., London, Series A, 110:709-737 (1926). 5. Gifford, Frank, Jr., "Statistical properties of a fluctuating plume dispersion model," in: Advances in Geophysics, Vol. 6, New York, Academic Press, PP. 117-138, 1958. 6. Hewson, E. Wendell, "Meteorological factors affecting causes and controls of air pollution," J. Air Pollution Control Assoc., 5:2-8 (1956). 7. Hewson, E. Wendell, "Atmospheric pollution," in: Compendium of Meteorology, s., The American Meteorological Society, pp. 119 8. Hewson, E. Wendell, "Industrial air poollution meteorology,," in: Vol. III., Industrial Hygiene and Toxicology., New York, Interscience Publishers, in press. 9. Brock, F. V., Analog Computing Techniques Applied to Atmospheric Diffusion: Infinite Line Source,1 Un'iv. of Mich. ORA Report 03632-2-T., Ann Arbor, 1961. 10. Hewson., E. Wendell, "The meteorological control of atmospheric pollution by heavy industry,," Quart. J. 1~oy Meteor. Soc., 711:266-282 (1945).

UNIVERSITY OF MICHIGAN 3 901 5 0252611110111 3 9015 02526 1 150