- Email: [email protected]
A numerical model for wind field and pollutant concentration calculations over complex terrain. Application to Athens, Greece
Journal of Wind Engineering and Industrial Aerodynamics 73 (1998) 285—306
A numerical model for wind field and pollutant concentration calculations over complex terrain. Application to Athens, Greece J.S. Anagnostopoulos, G. Bergeles* Department of Mechanical Engineering, Laboratory of Aerodynamics, National Technical University of Athens, 9 Heroon Polytechniou Avenue, 15 773 Zografou, Athens, Greece Received 25 September 1995; revised 26 March 1997; accepted 10 September 1997
Abstract A non-hydrostatic, three-dimensional, prognostic, mesoscale model is developed for the calculation of the flow field and pollutant concentrations over complex terrain. The model solves numerically in an implicit, time-dependent manner the full Reynolds equations, in Cartesian, collocated grid arrangement. The differential equations are discretized by the finite volume method and solved iteratively with the SIMPLE algorithm. Turbulence closure is obtained by the two-equation (k—e) turbulence model. A technique of partially covered cells is introduced for an accurate representation of complex terrain features. The greater area of Athens, Greece, which combines significant topographic and sea-breeze aerodynamic complexities, is chosen for the model performance evaluation. The results reproduce reasonably well the diurnal variation of the wind and pollutant concentration fields, and the comparison with available observations is satisfactory. ( 1998 Elsevier Science B.V. All rights reserved. Keywords: Mesoscale model; Finite volume; Cartesian grid; Atmospheric pollution; Sea-breeze
1. Introduction The continuous growth of industrial and transportation activities during the last decades has resulted in a severe environmental degradation in many urban areas, as in the city of Athens where the dispersion of atmospheric pollutants is not favoured by the topographic or climate features. The three-dimensional and time-dependent wind field simulation is a prerequisite in order to develop a systematic strategy for air pollution control in such complex terrain.
* Corresponding author. 0167-6105/98/$19.00 ( 1998 Elsevier Science B.V. All rights reserved. PII S 0 1 6 7 - 6 1 0 5 ( 9 7 ) 0 0 2 8 6 - 9
286
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
Most of the models used today for complex mesoscale terrain analysis are based on the hydrostatic approximation for the pressure field [1—3]; however, this approach cannot handle complex terrain with slopes greater than 15%, and a smoothing is required in case of higher slopes. The number of non-hydrostatic models is limited in the literature. Spatial discretization is performed by the use of finite elements in the model of Lee et al. [4], whilst the finite differences or volume discretization are used by other models [5—9]. The time-dependent terms of the equations are introduced explicitly in the model of Moussiopoulos et al. [6]. However, this formulation restricts the maximum-allowed integration time step to few seconds, and moreover, the higher the grid resolution, the shorter the time step which is required [10]. A usual way of describing a complex terrain is by the use of non-orthogonal co-ordinate transformation techniques [11,12], which however increase the complexity of the equations and in some cases affect the numerical stability. Alternatively, Bartzis et al. [5] use the simple Cartesian co-ordinate system with parallel treatment of the ground as a porous medium. The aim of this paper is to present and evaluate a new prognostic numerical model capable of predicting the wind and pollutant concentration field at both the mesoscale and microscale level over terrain of any degree of complexity. Special care has been given during the model development, towards a relatively simple but nonetheless accurate computational procedure. For this purpose, the collocated grid approach [13] is employed instead of the staggered grid arrangement, that has been used in all previous models; also, the equations are solved in Cartesian co-ordinates, while an effective method of partially-covered cells is constructed to allow for an accurate representation of the complex topography. For the same reason, so as to obtain details of the wind field over certain areas, the model can employ a local grid refinement technique. All the previous features have been embodied in a innovative way, and it must be pointed that the final source code is smaller in length and faster than the conventional TEACH code computational procedure [14], on which it is based. The greater area of Athens was chosen as a test case, since it combines all the features characterizing a complex terrain: The wind field over the Athens basin is strongly influenced by the presence of surrounding mountains, hills and large water basins (Fig. 1). The local circulation is mainly driven by the sea breeze, which appears frequently, and during the summer months occurs on nearly 50% of all days [15]. The present code is used to simulate the wind field and inert pollutant concentrations for the day of 25 May 1991, for which some measurements are available and quasi-constant geostrophic wind is reported. Following a brief description of the numerical model, the predictions of the observed diurnal variation of the flow field and its effect on the pollutant dispersion are discussed.
2. Mathematical formulation The final form of the conservation equations describing the time-dependent, mesoscale wind flow are obtained following simplification and averaging techniques
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
287
Fig. 1. The modelling domain of the greater area of Athens. Residential areas are stippled, industrial areas are solid. Measuring stations: 1 — Philadelphia (altitude 150 m), 2 — Tatoi (210 m), 3 — Helliniko (20 m), 4 — Piraeus (80 m), 5 — Eleysina (30 m), 6 — Peristeri (80 m), 7 — Marousi (80 m), 8 — Geoponiki (60 m), 9 — Nea Smirni (50 m), 10 — Patision (100 m), 11 — Athinas (100 m).
presented in Ref. [11]. The solution provides the “layer-domain-averaged” values of a dependent variable U: U"U #UI , (1) 0 where U is assumed to represent the synoptic scale atmospheric conditions at a given 0 altitude, and UI are the mesoscale variations from this larger scale [11]. The conservation equations are expressed below, in Cartesian tensorial form. The anelastic form of mass continuity equation is adopted [16]: L (o u )"0, j"1,2,3. Lx 0 j j Momentum equations:
CA
(2)
B
D
Lu D(o u ) LpJ L Lu hI i # j !ou@ u@ #o g #f · (d u !d u ) 0 i "! # k i j 0 i 1 i1 2 i2 1 Lx Lx Dt Lx Lx h j i i j 0 !f d u #f d u , (3) 2 i1 3 2 i3 1 where f and f are the coriolis parameters: f "2X sin /, f "2X cos /; u is the 1 2 1 2 i velocity component at the i-direction, p is the pressure, o the density, k the dynamic
288
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
viscosity, h the potential temperature, g the gravity acceleration, d the Kronecker’s ij delta and D/Dt represents the substantial derivative. X is the angular velocity of the earth’s rotation and / the latitude. The tilde denotes fluctuating components. The potential temperature equation is written as follows — a similar equation describes the concentration field:
A
B
D(o h) L k Lh 0 " !ou@ h@ #Q, (4) j Dt Lx p Lx j T j where Q represents the source-sink term (surface heat flux for temperature, and mass flow rate of pollutant emissions for concentration equation, respectively). The potential temperature is solved instead of the thermodynamic one ¹, because it is a conservable quantity and its integration is independent of the path. A closed set of the above differential equations is obtained by the incorporation of the well-known turbulence model k—e [17] using the constant values of Table 1. The kinematic heat flux is modelled as follows: k Lh !ou@ h@" t . (5) j p Lx T j The transport equations for the turbulent kinetic energy k and its dissipation rate e take the form
A B A B
D(o k) L k Lk t 0 " #P#G!o e, (6) 0 Dt Lx p Lx j k j L k Le e e2 D(o e) t 0 " #C [P#G (1!C )]!C o , (7) 1k 3 2 0k Dt Lx p Lx j e j where P is the turbulent shear stress production and the term G represents the buoyant contribution to the production term of turbulent kinetic energy [18]: k Lh t . (8) i ¹ p Lx 0 T i In stable stratification G becomes a sink term, reducing the turbulence intensity, whereas in unstable stratification the buoyancy enhances turbulence, since G is positive. The parameter C controls the buoyancy effect on the e production. Rodi [18] used 3 the standard value of 0.8 for heated surface jets. For sea-breeze simulation, Kitada and G"!g
Table 1 Model parameters C k
C 1
C 2
p k
p e
p T
0.09
1.44
1.92
1
1, 3
0.9
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
289
Takagi [19] suggest two different values of C according to local thermal stratifica3 tion: 0 for an unstable and 1 for a stable layer. These values were adopted in the present work. The values of all the other constant coefficients involved in the above equations are summarized in Table 1.
3. Computational details The numerical code is used to predict the wind field and inert pollutant concentration variation during the day of 25 May 1991, for which measurements are available and a quasi-constant geostrophic wind was present over the greater area of the city of Athens. The 70]70 km modelling domain is shown in Fig. 1, taken from Ref. [5]. The computational space is divided in orthogonal cells by the use of a 39]39]30 Cartesian grid, with uniform horizontal resolution of 2 km. This relatively coarse resolution cannot secure the grid independence of the results; however, it was selected in order to fit the corresponding resolution of the available topographic and emissions data. The use of finer grids over the entire domain would result in an enormous increase of computational demands and for this reason the development and incorporation of local grid refinement techniques is planned for the future. The vertical grid-line spacing equals 30 m from zero to 400 m altitude, and then is gradually expanded up to 6 km. 3.1. The iterative solution The governing equations are integrated over finite volumes surrounding the grid nodes, and then discretized with the hybrid scheme approach [20]. All variables are stored at the centre of the control volumes, following the collocated grid arrangement [13]. After discretization, all equations are cast into the same general form, and the system is solved iteratively using the SIMPLE algorithm [21]. The discretization procedure is implicit in time; this formulation provides stability in the solution procedure and allows large integration time steps (several hours), independently of the grid resolution. A time step of 1 h is used in the present calculations, and thus 24 “runs” are performed to reproduce the flow field variation during the examined day. The resulting field of each “run” forms the initial conditions for the next one. The hourly interval was found adequate to produce reliable results: the use of a 3 h time step gave similar results under constant (stable) atmospheric stratification. Moreover, since the available emission data are averaged over hourly intervals, the use of smaller time steps would not present any advantage. Due to the strong stratification in the examined buoyant flows, some stability problems arose, which were handled by using small values of the pressure underelaxation factor (0.1) and a smooth introduction of the buoyancy term included in the vertical momentum equation.
290
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
Each “run” requires on the average 1500 iterations to converge, which corresponds to about 5 h execution time in a 586 Pentium PC (100 MHz), and needs about 10 Mbytes CPU memory. 3.2. Inlet conditions At the inlet boundaries, the vertical distribution of the wind velocity and turbulent properties must be prescribed. In the surface layer this was done using the relationships of Ref. [11], in which a stability function t is introduced into the logarithmic law: u"u /i[1n(z/z )!W(z/¸)], * 0 h u3 ¸"! w * , gw@h@i
(9) (10)
where h is the surface potential temperature, z the surface roughness length, u the 8 0 * friction velocity, g the acceleration due to gravity and i the Von Karman constant ("0.4187). The surface heat flux w@h@ was calculated as described in Section 3.3 below. In the transition layer the wind velocity and direction was obtained from the Ekman equations [11]. The heights of the surface and the planetary layer were taken from Refs. [11,22], respectively, as function of the stratification conditions and the time of the day. The vertical distribution of the turbulence parameters (k, e) throughout the whole planetary boundary layer was calculated using the expressions given in Ref. [23], depending on the prevailing stability conditions. Above the planetary layer and at the upper boundary of the computational domain (6 km), the vertical velocity component is assumed zero, while the geostrophic northwest wind of 3.5 m s~1 is imposed for the horizontal component. Based on the observations, the synoptic field is kept constant during the day. Entrainment conditions are applied to all the other free boundaries of the domain. This is a simple numerical technique to account for situations when the developed pressure field inside the computational domain causes inflow conditions at a free boundary. According to this technique, the amount of the entering mass and the corresponding inlet velocities are calculated from a mass balance over the faces of each cell adjacent to a free boundary. The entering flow mass and momentum as well as the rest ambient fluid characteristics (temperature, turbulence quantities, etc.) are then included as extra source terms into the corresponding equations. The inlet temperature profile at a given time of the day is estimated via interpolation between the two available radiosonde measurements at 02:00 and 14:00 LST. 3.3. Ground boundary conditions The inability of the Cartesian grid to follow the exact surface configuration is eliminated with the development of a partially-covered cell technique. According to
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
291
Fig. 2. Terrain simulation in the computational grid.
this, a computational cell that is partly under the ground surface, is solved after multiplying its areas and volume by proper geometrical coefficients (y and y ) i 7 that represent their free (above the ground) portion. For example, in the simplified two-dimensional sketch of Fig. 2, the east surface of the (I,J) cell should be multiplied by the factor y "dy/D½, while the south surface by the factor y " E S dx/Dx. Also, the volume of the cell is multiplied by a factor y , equal to the free volume 7 portion. Thus, y"1 denotes a completely free and y"0 a totally blocked cell face or volume. The significant improvement of the ground boundary description achieved by the above technique is demonstrated in Fig. 2: It is obvious that the standard Cartesian step approximation may lead to substantial deviations from the actual boundary, which in contrast, can be simulated quite accurately using the new technique. Using these new coefficients, the general form of the discretized equations becomes (A !c S )U "+ c A U #c S , 1 7 1 P i i i 7 U i A "+ c A , i"East, West, North, South, Up, and Down face, (11) P i i i where A are the coefficients linking a dependent variable of a cell to its neighbours, on i the adjacent grid volumes. All the necessary geometrical coefficients are calculated and stored once at the beginning of the iterative procedure, and only for the cells adjacent to the ground. This practice reduces the additional storage requirements, without increasing the total execution time. Also, the wall function boundary conditions can be applied in a similar manner to all these cells, thus, considerably simplifying the algorithm.
292
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
Moreover, the actual slope and the exact area of the bounding surfaces are both calculated accurately. The surface heat flux is calculated from an energy balance on the ground surface using algebraic expressions to obtain the solar radiation as function of the incident radiation angle, which in turn depends mainly on the time of the day, as also on several other parameters (latitude, longitude, day of the year, cloudiness, surface temperature [24]). Typical values for the various constants involved were taken [8,25]. A typical daily variation (based on the measured data) was assumed for the surface temperature to enhance stability, whereas the actual slope of the ground at each location was not considered. Therefore, the calculated heat flux is uniform over the ground surface at any specific time of the day. A uniform ground roughness was also assumed for the land surface, whereas the sea surface was considered adiabatic [8]. 3.4. Code starting The time-dependent solution of the differential equations requires the knowledge of an initial flow field over the entire calculation domain. This field is obtained at 00:00 by the following practice: The code starts at 21:00 of 24 May 91, namely 3 h earlier, assuming a neutral atmosphere, same synoptic wind and constant temperature at levels parallel to the sea surface, which is the most convenient solution in case of unknown temperature distribution [26]. For the next 3 h the time-dependent solution is performed using a typical night-time heat flux at the ground. Thus, the resulting field approaches the real conditions prevailing at the beginning of the new day, and can be used as the initial state of the atmosphere at 00:00.
4. Results and discussion Meteorological measurements available for the examined day of 25 May 92 consist of wind and temperature vertical profiles at Helliniko airport station every 12 h, and surface (10 m) temperature and wind speeds at five stations every 3 h. The pollutant concentration pattern is available at seven different surface locations, at hourly intervals. The emission database used for the present investigation has been created during the APSIS exercise [27,28] for model evaluation purposes. These data include mean hourly values for various pollutants (CO, SO , NO, NO , HCHO and other hydro2 2 carbons), averaged in a 2]2 km2 space resolution over the entire terrain, and they are given in two layers: 0—30 m and 30—60 m above the ground. For this reason, buoyancy effects or plume rising phenomena (i.e. from stack emissions) were not considered. The numerical investigation concerns the dispersion characteristics of the relatively inert pollutants (CO and SO ), to avoid at this stage the modelling of the complex 2 photochemical kinetics. Also, the calculations were restricted to one day in order to check the reliability of the results and to spot any numerically imposed limitations affecting the validity of the predictions, before trying to assess the overall performance by applying the algorithm to different days or situations.
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
293
In order to make clear the various topographic and thermal influences, the general characteristics of the diurnal variation of the flow field are first demonstrated, while the comparison with observations follows. 4.1. Diurnal flow field variation The predicted diurnal variation of the flow field and pollutant concentration on the ground surface is illustrated in Figs. 3—10. Calculations are made for every hour of the day, but only the most representative results are reported and discussed here. The atmosphere is clearly stable during the nighttime (Fig. 3). The cooling of the air near the mountain surfaces generates downward breezes. This phenomenon is evident in the predictions of Fig. 4a. Air masses from the surrounding mountains of Parnis, Pendeli and Hymettus are accumulated in a main land breeze directed from the interior of the basin offshore. This northeast breeze protects the basin against the SO 2 pollution which is emitted mainly at the industrial zones of Eleysina and Aspropirgos (Fig. 4b). The CO concentration is also kept low during the night (not exceeding 2 mg m~3, Fig. 4c), since car traffic is light. The wind field and pollution over the Athens basin remains with minor changes as described above, until the first daylight. At 08:00 the land breeze is still prevailing (Fig. 5a), but the pollutant concentration exhibits appreciable changes and both SO 2 and CO values are considerably increased, especially at the centre of Athens and Piraeus, where traffic is now heavy (Fig. 5b and Fig. 5c). SO and CO concentration 2 exceeds 100 lg m~3 and 5 mg m~3 respectively, over these areas, and it is noticeable that the pollutants are transferred far out to the sea, thus substantial surface concentrations are predicted up to the nearby Aegina island. After sunrise, the temperature field changes radically due to the rapid heating of the land surface (Fig. 3). The stable stratification is partially destroyed, and transition to
Fig. 3. Predicted diurnal variation of the temperature profile at the Helliniko station.
294
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
Fig. 4. Predicted surface flow at 04:00 LST: (a) horizontal velocity vectors (scale: 1 mm"2 m/s); (b) SO 2 concentration (lg m~3); (c) CO concentration (mg m~3).
an unstable atmosphere takes place. As a consequence, the land breeze ceases to exist and the synoptic northwest wind field initially dominates (Fig. 6a). At 10:00, the air stagnation above the basin causes the CO and SO concentrations to reach their 2 maximum values, which are above 9 mg m~3 and 200 mg m~3, respectively, over the centre of Athens and Piraeus, while the surface distribution patterns remain similar to these at 08:00 (Fig. 5b and Fig. 5c). The continuation of land surface heating leads soon to the establishment of upward air streams along the mountainsides. These streams, in conjunction with the seabreeze formation, due to the relatively slower rate of temperature rising over the sea surface, result in the development of a completely different aerodynamic field over the basin at 14:00, as can be seen in Fig. 7a. The most interesting characteristic of this complicated field is the formation of a distinct line extending along the peninsula, where the fronts of two opposite sea-breeze circulation systems meet. The existence of the northern part of this line, between the mountains of Pendeli and Hymettus, has been observed as a common situation at this time of the day, and it seems to block the physical channel for the pollutants to be transported away from the basin. This
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
295
Fig. 5. Predicted surface flow at 08:00 LST: (a) horizontal velocity vectors (scale: 1 mm"2 m/s); (b) SO 2 concentration (lg m~3); (c) CO concentration (mg m~3).
Fig. 6. Predicted surface horizontal velocity field at 10:00 LST (scale: 1 mm"2 m/s).
296
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
Fig. 7. Predicted surface flow at 14:00 LST: (a) horizontal velocity vectors (scale: 1 mm"2 m/s); (b) SO 2 concentration (lg m~3); (c) CO concentration (mg m~3).
Fig. 8. Predicted surface horizontal velocity field at 18:00 LST (scale: 1 mm"2 m/s).
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
297
Fig. 9. Predicted surface horizontal velocity field at 24:00 LST (scale: 1 mm"2 m/s).
Fig. 10. Predicted surface pollutant concentration levels at 22:00 LST: (a) SO in lg m~3; (b) CO in 2 mg m~3.
blockage may be responsible for the dramatic increase of pollutant concentrations that occur when strong temperature inversion conditions exist over the basin and the upward streams are weak. This is not the case in the examined situation, and therefore, the predicted values of SO and CO shown in Fig. 7b and Fig. 7c are both 2 decreased due to unstable conditions. The polluted area is now split in two parts, with the one over the sea diluting due to diffusion, since no more pollutants are transferred there. The relation between land and sea surface temperature is reversed in the evening, due to the faster cooling of the ground (Fig. 3). A weak sea breeze still exists at 18:00 (Fig. 8), but the synoptic wind field dominates. A down-slope breeze has already been established at 22:00, while the aerodynamic field at 24:00 (Fig. 9) is almost the same as
298
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
the one at the beginning of the day (Fig. 4). This is an expected result, since the synoptic field is kept constant during the whole day, and yet confirms the accurate operation of the entrainment boundary conditions used in the computer code. The wind direction is uncertain over the basin during the re-transition from unstable to stable stratification, while traffic is still dense. Consequently, the pollutant concentrations start to increase again, and obtain high maximum values at about 22:00 (Fig. 10a and Fig. 10b), approaching the morning maxima. 4.2. Comparison with available measurements The daytime temperature vertical profile is compared with measurements at Helliniko in Fig. 11a. The agreement is satisfactory, except in the region between 300
Fig. 11. Comparison between measurements and predictions at Helliniko station, at 14:00 LST: (a) Temperature profile; (b) west wind component; (c) south wind component.
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
299
Fig. 12. Diurnal variation of surface temperature at five measuring stations.
and 600 m above the ground, where an inversion layer is indicated but not reproduced by the numerical model. A possible reason for this discrepancy is that the temperature profile for a sea-shore location as Helliniko is influenced by the sea surface temperature spatial variability [27], that is not discernible by the horizontal resolution of 2 km used. Unfortunately, no additional measurements of the vertical temperature profile are available at more surface stations in order to estimate the extent of this inversion and the influence of the above discrepancy on the air quality predictions. Similar are the predictions of Varvayanni et al. [28] that are also plotted in Fig. 11a, whereas an inversion layer like the measured one can be observed in the results of Thunis et al. [27], but temperatures are strongly underpredicted.
300
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
Fig. 13. Diurnal variation of west wind component at five measuring stations.
The calculated versus the measured profiles of mean wind speed at the same station and 14:00 LST, are shown in Fig. 11b and Fig. 11c. The results are considered acceptable, mainly in their qualitative behaviour concerning the wind direction. The predicted surface temperatures are in reasonable agreement with the measured ones at the low-altitude stations (Helliniko, Piraeus and Eleysina), as can be seen in Fig. 12. The discrepancies can be attributed to the approximate estimation of the thermal air—ground interaction (constant surface roughness length and a typical diurnal variation of surface heat flux), adopted in the model. More pronounced are the disagreements at the stations of Philadelphia and Tatoi at the foot of the mountain
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
301
Fig. 14. Diurnal variation of south wind component at five measuring stations.
where the dynamic effects are strong, indicating a difference between calculated and actual time variation of the convergence line location. Unfortunately, no measurements are available on the mountain slopes, thus making difficult the evaluation of the calculations concerning the up-slope and down-slope streams developed there. The percentage error, (¹ -¹ )/*¹ , ranges from almost zero to less than 45% in 13%$. .%!4. .%!4. all the above cases, except at 11:00 at Helliniko, where it exceeds 100%. The surface wind field is reported in Figs. 13 and 14 at the same five meteorological stations. The calculated velocity components agree fairly well with the observed ones. The model reproduces the low nighttime wind speeds, while non-systematic discrepancies can be observed at various stations and times of the day because the
302
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
Fig. 15. Diurnal variation of SO concentration level at seven measuring stations. 2
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
Fig. 16. Diurnal variation of CO concentration level at seven measuring stations.
303
304
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
2]2 km2 grid used cannot simulate the smaller-scale surface characteristics which influence the local measurements. Some more systematic deviations occurring in regions near the boundaries of the calculation domain (i.e. Eleysina, Fig. 14e), can be assigned to inaccuracies in the inlet conditions used. It must be also noted that the velocities are calculated at the centre of the first solvable grid volume above the ground, where the exact altitude cannot be calculated from the available data for the four neighbouring corners of the grid volume, and may not correspond well to the measured position (10 m). The most decisive evaluation point concerns the ability of the model to predict the diurnal pollutant concentration at ground level. The numerical results are plotted in Figs. 15 and 16 for SO and CO, respectively, and show a reasonable agreement with 2 the observed values, at least in the shape of the distribution and the timing of the concentration maxima. As for the quantitative differences, it must be noticed that the numerical results represent averaged values inside relatively large computational cells (2 km]2 km]30 m), whereas the measurements are taken at certain fixed points. It is obvious, for example, that the significant underprediction of pollutant concentrations at Patision station, is due to this station located just over one of the busiest streets in Athens. In general, all the diurnal concentration curves exhibit two maxima (Figs. 15 and 16): the first appears between 08:00 and 10:00 LST, because traffic is enhanced at that period, whilst stagnant flow is created over the basin, due to the change from a landbreeze to sea-breeze (Fig. 6). The unstable atmosphere which is developed and established afterwards, contributes to a significant reduction of the pollution at all measuring stations. This reduction continues till about 18:00. From that time, the wind above the basin starts slowing down due to the change over from sea-breeze to land-breeze with the appearance of the first stable layer aloft, allowing thus again the accumulation of the produced pollutants. A second peak is therefore formed between 20:00 and 22:00 LST, but with a smaller intensity in most cases. Some corresponding results, obtained by the explicit solution of Ref. [6], are also included in Fig. 16, and exhibit analogous variability and similar average degree of agreement with the measurements, although the integration time-step was very small (several seconds). Therefore, the hourly integration used for the present calculations can be considered adequate to reproduce the diurnal variation of the pollutant concentrations.
5. Conclusions A new mesoscale numerical model has been used to predict the wind field and inert pollutant concentrations over the complex terrain of the greater area of Athens, during the day of 25 May 91. A technique for the solution of partially-covered cells, in conjunction with the collocated Cartesian grid arrangement, achieved a close representation of the complex topography. The results reveal the great influence of the surrounding mountains and sea on the flow field characteristics. Thermally driven down-slope and up-slope air streams are
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
305
developed during the whole day along the mountainsides, while two strong sea breeze systems are established during the day-hours from the opposite sea-shores of the peninsula. The fronts meet and block the natural ventilation of the basin. The resulting aerodynamic field is extremely complex and dominates over the weak synoptic wind, thus controlling the dispersion of pollutants. Pollutant concentrations exhibit at all stations two maxima during the day: the first and higher around 10:00 and the second at about 22:00, both due to a combination of increased traffic and the change over from sea breeze to land breeze and vice versa, associated with the change in the near-ground atmospheric stability. Considering this high degree of complexity, the comparison of the model predictions with available observations is relatively satisfactory and the general diurnal variation follows in most cases the corresponding measurements. The consideration of smaller-scale ground surface characteristics (local roughness length, emissivity, surface slope, etc.), along with the planned incorporation of local grid refining techniques, higher-order differencing schemes and a photochemical reaction model, should provide more detailed and accurate information in order to achieve a better understanding of the air pollution problem. Additional surface station measurements are however necessary, especially in the northern mountainous region of the area examined, in order to evaluate the model performance at the microscale level.
References [1] P.C. Patnaik, A preliminary users guide for the SIGMET Mesoscale Meteorology Code, SAI Report SAI-78-768-LJ, RLO/2440-77-9, San Diego, CA, 1979. [2] N.L. Seaman, R.A. Anthes, A mesoscale semi-implicit numerical model, Quart. J. Roy. Meteorol. Soc. 107 (1981) 167—190. [3] M. Tjernstrom, Numerical simulations of stratiform boundary layer clouds on the Meso-c-Scale, Part I: the influence of terrain height differences, Boundary Layer Meteorol. 44 (1988) 33. [4] R.L. Lee, J.H. Leone, D.M. Gresho, S.T. Chan, A finite element model for environmental problems involving complex terrain, in: Proc. 4th Int. Conf. of Finite Elements in Flow Problems, Tokyo, Japan, 1982. [5] J.G. Bartzis, A.G. Venetsanos, M. Varvayanni, N. Catsaros, A. Megaritou, ADREA-I: A threedimensional transient transport code for complex terrain and other applications, Nucl. Technol. 94 (1991) 135—148. [6] N. Moussiopoulos, Th. Flassak, P. Sahm, D. Berlowitz, Simulations of the wind field in Athens with the nonhydrostatic mesoscale model MEMO, Environ. Software 8 (1993). [7] R. Kunz, N. Moussiopoulos, Simulation of the wind field in Athens using refined boundary conditions, Atmos. Environ. (1995), in press. [8] D. Tryfonopoulos, G. Bergeles, Temperature fields and air pollution build-up over the Athens basin, Environ. Software 9 (1994) 269—283. [9] J. Glekas, G. Bergeles, Dispersion under neutral atmospheric conditions, Int. J. Numer. Meth. Fluids 19 (1994) 237—257. [10] T. Flassak, M. Wortmann, N. Moussiopoulos, Influence of horizontal grid resolution on simulated wind and pollution concentration field in the greater Athens area, in: N. Moussiopoulos, G. Kaiser (Ed.), Proc. Seminar on Monitoring and Modelling in the Mesoscale, Thessaloniki, 1991, pp. 155—164. [11] R.A. Pielke, Mesoscale Meteorological Modeling, Academic Press, London, 1984.
306
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
[12] J. Glekas J, G. Bergeles, N. Athanasiadis, Numerical solution of the transport equation for passive contaminants in three-dimensional complex terrain, Int. J. Numer. Meth. Fluids 7 (1987) 319—335. [13] C. Rhie, W. Chow, Numerical study of the turbulent flow past an airfoil with trailing edge separation, A.I.A.A. J. 21 (1983) 1525—1532. [14] A.D. Gosman, F.J.K. Ideriah, TEACH-2E: A general computer program for two-dimensional turbulent recirculating flows, Report, Department of Mechanical Engineering Imperial College, 1976. [15] D.P. Lalas, D.N. Asimakopoulos, D.G. Deligiorgi, C.G. Helmis, Sea-breeze circulation and photochemical pollution in Athens, Greece, Atmos. Environ. 17 (1983) 1621—1632. [16] F.B. Lipps, R.S. Hemler, A scale analysis of deep moist convection and some related numerical calculations, J. Atmos. Sci. 39 (1982) 2192—2210. [17] B.E. Launder, D.B. Spalding, The numerical computation of turbulent flows, Comput. Meth. Appl. Mech. Eng. 3 (1974) 269—289. [18] W. Rodi, Turbulence Models and their Application in Hydraulics. A State of the Art Review, Int. Assoc. Hydraulic Research Publ., 1980. [19] T. Kitada, H. Takagi, Some remarks on the k—e turbulence model applied to sea breeze simulation: Buoyancy effect of the e-equation and horizontal eddy diffusivity, in: N. Moussiopoulos, G. Kaiser (Ed.), Proc. Seminar on Monitoring and Modelling in the Mesoscale, Thessaloniki, 1991 pp. 135—146. [20] A.D. Gosman, W.M. Pun, A.K. Ruchal, D.B. Spalding, Heat and Mass Transfer in Recirculating Flows, Academic Press, London, 1969. [21] S.V. Patankar, D.B. Spalding, A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows, Int. J. Heat Mass Transfer 15 (1972) 1787—1792. [22] F.B. Smith, Atmospheric turbulent transport, in: S. Sandroni (Ed.), Course on Regional and Longrange Transport of Air Pollution, Ispra, Italy, Elsevier, Amsterdam, 1986, pp. 43—69. [23] S.R. Hanna, G.A. Briggs, R.P. Hosker, Handbook on Atmospheric Diffusion, Technical Information Center, U.S. Department of Energy, 1982. [24] A.A.M. Holtslag, A.P. Van Ulden, A simple scheme for daytime estimates of the surface fluxes from routine weather data, J. Clim. Appl. Meteorol. 22 (1983) 517—529. [25] D. Tryfonopoulos, Numerical solution of atmospheric flows over complex terrain, Ph.D. Thesis, National Technical University of Athens, 1992. [26] R.T. McNider, R.A. Pielke, Diurnal boundary layer development over sloping terrain, J. Atmos. Sci. 38 (1981) 2198—2212. [27] P. Thunis, P. Grossi, G. Graziani, H. Gallee, B. Moyaux, G. Schayes, Preliminary simulation of the flow field over the Attic Peninsula, Environ. Software 8 (1993). [28] M. Varvayanni, J.G. Bartzis, C.G. Helmis, D.N. Asimakopoulos, Simulation of the sea breeze under opposing synoptic conditions, Environ. Software 8 (1993).
A numerical model for wind field and pollutant concentration calculations over complex terrain. Application to Athens, Greece J.S. Anagnostopoulos, G. Bergeles* Department of Mechanical Engineering, Laboratory of Aerodynamics, National Technical University of Athens, 9 Heroon Polytechniou Avenue, 15 773 Zografou, Athens, Greece Received 25 September 1995; revised 26 March 1997; accepted 10 September 1997
Abstract A non-hydrostatic, three-dimensional, prognostic, mesoscale model is developed for the calculation of the flow field and pollutant concentrations over complex terrain. The model solves numerically in an implicit, time-dependent manner the full Reynolds equations, in Cartesian, collocated grid arrangement. The differential equations are discretized by the finite volume method and solved iteratively with the SIMPLE algorithm. Turbulence closure is obtained by the two-equation (k—e) turbulence model. A technique of partially covered cells is introduced for an accurate representation of complex terrain features. The greater area of Athens, Greece, which combines significant topographic and sea-breeze aerodynamic complexities, is chosen for the model performance evaluation. The results reproduce reasonably well the diurnal variation of the wind and pollutant concentration fields, and the comparison with available observations is satisfactory. ( 1998 Elsevier Science B.V. All rights reserved. Keywords: Mesoscale model; Finite volume; Cartesian grid; Atmospheric pollution; Sea-breeze
1. Introduction The continuous growth of industrial and transportation activities during the last decades has resulted in a severe environmental degradation in many urban areas, as in the city of Athens where the dispersion of atmospheric pollutants is not favoured by the topographic or climate features. The three-dimensional and time-dependent wind field simulation is a prerequisite in order to develop a systematic strategy for air pollution control in such complex terrain.
* Corresponding author. 0167-6105/98/$19.00 ( 1998 Elsevier Science B.V. All rights reserved. PII S 0 1 6 7 - 6 1 0 5 ( 9 7 ) 0 0 2 8 6 - 9
286
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
Most of the models used today for complex mesoscale terrain analysis are based on the hydrostatic approximation for the pressure field [1—3]; however, this approach cannot handle complex terrain with slopes greater than 15%, and a smoothing is required in case of higher slopes. The number of non-hydrostatic models is limited in the literature. Spatial discretization is performed by the use of finite elements in the model of Lee et al. [4], whilst the finite differences or volume discretization are used by other models [5—9]. The time-dependent terms of the equations are introduced explicitly in the model of Moussiopoulos et al. [6]. However, this formulation restricts the maximum-allowed integration time step to few seconds, and moreover, the higher the grid resolution, the shorter the time step which is required [10]. A usual way of describing a complex terrain is by the use of non-orthogonal co-ordinate transformation techniques [11,12], which however increase the complexity of the equations and in some cases affect the numerical stability. Alternatively, Bartzis et al. [5] use the simple Cartesian co-ordinate system with parallel treatment of the ground as a porous medium. The aim of this paper is to present and evaluate a new prognostic numerical model capable of predicting the wind and pollutant concentration field at both the mesoscale and microscale level over terrain of any degree of complexity. Special care has been given during the model development, towards a relatively simple but nonetheless accurate computational procedure. For this purpose, the collocated grid approach [13] is employed instead of the staggered grid arrangement, that has been used in all previous models; also, the equations are solved in Cartesian co-ordinates, while an effective method of partially-covered cells is constructed to allow for an accurate representation of the complex topography. For the same reason, so as to obtain details of the wind field over certain areas, the model can employ a local grid refinement technique. All the previous features have been embodied in a innovative way, and it must be pointed that the final source code is smaller in length and faster than the conventional TEACH code computational procedure [14], on which it is based. The greater area of Athens was chosen as a test case, since it combines all the features characterizing a complex terrain: The wind field over the Athens basin is strongly influenced by the presence of surrounding mountains, hills and large water basins (Fig. 1). The local circulation is mainly driven by the sea breeze, which appears frequently, and during the summer months occurs on nearly 50% of all days [15]. The present code is used to simulate the wind field and inert pollutant concentrations for the day of 25 May 1991, for which some measurements are available and quasi-constant geostrophic wind is reported. Following a brief description of the numerical model, the predictions of the observed diurnal variation of the flow field and its effect on the pollutant dispersion are discussed.
2. Mathematical formulation The final form of the conservation equations describing the time-dependent, mesoscale wind flow are obtained following simplification and averaging techniques
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
287
Fig. 1. The modelling domain of the greater area of Athens. Residential areas are stippled, industrial areas are solid. Measuring stations: 1 — Philadelphia (altitude 150 m), 2 — Tatoi (210 m), 3 — Helliniko (20 m), 4 — Piraeus (80 m), 5 — Eleysina (30 m), 6 — Peristeri (80 m), 7 — Marousi (80 m), 8 — Geoponiki (60 m), 9 — Nea Smirni (50 m), 10 — Patision (100 m), 11 — Athinas (100 m).
presented in Ref. [11]. The solution provides the “layer-domain-averaged” values of a dependent variable U: U"U #UI , (1) 0 where U is assumed to represent the synoptic scale atmospheric conditions at a given 0 altitude, and UI are the mesoscale variations from this larger scale [11]. The conservation equations are expressed below, in Cartesian tensorial form. The anelastic form of mass continuity equation is adopted [16]: L (o u )"0, j"1,2,3. Lx 0 j j Momentum equations:
CA
(2)
B
D
Lu D(o u ) LpJ L Lu hI i # j !ou@ u@ #o g #f · (d u !d u ) 0 i "! # k i j 0 i 1 i1 2 i2 1 Lx Lx Dt Lx Lx h j i i j 0 !f d u #f d u , (3) 2 i1 3 2 i3 1 where f and f are the coriolis parameters: f "2X sin /, f "2X cos /; u is the 1 2 1 2 i velocity component at the i-direction, p is the pressure, o the density, k the dynamic
288
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
viscosity, h the potential temperature, g the gravity acceleration, d the Kronecker’s ij delta and D/Dt represents the substantial derivative. X is the angular velocity of the earth’s rotation and / the latitude. The tilde denotes fluctuating components. The potential temperature equation is written as follows — a similar equation describes the concentration field:
A
B
D(o h) L k Lh 0 " !ou@ h@ #Q, (4) j Dt Lx p Lx j T j where Q represents the source-sink term (surface heat flux for temperature, and mass flow rate of pollutant emissions for concentration equation, respectively). The potential temperature is solved instead of the thermodynamic one ¹, because it is a conservable quantity and its integration is independent of the path. A closed set of the above differential equations is obtained by the incorporation of the well-known turbulence model k—e [17] using the constant values of Table 1. The kinematic heat flux is modelled as follows: k Lh !ou@ h@" t . (5) j p Lx T j The transport equations for the turbulent kinetic energy k and its dissipation rate e take the form
A B A B
D(o k) L k Lk t 0 " #P#G!o e, (6) 0 Dt Lx p Lx j k j L k Le e e2 D(o e) t 0 " #C [P#G (1!C )]!C o , (7) 1k 3 2 0k Dt Lx p Lx j e j where P is the turbulent shear stress production and the term G represents the buoyant contribution to the production term of turbulent kinetic energy [18]: k Lh t . (8) i ¹ p Lx 0 T i In stable stratification G becomes a sink term, reducing the turbulence intensity, whereas in unstable stratification the buoyancy enhances turbulence, since G is positive. The parameter C controls the buoyancy effect on the e production. Rodi [18] used 3 the standard value of 0.8 for heated surface jets. For sea-breeze simulation, Kitada and G"!g
Table 1 Model parameters C k
C 1
C 2
p k
p e
p T
0.09
1.44
1.92
1
1, 3
0.9
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
289
Takagi [19] suggest two different values of C according to local thermal stratifica3 tion: 0 for an unstable and 1 for a stable layer. These values were adopted in the present work. The values of all the other constant coefficients involved in the above equations are summarized in Table 1.
3. Computational details The numerical code is used to predict the wind field and inert pollutant concentration variation during the day of 25 May 1991, for which measurements are available and a quasi-constant geostrophic wind was present over the greater area of the city of Athens. The 70]70 km modelling domain is shown in Fig. 1, taken from Ref. [5]. The computational space is divided in orthogonal cells by the use of a 39]39]30 Cartesian grid, with uniform horizontal resolution of 2 km. This relatively coarse resolution cannot secure the grid independence of the results; however, it was selected in order to fit the corresponding resolution of the available topographic and emissions data. The use of finer grids over the entire domain would result in an enormous increase of computational demands and for this reason the development and incorporation of local grid refinement techniques is planned for the future. The vertical grid-line spacing equals 30 m from zero to 400 m altitude, and then is gradually expanded up to 6 km. 3.1. The iterative solution The governing equations are integrated over finite volumes surrounding the grid nodes, and then discretized with the hybrid scheme approach [20]. All variables are stored at the centre of the control volumes, following the collocated grid arrangement [13]. After discretization, all equations are cast into the same general form, and the system is solved iteratively using the SIMPLE algorithm [21]. The discretization procedure is implicit in time; this formulation provides stability in the solution procedure and allows large integration time steps (several hours), independently of the grid resolution. A time step of 1 h is used in the present calculations, and thus 24 “runs” are performed to reproduce the flow field variation during the examined day. The resulting field of each “run” forms the initial conditions for the next one. The hourly interval was found adequate to produce reliable results: the use of a 3 h time step gave similar results under constant (stable) atmospheric stratification. Moreover, since the available emission data are averaged over hourly intervals, the use of smaller time steps would not present any advantage. Due to the strong stratification in the examined buoyant flows, some stability problems arose, which were handled by using small values of the pressure underelaxation factor (0.1) and a smooth introduction of the buoyancy term included in the vertical momentum equation.
290
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
Each “run” requires on the average 1500 iterations to converge, which corresponds to about 5 h execution time in a 586 Pentium PC (100 MHz), and needs about 10 Mbytes CPU memory. 3.2. Inlet conditions At the inlet boundaries, the vertical distribution of the wind velocity and turbulent properties must be prescribed. In the surface layer this was done using the relationships of Ref. [11], in which a stability function t is introduced into the logarithmic law: u"u /i[1n(z/z )!W(z/¸)], * 0 h u3 ¸"! w * , gw@h@i
(9) (10)
where h is the surface potential temperature, z the surface roughness length, u the 8 0 * friction velocity, g the acceleration due to gravity and i the Von Karman constant ("0.4187). The surface heat flux w@h@ was calculated as described in Section 3.3 below. In the transition layer the wind velocity and direction was obtained from the Ekman equations [11]. The heights of the surface and the planetary layer were taken from Refs. [11,22], respectively, as function of the stratification conditions and the time of the day. The vertical distribution of the turbulence parameters (k, e) throughout the whole planetary boundary layer was calculated using the expressions given in Ref. [23], depending on the prevailing stability conditions. Above the planetary layer and at the upper boundary of the computational domain (6 km), the vertical velocity component is assumed zero, while the geostrophic northwest wind of 3.5 m s~1 is imposed for the horizontal component. Based on the observations, the synoptic field is kept constant during the day. Entrainment conditions are applied to all the other free boundaries of the domain. This is a simple numerical technique to account for situations when the developed pressure field inside the computational domain causes inflow conditions at a free boundary. According to this technique, the amount of the entering mass and the corresponding inlet velocities are calculated from a mass balance over the faces of each cell adjacent to a free boundary. The entering flow mass and momentum as well as the rest ambient fluid characteristics (temperature, turbulence quantities, etc.) are then included as extra source terms into the corresponding equations. The inlet temperature profile at a given time of the day is estimated via interpolation between the two available radiosonde measurements at 02:00 and 14:00 LST. 3.3. Ground boundary conditions The inability of the Cartesian grid to follow the exact surface configuration is eliminated with the development of a partially-covered cell technique. According to
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
291
Fig. 2. Terrain simulation in the computational grid.
this, a computational cell that is partly under the ground surface, is solved after multiplying its areas and volume by proper geometrical coefficients (y and y ) i 7 that represent their free (above the ground) portion. For example, in the simplified two-dimensional sketch of Fig. 2, the east surface of the (I,J) cell should be multiplied by the factor y "dy/D½, while the south surface by the factor y " E S dx/Dx. Also, the volume of the cell is multiplied by a factor y , equal to the free volume 7 portion. Thus, y"1 denotes a completely free and y"0 a totally blocked cell face or volume. The significant improvement of the ground boundary description achieved by the above technique is demonstrated in Fig. 2: It is obvious that the standard Cartesian step approximation may lead to substantial deviations from the actual boundary, which in contrast, can be simulated quite accurately using the new technique. Using these new coefficients, the general form of the discretized equations becomes (A !c S )U "+ c A U #c S , 1 7 1 P i i i 7 U i A "+ c A , i"East, West, North, South, Up, and Down face, (11) P i i i where A are the coefficients linking a dependent variable of a cell to its neighbours, on i the adjacent grid volumes. All the necessary geometrical coefficients are calculated and stored once at the beginning of the iterative procedure, and only for the cells adjacent to the ground. This practice reduces the additional storage requirements, without increasing the total execution time. Also, the wall function boundary conditions can be applied in a similar manner to all these cells, thus, considerably simplifying the algorithm.
292
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
Moreover, the actual slope and the exact area of the bounding surfaces are both calculated accurately. The surface heat flux is calculated from an energy balance on the ground surface using algebraic expressions to obtain the solar radiation as function of the incident radiation angle, which in turn depends mainly on the time of the day, as also on several other parameters (latitude, longitude, day of the year, cloudiness, surface temperature [24]). Typical values for the various constants involved were taken [8,25]. A typical daily variation (based on the measured data) was assumed for the surface temperature to enhance stability, whereas the actual slope of the ground at each location was not considered. Therefore, the calculated heat flux is uniform over the ground surface at any specific time of the day. A uniform ground roughness was also assumed for the land surface, whereas the sea surface was considered adiabatic [8]. 3.4. Code starting The time-dependent solution of the differential equations requires the knowledge of an initial flow field over the entire calculation domain. This field is obtained at 00:00 by the following practice: The code starts at 21:00 of 24 May 91, namely 3 h earlier, assuming a neutral atmosphere, same synoptic wind and constant temperature at levels parallel to the sea surface, which is the most convenient solution in case of unknown temperature distribution [26]. For the next 3 h the time-dependent solution is performed using a typical night-time heat flux at the ground. Thus, the resulting field approaches the real conditions prevailing at the beginning of the new day, and can be used as the initial state of the atmosphere at 00:00.
4. Results and discussion Meteorological measurements available for the examined day of 25 May 92 consist of wind and temperature vertical profiles at Helliniko airport station every 12 h, and surface (10 m) temperature and wind speeds at five stations every 3 h. The pollutant concentration pattern is available at seven different surface locations, at hourly intervals. The emission database used for the present investigation has been created during the APSIS exercise [27,28] for model evaluation purposes. These data include mean hourly values for various pollutants (CO, SO , NO, NO , HCHO and other hydro2 2 carbons), averaged in a 2]2 km2 space resolution over the entire terrain, and they are given in two layers: 0—30 m and 30—60 m above the ground. For this reason, buoyancy effects or plume rising phenomena (i.e. from stack emissions) were not considered. The numerical investigation concerns the dispersion characteristics of the relatively inert pollutants (CO and SO ), to avoid at this stage the modelling of the complex 2 photochemical kinetics. Also, the calculations were restricted to one day in order to check the reliability of the results and to spot any numerically imposed limitations affecting the validity of the predictions, before trying to assess the overall performance by applying the algorithm to different days or situations.
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
293
In order to make clear the various topographic and thermal influences, the general characteristics of the diurnal variation of the flow field are first demonstrated, while the comparison with observations follows. 4.1. Diurnal flow field variation The predicted diurnal variation of the flow field and pollutant concentration on the ground surface is illustrated in Figs. 3—10. Calculations are made for every hour of the day, but only the most representative results are reported and discussed here. The atmosphere is clearly stable during the nighttime (Fig. 3). The cooling of the air near the mountain surfaces generates downward breezes. This phenomenon is evident in the predictions of Fig. 4a. Air masses from the surrounding mountains of Parnis, Pendeli and Hymettus are accumulated in a main land breeze directed from the interior of the basin offshore. This northeast breeze protects the basin against the SO 2 pollution which is emitted mainly at the industrial zones of Eleysina and Aspropirgos (Fig. 4b). The CO concentration is also kept low during the night (not exceeding 2 mg m~3, Fig. 4c), since car traffic is light. The wind field and pollution over the Athens basin remains with minor changes as described above, until the first daylight. At 08:00 the land breeze is still prevailing (Fig. 5a), but the pollutant concentration exhibits appreciable changes and both SO 2 and CO values are considerably increased, especially at the centre of Athens and Piraeus, where traffic is now heavy (Fig. 5b and Fig. 5c). SO and CO concentration 2 exceeds 100 lg m~3 and 5 mg m~3 respectively, over these areas, and it is noticeable that the pollutants are transferred far out to the sea, thus substantial surface concentrations are predicted up to the nearby Aegina island. After sunrise, the temperature field changes radically due to the rapid heating of the land surface (Fig. 3). The stable stratification is partially destroyed, and transition to
Fig. 3. Predicted diurnal variation of the temperature profile at the Helliniko station.
294
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
Fig. 4. Predicted surface flow at 04:00 LST: (a) horizontal velocity vectors (scale: 1 mm"2 m/s); (b) SO 2 concentration (lg m~3); (c) CO concentration (mg m~3).
an unstable atmosphere takes place. As a consequence, the land breeze ceases to exist and the synoptic northwest wind field initially dominates (Fig. 6a). At 10:00, the air stagnation above the basin causes the CO and SO concentrations to reach their 2 maximum values, which are above 9 mg m~3 and 200 mg m~3, respectively, over the centre of Athens and Piraeus, while the surface distribution patterns remain similar to these at 08:00 (Fig. 5b and Fig. 5c). The continuation of land surface heating leads soon to the establishment of upward air streams along the mountainsides. These streams, in conjunction with the seabreeze formation, due to the relatively slower rate of temperature rising over the sea surface, result in the development of a completely different aerodynamic field over the basin at 14:00, as can be seen in Fig. 7a. The most interesting characteristic of this complicated field is the formation of a distinct line extending along the peninsula, where the fronts of two opposite sea-breeze circulation systems meet. The existence of the northern part of this line, between the mountains of Pendeli and Hymettus, has been observed as a common situation at this time of the day, and it seems to block the physical channel for the pollutants to be transported away from the basin. This
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
295
Fig. 5. Predicted surface flow at 08:00 LST: (a) horizontal velocity vectors (scale: 1 mm"2 m/s); (b) SO 2 concentration (lg m~3); (c) CO concentration (mg m~3).
Fig. 6. Predicted surface horizontal velocity field at 10:00 LST (scale: 1 mm"2 m/s).
296
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
Fig. 7. Predicted surface flow at 14:00 LST: (a) horizontal velocity vectors (scale: 1 mm"2 m/s); (b) SO 2 concentration (lg m~3); (c) CO concentration (mg m~3).
Fig. 8. Predicted surface horizontal velocity field at 18:00 LST (scale: 1 mm"2 m/s).
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
297
Fig. 9. Predicted surface horizontal velocity field at 24:00 LST (scale: 1 mm"2 m/s).
Fig. 10. Predicted surface pollutant concentration levels at 22:00 LST: (a) SO in lg m~3; (b) CO in 2 mg m~3.
blockage may be responsible for the dramatic increase of pollutant concentrations that occur when strong temperature inversion conditions exist over the basin and the upward streams are weak. This is not the case in the examined situation, and therefore, the predicted values of SO and CO shown in Fig. 7b and Fig. 7c are both 2 decreased due to unstable conditions. The polluted area is now split in two parts, with the one over the sea diluting due to diffusion, since no more pollutants are transferred there. The relation between land and sea surface temperature is reversed in the evening, due to the faster cooling of the ground (Fig. 3). A weak sea breeze still exists at 18:00 (Fig. 8), but the synoptic wind field dominates. A down-slope breeze has already been established at 22:00, while the aerodynamic field at 24:00 (Fig. 9) is almost the same as
298
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
the one at the beginning of the day (Fig. 4). This is an expected result, since the synoptic field is kept constant during the whole day, and yet confirms the accurate operation of the entrainment boundary conditions used in the computer code. The wind direction is uncertain over the basin during the re-transition from unstable to stable stratification, while traffic is still dense. Consequently, the pollutant concentrations start to increase again, and obtain high maximum values at about 22:00 (Fig. 10a and Fig. 10b), approaching the morning maxima. 4.2. Comparison with available measurements The daytime temperature vertical profile is compared with measurements at Helliniko in Fig. 11a. The agreement is satisfactory, except in the region between 300
Fig. 11. Comparison between measurements and predictions at Helliniko station, at 14:00 LST: (a) Temperature profile; (b) west wind component; (c) south wind component.
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
299
Fig. 12. Diurnal variation of surface temperature at five measuring stations.
and 600 m above the ground, where an inversion layer is indicated but not reproduced by the numerical model. A possible reason for this discrepancy is that the temperature profile for a sea-shore location as Helliniko is influenced by the sea surface temperature spatial variability [27], that is not discernible by the horizontal resolution of 2 km used. Unfortunately, no additional measurements of the vertical temperature profile are available at more surface stations in order to estimate the extent of this inversion and the influence of the above discrepancy on the air quality predictions. Similar are the predictions of Varvayanni et al. [28] that are also plotted in Fig. 11a, whereas an inversion layer like the measured one can be observed in the results of Thunis et al. [27], but temperatures are strongly underpredicted.
300
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
Fig. 13. Diurnal variation of west wind component at five measuring stations.
The calculated versus the measured profiles of mean wind speed at the same station and 14:00 LST, are shown in Fig. 11b and Fig. 11c. The results are considered acceptable, mainly in their qualitative behaviour concerning the wind direction. The predicted surface temperatures are in reasonable agreement with the measured ones at the low-altitude stations (Helliniko, Piraeus and Eleysina), as can be seen in Fig. 12. The discrepancies can be attributed to the approximate estimation of the thermal air—ground interaction (constant surface roughness length and a typical diurnal variation of surface heat flux), adopted in the model. More pronounced are the disagreements at the stations of Philadelphia and Tatoi at the foot of the mountain
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
301
Fig. 14. Diurnal variation of south wind component at five measuring stations.
where the dynamic effects are strong, indicating a difference between calculated and actual time variation of the convergence line location. Unfortunately, no measurements are available on the mountain slopes, thus making difficult the evaluation of the calculations concerning the up-slope and down-slope streams developed there. The percentage error, (¹ -¹ )/*¹ , ranges from almost zero to less than 45% in 13%$. .%!4. .%!4. all the above cases, except at 11:00 at Helliniko, where it exceeds 100%. The surface wind field is reported in Figs. 13 and 14 at the same five meteorological stations. The calculated velocity components agree fairly well with the observed ones. The model reproduces the low nighttime wind speeds, while non-systematic discrepancies can be observed at various stations and times of the day because the
302
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
Fig. 15. Diurnal variation of SO concentration level at seven measuring stations. 2
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
Fig. 16. Diurnal variation of CO concentration level at seven measuring stations.
303
304
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
2]2 km2 grid used cannot simulate the smaller-scale surface characteristics which influence the local measurements. Some more systematic deviations occurring in regions near the boundaries of the calculation domain (i.e. Eleysina, Fig. 14e), can be assigned to inaccuracies in the inlet conditions used. It must be also noted that the velocities are calculated at the centre of the first solvable grid volume above the ground, where the exact altitude cannot be calculated from the available data for the four neighbouring corners of the grid volume, and may not correspond well to the measured position (10 m). The most decisive evaluation point concerns the ability of the model to predict the diurnal pollutant concentration at ground level. The numerical results are plotted in Figs. 15 and 16 for SO and CO, respectively, and show a reasonable agreement with 2 the observed values, at least in the shape of the distribution and the timing of the concentration maxima. As for the quantitative differences, it must be noticed that the numerical results represent averaged values inside relatively large computational cells (2 km]2 km]30 m), whereas the measurements are taken at certain fixed points. It is obvious, for example, that the significant underprediction of pollutant concentrations at Patision station, is due to this station located just over one of the busiest streets in Athens. In general, all the diurnal concentration curves exhibit two maxima (Figs. 15 and 16): the first appears between 08:00 and 10:00 LST, because traffic is enhanced at that period, whilst stagnant flow is created over the basin, due to the change from a landbreeze to sea-breeze (Fig. 6). The unstable atmosphere which is developed and established afterwards, contributes to a significant reduction of the pollution at all measuring stations. This reduction continues till about 18:00. From that time, the wind above the basin starts slowing down due to the change over from sea-breeze to land-breeze with the appearance of the first stable layer aloft, allowing thus again the accumulation of the produced pollutants. A second peak is therefore formed between 20:00 and 22:00 LST, but with a smaller intensity in most cases. Some corresponding results, obtained by the explicit solution of Ref. [6], are also included in Fig. 16, and exhibit analogous variability and similar average degree of agreement with the measurements, although the integration time-step was very small (several seconds). Therefore, the hourly integration used for the present calculations can be considered adequate to reproduce the diurnal variation of the pollutant concentrations.
5. Conclusions A new mesoscale numerical model has been used to predict the wind field and inert pollutant concentrations over the complex terrain of the greater area of Athens, during the day of 25 May 91. A technique for the solution of partially-covered cells, in conjunction with the collocated Cartesian grid arrangement, achieved a close representation of the complex topography. The results reveal the great influence of the surrounding mountains and sea on the flow field characteristics. Thermally driven down-slope and up-slope air streams are
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
305
developed during the whole day along the mountainsides, while two strong sea breeze systems are established during the day-hours from the opposite sea-shores of the peninsula. The fronts meet and block the natural ventilation of the basin. The resulting aerodynamic field is extremely complex and dominates over the weak synoptic wind, thus controlling the dispersion of pollutants. Pollutant concentrations exhibit at all stations two maxima during the day: the first and higher around 10:00 and the second at about 22:00, both due to a combination of increased traffic and the change over from sea breeze to land breeze and vice versa, associated with the change in the near-ground atmospheric stability. Considering this high degree of complexity, the comparison of the model predictions with available observations is relatively satisfactory and the general diurnal variation follows in most cases the corresponding measurements. The consideration of smaller-scale ground surface characteristics (local roughness length, emissivity, surface slope, etc.), along with the planned incorporation of local grid refining techniques, higher-order differencing schemes and a photochemical reaction model, should provide more detailed and accurate information in order to achieve a better understanding of the air pollution problem. Additional surface station measurements are however necessary, especially in the northern mountainous region of the area examined, in order to evaluate the model performance at the microscale level.
References [1] P.C. Patnaik, A preliminary users guide for the SIGMET Mesoscale Meteorology Code, SAI Report SAI-78-768-LJ, RLO/2440-77-9, San Diego, CA, 1979. [2] N.L. Seaman, R.A. Anthes, A mesoscale semi-implicit numerical model, Quart. J. Roy. Meteorol. Soc. 107 (1981) 167—190. [3] M. Tjernstrom, Numerical simulations of stratiform boundary layer clouds on the Meso-c-Scale, Part I: the influence of terrain height differences, Boundary Layer Meteorol. 44 (1988) 33. [4] R.L. Lee, J.H. Leone, D.M. Gresho, S.T. Chan, A finite element model for environmental problems involving complex terrain, in: Proc. 4th Int. Conf. of Finite Elements in Flow Problems, Tokyo, Japan, 1982. [5] J.G. Bartzis, A.G. Venetsanos, M. Varvayanni, N. Catsaros, A. Megaritou, ADREA-I: A threedimensional transient transport code for complex terrain and other applications, Nucl. Technol. 94 (1991) 135—148. [6] N. Moussiopoulos, Th. Flassak, P. Sahm, D. Berlowitz, Simulations of the wind field in Athens with the nonhydrostatic mesoscale model MEMO, Environ. Software 8 (1993). [7] R. Kunz, N. Moussiopoulos, Simulation of the wind field in Athens using refined boundary conditions, Atmos. Environ. (1995), in press. [8] D. Tryfonopoulos, G. Bergeles, Temperature fields and air pollution build-up over the Athens basin, Environ. Software 9 (1994) 269—283. [9] J. Glekas, G. Bergeles, Dispersion under neutral atmospheric conditions, Int. J. Numer. Meth. Fluids 19 (1994) 237—257. [10] T. Flassak, M. Wortmann, N. Moussiopoulos, Influence of horizontal grid resolution on simulated wind and pollution concentration field in the greater Athens area, in: N. Moussiopoulos, G. Kaiser (Ed.), Proc. Seminar on Monitoring and Modelling in the Mesoscale, Thessaloniki, 1991, pp. 155—164. [11] R.A. Pielke, Mesoscale Meteorological Modeling, Academic Press, London, 1984.
306
J.S. Anagnostopoulos, G. Bergeles/J. Wind Eng. Ind. Aerodyn. 73 (1998) 285–306
[12] J. Glekas J, G. Bergeles, N. Athanasiadis, Numerical solution of the transport equation for passive contaminants in three-dimensional complex terrain, Int. J. Numer. Meth. Fluids 7 (1987) 319—335. [13] C. Rhie, W. Chow, Numerical study of the turbulent flow past an airfoil with trailing edge separation, A.I.A.A. J. 21 (1983) 1525—1532. [14] A.D. Gosman, F.J.K. Ideriah, TEACH-2E: A general computer program for two-dimensional turbulent recirculating flows, Report, Department of Mechanical Engineering Imperial College, 1976. [15] D.P. Lalas, D.N. Asimakopoulos, D.G. Deligiorgi, C.G. Helmis, Sea-breeze circulation and photochemical pollution in Athens, Greece, Atmos. Environ. 17 (1983) 1621—1632. [16] F.B. Lipps, R.S. Hemler, A scale analysis of deep moist convection and some related numerical calculations, J. Atmos. Sci. 39 (1982) 2192—2210. [17] B.E. Launder, D.B. Spalding, The numerical computation of turbulent flows, Comput. Meth. Appl. Mech. Eng. 3 (1974) 269—289. [18] W. Rodi, Turbulence Models and their Application in Hydraulics. A State of the Art Review, Int. Assoc. Hydraulic Research Publ., 1980. [19] T. Kitada, H. Takagi, Some remarks on the k—e turbulence model applied to sea breeze simulation: Buoyancy effect of the e-equation and horizontal eddy diffusivity, in: N. Moussiopoulos, G. Kaiser (Ed.), Proc. Seminar on Monitoring and Modelling in the Mesoscale, Thessaloniki, 1991 pp. 135—146. [20] A.D. Gosman, W.M. Pun, A.K. Ruchal, D.B. Spalding, Heat and Mass Transfer in Recirculating Flows, Academic Press, London, 1969. [21] S.V. Patankar, D.B. Spalding, A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows, Int. J. Heat Mass Transfer 15 (1972) 1787—1792. [22] F.B. Smith, Atmospheric turbulent transport, in: S. Sandroni (Ed.), Course on Regional and Longrange Transport of Air Pollution, Ispra, Italy, Elsevier, Amsterdam, 1986, pp. 43—69. [23] S.R. Hanna, G.A. Briggs, R.P. Hosker, Handbook on Atmospheric Diffusion, Technical Information Center, U.S. Department of Energy, 1982. [24] A.A.M. Holtslag, A.P. Van Ulden, A simple scheme for daytime estimates of the surface fluxes from routine weather data, J. Clim. Appl. Meteorol. 22 (1983) 517—529. [25] D. Tryfonopoulos, Numerical solution of atmospheric flows over complex terrain, Ph.D. Thesis, National Technical University of Athens, 1992. [26] R.T. McNider, R.A. Pielke, Diurnal boundary layer development over sloping terrain, J. Atmos. Sci. 38 (1981) 2198—2212. [27] P. Thunis, P. Grossi, G. Graziani, H. Gallee, B. Moyaux, G. Schayes, Preliminary simulation of the flow field over the Attic Peninsula, Environ. Software 8 (1993). [28] M. Varvayanni, J.G. Bartzis, C.G. Helmis, D.N. Asimakopoulos, Simulation of the sea breeze under opposing synoptic conditions, Environ. Software 8 (1993).