ANALOG COMPUTING TECHNIQUES APPLIED TO ATMOSPHERIC DIFFUSION: CONTINUOUS LINE SOURCE by Fred V. Brock Technical Report No. 2 ORA Project 03632.NATIONAL SCIENCE FOUNDATI)ON GRANT G-11404 College of Engineering Institute of Science & Technology Department of Engineering Mechanics Great Lakes Research Division Meteorological Laboratories Special Report No. 11 THE UNIVERSITY OF MICHIGAN ANN ARBOR, MICHIGAN June 1961

-e r>) C tAi U - ~10r9 )Zn~,

ACKNOWLEDGMENTS The author wishes to acknowledge the many helpful suggestions and encouragement given by Professor E. Wendell Hewson. The wholehearted cooperation and assistance of Mrs. Anne C. Rivette is also appreciated. iii

TABLE OF CONTENTS Page LIST OF FIGURES vii ABSTRACT ix 1. INTRODUCTION 1 2, THE PARABOLIC DIFFUSION EQUATION 3 3. VARIOUS MATHEMATICAL-PHYSICAL MODELS 7 4. PREPARATION FOR ANALOG COMPUTATION 9 5. MINIMIZATION OF ERROR DUE TO FINITE DIFFERENCE APPROXIMATION 13 6. EQUATIONS FOR VARIOUS MATHEMATICAL-PHYSICAL MODELS 15 7. ANALOG MODEL 21 8. COMPARISON OF ANALOG AND ANALYTICAL SOLUTIONS 27 9. ADDITIONAL RESULTS OF COMPUTATION 31 10o SUMMARY AND CONCLUSIONS 45 REFERENCES 47 APPENDIX 49

LIST OF FIGURES Fig. Page 1 Physical model of diffusion from an infinite line source showing coordinate axes and cloud outlines. 2 2 Finite difference model showing the vertical coordinate axis relations. 10 3 Typical station in the passive network analog showing the analogy between circuit components and the physical properties of the atmosphere pertinent to the diffusion problem. 22 4 Electronic analog, passive network analog, and physical, finite difference model of the atmosphere. Note that each group of components represents a layer of the atmosphere. 23 5 Typical station in the active or electronic analog circuit. The concentration is represented as a voltage, and equation coefficients as potentiometer settings. 25 6 Plot of the concentration as a function of height for various values of the distance downstream for the case with constant wind speed profile. 32 7 Histogram of concentration versus height. Constant wind profile at X = 50. This figure indicates the real form of the computer solutions. 33 8 Plot of concentration versus height. Constant wind profile case with ground absorption. 34 9 Plot of concentration versus height. Constant wind profile with gravitational settling. 35 10 Plot of concentration versus height. Constant wind profile with gravitational settling and ground absorption~ 36 11 Plot of concentration versus height. Seventh root power law wind profile. 37 12 Plot of concentration versus height. Logarithmic wind profile. 38 vii

LIST OF FIGURES (Concluded) Fig. Page 13 Lines of constant concentration in the X,Z space for the constant velocity profile case. 39 14 Lines of constant concentration in the X,Z space for the constant velocity profile with ground absorption. 40 15 Lines of constant concentration in the X,Z space for the constant velocity profile with gravitational settling. 41 16 Lines of constant concentration in the X,Z space for the constant velocity profile with ground absorption and gravitational settling. 42 17 Lines of constant concentration in the X,Z space for the power law wind profile case. 43 18 Lines of constant concentration in the X,Z space for the logarithmic law wind profile case, 44 vii

ABSTRACT The partial differential equation which describes steady-state diffusion from an infinite line source has been replaced with a set of simultaneous ordinary differential equations solved on an electronic analog computer. One space dimension, distance downwind, was represented by computer time; the other, height, was replaced with finite differences. Solutions were obtained for constant, power law, and logarithmic wind profiles, and for diffusion of particulates which can settle out and deposit on the ground. All solutions were obtained with one basic computing circuit. Each problem required only a particular setting of the coefficient potentiometers in the circuit. Tmplementation of this circuit required only 9 integrating amplifiers and 26 coefficient potentiometers, available in any medium sized computer. The solution's accuracy was measured by comparing the computer plots with the analytical solution for constant wind profile. This measured the total error due to the finite difference approximation and to computer errors. The solution's accuracy was found to be 5% or better over most of the field. ix

1. INTRODUCTION Diffusion due to atmospheric turbulence can be described by the parabolic diffusion equation. This is the direct physical approach in contrast to that which invokes statistical concepts. Through the use of such an equation and its boundary conditions, diffusion problems are conveniently specified, but analytical solutions are not so readily obtained, So great is the difficulty, indeed, that it is unusual to find a solution that fits problems of interest. With the aid of various computing devices, such solutions can be found. Given such a device, in this case an electronic analog computer, it is of interest to explore the form of the solutions as a function of the various equation parameters. The physical problem considered here is that of steady-state diffusion from an infinite line source, oriented normal to the wind, The solution will specify the concentration of the diffusing substance as a function of position in the region after the source has been emitting at a steady rate for a long time. It will be obtained as a plot of concentration as a function of distance downwind for discrete height intervals, As examples of this process, consider diffusion of exhaust gases from automobiles on a busy road, or a uniform, slowly advancing grass fire. The physical model, that is, the idealized case, is shown in Fig. 1 and consists of a line source of infinite extent along the y axis. The x axis is oriented downwind which is normal to the source and the z axis extends vertically. The wind vector is assumed to have no components in the y or z directions, Since the source has no height or width, the concentration at the source must be infinite,

