- Email: [email protected]
The adaptive-filtering transport model for prediction and control of pollutant concentration in an urban airshed
Atmosphrric
Ensironmmt Vol 9. pp. 793-808. Pergamon Press 1975. Printed in Great Britain.
THE ADAPTIVE-FILTERING TRANSPORT MODEL FOR PREDICTION AND CONTROL OF POLLUTANT CONCENTRATION IN AN URBAN AIRSHED S. G. BANKOFF
and E. L. HANZEVACK*
Chemical Engineering Department, Northwestern
University, Evanston, Illinois 60201, U.S.A.
(First received 22 May 1974 and in jinal form 28 January 1975)
Abstract-An adaptive-filtering transport model for estimation and prediction of pollutant concentrations in an urban or regional airshed is developed. The program uses the method of fractional steps for direct numerical solution of the turbulent transport partial differential equation, coupled with a regressive least-squares filter which estimates both the concentrations and a spatially-distributed parameter for the adaptive feature at each time interval. The model is applied to historical data on SO,-concentrations in the Chicago area, and gives very good 4-h predictions.
1. INTRODUCTION
In this work we present the AFT (adaptive-filtering transport) model for estimation, prediction and control of sulfur dioxide concentrations in an urban airshed. Following the work of Roth et al. (1971) the pollutant transport is assumed to be governed by a parabolic partial differential equation in four independent variables. The solution method, called the method of fractional steps, relates this problem to the simpler problem of solving a sequence of three two-dimensional differential equations, which are solved successively at each time step. The new feature, introduced to air pollution modeling for the first time (Bankoff and Hanzevack, 1973) is the addition of an adaptive, or learning, algorithm to the overall model for on-line corrections to the concentration predictions and selected model parameters, based on incoming measurement data. This on-line self-correction scheme is a modification of Kalman filtering techniques which were initially developed for use in satellite tracking and other aerospace applications. This gives emphasis to the short term or incident control of urban air pollution, while most previous studies have emphasized long term aspects of air pollution problems. Recently, Desalu et al. (1974) have also employed Kalman filtering for dynamic estimation of air pollution. It differs in several important details from the one given here, including a direct Crank-Nicolson integration scheme, rather than the method of fractional steps found here to be stable and rapid; absence of a divergence-prevention scheme for the Kalman filter; a smoothing maximum likelihood algorithm for parameter estimation, instead of the augmented state vector used herein; and a more complex approach towards partitioning the filter and state equations. The predictions have been compared only with the theoretical plume solution from a single stack. One of the major difficulties of any air pollution model is that the primary meteorological parameters are incompletely measured, making the calculation of secondary parameters extremely difficult. The inclusion of a self-correction algorithm provides a means of overcoming this difficulty by continual re-estimation of one or more model parameters, For example, one of the explored parameters with the greatest sensitivity is the low level source term. This may be looked upon as a bias to the assumed zero-mean Gaussian additive noise in the transport equation, and is used to characterize parameter uncertainties and neglected nonlinearities. The identification of the parameter significantly improves model performance. * Present address: Exxon Research and Engineering Co.. Florham Park, New Jersey, U.S.A. 793
S. G. BANKOFFand E. L. HANZEVACK
794
The set of differential equations derived from the time-splitting method of fractional steps is considered to be a difference operator, advancing the state in time over the entire three-dimensional grid. These differential equations are normally deterministic, but the addition of the previously mentioned noise term makes them stochastic. Coupled to this is a vector measurement equation for SO2 concentrations at actual pollutant monitoring stations which is also stochastic because of measurement errors. Techniques to handle this problem as formulated above are well established by Kalman and others Aoki (1967), Jazwinski (1970), Sage and Melsa (1971). Kalman filtering represents a recursive on-line estimation procedure which can be viewed as an extension of classical least square regression techniques. It is basically a two-step process. Given the current estimate of pollution concentration profile and model parameters, a prediction is made of the same quantities at the next time step. This prediction is then compared with concentration measurements for that time as they become available, which results in a corrected posteriori estimate for all variables. An estimate of the error covariance is simultaneously computed by the matrix covariance equations, which constitutes the principal part of the computational burden. In the current work, Kalman filtering was applied to the estimation of SO2 concentrations at eight monitoring stations in the Chicago area during a short term pollution incident. Concentrations at all other points were obtained from the numerical integration of the model partial differential equation. Input data obtained from a study conducted at Argonne National Laboratory included point and area emission source information, hourly meteorological data, and hourly monitoring station pollution readings. 2. THEORY
2.1 Prediction
model
The starting point for development of a prediction model is the time-averaged tion for conservation of mass in turbulent flow (Bird et al., 1967):
ac t+jX
3 f3(v,c)+ a($?)
(j dx
F)
=il
(Djg)
+
R(cJ+
S(X1~X23X39t)
equa-
(l)
where c and 5j are the pollutant concentration and jth velocity component at the point (x1,xz,x3), averaged over a time interval which is large compared to the dominant time scale of the turbulent fluctuations, but small compared to the time scale for variations in the mean concentrations. Similarly, C’ and U’j are the instantaneous deviations of the concentration and velocity component at any point from their mean values. R(c) and S(x1,x2,x3,t) are distributed-parameter reaction rate and source terms, respectively. The initial condition is ChJ2,X3r&J)
and boundary (a) -K3g (b) (c)
= fh-bx,)
conditions are = Q(xI,x2,t)
= 0 C = g(q,x2,x3,t)
at x3 = h(x,, x2) (ground level) at x3 = H(xi, x2, t) (inversion base) at x1 = xW or xE (west or east) = xS or xN (south or north) (wtkhever boundaries are points of bulk inflow)
wheref(x,,x2,x3) and g(xI,x2,x3,t) are assumed known functions over the urban airshed domain. It is assumed that the flow is incompressible; molecular diffusion is negligible compared to turbulent diffusion, and only the vertical diffusion component need be considered; and the surface flux is negligible, except during precipitation.
795
Prediction and control of pollutant concentration
The following transformation of coordinates is applied to simplify the application of the boundary conditions (Roth et al., 1971): x=
x1 x.E - XIV
# = u1
))=-%XN- xs z
=
x3
-
v = vz
NXl,XZ)
w = v3
AH
where AH = H(x1,x2, t) - h(x,,x2). With these assumptions, the transformed equation becomes
with initial condition C(x, Y,z, t) = f(x, Y, 2)
at t = to
(3)
atz=O
(4)
and boundary conditions
-KC&J =az = 0
atz=
(5)
at inflow points on the sides of the unit square
C = Qk Y, z, t) Olx,yl
1
(6)
1.
Classical purely explicit methods for solving partial differential equations are avoided because of their traditionally poor stability characteristics. Similarly, purely implicit methods involve inversion of very large matrices for three dimensional grids, and therefore require excessive computer time and storage. However, variations of the method of alternating directions can avoid both of these difficulties. The method of fractional steps, which is a time-splitting method, is used here. Essentially, a four-dimensional PDE is broken into three two-dimensional PDEs which are then solved successively at each time step. This retains stability advantages of classical implicit methods, but greatly reduces the computational burden by requiring inversion of much smaller matrices. Even so, for the large scale system encountered in air pollution modeling, the computer storage requirements are excessive if each segment is solved implicitly because the, entire three-dimensional grid distribution must be stored for three different points in time (current and two previous). Therefore, a compromise is reached in which the vertical direction parabolic equation is solved implicitly, and the horizontal direction hyperbolic equations are then solved in turn explicitly. The product form of the fractional-steps recursion matrix insures that the entire process is stable if each fractional step is stable. The Crank-Nicolson finite-difference equation, which is unconditionally stable, is used for the implicit portion, and the explicit steps can be shown to be stable if the Courant number, aj 3 [(ojAt)/(Axj)] I 1 (Roache, 1972). Therefore, stability can be assured by appropriate choice of the time increment. The following series of finite-difference equations, which are second-order correct in both time and space variables, is used: _ C$k - Cyj, + At =
&tc~j,L'
2~~j,~,(C~j.t-l
+
Ryjk+ Syj,,
1 +
+
C,j,k+ 1 -
C,j,k- 1 -
2G$
Ctj,k- 1 -
-
2cYjk +
CY,j,k- 1)
Czj,k+l
+
G,j,k+
1) (7)
S. G. BANKOFF and E. L. HANZEVACK
796
with boundary conditions
-&(G’,j,l - Ci’,j,o)= QYj C,j,K
=
1,
G,j,K+
where K is the mesh point just under the inversion base;
where
_tw)~XF?-~(C:_ ll2,j.k
txE
f
l,j,k
c;k) -
+
%(c$k- CT-
l,j,k)
and uAt ’ = (xE - xW) AX
with boundary conditions At txE
q,;,k
= j(c?,j-
_
xw)~xF:!2,j.k = al12gfPd,k
1,k + c!,j+ l,k) + aI- l/2 CC,*_ l,j,k - 3(c?,j-
l,k
+
c?,j+
l,k)l
for u 2 Oand
_tw)AxFF+ 112,j.k=aI+ lpgr+
txE
c:$,k
= 3cc?,j+
1,k + c:,j-
I,k) + al,Z[+(c:,j-
112,j.k
1,k + c?,j+ l,k) -
cy,j.d
for 1.4< 0, where Z is the mesh point just west of xE; C””
= c$,*[l
+ /jj+1,2
- pj-l,2]
+
(xN
_A;s)&(G?_fG::+l.k)’ l,k
-
(9)
where At (xiv -
xs)AY
GF$_1 k s ---+ Bj-1j2 CW r,,’
1,k
+
c$:)
-2
and p=
vAt
h -
XSWY'
with boundary conditions analogous to those for equation (8). Equation (7) is a Crank-Nicolson formulation for the (z,t) fragment, and is solved implicitly. Equations (8 and 9) are second-order formulations for the horizontal advection ((x,t) and (y,t), respectively) written in conservative form and are solved explicitly. The three sets of equations solved consecutively constitute one time step. Due to the coarse grid spacing, the concentration distribution is often not smooth. The following additional conditions are imposed upon the system to prevent truncation error from increasing: F*_ = i
F*, =
l,j,k
ifa- IO, C,k < -FE ifa_ > 0, CL_1,j.k < F*_
ci+ l,j,k
ifa+ >O,C,k
cijk
ciCijk
Prediction and control of pollutant concentration
197
A similar condition is imposed in the y-direction. The treatment of physical parameters which are only incompletely measured is now detailed and the resulting approximations listed. The physical parameters necessary for model prediction are wind field (aj), inversion height (H), vertical turbulent eddy diffusivity (KS), reaction rate (R), and source emission inventory (s). The wind field must be specified at each of the three-dimensional grid points. As an idealized extreme, it would therefore be desirable to have an hourly measurement at each grid point. Obviously, this is never the case in real situations. At the other extreme, the minimum allowable information is one average bulk wind velocity vector for the entire urban area. Usually, however, an intermediate amount of information is available. Often 2-10 hourly ground level measurements are taken at various points throughout the city. These can either be interpolated spatially to obtain a ground level wind profile, or (preferably) a subjective analysis map of the data by an experienced meteorologist could be employed. In either case, the upper wind profile must still be specified. Roth et al. (1971) suggest using c?u~/~.x,= 0, i.e. the upper wind field is the same as the ground level profile. Alternatively, they suggest a two-step wind field with a discontinuity. Either of these methods may suffice for a rough approximation, but most authors favor a logarithmic or power law distribution. Shir and Shieh (1973) present a typical expression, IuJ = IuO&z/ z$ where the subscript 0 refers to ground level and the exponent p - In (Iuzl/lull)/ln (zJzJ where the subscripts refer to two heights where wind measurements are available (e.g. a television tower). The further restriction 0.15 c p < 0.65, is imposed. This latter formulation is used in validation studies when possible; one of the methods listed by Roth et al., is used when necessary. Furthermore, it is well known that there is an angle shift in the upper wind field, to the right of the surface wind in the northern hemisphere. However, since previous studies proved relatively insensitive to small angle variations (Roth et al., 1971; Shir and Shieh, 1971X this is considered a minor effect and is neglected here. The mixing depth is an important prediction model parameter, but there is little quantitative information about it available. A mixing height value is needed for each vertical column of grid volumes at each time interval, but, typically, measurements are taken at only two or three points, and only once or twice per day. Furthermore, the measured values are for temperature inversion-base heights, which may only roughly correspond to mixing heights. With this sparse data base, most studies (Roberts, 1970; Roth et al., 1971; Shir and Shieh, 1973) use a crude interpolation scheme in both space and time, constrained by assumed maximum and minimum values. The vertical turbulent eddy diffisivity cannot be directly measured, and so it must be calculated from other measurable quantities. One approach describes KZ as a component of a second-order tensor K. It is assumed that off-diagonal terms of this tensor are equal to zero, and if the mean velocity is parallel to the ground, this further implies that all other components can be set equal to zero, by virtue of the neglect of horizontal diffusion compared to horizontal advection. It is generally agreed that KZ is a function of wind velocity, shear stress and the change of temperature with height (lapse rate). However, the functional relationship is not well understood, and a variety of expressions have been proposed. A trapezoidal function at least qualitatively describes the expected behavior of K, with respect to height, and since fairly simple calculations are involved, a relation of that form is sought for specific situations. Further work in this area is also indicated. The reactions of SO2 in the atmosphere include catalytic conversion to H2S04, photochemical conversion to HzS04, and adsorption onto particulates. Of these the first reaction is considered dominant (Turner, 1964). Various researchers (Kontnik, 1973; Shir and Shieh, 1973; Stern, 1968) have listed the half life of atmospheric SO2 due to reaction between 3 and 4.5 h, with 4 h being used here as an intermediate value. An emission source inventory is one of the most basic and important components of a pollution dispersion model. Without a comprehensive picture of emission patterns, accurate short-term pollutant concentration predictions are nearly impossible. Therefore,
798
S. G.
BANKOFFand
E. L.
HANZEVACK
serious attention must be given to this rather mundane aspect of the problem. Most urban inventories are divided into point and area source strengths. Electric power plants and major industrial emitters are considered to be point sources. Hourly or bihourly emission values are computed based on regular process emissions plus a term calculated from the degree day concept for emissions due to heating. These values are additionally corrected to account for varying normal patterns in diurnal and weekly levels of activity. Smaller industrial plants, commercial establishments and residential heating are computed as area sources. Most empirical correlations are of the form where
S = YICI& + Y2C2 ww? $7, S = area source emission (lb h- ’ miles2) S, = annual emission total (tons y- ‘) T = temperature DD = degree days = T - 65°F if T < 65” = 0 if T 2 65” DD, = total annual degree days yl, y2 = hot water heating fraction, space heating fraction, y1 + y2 = 1 cl, c2 = unit conversion factors.
Typical examples can be found in Shir (1972) and Turner (1968). The foregoing represents the urban airshed pollution dispersion model, which now must be reformulated slightly in the next section in order to be amenable to attack by Kalman filtering, techniques. 2.2 Filter algorithm The Kalman filter is related to classical least squares techniques, the main difference lying in the recursive nature of the filter, while least squares methods are for batch processing of data. This difference is crucial, however, for this application. With the filter, new state and parameter estimates can be obtained on-line, as each new measurement set becomes available, by considering only the effect of this most recent observation, instead of re-solving the entire problem from t = 0 each time with a steadily increasing data base and corresponding computational load. In addition, the error covariance matrix, P(k), at the kth time interval, is continually updated, giving a quantitative estimate of model and filter performance. Consider a discrete dynamic system equation and associated measurement equation: x(k + 1) = cp[x(k), k] + T(k) w(k)
(10)
y(k) = H(k) x(k) + o(k),
(11)
where k0 4 k I k, represents the time interval under consideration. In these equations x(k) is a vector of n state variables at the kth time increment, which are, in fact, concentrations at each of the grid points throughout the airshed. Similarly, y(k) is a much smaller vector of concentration measurements at the same time. cp[x(k),k] is a vector of nonlinear functions resulting from the spatial discretization of the transport equation, which may also include time-dependent inputs, such as source terms. The stochastic inputs are the terms w(k), which represent the uncertainties due to neglected non-linearity, poorly-known parameters, model inaccuracies, and physical inputs of a random nature, while v(k) is the vector of measurement errors. H(k) and T(k) are appropriate coefficient matrices, here taken to be constant. We assume that w(k) and v(k) are uncorrelated Gaussian white noises, such that the expected value of w(k) is E(w(k)) = A(k). This represents a bias term in the process equation, which may be re-estimated at each time for the adaptive feature. Similarly, it is assumed that,
~C-+dl = A -QIdk)l= 0 ECx(kcJ xTWl = J=o
(12)
799
Prediction and control of pollutant concentration
E[w(k) wT(k)]= Q(k)G(k - j)
vT(j)]=
E[u(k)
R6(k
- j)
k
k 2 j
2 j,
where ~l~p,,,,P,, and R are assumed known prior statistics, and 6(k - j)is the Kronecker delta. Q(k) is also assumed to be known, but its norm is artificially increased at each time in order to prevent divergence of the filter. It is desired to obtain a minimum-variance recursive estimator for x(k) and to identify a parameter vector p(k). The state is augmented by defining an extended state vector X(k) =
in which the parameter
[ 1, ;;“kl
(13)
vector is assumed to evolve according to p(k + 1) = A(k) p(k) + B(k) w(k).
(14)
Then the extended state equation becomes X(k + 1) = @[X(k)&]
+ G(k) w(k),
where
kl A(k)p(k)
cp[x(k), p(k),
@[X(k), k] =
G(k) =[W BtkI 1a
(15)
1
Thus it is desired to minimize the cost function (Jazwinski, 1970) J = ;
+
i
IlWko) - Ako)ll;,-1 k,-
1
,zk
0
ClMk
+
1)
-
W
+
1)-W
+
1,llih
+
IlWll&~l 7
(16)
subject to the constraints implied by the extended state equation and prior statistics. The modtied extended Kalman filter with an extended state vector is the solution to this problem, yielding the following algorithm: i?(k + l/k) = @[2(k), M(k + 1) = GQ(k + l)GT + “;g;)]
k]
(17)
p(k+‘;%$)],
(18)
P(k + 1) = M(k + 1) - M(k $ 1) HT(k + 1) [H(k + l)M(k + l)HT(k + 1) + RI- ‘. .H(k + l)M(k + 1) (19) i(k
+ 1) = i(k
+ l/k) + M(k + l)HT(k + l)[H(k + l)M(k + l)HT(k + l)]- I.
. [y(k + 1) - H(k + 1) k(k + l/k)],
(20)
where k is the running discrete time index. The updating algorithm (I) is then: 1. 2. 3. 4. 5.
Use initial conditions P(k,) = PO, _$(k,) = px, Predict from (17). Compute prior variance from (18). Compute error variance from (19). Filter to get updated state and parameter estimates from (20).
This constitutes the extended Kalman filter with an augmented state vector, derived from the viewpoint of invariant imbedding (Sage and Melsa, 1971). Note that X(k/k - 1)
S. G.
800
BANKOFF and
E. L.
HANZEVACK
and _f(k) are state vectors consisting of a priori and a posteriori, respectively, concentration and parameter estimates at each of the n,nynz three-dimensional grid points at time k. Equation (20) comprises the updating step; equation (17) is the prediction step. An essential further step, in order to prevent divergence of the Kalman filter, is the addition of an algorithm for articially increasing the norm of the process noise matrix, Q(k), as outlined in Jazwinski. The full problem can lead to manipulation of filter matrices on the order of 4000 x 4000. This is clearly not desirable, so that the problem must be partitioned into related subproblems which are then connected in a second-level solution. Aoki (1967) presents the appropriate state vector partition theory. Let X(k) be an n-dimensional state vector partitioned into 1 subvectors Z(k)‘, . . . ,2(k)’ where Z(k)’ has dimension trj.
This leads to 1 relations of the form Z(ky’ = DjX(k) 1 5 j < 1. The submatrices Dj are taken to be time invariant. Thus, the same partitioning throughout, although this is not absolutely necessary. The prior and posterior state estimates are then given by k(k
+ l/k) = i @(k j= 1
is used
+ l/k)’
_$(k + 1) = i Fjk(k j= 1
+ ly
where
so that the state vector is reconstructed from the partitioned subvectors. Similarly, a matrix Bj is defined which picks out the subvector of the measurement vector Y(k) that is used in updating Z(ky’. The filter algorithm is now applied in turn to each of the subvectors Z(k) and the complete state vector X(k) is then reconstructed using the above relations. This process yields the filter i(k
-t l/k) = @[k(k),k],
(21)
M(k + 1) = GQ(k + l)GT + a;$$)b(k)a@T!2(k)1
(22)
8X(k)
i
j= 1
)
.M(k + 1) +
jil
F,K(k
1 T
E F’K(k + 1yGj H(k + 1) j= 1 ) + l)‘Gj Rj i f’jK(k + l)‘Gj j= 1 >(
T?(k + 1) = i?(k + l/k) + i j=
where
1
F’K(k + 1yGj H(k + 1)
1
F’K(k + ly’Bj[y(k + 1) - H(k + l)Z(k + l/k)],
K(k)j = filter gain matrix = M(kY’H(kY’[H(kY’M(kY’H(kY” + Rj]- ’ @j s Dj~Di’ where + indicates pseudoinverse Q(ky’ ~ DjQ(k)Dj H(ky E BjH(k)Dj+ M(k + ly’ = GQj(k)GT + @jJ’(k)@jT P(Oy z DjP(O)D,’ P(kY’ = [Z - K(ky’H(k)‘]M(ky’[Z - K(k)‘H(k)‘lT + K(k)‘RjK(kY’
and g(O). f’(O)are assumed known.
T,
(23)
(24)
Prediction and control of pollutant concentration
801
This process reduces the size of the matrices which must be manipulated. If the subvectors are not coupled through the state or measurement equations, then the above scheme is optimal. If the system is only weakly coupled, performance may still be tolerable. The specific partitioning scheme used in the present study, which yields a weakly coupled system, will now be described. A three-dimensional grid is superimposed on the urban area in order to solve the model equations by finite difference techniques. Eight x 16 x 4 nodes are used for the x-(west-east), y-(south-north) and z-(vertical) directions, respectively. This yields constant, uniform two-mile spacings in both horizontal directions; however, the vertical node spacing is neither constant nor uniform. The first two vertical nodes from the ground are assigned constant spacing between 100 and 200 ft, and the remaining distance to the mixing height is apportioned between the other nodes, yielding spacings typically between 100 and 1000 ft. Since most source inputs occur in the bottom two grid volumes, it is considered desirable to keep these volumes constant for simplicity. Certain ground level grid volumes contain a pollutant monitoring station. Consider a relatively small volume around one of these monitoring station grid volumes m, defined by a dia. r in the xy-(horizontal) plane and a height h in the z-(vertical) direction. r and h are chosen to be integers in units corresponding to grid dimensions, and r is odd for symmetry. This defines a semi-ellipsoidal subregion centered at m. Over a time arc of At (typically l-h), it is assumed that the state at points outside this subregion does not affect the state at m, in particular, or at other points inside the subregion. This is the basis for a series of partitioning schemes, depending on the value of r and h chosen. If r 2 3, h 2 2 the state at points outside the subregion affects the state estimate at m only to the second-order over a l-h time interval (i.e. for one filter application). If this process is repeated at each of 1 monitoring stations, then 1 subvectors zj, each of dimension r’h, are defined. The 1 subfilters operate on square matrices of dimension r2h, instead of dimension ~l,n,,~, resulting in a considerable computational saving. Note, however, that depending on the size of r and h, all of the grid points may not be included in one of the 1 subvectors by this particular partition. This is equivalent to considering all such excluded points to be in one z’+ ’ subvector in which the subfilter gain matrix K’+ 1 is always the null matrix. Therefore, between subregions, the parameter estimates at excluded points are spatially interpolated after each filter application in order to smooth the surfaces, primarily for computational reasons, since prediction estimates are relatively insensitive to small changes in these parameter values. Each monitoring station grid node and the nodes immediately adjacent (totaling 18 nodes) are considered to be the state subvector for one of the matrix subfilters. One or two parameters are added to this to give the 19- or 20-component, augmented state subvector. The constant (or slowly drifting) parameters represent: (a) A multiplicative source strength adjustment to account for the seasonal or diurnal source strength variations, (b) An additive bias or fictitious source parameter to identify the effect of neglected nonlinearities yielding a background pollutant level. For simplicity and computational savings, these parameters are assumed to be uniform throughout the filter subregion, although this is of course not necessary. We have thus developed a multistage solution intended to make maximum use of available information at each stage individually while maintaining a tolerable computational load. Note that there are several possible formulations to the predictorcorrector approach. The most exact method is to predict by solving the original PDE, and then filter using the entire state vector. However, as previously discussed, this leads to excessive computational requirements. Alternatively, the system may be partitioned into subregions in which a PDE is solved and a filter is applied to each of these subregions. This problem is computationally feasible, but is equivalent to assuming that the subregions are uncoupled over all time.
802
S. G. BANKOFFand E. L. HANZEVACK
Finally, the method actually used in this study employs partial uncoupling. The mass transport PDE is solved numerically for the entire system, and then the system is partitioned into subregions for filtering. This requires the much less stringent assumption that the subregions are approximately uncoupled over l-h time intervals. The same basic PDE is considered from two different viewpoints, over the whole urban region for prediction and locally in space and time for correction. 3.
MODEL
VALIDATION
Most raw data used in validation studies, meteorological measurements, source strength inventories and pollutant concentration monitoring station data, were obtained from studies conducted at Argonne National Laboratory for SO2 levels in the city of Chicago, Illinois during January 1967 (Roberts, 1970, 1971). Other pertinent references concerning source and meteorological correlations are listed. An hourly piecewise-constant ground level wind velocity field was calculated, based on an average of wind speed and direction measurements taken at Midway Airport and monitoring station No. 4. These continuous measurements are integrated over 5-min intervals every 15 min, to give an hourly average. The Norco objective mixing depth calculations (Roberts, 1971) were used for lid height values. These values are calculated from rawinsonde vertical temperature soundings taken at Peoria, Illinois and Green Bay, Wisconsin, at 600 and 1800 CST. The Chicago rural temperature profile is formed by linear interpolation. Then hourly rural profiles are calculated using temporal weighting for interpolation between the sounding times. The hourly mixing layer height is defined by the intersection of a dry adiabat from the urban surface measured at Midway Airport with the constructed rural (Argonne) temperature profile. Residential heating and industrial area source and power plant point source strength information was obtained from comprehensive source inventory questionnaires and calculations based on fuel consumption and ambient temperature. These values were adjusted for varying activity levels during different seasons of the year, day of the week and time of day or night. Pollutant concentration data were obtained from a network of eight telemetered air monitoring stations, which continuously and automatically record 5-min average SO2 concentrations at 15 min intervals throughout the day. Therefore, four 5-min averages are used to compute the hourly average at each location. Integration time increments of 5-6 min are found to be most satisfactory to meet stability, accuracy, and computational burden requirements. The initial concentration distribution necessary to start a simulation run was based on the monitoring station measurements taken at the time interval immediately preceding the run. These values were smoothed visually to yield a ground-level initial concentration value for each grid node. The ground-level initial concentration values were then multiplied by an arbitrary factor of 0.50-0.75 to establish estimates for the upper grid nodes. It was found that these estimates, although somewhat crude and arbitrary, were sufficiently accurate to yield good model performance because the sensitivity of model predictions to changes in initial concentration distribution was small after 3-4 simulation h. Other authors (Shir and Shieh, 1973) have also noted this phenomenon and have started simulation with their model assuming a uniform zero initial concentration distribution. Additional empirical relations for the upper wind field, vertical turbulent eddy diffisivity. and SO, reaction rate are necessary model inputs. A power law formula was used for upper wind field calculations (Shir and Shieh, 1973) IUI= Iv01@/%,‘?
(25)
where z. is the height at which the reference velocity u. is measured, with the value of the exponent p being on the order of 0.5. In addition, no angle shift is assumed in the calculations.
803
Prediction and control of pollutant concentration
The vertical turbulent eddy diffusivity is generally believed to be a function of temperature, height, atmospheric stability, and other variables. However, the exact nature of the functional dependence is largely unknown. Therefore, based on qualitative information (Roth et al., 1971), a trapezoidal function was used, approx. 0.039 km2 h- ’ at ground level and inversion height, and limited to twice that value at intermediate points. To avoid complicated atmospheric SO2 reaction mechanism expressions which would have only a relatively small effect compared to other terms in the model equations, a first-order decay rate was used for the reaction term. The half life for SO2 in the atmosphere was assumed to be approx. 4 h (Turner, 1964). 4. RESULTS
(a) Sensitivity. To determine the approximate sensitivity of model predictions to variations in input parameters, simulated data were computed by assuming that the transport equation models the real system exactly and that there are no errors associated with the finite difference solution of the equation. Twelve-hour simulation runs were then made with IO, 20 and 40 per cent average distortion factors built into each of the following inputs in turn: wind speed, wind direction, source strength, initial concentration profile, inversion height or mixing depth, and vertical eddy diffusivity. These inputs were considered to be the most weakly known and the most likely to have systematic errors. Table 1. Approximate relative sensitivity of 1 h concentration tions to input parameters
Parameter
Wind speed Wind direction Lid height Source strength Vertical eddy diffusivity Concentration at previous time step
predic-
Parameter variation (%) 10 20 40 Change in average concentration (%) 28 4 8 21 3 6 14 6 10 16 4 8 6 2 4 37 9 19
Table 1 indicates the approximate magnitude of the fractional change in calculated pollutant concentration for a specified change in input parameter. The model is most sensitive to the source strength and lid height with moderate parameter changes, although at high distortion levels it is more sensitive to both wind speed and direction. The model is least sensitive to vertical eddy diffusivity changes, It should be noted that an error in the pollutant concentration at the previous time step leads to approximately the same fractional error in the calculated concentration, as would be expected. (b) Filter efictiueness. To insure that the proposed suboptimal filtering scheme could perform efficiently, simulation runs were made attempting to filter out the previously mentioned simulated data distortions by means of a simple additive bias term. These studies can be briefly summarized by stating that in every case the filtered state estimates were consistently more accurate than the model predictions alone. When estimates strayed far from measured values, the resulting increase in noise covariance from the divergence prevention algorithm quickly (within l-2 h) brought estimates back into line. In addition, a slowly drifting bias parameter in many cases could be identified within a relatively short (12 h) run. In particular, initial concentration profile distortions of any size could quickly (within 3-5 h) be identified, involving a substantial bias term, as might be expected. Source strength distortions of 20 per cent could be filtered out of state estimates, although bias term identification took much longer.
S. G. BANKOFFand E. L. HANZEVACK
804
Measured
So2 Cone
(wm)
Fig. 1. Scatter plot of measured vs. predicted concentration
during high pollution
incidents.
For low velocities (less than approx. 8 miles h- ‘), moderate errors in wind speed and/or direction could be overcome in the state estimates, but bias terms tended to oscillate. At high wind speeds, accurate estimation using distorted values was more difficult. Fortunately, however, although current wind measuring devices are notoriously poor for low wind speeds, they are accurate at higher speeds. Tests with the other parameters resulted in similar recognizable patterns being observed. (c) Modellfilter studies. The combined prediction model and filtering scheme was examined by attempting l-, 2- and 4-h predictions using historical data records. In this context, 4-h prediction is taken to mean concentration predictions up to and including time t + 4 given only concentration measurement (and associated filtering) information up to time t. See Fig. 1 for a scatter plot of observed vs predicted concentration levels, and Figs. 2 and 3 for predicted and measured pollution values vs time. These figures show consistent agreement with measured pollutant concentrations at both high and low concentration levels, correctly predicting both trends and approximate magnitudes of peaks.
J 4
8
I2
TIME
Fig. 2. Comparison
16
20
z4
(h)
of measured and predicted SO2 concentrations incident.
during severe pollution
Prediction
0’0
and control
of pollutant
I 4
8
12
TIME
Fig.
3. Comparison
805
concentration
of measured
and
16
20
24
(h)
predicted SO2 incident.
concentrations
during
mild
pollution
Additional studies compared model performance using only current meteorological information for 4-h predictions compared to using actual meteorological (but not pollution) information from the historical data record. One can thus develop upper and lower bounds for the prediction accuracy, using at one extreme the actual meteorological records (corresponding to perfect weather prediction) and at the other extreme the procedure used herein, assuming no changes in wind or inversion height during the 4-h prediction interval. However, except in a few cases (Fig. 4) often it was found that model prediction was relatively insensitive to which set of meteorological information was used, not yielding a well-defined “envelope” (Fig. 5). This is understandable in view of the observed fact that most severe pollution incidents occur during periods of atmospheric stability, during which there are only relatively slow changes in wind velocity and mixing height. Longer-range predictions were then attempted in order to determine the prediction horizon. Eight-h predictions seem feasible, but 12-h predictions are often erratic. In practice, much longer data records would be available to the predictive model than in these simulations, so that even better performance can be anticipated.
TIME
(h)
Fig. 4. Comparison of predictions at time t with measured values of time t + 4; assuming (a) only current meterological information available. (b) meterological measurement data available from historical records during severe pollution incident.
S. G. BANKOFFand E. L. HANZEVACK
806
4
8
TIME
I*
12 (h)
Fig. 5. Comparison of predictions at time c with measured values of time t + 4; assuming (a) only current meterological information available, (b) meterological measurement data available from historical records during mild pollution incident.
All simulations were performed on a C.D.C. 6400 computer. For the forward integration step and update, less than 110 K of central storage and less than 30 s central processor time/stimulation hour were required. These figures compare favorably with those of models developed for other cities or pollutants. The filtering algorithm was responsible for approx. 20 per cent of the total simulation time, which appears to be a small price to pay for the dramatically improved results. (d) Model comparison. Table 2 summarizes the results of validation studies. It also indicates the relative performance, using the same data, of the Argonne National Laboratory Integrated Puff Model, previously the most sophisticated SO2 prediction model developed for Chicago. This table clearly shows the superiority of the APT model, even at the 4-h prediction level, compared to the one-hour ANL predictions. Figure 6 shows a typical comparison for a particular monitoring station. (e) Control strategies. Figure 7 shows the predicted concentration patterns assuming a 50 per cent source reduction at hour 11. This results in partial abatement of the incident, but points out the need for early control action. Other simulation runs were made using other feasible strategies, such as power plant source reduction only, source reduction only in areas of the city with predicted unacceptable concentrations, etc. However, such predictions cannot be validated by historical data but only by actual field tests. Nevertheless, since the model performs well with the historical data, it seems reasonable to assume that similar accuracy could be expected from control strategy predictions, and the implications for possible savings in money and improved air quality for short term pollution incident control are obvious.
Table 2. Comparison of model performance* Prediction time (h) AFT Model
1 2 4
Argonne National Lab. Integrated Puff Model
1
* 192 Observations (Roberts, 1970, 1971). I-Per cent of predictions within a factor of 2 of measured values. 2-Per cent of predictions within 0.05 ppm of measured vaiues.
Prediction error band 1 2 99 93 85 74
98 87 77 65
Prediction and control of pollutant concentration
. 0
4
? ‘:h,
8
16
20
24
TIME
Fig. 6. Comparison of 4 h AFT model predictions with 1 h integrated-puff
model predictions.
5. DISCUSSION
Sparse meteorological information and a rough source inventory used as inputs to an approximate (fractional steps) solution to the continuity equation and followed by a suboptimal filter have been shown to be sufficient to produce reasonably accurate and potentially useful short-term pollutant concentration predictions with only moderate computer storage and computation time requirements. The appropriateness and potential utility of short term prediction is indicated by historical data records (Roberts, 1971) which show that most serious air pollution episodes in Chicago range from about 12 to 72 h in duration, precisely the time scale in which the proposed model and filter scheme is most effective. Here an SO2 concentration of 0.3 ppm for a duration of four or more consecutive hours was arbitrarily used as a minimum standard for the existence of an air pollution episode. The most important variables in air pollution modeling are wind direction and speed mixing depth, ambient temperature, atmospheric stability (vertical turbulent eddy diffusivity), and hour of day. The relationships between these variables and pollution concentrations are generally not linear. Some of the functions are monotonic while others 57
4-
3-
2-
I -
O-
0
4
16
8
20
TIME I2 (h). SO, CONC (PPM) (X10)
Fig. 7. Predicted effect of source reduction on pollution concentrations.
808
S. G. BANKOFFand E. L. HANZEVACK
are cyclic. It is clear that improved meteorological measurement information would be useful in pollution prediction modeling. Specifically, more urban vertical temperature soundings (for mixing depth estimates) and wind measurements (for upper wind field estimates, including heat island and lake breeze effects) would be helpful. Some SO2 concentration measurements at selected elevated points throughout the city would aid in verifying the predicted vertical concentration profile. Also, better correlations for relating turbulent eddy diffusivity to measurable quantities are needed. Different filtering methods could be devised for specific uses. For example, different partitioning schemes could be tested or different parameters updated and identified. Alternatively, a sequential estimation algorithm could be developed in which certain parameters are identified for a specified period (say 12-24 h), then a different parameter set is identified for the next period, and the process could be repeated indefinitely. Additional 24-h average pollutant concentration measurements are taken at 23 points throughout the city of Chicago. An intriguing possibility would be to use these values to develop a two-frequency filter in which partitioning and filtering previously described is used to update the model every hour and a different partition and filter are used to update the model once every 24 h. This “interlocking” of the urban modeling area by means of two different partitioning schemes applied at different frequencies deserves theoretical as well as experimental attention. REFERENCES Aoki M. (1967) Optimization of Stochastic Systems. Academic Press, New York. Bankoff S. G. and Hanzevack E. L., Jr. (1973) Parameter updating for air pollution dispersion model. Symp. on Air Pollution Modeling, AIChE National Meeting, New Orleans, La. Bankoff S. G. and Hanzevack E. L., Jr. (1974) Real-time adaptive air pollution modeling. Presented at GVC/ AIChE Meeting, Munich, Germany. Bird R. B., Stewart W. E. and Lightfoot E. N. (1967) Transport Phenomena. Wiley, New York. De&u A. A., Gould L. A. and Schweppe F. C. (1974) Dynamic estimation of air pollution. JACC, Atlanta, Ga. Hanzevack E. L., Jr. (1974) A self-correcting model for prediction, estimation, and control of sulfur dioxide concentrations in an urban airshed. Ph.D. Dissertation, Northwestern University, Evanston, Ill. Jazwinski A. H. (1970) Stochastic Processes and Filtering Theory. Academic Press, New York. Kontnik L. T. (1973) Personal communication. Environmental Protection Agency. Lamb R. G. and Neiburger M. (1971) An interim version of a generalized urban air pollution model. Atmospheric Environment 5, 23%264.
Peterson J. T. (1972) Calculations of SO, concentrations
over metropolitan
St. Louis. Atmospheric Environ-
ment 6.433-442.
Randerson D. (1970) A numeridal experiment in simulating the transport of sulfur dioxide through the atmosphere. Atmospheric Environment 4, 615-632. Roache P. J. (1972) Computational Fluid Dynamics. Hermosa Press, Albuquerque. Roberts J. J. et al. (1970) A multiple source atmospheric dispersion model. Argonne National Laboratory Report ANLIES-CC-007. Roberts J. J. et al. (1971) Chicago air pollution systems analysis program-Final report. Argonne National Laboratory Report ANL/ES-CC-009. Roth P. M., Reynolds S. D. and Seinfeld J. H. (1971) Development of a simulation model for estimating ground level concentrations of photochemical pollutants. Systems Applications, Inc., Report 71SIA-21. Sage A. P. and Melsa J. L. (1971) System Identijcation. Academic Press, New York. Shieh L. J. et al. (1972) Air quality diffusion model; application to New York City. IBM J. Res. Develop. 16 (2). 162-170. Shir C. C. (1972) Numerical investigation of the atmospheric dispersion of stack effluents. IBM J. Res. Develop. 16 (2), 172-179.
Shir C. C. and Shieh L. J. (1973) A generalized urban air pollution model and its application to the study of SO2 distributions in the St. Louis metropolitan area. IBM Research Report RJ 1227. Stern A. C. (Editor) (1968) Air Pollution, 2nd Rev. Edn., Volume 3. Academic Press, New York. Turner D. B. (1964) A diffusion model for an urban area. J. Appl. Meteor. 3 (1). 8391. Turner D. B. (1968) The diurnal and day-to-day variations of fuel usage for space heating in St. Louis, Missouri. Atmospheric Environment 2, 33%351.
Ensironmmt Vol 9. pp. 793-808. Pergamon Press 1975. Printed in Great Britain.
THE ADAPTIVE-FILTERING TRANSPORT MODEL FOR PREDICTION AND CONTROL OF POLLUTANT CONCENTRATION IN AN URBAN AIRSHED S. G. BANKOFF
and E. L. HANZEVACK*
Chemical Engineering Department, Northwestern
University, Evanston, Illinois 60201, U.S.A.
(First received 22 May 1974 and in jinal form 28 January 1975)
Abstract-An adaptive-filtering transport model for estimation and prediction of pollutant concentrations in an urban or regional airshed is developed. The program uses the method of fractional steps for direct numerical solution of the turbulent transport partial differential equation, coupled with a regressive least-squares filter which estimates both the concentrations and a spatially-distributed parameter for the adaptive feature at each time interval. The model is applied to historical data on SO,-concentrations in the Chicago area, and gives very good 4-h predictions.
1. INTRODUCTION
In this work we present the AFT (adaptive-filtering transport) model for estimation, prediction and control of sulfur dioxide concentrations in an urban airshed. Following the work of Roth et al. (1971) the pollutant transport is assumed to be governed by a parabolic partial differential equation in four independent variables. The solution method, called the method of fractional steps, relates this problem to the simpler problem of solving a sequence of three two-dimensional differential equations, which are solved successively at each time step. The new feature, introduced to air pollution modeling for the first time (Bankoff and Hanzevack, 1973) is the addition of an adaptive, or learning, algorithm to the overall model for on-line corrections to the concentration predictions and selected model parameters, based on incoming measurement data. This on-line self-correction scheme is a modification of Kalman filtering techniques which were initially developed for use in satellite tracking and other aerospace applications. This gives emphasis to the short term or incident control of urban air pollution, while most previous studies have emphasized long term aspects of air pollution problems. Recently, Desalu et al. (1974) have also employed Kalman filtering for dynamic estimation of air pollution. It differs in several important details from the one given here, including a direct Crank-Nicolson integration scheme, rather than the method of fractional steps found here to be stable and rapid; absence of a divergence-prevention scheme for the Kalman filter; a smoothing maximum likelihood algorithm for parameter estimation, instead of the augmented state vector used herein; and a more complex approach towards partitioning the filter and state equations. The predictions have been compared only with the theoretical plume solution from a single stack. One of the major difficulties of any air pollution model is that the primary meteorological parameters are incompletely measured, making the calculation of secondary parameters extremely difficult. The inclusion of a self-correction algorithm provides a means of overcoming this difficulty by continual re-estimation of one or more model parameters, For example, one of the explored parameters with the greatest sensitivity is the low level source term. This may be looked upon as a bias to the assumed zero-mean Gaussian additive noise in the transport equation, and is used to characterize parameter uncertainties and neglected nonlinearities. The identification of the parameter significantly improves model performance. * Present address: Exxon Research and Engineering Co.. Florham Park, New Jersey, U.S.A. 793
S. G. BANKOFFand E. L. HANZEVACK
794
The set of differential equations derived from the time-splitting method of fractional steps is considered to be a difference operator, advancing the state in time over the entire three-dimensional grid. These differential equations are normally deterministic, but the addition of the previously mentioned noise term makes them stochastic. Coupled to this is a vector measurement equation for SO2 concentrations at actual pollutant monitoring stations which is also stochastic because of measurement errors. Techniques to handle this problem as formulated above are well established by Kalman and others Aoki (1967), Jazwinski (1970), Sage and Melsa (1971). Kalman filtering represents a recursive on-line estimation procedure which can be viewed as an extension of classical least square regression techniques. It is basically a two-step process. Given the current estimate of pollution concentration profile and model parameters, a prediction is made of the same quantities at the next time step. This prediction is then compared with concentration measurements for that time as they become available, which results in a corrected posteriori estimate for all variables. An estimate of the error covariance is simultaneously computed by the matrix covariance equations, which constitutes the principal part of the computational burden. In the current work, Kalman filtering was applied to the estimation of SO2 concentrations at eight monitoring stations in the Chicago area during a short term pollution incident. Concentrations at all other points were obtained from the numerical integration of the model partial differential equation. Input data obtained from a study conducted at Argonne National Laboratory included point and area emission source information, hourly meteorological data, and hourly monitoring station pollution readings. 2. THEORY
2.1 Prediction
model
The starting point for development of a prediction model is the time-averaged tion for conservation of mass in turbulent flow (Bird et al., 1967):
ac t+jX
3 f3(v,c)+ a($?)
(j dx
F)
=il
(Djg)
+
R(cJ+
S(X1~X23X39t)
equa-
(l)
where c and 5j are the pollutant concentration and jth velocity component at the point (x1,xz,x3), averaged over a time interval which is large compared to the dominant time scale of the turbulent fluctuations, but small compared to the time scale for variations in the mean concentrations. Similarly, C’ and U’j are the instantaneous deviations of the concentration and velocity component at any point from their mean values. R(c) and S(x1,x2,x3,t) are distributed-parameter reaction rate and source terms, respectively. The initial condition is ChJ2,X3r&J)
and boundary (a) -K3g (b) (c)
= fh-bx,)
conditions are = Q(xI,x2,t)
= 0 C = g(q,x2,x3,t)
at x3 = h(x,, x2) (ground level) at x3 = H(xi, x2, t) (inversion base) at x1 = xW or xE (west or east) = xS or xN (south or north) (wtkhever boundaries are points of bulk inflow)
wheref(x,,x2,x3) and g(xI,x2,x3,t) are assumed known functions over the urban airshed domain. It is assumed that the flow is incompressible; molecular diffusion is negligible compared to turbulent diffusion, and only the vertical diffusion component need be considered; and the surface flux is negligible, except during precipitation.
795
Prediction and control of pollutant concentration
The following transformation of coordinates is applied to simplify the application of the boundary conditions (Roth et al., 1971): x=
x1 x.E - XIV
# = u1
))=-%XN- xs z
=
x3
-
v = vz
NXl,XZ)
w = v3
AH
where AH = H(x1,x2, t) - h(x,,x2). With these assumptions, the transformed equation becomes
with initial condition C(x, Y,z, t) = f(x, Y, 2)
at t = to
(3)
atz=O
(4)
and boundary conditions
-KC&J =az = 0
atz=
(5)
at inflow points on the sides of the unit square
C = Qk Y, z, t) Olx,yl
1
(6)
1.
Classical purely explicit methods for solving partial differential equations are avoided because of their traditionally poor stability characteristics. Similarly, purely implicit methods involve inversion of very large matrices for three dimensional grids, and therefore require excessive computer time and storage. However, variations of the method of alternating directions can avoid both of these difficulties. The method of fractional steps, which is a time-splitting method, is used here. Essentially, a four-dimensional PDE is broken into three two-dimensional PDEs which are then solved successively at each time step. This retains stability advantages of classical implicit methods, but greatly reduces the computational burden by requiring inversion of much smaller matrices. Even so, for the large scale system encountered in air pollution modeling, the computer storage requirements are excessive if each segment is solved implicitly because the, entire three-dimensional grid distribution must be stored for three different points in time (current and two previous). Therefore, a compromise is reached in which the vertical direction parabolic equation is solved implicitly, and the horizontal direction hyperbolic equations are then solved in turn explicitly. The product form of the fractional-steps recursion matrix insures that the entire process is stable if each fractional step is stable. The Crank-Nicolson finite-difference equation, which is unconditionally stable, is used for the implicit portion, and the explicit steps can be shown to be stable if the Courant number, aj 3 [(ojAt)/(Axj)] I 1 (Roache, 1972). Therefore, stability can be assured by appropriate choice of the time increment. The following series of finite-difference equations, which are second-order correct in both time and space variables, is used: _ C$k - Cyj, + At =
&tc~j,L'
2~~j,~,(C~j.t-l
+
Ryjk+ Syj,,
1 +
+
C,j,k+ 1 -
C,j,k- 1 -
2G$
Ctj,k- 1 -
-
2cYjk +
CY,j,k- 1)
Czj,k+l
+
G,j,k+
1) (7)
S. G. BANKOFF and E. L. HANZEVACK
796
with boundary conditions
-&(G’,j,l - Ci’,j,o)= QYj C,j,K
=
1,
G,j,K+
where K is the mesh point just under the inversion base;
where
_tw)~XF?-~(C:_ ll2,j.k
txE
f
l,j,k
c;k) -
+
%(c$k- CT-
l,j,k)
and uAt ’ = (xE - xW) AX
with boundary conditions At txE
q,;,k
= j(c?,j-
_
xw)~xF:!2,j.k = al12gfPd,k
1,k + c!,j+ l,k) + aI- l/2 CC,*_ l,j,k - 3(c?,j-
l,k
+
c?,j+
l,k)l
for u 2 Oand
_tw)AxFF+ 112,j.k=aI+ lpgr+
txE
c:$,k
= 3cc?,j+
1,k + c:,j-
I,k) + al,Z[+(c:,j-
112,j.k
1,k + c?,j+ l,k) -
cy,j.d
for 1.4< 0, where Z is the mesh point just west of xE; C””
= c$,*[l
+ /jj+1,2
- pj-l,2]
+
(xN
_A;s)&(G?_fG::+l.k)’ l,k
-
(9)
where At (xiv -
xs)AY
GF$_1 k s ---+ Bj-1j2 CW r,,’
1,k
+
c$:)
-2
and p=
vAt
h -
XSWY'
with boundary conditions analogous to those for equation (8). Equation (7) is a Crank-Nicolson formulation for the (z,t) fragment, and is solved implicitly. Equations (8 and 9) are second-order formulations for the horizontal advection ((x,t) and (y,t), respectively) written in conservative form and are solved explicitly. The three sets of equations solved consecutively constitute one time step. Due to the coarse grid spacing, the concentration distribution is often not smooth. The following additional conditions are imposed upon the system to prevent truncation error from increasing: F*_ = i
F*, =
l,j,k
ifa- IO, C,k < -FE ifa_ > 0, CL_1,j.k < F*_
ci+ l,j,k
ifa+ >O,C,k
cijk
ciCijk
Prediction and control of pollutant concentration
197
A similar condition is imposed in the y-direction. The treatment of physical parameters which are only incompletely measured is now detailed and the resulting approximations listed. The physical parameters necessary for model prediction are wind field (aj), inversion height (H), vertical turbulent eddy diffusivity (KS), reaction rate (R), and source emission inventory (s). The wind field must be specified at each of the three-dimensional grid points. As an idealized extreme, it would therefore be desirable to have an hourly measurement at each grid point. Obviously, this is never the case in real situations. At the other extreme, the minimum allowable information is one average bulk wind velocity vector for the entire urban area. Usually, however, an intermediate amount of information is available. Often 2-10 hourly ground level measurements are taken at various points throughout the city. These can either be interpolated spatially to obtain a ground level wind profile, or (preferably) a subjective analysis map of the data by an experienced meteorologist could be employed. In either case, the upper wind profile must still be specified. Roth et al. (1971) suggest using c?u~/~.x,= 0, i.e. the upper wind field is the same as the ground level profile. Alternatively, they suggest a two-step wind field with a discontinuity. Either of these methods may suffice for a rough approximation, but most authors favor a logarithmic or power law distribution. Shir and Shieh (1973) present a typical expression, IuJ = IuO&z/ z$ where the subscript 0 refers to ground level and the exponent p - In (Iuzl/lull)/ln (zJzJ where the subscripts refer to two heights where wind measurements are available (e.g. a television tower). The further restriction 0.15 c p < 0.65, is imposed. This latter formulation is used in validation studies when possible; one of the methods listed by Roth et al., is used when necessary. Furthermore, it is well known that there is an angle shift in the upper wind field, to the right of the surface wind in the northern hemisphere. However, since previous studies proved relatively insensitive to small angle variations (Roth et al., 1971; Shir and Shieh, 1971X this is considered a minor effect and is neglected here. The mixing depth is an important prediction model parameter, but there is little quantitative information about it available. A mixing height value is needed for each vertical column of grid volumes at each time interval, but, typically, measurements are taken at only two or three points, and only once or twice per day. Furthermore, the measured values are for temperature inversion-base heights, which may only roughly correspond to mixing heights. With this sparse data base, most studies (Roberts, 1970; Roth et al., 1971; Shir and Shieh, 1973) use a crude interpolation scheme in both space and time, constrained by assumed maximum and minimum values. The vertical turbulent eddy diffisivity cannot be directly measured, and so it must be calculated from other measurable quantities. One approach describes KZ as a component of a second-order tensor K. It is assumed that off-diagonal terms of this tensor are equal to zero, and if the mean velocity is parallel to the ground, this further implies that all other components can be set equal to zero, by virtue of the neglect of horizontal diffusion compared to horizontal advection. It is generally agreed that KZ is a function of wind velocity, shear stress and the change of temperature with height (lapse rate). However, the functional relationship is not well understood, and a variety of expressions have been proposed. A trapezoidal function at least qualitatively describes the expected behavior of K, with respect to height, and since fairly simple calculations are involved, a relation of that form is sought for specific situations. Further work in this area is also indicated. The reactions of SO2 in the atmosphere include catalytic conversion to H2S04, photochemical conversion to HzS04, and adsorption onto particulates. Of these the first reaction is considered dominant (Turner, 1964). Various researchers (Kontnik, 1973; Shir and Shieh, 1973; Stern, 1968) have listed the half life of atmospheric SO2 due to reaction between 3 and 4.5 h, with 4 h being used here as an intermediate value. An emission source inventory is one of the most basic and important components of a pollution dispersion model. Without a comprehensive picture of emission patterns, accurate short-term pollutant concentration predictions are nearly impossible. Therefore,
798
S. G.
BANKOFFand
E. L.
HANZEVACK
serious attention must be given to this rather mundane aspect of the problem. Most urban inventories are divided into point and area source strengths. Electric power plants and major industrial emitters are considered to be point sources. Hourly or bihourly emission values are computed based on regular process emissions plus a term calculated from the degree day concept for emissions due to heating. These values are additionally corrected to account for varying normal patterns in diurnal and weekly levels of activity. Smaller industrial plants, commercial establishments and residential heating are computed as area sources. Most empirical correlations are of the form where
S = YICI& + Y2C2 ww? $7, S = area source emission (lb h- ’ miles2) S, = annual emission total (tons y- ‘) T = temperature DD = degree days = T - 65°F if T < 65” = 0 if T 2 65” DD, = total annual degree days yl, y2 = hot water heating fraction, space heating fraction, y1 + y2 = 1 cl, c2 = unit conversion factors.
Typical examples can be found in Shir (1972) and Turner (1968). The foregoing represents the urban airshed pollution dispersion model, which now must be reformulated slightly in the next section in order to be amenable to attack by Kalman filtering, techniques. 2.2 Filter algorithm The Kalman filter is related to classical least squares techniques, the main difference lying in the recursive nature of the filter, while least squares methods are for batch processing of data. This difference is crucial, however, for this application. With the filter, new state and parameter estimates can be obtained on-line, as each new measurement set becomes available, by considering only the effect of this most recent observation, instead of re-solving the entire problem from t = 0 each time with a steadily increasing data base and corresponding computational load. In addition, the error covariance matrix, P(k), at the kth time interval, is continually updated, giving a quantitative estimate of model and filter performance. Consider a discrete dynamic system equation and associated measurement equation: x(k + 1) = cp[x(k), k] + T(k) w(k)
(10)
y(k) = H(k) x(k) + o(k),
(11)
where k0 4 k I k, represents the time interval under consideration. In these equations x(k) is a vector of n state variables at the kth time increment, which are, in fact, concentrations at each of the grid points throughout the airshed. Similarly, y(k) is a much smaller vector of concentration measurements at the same time. cp[x(k),k] is a vector of nonlinear functions resulting from the spatial discretization of the transport equation, which may also include time-dependent inputs, such as source terms. The stochastic inputs are the terms w(k), which represent the uncertainties due to neglected non-linearity, poorly-known parameters, model inaccuracies, and physical inputs of a random nature, while v(k) is the vector of measurement errors. H(k) and T(k) are appropriate coefficient matrices, here taken to be constant. We assume that w(k) and v(k) are uncorrelated Gaussian white noises, such that the expected value of w(k) is E(w(k)) = A(k). This represents a bias term in the process equation, which may be re-estimated at each time for the adaptive feature. Similarly, it is assumed that,
~C-+dl = A -QIdk)l= 0 ECx(kcJ xTWl = J=o
(12)
799
Prediction and control of pollutant concentration
E[w(k) wT(k)]= Q(k)G(k - j)
vT(j)]=
E[u(k)
R6(k
- j)
k
k 2 j
2 j,
where ~l~p,,,,P,, and R are assumed known prior statistics, and 6(k - j)is the Kronecker delta. Q(k) is also assumed to be known, but its norm is artificially increased at each time in order to prevent divergence of the filter. It is desired to obtain a minimum-variance recursive estimator for x(k) and to identify a parameter vector p(k). The state is augmented by defining an extended state vector X(k) =
in which the parameter
[ 1, ;;“kl
(13)
vector is assumed to evolve according to p(k + 1) = A(k) p(k) + B(k) w(k).
(14)
Then the extended state equation becomes X(k + 1) = @[X(k)&]
+ G(k) w(k),
where
kl A(k)p(k)
cp[x(k), p(k),
@[X(k), k] =
G(k) =[W BtkI 1a
(15)
1
Thus it is desired to minimize the cost function (Jazwinski, 1970) J = ;
+
i
IlWko) - Ako)ll;,-1 k,-
1
,zk
0
ClMk
+
1)
-
W
+
1)-W
+
1,llih
+
IlWll&~l 7
(16)
subject to the constraints implied by the extended state equation and prior statistics. The modtied extended Kalman filter with an extended state vector is the solution to this problem, yielding the following algorithm: i?(k + l/k) = @[2(k), M(k + 1) = GQ(k + l)GT + “;g;)]
k]
(17)
p(k+‘;%$)],
(18)
P(k + 1) = M(k + 1) - M(k $ 1) HT(k + 1) [H(k + l)M(k + l)HT(k + 1) + RI- ‘. .H(k + l)M(k + 1) (19) i(k
+ 1) = i(k
+ l/k) + M(k + l)HT(k + l)[H(k + l)M(k + l)HT(k + l)]- I.
. [y(k + 1) - H(k + 1) k(k + l/k)],
(20)
where k is the running discrete time index. The updating algorithm (I) is then: 1. 2. 3. 4. 5.
Use initial conditions P(k,) = PO, _$(k,) = px, Predict from (17). Compute prior variance from (18). Compute error variance from (19). Filter to get updated state and parameter estimates from (20).
This constitutes the extended Kalman filter with an augmented state vector, derived from the viewpoint of invariant imbedding (Sage and Melsa, 1971). Note that X(k/k - 1)
S. G.
800
BANKOFF and
E. L.
HANZEVACK
and _f(k) are state vectors consisting of a priori and a posteriori, respectively, concentration and parameter estimates at each of the n,nynz three-dimensional grid points at time k. Equation (20) comprises the updating step; equation (17) is the prediction step. An essential further step, in order to prevent divergence of the Kalman filter, is the addition of an algorithm for articially increasing the norm of the process noise matrix, Q(k), as outlined in Jazwinski. The full problem can lead to manipulation of filter matrices on the order of 4000 x 4000. This is clearly not desirable, so that the problem must be partitioned into related subproblems which are then connected in a second-level solution. Aoki (1967) presents the appropriate state vector partition theory. Let X(k) be an n-dimensional state vector partitioned into 1 subvectors Z(k)‘, . . . ,2(k)’ where Z(k)’ has dimension trj.
This leads to 1 relations of the form Z(ky’ = DjX(k) 1 5 j < 1. The submatrices Dj are taken to be time invariant. Thus, the same partitioning throughout, although this is not absolutely necessary. The prior and posterior state estimates are then given by k(k
+ l/k) = i @(k j= 1
is used
+ l/k)’
_$(k + 1) = i Fjk(k j= 1
+ ly
where
so that the state vector is reconstructed from the partitioned subvectors. Similarly, a matrix Bj is defined which picks out the subvector of the measurement vector Y(k) that is used in updating Z(ky’. The filter algorithm is now applied in turn to each of the subvectors Z(k) and the complete state vector X(k) is then reconstructed using the above relations. This process yields the filter i(k
-t l/k) = @[k(k),k],
(21)
M(k + 1) = GQ(k + l)GT + a;$$)b(k)a@T!2(k)1
(22)
8X(k)
i
j= 1
)
.M(k + 1) +
jil
F,K(k
1 T
E F’K(k + 1yGj H(k + 1) j= 1 ) + l)‘Gj Rj i f’jK(k + l)‘Gj j= 1 >(
T?(k + 1) = i?(k + l/k) + i j=
where
1
F’K(k + 1yGj H(k + 1)
1
F’K(k + ly’Bj[y(k + 1) - H(k + l)Z(k + l/k)],
K(k)j = filter gain matrix = M(kY’H(kY’[H(kY’M(kY’H(kY” + Rj]- ’ @j s Dj~Di’ where + indicates pseudoinverse Q(ky’ ~ DjQ(k)Dj H(ky E BjH(k)Dj+ M(k + ly’ = GQj(k)GT + @jJ’(k)@jT P(Oy z DjP(O)D,’ P(kY’ = [Z - K(ky’H(k)‘]M(ky’[Z - K(k)‘H(k)‘lT + K(k)‘RjK(kY’
and g(O). f’(O)are assumed known.
T,
(23)
(24)
Prediction and control of pollutant concentration
801
This process reduces the size of the matrices which must be manipulated. If the subvectors are not coupled through the state or measurement equations, then the above scheme is optimal. If the system is only weakly coupled, performance may still be tolerable. The specific partitioning scheme used in the present study, which yields a weakly coupled system, will now be described. A three-dimensional grid is superimposed on the urban area in order to solve the model equations by finite difference techniques. Eight x 16 x 4 nodes are used for the x-(west-east), y-(south-north) and z-(vertical) directions, respectively. This yields constant, uniform two-mile spacings in both horizontal directions; however, the vertical node spacing is neither constant nor uniform. The first two vertical nodes from the ground are assigned constant spacing between 100 and 200 ft, and the remaining distance to the mixing height is apportioned between the other nodes, yielding spacings typically between 100 and 1000 ft. Since most source inputs occur in the bottom two grid volumes, it is considered desirable to keep these volumes constant for simplicity. Certain ground level grid volumes contain a pollutant monitoring station. Consider a relatively small volume around one of these monitoring station grid volumes m, defined by a dia. r in the xy-(horizontal) plane and a height h in the z-(vertical) direction. r and h are chosen to be integers in units corresponding to grid dimensions, and r is odd for symmetry. This defines a semi-ellipsoidal subregion centered at m. Over a time arc of At (typically l-h), it is assumed that the state at points outside this subregion does not affect the state at m, in particular, or at other points inside the subregion. This is the basis for a series of partitioning schemes, depending on the value of r and h chosen. If r 2 3, h 2 2 the state at points outside the subregion affects the state estimate at m only to the second-order over a l-h time interval (i.e. for one filter application). If this process is repeated at each of 1 monitoring stations, then 1 subvectors zj, each of dimension r’h, are defined. The 1 subfilters operate on square matrices of dimension r2h, instead of dimension ~l,n,,~, resulting in a considerable computational saving. Note, however, that depending on the size of r and h, all of the grid points may not be included in one of the 1 subvectors by this particular partition. This is equivalent to considering all such excluded points to be in one z’+ ’ subvector in which the subfilter gain matrix K’+ 1 is always the null matrix. Therefore, between subregions, the parameter estimates at excluded points are spatially interpolated after each filter application in order to smooth the surfaces, primarily for computational reasons, since prediction estimates are relatively insensitive to small changes in these parameter values. Each monitoring station grid node and the nodes immediately adjacent (totaling 18 nodes) are considered to be the state subvector for one of the matrix subfilters. One or two parameters are added to this to give the 19- or 20-component, augmented state subvector. The constant (or slowly drifting) parameters represent: (a) A multiplicative source strength adjustment to account for the seasonal or diurnal source strength variations, (b) An additive bias or fictitious source parameter to identify the effect of neglected nonlinearities yielding a background pollutant level. For simplicity and computational savings, these parameters are assumed to be uniform throughout the filter subregion, although this is of course not necessary. We have thus developed a multistage solution intended to make maximum use of available information at each stage individually while maintaining a tolerable computational load. Note that there are several possible formulations to the predictorcorrector approach. The most exact method is to predict by solving the original PDE, and then filter using the entire state vector. However, as previously discussed, this leads to excessive computational requirements. Alternatively, the system may be partitioned into subregions in which a PDE is solved and a filter is applied to each of these subregions. This problem is computationally feasible, but is equivalent to assuming that the subregions are uncoupled over all time.
802
S. G. BANKOFFand E. L. HANZEVACK
Finally, the method actually used in this study employs partial uncoupling. The mass transport PDE is solved numerically for the entire system, and then the system is partitioned into subregions for filtering. This requires the much less stringent assumption that the subregions are approximately uncoupled over l-h time intervals. The same basic PDE is considered from two different viewpoints, over the whole urban region for prediction and locally in space and time for correction. 3.
MODEL
VALIDATION
Most raw data used in validation studies, meteorological measurements, source strength inventories and pollutant concentration monitoring station data, were obtained from studies conducted at Argonne National Laboratory for SO2 levels in the city of Chicago, Illinois during January 1967 (Roberts, 1970, 1971). Other pertinent references concerning source and meteorological correlations are listed. An hourly piecewise-constant ground level wind velocity field was calculated, based on an average of wind speed and direction measurements taken at Midway Airport and monitoring station No. 4. These continuous measurements are integrated over 5-min intervals every 15 min, to give an hourly average. The Norco objective mixing depth calculations (Roberts, 1971) were used for lid height values. These values are calculated from rawinsonde vertical temperature soundings taken at Peoria, Illinois and Green Bay, Wisconsin, at 600 and 1800 CST. The Chicago rural temperature profile is formed by linear interpolation. Then hourly rural profiles are calculated using temporal weighting for interpolation between the sounding times. The hourly mixing layer height is defined by the intersection of a dry adiabat from the urban surface measured at Midway Airport with the constructed rural (Argonne) temperature profile. Residential heating and industrial area source and power plant point source strength information was obtained from comprehensive source inventory questionnaires and calculations based on fuel consumption and ambient temperature. These values were adjusted for varying activity levels during different seasons of the year, day of the week and time of day or night. Pollutant concentration data were obtained from a network of eight telemetered air monitoring stations, which continuously and automatically record 5-min average SO2 concentrations at 15 min intervals throughout the day. Therefore, four 5-min averages are used to compute the hourly average at each location. Integration time increments of 5-6 min are found to be most satisfactory to meet stability, accuracy, and computational burden requirements. The initial concentration distribution necessary to start a simulation run was based on the monitoring station measurements taken at the time interval immediately preceding the run. These values were smoothed visually to yield a ground-level initial concentration value for each grid node. The ground-level initial concentration values were then multiplied by an arbitrary factor of 0.50-0.75 to establish estimates for the upper grid nodes. It was found that these estimates, although somewhat crude and arbitrary, were sufficiently accurate to yield good model performance because the sensitivity of model predictions to changes in initial concentration distribution was small after 3-4 simulation h. Other authors (Shir and Shieh, 1973) have also noted this phenomenon and have started simulation with their model assuming a uniform zero initial concentration distribution. Additional empirical relations for the upper wind field, vertical turbulent eddy diffisivity. and SO, reaction rate are necessary model inputs. A power law formula was used for upper wind field calculations (Shir and Shieh, 1973) IUI= Iv01@/%,‘?
(25)
where z. is the height at which the reference velocity u. is measured, with the value of the exponent p being on the order of 0.5. In addition, no angle shift is assumed in the calculations.
803
Prediction and control of pollutant concentration
The vertical turbulent eddy diffusivity is generally believed to be a function of temperature, height, atmospheric stability, and other variables. However, the exact nature of the functional dependence is largely unknown. Therefore, based on qualitative information (Roth et al., 1971), a trapezoidal function was used, approx. 0.039 km2 h- ’ at ground level and inversion height, and limited to twice that value at intermediate points. To avoid complicated atmospheric SO2 reaction mechanism expressions which would have only a relatively small effect compared to other terms in the model equations, a first-order decay rate was used for the reaction term. The half life for SO2 in the atmosphere was assumed to be approx. 4 h (Turner, 1964). 4. RESULTS
(a) Sensitivity. To determine the approximate sensitivity of model predictions to variations in input parameters, simulated data were computed by assuming that the transport equation models the real system exactly and that there are no errors associated with the finite difference solution of the equation. Twelve-hour simulation runs were then made with IO, 20 and 40 per cent average distortion factors built into each of the following inputs in turn: wind speed, wind direction, source strength, initial concentration profile, inversion height or mixing depth, and vertical eddy diffusivity. These inputs were considered to be the most weakly known and the most likely to have systematic errors. Table 1. Approximate relative sensitivity of 1 h concentration tions to input parameters
Parameter
Wind speed Wind direction Lid height Source strength Vertical eddy diffusivity Concentration at previous time step
predic-
Parameter variation (%) 10 20 40 Change in average concentration (%) 28 4 8 21 3 6 14 6 10 16 4 8 6 2 4 37 9 19
Table 1 indicates the approximate magnitude of the fractional change in calculated pollutant concentration for a specified change in input parameter. The model is most sensitive to the source strength and lid height with moderate parameter changes, although at high distortion levels it is more sensitive to both wind speed and direction. The model is least sensitive to vertical eddy diffusivity changes, It should be noted that an error in the pollutant concentration at the previous time step leads to approximately the same fractional error in the calculated concentration, as would be expected. (b) Filter efictiueness. To insure that the proposed suboptimal filtering scheme could perform efficiently, simulation runs were made attempting to filter out the previously mentioned simulated data distortions by means of a simple additive bias term. These studies can be briefly summarized by stating that in every case the filtered state estimates were consistently more accurate than the model predictions alone. When estimates strayed far from measured values, the resulting increase in noise covariance from the divergence prevention algorithm quickly (within l-2 h) brought estimates back into line. In addition, a slowly drifting bias parameter in many cases could be identified within a relatively short (12 h) run. In particular, initial concentration profile distortions of any size could quickly (within 3-5 h) be identified, involving a substantial bias term, as might be expected. Source strength distortions of 20 per cent could be filtered out of state estimates, although bias term identification took much longer.
S. G. BANKOFFand E. L. HANZEVACK
804
Measured
So2 Cone
(wm)
Fig. 1. Scatter plot of measured vs. predicted concentration
during high pollution
incidents.
For low velocities (less than approx. 8 miles h- ‘), moderate errors in wind speed and/or direction could be overcome in the state estimates, but bias terms tended to oscillate. At high wind speeds, accurate estimation using distorted values was more difficult. Fortunately, however, although current wind measuring devices are notoriously poor for low wind speeds, they are accurate at higher speeds. Tests with the other parameters resulted in similar recognizable patterns being observed. (c) Modellfilter studies. The combined prediction model and filtering scheme was examined by attempting l-, 2- and 4-h predictions using historical data records. In this context, 4-h prediction is taken to mean concentration predictions up to and including time t + 4 given only concentration measurement (and associated filtering) information up to time t. See Fig. 1 for a scatter plot of observed vs predicted concentration levels, and Figs. 2 and 3 for predicted and measured pollution values vs time. These figures show consistent agreement with measured pollutant concentrations at both high and low concentration levels, correctly predicting both trends and approximate magnitudes of peaks.
J 4
8
I2
TIME
Fig. 2. Comparison
16
20
z4
(h)
of measured and predicted SO2 concentrations incident.
during severe pollution
Prediction
0’0
and control
of pollutant
I 4
8
12
TIME
Fig.
3. Comparison
805
concentration
of measured
and
16
20
24
(h)
predicted SO2 incident.
concentrations
during
mild
pollution
Additional studies compared model performance using only current meteorological information for 4-h predictions compared to using actual meteorological (but not pollution) information from the historical data record. One can thus develop upper and lower bounds for the prediction accuracy, using at one extreme the actual meteorological records (corresponding to perfect weather prediction) and at the other extreme the procedure used herein, assuming no changes in wind or inversion height during the 4-h prediction interval. However, except in a few cases (Fig. 4) often it was found that model prediction was relatively insensitive to which set of meteorological information was used, not yielding a well-defined “envelope” (Fig. 5). This is understandable in view of the observed fact that most severe pollution incidents occur during periods of atmospheric stability, during which there are only relatively slow changes in wind velocity and mixing height. Longer-range predictions were then attempted in order to determine the prediction horizon. Eight-h predictions seem feasible, but 12-h predictions are often erratic. In practice, much longer data records would be available to the predictive model than in these simulations, so that even better performance can be anticipated.
TIME
(h)
Fig. 4. Comparison of predictions at time t with measured values of time t + 4; assuming (a) only current meterological information available. (b) meterological measurement data available from historical records during severe pollution incident.
S. G. BANKOFFand E. L. HANZEVACK
806
4
8
TIME
I*
12 (h)
Fig. 5. Comparison of predictions at time c with measured values of time t + 4; assuming (a) only current meterological information available, (b) meterological measurement data available from historical records during mild pollution incident.
All simulations were performed on a C.D.C. 6400 computer. For the forward integration step and update, less than 110 K of central storage and less than 30 s central processor time/stimulation hour were required. These figures compare favorably with those of models developed for other cities or pollutants. The filtering algorithm was responsible for approx. 20 per cent of the total simulation time, which appears to be a small price to pay for the dramatically improved results. (d) Model comparison. Table 2 summarizes the results of validation studies. It also indicates the relative performance, using the same data, of the Argonne National Laboratory Integrated Puff Model, previously the most sophisticated SO2 prediction model developed for Chicago. This table clearly shows the superiority of the APT model, even at the 4-h prediction level, compared to the one-hour ANL predictions. Figure 6 shows a typical comparison for a particular monitoring station. (e) Control strategies. Figure 7 shows the predicted concentration patterns assuming a 50 per cent source reduction at hour 11. This results in partial abatement of the incident, but points out the need for early control action. Other simulation runs were made using other feasible strategies, such as power plant source reduction only, source reduction only in areas of the city with predicted unacceptable concentrations, etc. However, such predictions cannot be validated by historical data but only by actual field tests. Nevertheless, since the model performs well with the historical data, it seems reasonable to assume that similar accuracy could be expected from control strategy predictions, and the implications for possible savings in money and improved air quality for short term pollution incident control are obvious.
Table 2. Comparison of model performance* Prediction time (h) AFT Model
1 2 4
Argonne National Lab. Integrated Puff Model
1
* 192 Observations (Roberts, 1970, 1971). I-Per cent of predictions within a factor of 2 of measured values. 2-Per cent of predictions within 0.05 ppm of measured vaiues.
Prediction error band 1 2 99 93 85 74
98 87 77 65
Prediction and control of pollutant concentration
. 0
4
? ‘:h,
8
16
20
24
TIME
Fig. 6. Comparison of 4 h AFT model predictions with 1 h integrated-puff
model predictions.
5. DISCUSSION
Sparse meteorological information and a rough source inventory used as inputs to an approximate (fractional steps) solution to the continuity equation and followed by a suboptimal filter have been shown to be sufficient to produce reasonably accurate and potentially useful short-term pollutant concentration predictions with only moderate computer storage and computation time requirements. The appropriateness and potential utility of short term prediction is indicated by historical data records (Roberts, 1971) which show that most serious air pollution episodes in Chicago range from about 12 to 72 h in duration, precisely the time scale in which the proposed model and filter scheme is most effective. Here an SO2 concentration of 0.3 ppm for a duration of four or more consecutive hours was arbitrarily used as a minimum standard for the existence of an air pollution episode. The most important variables in air pollution modeling are wind direction and speed mixing depth, ambient temperature, atmospheric stability (vertical turbulent eddy diffusivity), and hour of day. The relationships between these variables and pollution concentrations are generally not linear. Some of the functions are monotonic while others 57
4-
3-
2-
I -
O-
0
4
16
8
20
TIME I2 (h). SO, CONC (PPM) (X10)
Fig. 7. Predicted effect of source reduction on pollution concentrations.
808
S. G. BANKOFFand E. L. HANZEVACK
are cyclic. It is clear that improved meteorological measurement information would be useful in pollution prediction modeling. Specifically, more urban vertical temperature soundings (for mixing depth estimates) and wind measurements (for upper wind field estimates, including heat island and lake breeze effects) would be helpful. Some SO2 concentration measurements at selected elevated points throughout the city would aid in verifying the predicted vertical concentration profile. Also, better correlations for relating turbulent eddy diffusivity to measurable quantities are needed. Different filtering methods could be devised for specific uses. For example, different partitioning schemes could be tested or different parameters updated and identified. Alternatively, a sequential estimation algorithm could be developed in which certain parameters are identified for a specified period (say 12-24 h), then a different parameter set is identified for the next period, and the process could be repeated indefinitely. Additional 24-h average pollutant concentration measurements are taken at 23 points throughout the city of Chicago. An intriguing possibility would be to use these values to develop a two-frequency filter in which partitioning and filtering previously described is used to update the model every hour and a different partition and filter are used to update the model once every 24 h. This “interlocking” of the urban modeling area by means of two different partitioning schemes applied at different frequencies deserves theoretical as well as experimental attention. REFERENCES Aoki M. (1967) Optimization of Stochastic Systems. Academic Press, New York. Bankoff S. G. and Hanzevack E. L., Jr. (1973) Parameter updating for air pollution dispersion model. Symp. on Air Pollution Modeling, AIChE National Meeting, New Orleans, La. Bankoff S. G. and Hanzevack E. L., Jr. (1974) Real-time adaptive air pollution modeling. Presented at GVC/ AIChE Meeting, Munich, Germany. Bird R. B., Stewart W. E. and Lightfoot E. N. (1967) Transport Phenomena. Wiley, New York. De&u A. A., Gould L. A. and Schweppe F. C. (1974) Dynamic estimation of air pollution. JACC, Atlanta, Ga. Hanzevack E. L., Jr. (1974) A self-correcting model for prediction, estimation, and control of sulfur dioxide concentrations in an urban airshed. Ph.D. Dissertation, Northwestern University, Evanston, Ill. Jazwinski A. H. (1970) Stochastic Processes and Filtering Theory. Academic Press, New York. Kontnik L. T. (1973) Personal communication. Environmental Protection Agency. Lamb R. G. and Neiburger M. (1971) An interim version of a generalized urban air pollution model. Atmospheric Environment 5, 23%264.
Peterson J. T. (1972) Calculations of SO, concentrations
over metropolitan
St. Louis. Atmospheric Environ-
ment 6.433-442.
Randerson D. (1970) A numeridal experiment in simulating the transport of sulfur dioxide through the atmosphere. Atmospheric Environment 4, 615-632. Roache P. J. (1972) Computational Fluid Dynamics. Hermosa Press, Albuquerque. Roberts J. J. et al. (1970) A multiple source atmospheric dispersion model. Argonne National Laboratory Report ANLIES-CC-007. Roberts J. J. et al. (1971) Chicago air pollution systems analysis program-Final report. Argonne National Laboratory Report ANL/ES-CC-009. Roth P. M., Reynolds S. D. and Seinfeld J. H. (1971) Development of a simulation model for estimating ground level concentrations of photochemical pollutants. Systems Applications, Inc., Report 71SIA-21. Sage A. P. and Melsa J. L. (1971) System Identijcation. Academic Press, New York. Shieh L. J. et al. (1972) Air quality diffusion model; application to New York City. IBM J. Res. Develop. 16 (2). 162-170. Shir C. C. (1972) Numerical investigation of the atmospheric dispersion of stack effluents. IBM J. Res. Develop. 16 (2), 172-179.
Shir C. C. and Shieh L. J. (1973) A generalized urban air pollution model and its application to the study of SO2 distributions in the St. Louis metropolitan area. IBM Research Report RJ 1227. Stern A. C. (Editor) (1968) Air Pollution, 2nd Rev. Edn., Volume 3. Academic Press, New York. Turner D. B. (1964) A diffusion model for an urban area. J. Appl. Meteor. 3 (1). 8391. Turner D. B. (1968) The diurnal and day-to-day variations of fuel usage for space heating in St. Louis, Missouri. Atmospheric Environment 2, 33%351.