EVALUATION OF ENVIRONMENTAL IMPACT OF AIR POLLUTION SOURCES

This paper addresses the problem of evaluation and comparison of environmental impact of emission sources in the case of a complex, multisource emission field. The analysis is based on the forecasts of a short-term, dynamic dispersion model. The aim is to get a quantitative evaluation of the contribution of the selected sources according to the predefined, environmental cost function. The approach utilizes the optimal control technique for distributed parameter systems. The adjoint equation, related to the main transport equation of the forecasting model, is applied to calculate the sensitivity of the cost function to the emission intensity of the specified sources. An example implementation of a regional-scale, multilayer dynamic model of SOx transport is discussed as the main forecasting tool. The test computations have been performed for a set of the major power plants in a selected industrial region of Poland.

makers and are not supported by the respective models. However, some optimization methods and environmental models give the possibility of implementation of air pollution control strategies (Ciechanowicz et al. 1996;Holnicki et al., to appear).
For example, the optimal strategy for emission reduction (Holnicki et al., to appear) or real-time emission control can be worked out based on the respective model solution. The algorithms that solve such problems usually need a certain procedure to evaluate the contribution of the controlled emission sources in the final environmental damage. This paper addresses some related problems. We consider here a class of air pollution dispersion models that can be applied as tools to support decisions in environmental quality control by evaluating the contribution of each source in the resulting pollution field. Ultimately, those data can be utilized in decision making concerning air quality control.
It is assumed that the pollution transport process can be considered a distributed parameter system governed by the transport equation. Implementation discussed in the sequel is sulfur-oriented, but the approach can be applied in a more general class of the forecasting models. The governing model generates short-term forecasts of air pollution related to the specified, complex emission field.
Computation of the transport of sulfur pollution is carried out by a Lagrangian-type, three-layer trajectory model (Holnicki et al. 1994(Holnicki et al. , 2001. The mass balance for the pollutants is calculated for air parcels following the wind trajectories. The model takes into account two basic polluting components: primary SO 2 and secondary SO ¼ 4 . Transport equations include chemical transformations SO 2 ) SO ¼ 4 , dry deposition, and scavenging by precipitation.
The main output constitutes the concentrations of SO 2 , averaged over the discretization element and the vertical layer height. The governing equations have the following general form: along with the boundary conditions K @c @ñ n ¼ 0 on S À ¼ f@O Â h0; Ti jũ u Áñ n ! 0g outflow of the domain ð1bÞ 596 P. HOLNICKI and the initial condition where O is the domain considered with the boundary @O ¼ S þ [ S À ; (0, T) is the time interval of the forecast; c is the pollution concentration;ũ u is the wind velocity vector;ñ n is the normal outward vector of the domain O; K is the horizontal diffusion coefficient; g is the pollution reduction coefficient (due to deposition and chemical transformation); and Q is the total emission field. The emission field on the right side of Eq. (1) can be expressed as follows: where qðx; y; tÞ is the background (uncontrolled) emission field, q i ðtÞ is the emission intensity of the controlled, ith source, and w i ðx; yÞ is the characteristic function of the ith source location. The numerical algorithm is based on the discrete in time, finite element spatial approximation, combined with the method of characteristics (Holnicki 1995;Holnicki et al. 1994). The uniform space discretization step, h ¼ Dx ¼ Dy, is applied in the computational algorithm. The mass balance for the pollutants is calculated for air parcels following the wind trajectories. Points along the trajectory are determined at discrete time points based on the predefined interval t.
The air quality damage is represented by an environmental cost function of the following form: where 0 wðx; yÞ 1 denotes the area weight function, which represents the sensitivity of the domain to the specific type of air pollution. The component c ad is the admissible level of concentration, which should not be violated. Environmental cost function (3), due to transport equation (1), implicity depends on the emission intensity of the controlled sources. Information about the quantitative contribution of emission sources in the total pollution field is necessary in some cases (for example, the real-time emission control algorithms). A direct method of evaluation of this impact can utilize the consecutive reduction of emission level of the sources under question; the impact is represented by the related change of environmental cost index (3). In this approach, however, the main transport equation must be consecutively solved for all the sources considered. This means that, in the case of emission control, the most time-consuming step of the analysis has to be repeated many times.
Another approach, discussed in the sequel, utilizes optimal control techniques, namely, the properties of the adjoint equation (Lyons, 1971;Marchuk, 1995) related to state equation (1). Letq qðtÞ ¼ ½q 1 ðtÞ; . . . ; q N ðtÞ denote the vector function representing emissions of the controlled sources. Environmental cost index (3) can be considered a function ofq q, which means where subintegral function F, according to Eq. (3), is expressed as Fðx; y; tÞ ¼ wðx; yÞ ½maxð0; cðx; y; tÞ À c ad ðx; y; tÞÞ 2 To evaluate sensitivity of index (4) to the emission of the sources, q i ðtÞ ði ¼ 1; . . . ; NÞ, the gradient of this function must be computed. Assume a small change of emission level and denote by Es s the linear part of this change. Let Ep be the respective change of concentration level related to emission by Eq. (1). Thus, the respected disturbed values arẽ q q c ¼q q þ Es s and Consequently, the linear part of the variation of cost function (4) can be expressed in the form where, according to Eq. (4), the derivative of the subintegral function is as follows: @F @c ðx; y; tÞ ¼ 2wðx; yÞ Á max½0; cðx; y; tÞ À c ad ðx; y; tÞ ð5aÞ 598 P. HOLNICKI To calculate variation (5), the value of function p, which is not explicitly available, must be known. The following procedure, based on the optimal control theory, is presented, which calculates this variation in one simulation run of the transport equation. It is known (Lyons, 1971;Marchuk, 1995) that the minimum of index (4) is characterized as the solution of state equation (1) and the solution p Ã of the following adjoint equation: along with the boundary conditions and the final condition (for the end of the time interval) It must be noted that Eq. (6) is solved for the negative time and the reversed direction of wind, whereas the right-hand side is the derivative of the subintegral function of environmental index (4). It can be shown (Lyons, 1971) that, due to the specific form of the boundary conditions in Eqs. (1) and (4), the increment of cost function (5) can be expressed in the following form: where The last relation makes it possible to effectively calculate the increment of the cost function related to variation of the specific emission source. Our goal is to evaluate the contribution of each emission source in air quality deterioration. This contribution, in the sense of the measure of Eq. (4), can be calculated as the respective component of the following gradient function: Thus, to calculate the contribution of the emission sources one must successively perform the following steps: Solve the transport equation, Eq. (1). Solve the adjoint equation, Eq. (6), for the reversed time. Substitute the adjoint variable p Ã into Eq. (8) and calculate gradient components of environmental index (4).
To get the final solution, the transport and adjoint equations must be solved only once. The aforementioned method has been applied for the real-data case study concerning the Upper Silesia region. Evaluation of the environmental impact of the major power plants has been performed by the adjoint variable algorithm.

TEST COMPUTATIONS
The test calculations have been performed for the set of 27 major power plants in the region of Upper Silesia and Krako´w in Poland. Locations of the sources are indicated in Figures 1 and 2. For computational purposes, the rectangle domain of 110 km Â 74 km is discretized by the uniform grid, with the space step h ¼ 2 km. Grid coordinates and the main parameters of the controlled emission sources are presented in Table 1.
The aim of the computational experiment was (i) to evaluate and compare the environmental impact of each source by the adjoint variable method, discussed in the previous section and (ii) to examine the accuracy 600 P. HOLNICKI

602
P. HOLNICKI of this technique by comparing results to some reference data. The area sensitivity function wðx; yÞ introduced in relation (3) was set to 1 for the Krako´w urban area and 0 outside this region (compare the area indicated by the dashed line in Figures 1 and 2). Thus, the environmental impact of the sources under consideration was computed in the sense of deterioration of this domain.
To estimate the accuracy of the results obtained by the adjoint variable algorithm, the reference influence of the emission sources was calculated. The relative contribution of a specific source has been directly obtained as the solution problem (1), for the emission of this source reduced by 50%, with respect to the nominal emission intensity. This procedure was repeated in the sequel for all the sources, giving the reference contribution of each of them. Test computations were performed for a selected representative year (1996). The meteorological conditions are characterized by the respective sequence of 12-h sets of data. A one-year interval was split down into four three-month periods, and calculations were performed for four quarters, respectively. Selected numerical results are presented in Table 2. For a selected quarter, the neighboring columns compare the relative impact of source 22 (Skawina power plant) and the intermediate contribution of sources 10, 20, and 21. On the other hand, there is a group of sources of minor or negligible influence, in the sense of the assumed criterion function.
Since the dimension of the adjoint variable depends of the form of Eq. (6) and has no a direct physical meaning, two columns in Table 2 cannot be directly compared. The accuracy of the computational method discussed was evaluated by correlation between computed and reference data. The correlation of two sets of results is above R ¼ 0:97 for all cases considered. Figure 1 presents the resulting maps of a long-term forecast (for the last quarter) of SO 2 concentration and the related map of the adjoint variable distribution, calculated according to Eq. (6), for the reversed direction of time. Both maps show distributions of SO 2 and the adjoint variable averaged over the simulations period. The high values of the adjoint variable indicate the area of potentially high influence in the sense of index (3). This interpretation of the adjoint variable is more directly seen in Figure 2, which shows the respective maps for the short-term (3-day) forecast. It illustrates the relation between the SO 2 concentration map and the distribution of the respective adjoint variable.