a) o0 CH 0 0.-i r.D rcl 0CH rc$ 0 Cd 0) ~~~~~~~~~~~~~~~~~~~~~~~~~4-3'~~~~~~~~~~~~~~~' Cd Cd m +~ Co )i~~a c1j 00'~~~~~~~~~~~'' —~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~.bOO N~~~~~~~~~~~~~~~~~~~~~4 =. m~

2. THE PARABOLIC DIFFUSION EQUATION The general diffusion equation is AX a Kx x + - Ky y -z L z (1)K where Lx = ax + u a2 + V 6x + w 2X dt 6t yx by az and X = concentration of the diffusing substance t = time Ki = coefficient of eddy diffusivity in the i direction u = mean wind speed in the x direction v = mean wind speed in the y direction w = mean wind speed in the z direction. In accordance with the physical model, the derivatives with respect to y are zero, and the wind components v and w are zero, The steady-state provision eliminates the partial derivative with respect to time. If we also neglect the diffusion term in x compared to the translation term, that is, if ua ax [X Kx x then the diffusion equation becomes u =a [ ax] (2)

With the equation in this form, the solution will depend upon the form of u and Kz and upon the boundary conditions. If u and Kz are given by m n u= -(, KZ(Z) = Kl(Z and the boundary and continuity conditions are KzaX/)z + O as z + O, x > 0O X + ~ along x = z = 0 00 f IX(x,z) dz = Q for all x > 0 (3) 0 where Q is the rate of emission of material per unit length of source, then the solution given by Sutton1 is Q _ Q [[ -(XY _) exp 4- () (x,z) = iulr(s) (m - n + 2) Klxl L (m - n + 2) ()Klx where m +1 S - m- n +2 X(x,Z) is the mean concentration of the diffusing substance. This is a very useful solution since it can be used for any power law wind profile provided that m - n + 2 > 0. For the very important case where the wind profile is given by u oC log z, Kz c z, no analytical solution has been found. An approximate solution for the logarithmic wind profile, u.C log z, Kz c z, was given by Karplus and Allder.2 They used the electric analog which is similar in principle to the technique reported here. It consisted of building

a ladder network of resistors and capacitors which satisfied a set of simultaneous differential equations which was an approximation to the partial differential equation (2). The electronic analog approach is different primarily in that a standard electronic analog computer is utilized instead of building a special circuit as was done by Karplus and Allder. This has the advantage of using standard, general purpose equipment, and is much more flexible.

3. VARIOUS MATHEMATICAL-PHYSICAL MODELS The basic mathematical model which expresses the physical system shown in Fig. 1 is Eq. (2). The special cases to be considered here are (a) Constant velocity profile with standard boundary conditions [those specified in (3) ]; u = const, Kz = const. (b) Constant velocity profile with standard boundary conditions and gravitational settling; u = const, Kz = const, w = fall speed. (c) Constant velocity profile with ground absorption; u = const, Kz = const, iz j OO (d) Constant velocity profile with gravitational settling and ground absorption; u = const, Kz = const, w = fall speed, X 0.O az z=O (e) Power law wind profile with standard boundary conditions;! 1/7 63/7 u(z) = 1ul/ K (z) = K- () 6/7 (f) Logarithmic wind profile with standard boundary conditions; u( z) = ul log( ) Kz(z) = K1 ( z The complete model is wfA(z) ax [B(z) + w X (5) where the functions A(z) and B(z) represent the wind speed and eddy duffusivity terms, respectively. Thus u(z) u1 K,(z) = K1B( U(Z) A Ul Kz(z) = K1B(z) ~ A( z) ae

earthward. The boundary conditions for (5) are X + oo along x = z = O f u(x,z)dz = Q for all x > O Kz + 0 as z + O, x > 0 az for cases a, b, e, and f. Kz # O at z = O, x > 0 Zz for cases c and d. To facilitate handling Eq. (5) and to increase the utility of the solutions, it will be expressed in nondimensional terms. Set S = X/Xo and Z = z/z where the reference parameters Xo and zo will be defined later. Then Eq. (5) may be written = A(Z) [(Z) S + (6 where K1 x - X zoul K1 and the symbol X is reserved for incorporation of some more constants which will appear later.

4, PREPARATION FOR ANALOG COMPUTATION Since this problem is being prepared for solution on an analog computer, we must observe the restriction that the analog computer can integrate only with respect to one independent variable, namely, computer time. As stated [Eq. (6)], the model has two independent variables x and Z so that one of them, in this case Z, must be eliminated. This is done by replacing derivatives with respect to Z by finite differences as follows: as | Sn+l - Sn-1 n (AZ) + (AZ) A..n_ n-I 2 2 a as 1 rSn+l - Sn Sn -Sn-1 az azj n ((AZ)n (AZ) n (AZ) n nImplementation of finite differences requires a new model to clarify the meaning of the position n and the incremeht AZ. Such a model is presented in Fig. 2 which shows the atmosphere dividea into 10 layers in the vertical, The thickness of the layers increases exponentially with height so that the thickness of the top layer is infinite, The reason for this exponential distribution is discussed more fully later. Each layer has a thickness AZ, is infinite in the y direction, and semi-infinite in the x direction, All properties of each layer are lumped in the center of the layer so that the model may be thought of as a stack of planes spaced exponentially. Each plane incorporates the properties of the layer of atmosphere which it replaces. For convenience, the planes are labeled 1 through 10 and represent the computing stations, Thus the

10 2.996 n Z 9.5 - 2.303 9 1.897 8.5 - 1.609 8 1.386 7.5 - 1.204 7 - 1.050 6.5 - 0.9163.6 0.7984 5.5 - 0.6932 5 0.5977 4.5- 0.51 10 4 0.4305 3.5 - 0.3570 3 0.2876 2.5 - 0.223 1.5- 8:1.5- 0.051630 Fig. 2. Finite difference model showing the vertical coordinate axis relations. 10

notation Sn(x) represents the concentration in the nth layer as a function of the distance downwind. The first partial derivative with respect to Z required for the settling term must be centered at the station n and must span from n - 1 to n + 1 so that the increment of Z is Zn+1 - Zn-lo If the layer spacing were uniform, this would be just 2(AZ). The notation (AZ) I indicates the span from Zn to e 2 2 Zn+l, and (AZ) n indicates Z - Z n. Now the meaning of Xo and zo can be defined in terms of the finite difference model. In (3), the continuity condition 00 Q u x(x,z)dz (x > 0) relates the source strength Q to the concentration and wind speed. In the finite difference model, this expression holds at x 2 O, so at x = O, a Q = f u X dz auX since the concentration is taken to be constant throughout the lowest layer, which extends to a height a = 0o1053zo, and is zero above this layer. Labeling the concentration at x = 0 as Xo, then Q xo Is (7) ~.1053zou If u is a function of height, the value used here is the average over the height interval 0 ~ z' 0o1055z0. In the model the source is a region of uniform emission through some height 11

increment and the spacing of the layers is set accordingly as a function of this height increment. Using finite difference expressions for the derivatives with respect to Z, Eq. (6) will be replaced with a set of 10 simultaneous differential equations, one for each layer. The form of these equations is dS An Sn+1 - Sn Sn -n Sn 1 n+- sn- dX( AZ) n n (AZ) n (AZ+ + (AZ) 1 Te2 nw 2 2 The new independent variable X incorporates any arbitrary constants needed. 12

5. MINIMIZATION OF ERROR DUE TO FINITE DIFFERENCE APPROXIMATION Since Eq. (8) is only an approximation to Eq. (6), it is important to analyze the inherent error and to find ways of minimizing it. Considering the case where the layer spacing in the vertical is uniform, the finite difference approximation for the second derivative of S with respect to Z is given by's Sn1 S S_-SnSn+ - 2Sn + Sn-i a2S 1 [ Z~]Sn+ll - SSn n Sn-1]2) z2 A-AZ Z - AZ = (AZ)2 nLA To evaluate the error involved, expand Sn+j and Sn-1 about the point n using a Taylor expansion,3 Sn+ = Sn+AzS + (AZ 2S +(AZ)3 i +* (9) 1' l z n 2' ZaZ 3Z Znn ~~n n Z 5~ S- = Sn AlZ S + (AZ)2 62s ( AZ)3 )S +... (10) n1l n 2' z n 3' z3- n Add the expressions (9) and (10) and solve for 62S/6Z2, )2S Sn+l1- 2Sn + Sn-1 a22'"~~v'~- )z2j n(AZ)2 where (Az)2 54s 12 )Z4 The error that we are seeking is e and it is a function of the station spacing and of the gradient at a given station. The error can be reduced by increasing 13

the number of stations which decreases the height increment. However, since the amount of equipment involved is proportional to the number of stations, we must keep the number to a minimum. Since the gradient is greatest near the source and falls off with distance from the source, the contribution of the gradient term to the error will decrease with height. Thus we space the stations closer together near the ground with small height increments and let the height increments increase exponentially with height. This reduces the error for a given number of stations and provides more data in the lower layers where the greatest changes occur.

6. EQUATIONS FOR VARIOUS MATHEMATICAL-PHYSICAL MODELS Equation (8) with the appropriate boundary conditions for each of the seven cases is given below. (a) Constant velocity profile. The simultaneous differential equations are dS l )1 Sn+l - Sn Sn- Sn-1 (12) n nn +2n where X = 100 x SI = 1.00 at X= O So = Si for all X Slo = 0 for all X An arbitrary factor, in this case lO0, is introduced in each case to facilitate handling the equations on the computer. Also, since the independent variable X corresponds to computer time, the arbitrary factor controls the solution speed which must be within the limits set by the frequency response of the two-coordinate plotter, used to record the solution, on the one hand, and by the patience of the operator on the other. The boundary conditions must be restated in terms of the finite differences. The source strength is no longer infinite at a line but finite at an area as given by (7). Thus, the concentration ratio S1 in the first layer is unity at X = 0. The requirement that the flux across the boundary, represented by 15

the ground at n = 2, be zero is satisfied by setting the concentration at 2 So, a virtual station below ground, equal to that at S1 for all values of X. The flux at n = is given by 2 as S Si-so _ az lBAZ 2 The concentration in the tenth layer must be zero because that layer is infinite in extent. Since Slo = 0 for all values of X, there are only 9 active cells and 9 simultaneous differential equations when (12) is expanded. (b) Constant velocity profile with gravitational settling. dS | 1 [Sn+l - n - Sn-1 + Sn+ - Sn-1 (13) dX In n0(Z) Z)n+ n — 100 (AZ) n+1 + (AZ) n _ 2 2' 2 2 X = 100 x S1 = 1.00 at X = 0 So = S for all X Slo = 0 for all X W = 1 for all X and Z. This equation is different from the previous one only in that the settling term has been added. Since the purpose here is only to show how to incorporate settling and the nature of the effect on the solutions, the settling coefficient was arbitrarily set at W = 1. If the eddy diffusivity K1 were 104 cm sec l and zo were 103 cm, the fall speed would be 10 cm sec-l since W = wZo/K1. (c) Constant velocity profile withground absorption. 16

IdS Sn +1 Sn S+l - Sn Sn Sn-1 (14) dX 100(A oZ)z (Az) (AZ) 2 2 where X = 100 x S1 = 1.00 at X = 0 SO = yS1 for all X (y = 0.5) Slo = 0 for all X. Equation (14) is just like (12) with the exception of one of the boundary conditions. If we define an absorption coefficient y such that when y = 0 there is complete absorption, the ground acts as a sink. When y = 1.00, there is complete reflection at the ground. This model cannot distinguish between passage through the boundary and deposition on it. For this case, y = 0.5, which means that the concentration in the virtual layer n = O is just one-half the concentration in the first layer for all X. In the case of deposition, enough material is deposited on the boundary to provide the required concentration in the virtual layer. (d) Constant velocity profile with gravitational settling and ground absorption. dS Sn+1 [ n+n - Sn Sn - Sn- 1 Sn+- Sn- (15) dX n 100(AZ)n (z (z) n100 (AZ)) -( + + (Az) n 2 22 2 where X = 100x S1 = 1.00 at X= O 17

SO = 7S1 for all X (7 = 0.5) S l = 0 for all X. This case provides the logical combination of the two previous ones. The material is assumed to have a finite fall speed and is allowed to deposit on the boundary surface. (e) Power law wind profile. 1 6 6 d (Z (Sn+- Sn) ) (Sn- Sn-1)| (16) d n n 86.6 /AZ A n+ n-1 2 2 where 5 Z 7 X= 86.36 ~ x Z1 So = 1.00 at X = O SO = S, for all X S'o = 0 for all X. In comparing the distribution of concentration in various wind profile regimes, differences occur due to a change in the average wind speed and to the wind speed distribution. If we set Zo = Z1 to eliminate the arbitrary constants raised to the 7 power in X, comparison of this case to, say case (a) or (f), at equal values of X will not show both effects but only that due to a different wind distribution. This is so because the number 86.36/100 compensates for the different average values of K and u. Comparison at equal values of x will show both effects. (f) Logarithmic wind profile. dSl n 65. 72 [ (2) ] [(iZ) a (Sn+l- Sn) ( Z (Sn - Sn l)0 nnl 2 2 (17) 18

where x = 63..72 x (Z1 = Zo) S1 = 1.00 at X = O So = S1 for all X Slo = O for all X log Z = natural logarithm of Z. As in the previous case, compensation has been introduced for the average u and K.

7. ANALOG MODEL The physical model now consists of a stack of planes spaced exponentially in the z direction. It is supposed that each plane has a capacity for containing some of the diffusing substance just equal to that of the layer of atmosphere that it replaces. At a given point in x, a traverse in the z direction shows the properties of the atmosphere lumped at discrete intervals. This model can be simulated with an electric analog, where each of the circuit components is analogous to some property of the atmosphere. A segment of such a circuit is shown in Fig. 3. The analogies shown were developed from the similarity of the node differential equations for a simple case, CndV = (Vn+- 2Vvn + V n-) dt n Rn (Az)ni. K (Az) n u d (Xn+l - 2Xn + Xn-1) To perform the simulation, set the resistor Rn proportional to Az/k, the capacitor Cn proportional to Az u and voltage V proportional to concentration. The problem time t will represent distance downstream. All these relations are stated as proportionalities so that scale factors can be used to provide reasonable values of resistance and capacitance. In general, one may relate electrical resistance to atmospheric resistance to diffusion and electrical capacitance to the capacity of a layer of the atmosphere to hold the diffusing substance. The complete network analog is shown in Fig. 4 along with the electronic analog and a representation of the physical model. 21

P.i ra o p D3 4 m 4l'o'A c fkJ.s - > C < ~t <'I~ -S I - -. X - 3 *H CIiI4 CC) T IJ Q) N1) p..1 22

so,~ Iz10 +S9 -SS.S -S8 + 57 S 7 -S4 S4 6 +53 5 s 3S L x x ELECTRONIC SO PHYSICAL MODEL ANALOG CIRCUIT PASSIVE NETWORK ANALOG Fig. 4. Electronic analog, passive network analog, and physical, finite difference model of the atmosphere. Note that each group of components represents a layer of the atmosphere.

The passive circuit analog was not used in this study, but it is helpful to construct it on paper because of the insight it provides. This is particularly useful in determining the method of applying the boundary conditions. To insure that S1(O) = 1.00, the S1 capacitor is disconnected from the circuit and charged to an arbitrary voltage which is taken to be unity, At the position (computer time) X = O, the capacitor is switched into the circuit. The requirement that SO = S1 for all X is simply provided by making the resistor Ro which connects the virtual station So to S1 infinite. Then the flux (current) between them is zero as required. If we want So = 0 (nonreflecting boundary), connect resistor Ro to ground as shown in Fig. 4. However, the requirement that So = 0.5 S1 is difficult to implement in the passive analog. Again, to insure that Slo O, make the position 10 at ground potential. The electronic analog circuit could be designed directly from the passive network or from the set of simultaneous differential equations. The latter method is much the easier course after Eqs. (12)-(17) are reduced to the form dX- = anSn+1 - bnSn + cnSn-I - (18) The principal element in the electronic analog is the integrating amplifier which can be used to represent one station as shown in Fig. 5. In the passive network analog, there was a direct correspondence between the circuit elements and the physical problem, while in the electronic analog, there exists a oneto-one correspondence between the circuit elements and the defining equation (18). The coefficients a, b, and c are set on the correspondingly designated coefficient

+Sn-i +Sn+: f I /L f AMPLIFIER Q —' COEFFICIENT POTENTIOMETER Fig. 5. Typical station in the active or electronic analog circuit. The concentration is represented as a voltage, and equation coefficients as potentiometer settings.

potentiometers. The integrating amplifier sums the three inputs, integrates the sum, and inverts the sign. Ordinarily, the gain of the integrating amplifier is unity, set by the 1 megohm input resistor and 1 microfarad feedback capacitor. If one of the coefficients is greater than unity, the corresponding input resistor used is 0.1 megohm which provides a gain of 10. This is necessary because the maximum setting of the coefficient potentiometers is 1. The electronic circuit is much more flexible than the passive network permitting almost any conceivable form of boundary conditions to be set up readily. 26

8. COMPARISON OF ANALOG AND ANALYTICAL SOLUTIONS The form of the solution is a voltage at the output of each amplifier which is proportional to the concentration at the appropriate level in Z and which is a continuous function of X. When recorded on a two-coordinate plotter, the result is a family of curves representing the concentration at discrete levels of Z. The evaluation of errors in this type of computation is especially important since there are two very different possible sources. One is the approximation involved in using finite differences and the other is due to the computer. Previously, we indicated how the error due to the use of finite differences could be minimized but said nothing about the magnitude of the errors. The procedure here will be to compare the computer solution of case (a) constant wind profile, with the analytical solution which may be obtained from Eq. (4) by setting m = n = 0. 2 2 lr KEquation4K (19) is the solution for diffusion from a line source of infinite Equation (19) is the solution for diffusion from a line source of infinite concentration while the computer solution is for a vertical area source of limited height and of finite concentration. There seems to be some problem in comparing these solutions. The actual source strengths can be related by substituting Eq. (7) into Eq. (19) and, for convenience, setting K1 = UZ = Z = 1. 27

s = X = 0.1053 exp - (20) Xo 2r( 2 xL 4x In terms of the nondimensional parameters X and Z, (20) is 1 25Z120 S = 0.59412X exp - - (20a) We can determine whether the finite difference solution is a good approximation to (20a) by comparing the form of the solutions with respect to X and Z and by comparing values of S at selected points in the X,Z space. A simple relation between S and X exists at Z = 0 such that = 0.5941x 2 (21) Z=O Then we can observe the relation of S to the exponential term by finding the slope of logarithm of S with respect to Z. Slope = log Sl - log S2 = _ 25 (22) z2 Z2 x Table I shows how well (21) and (22) are satisfied. The exponential term seems to fit well out to X > 100 where the slope was very small. At X 2 200, all the concentrations are small, S =<.04 so that a constant error which is small for S becomes a large percentage error when S is small. This is shown in Table II where computer values of concentration are compared with analytical values. The most significant feature of Table II is that the magnitude of the error is fairly constant. The percentage error increases as the concentration 28

TABLE I DETERMINATION OF THE FIT OF THE FINITE DIFFERENCE SOLUTION TO THE ANALYTICAL SOLUTION x x2S z0 -X [slope] 1 0.549 24.8 5 0.586 26.3 10 0.573 25.6 25 0.572 25.3 50 0.566 25.5 100 0.573 25.4 200 0.572 32.1 300 0.532 43.2 Analytical Values 0.5941 25.0 TABLE II CONPARISON,QF ANALYTICAL AND FINITE DIFFERENCE SOLUTIONS AT SELECTED POINTS n x Concentration Magnitude Percentage Analytical of Error Error 1 25 0.119 0. 0050 4.2 2.116 o.oo0040 3.5 3 0.109 0.04o 3.7 4 0.0987 0.0045 4.6 5 0.0831 0.0038 4.6 6 O. 0628 0.0027 4.3 7 0.0394 0.0019 4.8 8 0.0174 0.0011 6.3 9 0.00325 0.00025 7.7 1 100 0.0594 0.0021 3.5 2 0.0590 0.0020 3.4 3 0.0582 0.0021 3.6 4 0.0567 0.0020 3.5 5 0.0543 0.0020 3.7 6 0.0502 0.0015 30 7 0. 0447 0 0013 2.9 8 0.0368 0o 0016 4.3 9 0.0242 0.0022 9.1 29

decreases indicating that there is a limit to useful computation. When the concentration ratio S falls to less than 0.001, computer error in the form of noise becomes significant. Virtually all the error shown in Table II is due to the limitations of the finite difference model. The error in this computation is less than 5% provided that S > 0.001 and except in layer 9. Greater errors should be expected in this layer since it is the largest active one and is next to the infinite layer for which zero concentration is assumed. 3o

9. ADDITIONAL RESULTS OF COMPUTATION To present the information more compactly, the solutions have been replotted on log-linear paper as shown in Fig. 6. This presents concentration as a continuous function of height although what was actually obtained was a histogram as shown in Fig. 7. Taking Fig. 6 as a standard for comparison, it is seen from Fig. 8 that the effect of ground absorption is to produce lower concentrations in the profiles. Figure 8 also shows that the maximum concentration found in the profile for X > O occurs above the ground, thus indicating that the plume axis lifts off the ground. Figure 9 shows that gravitational settling produces profiles of somewhat lower concentrations than in Fig. 6 and with more slope. In this case, of course, the plume does not rise. The combination of settling and absorption or deposition is shown in Fig. 10. Solutions due to power law wind profile and logarithmic profile are shown in Figs. 11 and 12. The level of concentration in the profiles is very nearly the same but the slopes are slightly different. In both cases, the profiles for small values of X approximate to straight lines. To give an intuitive appreciation for the nature of the solutions obtained, they are presented in Figs. 13-18 as lines of constant concentration in the X,Z space. 31

1.00 X=X=5 ) 0.10 0 z X=25 _ Iz IX= 300 i = I 0 z 0 0.01 - 2 3 4 5 6 7 n 8 9 0 0.3 0,6 0.9 1.2 1.5 1.8 2.1 HEIGHT Z Fig. 6. Plot of the concentration as a function of height for various values of the distance downstream for the case with constant wind speed profile. 32

1.00 I I X 50 0.10 z) 0 I:0.01 12 31 4 5 6 17 8 9 00010 0.3 0.6 0.9 1.2 1.5 1.8 2.1 HEIGHT Z Fig. 7. Histogram of concentration versus height. Constant wind profile at X = 50. This figure indicates the real form of the computer solutions. 33

1.00 X=I \ = 0.50 0.10 X=5 z 1=10 z / PLUME AXIS 0 X=25 0.01 X=70 I 2 $ 4 5 6 7 n 8 9 0 0.3 0.6 0.9 1.2 1.5 1.9 2.1 HEIGHT Z Fig. 8. Plot of concentration versus height. Constant wind profile case with ground absorption. 34

1.00 X=I, W=1.00 X 5 0.10 ol t-4 \ X 5= IN. 30X=70 0.01 1 2 3 4 5 6 n 7 8 9 0 0.3 0.6 0.9 1.2 1.5 1.8 2 HEIGHT Z Fig. 9. Plot of concentration versus height. Constant wind profile with gravitational settling.

1.00 X=I 7= 0.50 W= 1. 00 0.10 = X=5 z Z - \P= A o ~0~? X 12 3 4 5 6,7 8 9. O l I I 0 0.3 0. 6 0.9 1.2 1.5 1.8 2.1 HEIGHT Z Fig. 10. Plot of concentration versus height. Constant wind profile with gravitational settling and ground absorption. \~~~3

1.00 0. 10 O) Z 0 -3 z X= 70 0.01 -I X= I X=5 X=I0 X=25 X= 50 234 4 5 6n7 8 9 L,,, I 0 0.3 0.6 0.9 1.2 1.5 1.8 2.1 HEIGHT Z Fig. 11. Plot of concentration versus height. Seventh root power law wind profile. 37

0.10 - U) z 0 ) — 0 0.01 - =70 X=I X=5 X=10 X =25 X =50 1 2 3 4 5 6 n 7 8 9 0 0.3 0.6 0.9 1.2 1.5 1.8 2.1 HEIGHT Z Fig. 12. Plot of concentration versus height. Logarithmic wind profile. 38

-1~FrI I I I 0 co 0 ODr) 0 (0 0 qt ~0.0 Z~ x _o D Z zl 0 C\Z 0 n/ 0 O Ci~ ttl~~~~~~~~~~~~~~~~~~C\ - 0 N N N~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~0 ci) d co O d~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~C 0z ci)C ~ 0 o 0 0' 0d 0 0 0 lkl~~~~~l4( ~1 0 d~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~r~C 0 0C) 0~~~4.I.~~~~~~~~ HH13 OD 9 - C~~ 00 ~~~n *H d~d CON C\j -n LO 0 a Z IHSI3H 59~~~~0

' -- I I I I'* I ~I (O cnU Cr)~~~~~~~~~~~~f x\rr 4oI0 N LIJ'~~~~~~~~~~~~~~P w 0 i x x'~~~~~~~~~~~~~~c - 0~~~~~ CD ~ ~~~~~~~~ 0 0 rl) 0~~~~~ I~~~~~~~~ i U) rOx \ ~ ~ ~ (0 N0 3:.: jw 0 0 z o If)~~~~~~~~~~~~~~~~~~~~~5 ~~~~~U~~U ~~~~~~~~o d/ ~~~~~~~~~~~~~~~~n o~~~~~~~~ cd.H O C~~~~~~~~~~~~(0H 0 14.Q~~~~~~~~~~~~~~~~~~~!1 Cja C_' 0 \ LO 0 \ O d~~~~~~~~daU \ O 4-3 ~ r -I O \ o\\ 40~~~~~~~~~~~~~~~~~~~~~ \~~(

I j I ) X (D w W: _ C\j I d-~~~~~~~~~~I aa o~~~~~~,, o 0 0) -r, Z6 ~+ co C.CC)~~~~~~~~~~~~~~~~~ C O'U~~ 0 0 (,03 3r.ao (:5 ~r p) 0 ~~~~~~~~~~~~~~~~~~~~~~ 0~Z= r..)L=) b.Ohl O d ~~~~oZ. o/ 0o o ~~~~~~~~~~ (0 0 ~A0 -- (D ( 4 o N~~~~~~~~\ 0 0 0 Nw\l * oCD ~ ~ D o o & 0 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~0 0 1 0 _0 Z\P ~~~~~~~~~~0 o N I-~~~~~~~~~~~~~~~~o a U) & oJ0~~o oJ ~ ~r) ro -- o~ ~- ~'~ roo - - U)0 ~~~~~t~o N 0 4-p d~~~~~~~~~~~U d H rl\HI 0 H+) (D~~~~0C ~~~~n re, a, F~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~D t~~~~~~~~~~~c) O~~~~~~b. - 41~~~~~~~~~~~~~ 0 0i( ~~~~~~~CJ N liii 1111.1~~~~~~~~~~~~( U) r~~~~~~) 0) U) ro - 0~~~~~~~~~~~C\

'I' 1~~~~~~~ II I1 I I 1 -o XI WI N OO o a.C 0 o31~~~~~C,, ~r. c~o I~~~~~~~~~P r( ~'r-I a 1~~~~~~~~~~~~~~~~~~~~~~~~~~a 00 0d ~~~~O o xH c i NrI %. Cd o> ca) H 4- 3'r-i 0 I~~~~~~~~~~~~~~()* 0 0.Hr-I I~~~~~~~~~~ cO. l'H rr~~~~o - - - 0 0 0 0 Z IHOI3H~~~~~~~~~~~~~~~~~~~~~~~~~~~~ZQ rj-t,,O0 C\ ed +).a ro L0 o Cj 0 0~~~~~~~~~~~Q ~r ci~~~~~~~~. 0 ~ \~~~~~~~~~~ 0.3.3 ~0,) d \Z. C\j Ln 0-~cr) c - d d d Z.LH~I3H

I I [ I Co a -0!0 a) 0D I0 It ro Co cu d F:~~ 0 (5 ~O Z~CU ~~~~ C) 4'\ I I Id CU O ICo C\j~~~~~~~l ~1a o ~o a~\ 00 ed0 ** o.~~~~~~~~~~~~~~~~~~~~~O dd Z o rq C)) O ~, 0 0D c~0 0 ~~~rd OD 1 CH O o~~~~~~~~~~~~~~' b.0 (D~~~~( pi( LO 6O6 6 6 I I, I I H H IZ3 rr -- 03 I u, O "=~~~~~~~L d c~ d d 6 Z iH91'qH

I i i~~~~~~~ I' I' ~~~~~~~~~~~~~i I'i LI (t 0 CH Pii p4 coo2 5 0 z~~~~~~~0 0 l 0 ci ci d 0 0 Zo c3. 3r, c0 ~~o~ o ) CMH4[ o o do old z ~~~~~~~~~~~~(.3 CC) Q - C C\j C. 4-:) 0) o C)) -r o: ~ 0 c or-I OD~~~ 0 C oO (D INr)~~Ci~r O -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~t - i i a Z IHOb3H 44 03 r,.- ~~~~~l rd')

10. SUMMARY AND CONCLUSIONS Solutions to the diffusion equation for the case of the infinite line source have been found for three types of wind speed profile and for the case of diffusion of particulates with an appreciable fall speed. These are sample problems that were chosen to demonstrate the method. In the problems selected, the functions of wind speed and eddy diffusivity with height were specified in analytical form and then translated to discrete steps in height. Since they are used in this form, one could as easily start with functions expressed graphically or in any arbitrary manner. Indeed, using the passive network analogy as a guide, it is possible to design the computing circuit directly from the physical problem. Even the layer spacing used need not follow some mathematical rule; it could be chosen to fit the problem. The wind speed and eddy diffusivity were assumed constant with distance downwind. Changing them with respect to x would involve changing resistors in time which is not difficult to implement with servo multipliers which are usually found in analog computer installations. In the same way, the boundary conditions could be made a function of x. Thus one could simulate diffusion over undulating ground. The accuracy obtained here was better than 5% over most of the field. This could be improved to 1% by increasing the number of cells inasmuch as analog computers are available with component quality of 0.01% of full scale. Extension beyond 1% solution accuracy is not practical with this technique, nor is there any theoretical or practical need for greater accuracy than this.

The advantages of the analog technique in problems of this type are the speed and convenience with which a solution may be obtained and the ease with which atmospheric parameters may be varied. In addition to its value as a research tool, the analog method should be a valuable teaching aid in the analysis of atmospheric diffusion. Taking advantage of the direct relation between mathematical operations and computer components, an instructor could demonstrate before a class the effects of changing parameters on the solution obtained. This method may be extended to problems of finite line sources, area sources both horizontal and vertical, and to elevated point sources. 46

REFERENCES 1. Sutton, 0. G., 1953: Micrometeorology, New York, McGraw-Hill Book Co., Inc., p, 281. 2. Karplus, W. J., and J. R. Allder, 1956: Atmospheric turbulent diffusion from infinite line sources: an electric analog solution. J. Meteor., 13, 583-586. 3. Howe, R. M., 1960: Solution of partial differential equations using the electronic differential analyzer. The University of Michigan, Department of Aeronautical and Astronautical Engineering.

APPENDIX The following tables list the equation coefficients in the form of Eq. (18). These are the coefficients that were applied to the computer in the form of coefficient potentiometer settings. Equation (18) is repeated for convenience: dX n = anSn+1 - bnSn + cnSn-1 (18) TABLE A- I CONSTANT WIND PROFILE n a b c 1 0.8571 0.8571 0 2 0.6780 1.444 0.7661 3 0.5226 1.089 0.5665 4 0.3887 0.8422 0.4545 5 0.2735 0.6011 0.3276 6 0.1782 0.4015 0.2233 7 0.1034 0.2416 0.1382 8 O. o4833 0.1218 0.07348 9 0 0o04131 0. 02820 TABLE A-II CONSTANT WIND PROFILE WITH GROUND ABSORPTION Same as Table A-I except bl = 1.308 for y = 0.5 49

TABLE A-III CONSTANT WIND PROFILE WITH GRAVITATIONAL SETTLING n a b c 1 0.9046 0.8571 0 2 0.7205 1.444 0.7236 3 0.5595 1.o89 0.5292 4 0o.4202 0o. 8422 0.4220 5 0.3009 o.6o011 0.3002 6 0.2006 0.4015 0.2007 7 0.1208 0.2416 0.1208 8 o. 6068 0.1218 o.o6113 9 0 0. 04131 0.02099 TABLE A- IV CONSTANT WIND PROFILE WITH GRAVITATIONAL SETTLING AND GROUND ABSORPTION Same as Table A-III except bl = 1.308 for y = 0.5 TABLE A-V POWER LAW WIND PROFILE n a b c 1 0.1956 0.1956 0 2 0.2499 0.3982 o.1483 3 0.2663 O.4690 0.2027 4 0.2534 0.4721 0.2187 5 0.2213 0.4255 0.2042 6 O.1759 O.3494 0.1735 7 0.1240 0.2551 0.1311 8 0.07137 0 1549 0.08355 9 0 0.o6504 0.03984 5o

TABLE A- VI LOGARITHMIC WIND PROFILE n a b c 1 0.1419 0.1419 o 2 0.2220 0.3403 0.1183 3 0.2550 o.4369 o.1819 4 0.2523 0.4589 002066 5 0.2246 o.4230 o.1984 6 0.1795 0.3497 0.1702 7 0.1264 0.2548 0.1284 8. 06782 0.1450 0.07719 9 0 0.06286 0.03774 51

UNIVERSITY OF MICHIGAN 3901111111115 02526 1 226111111 3 9015 02526 1226