Numerical Modeling for the Accidental Dispersion of Hazardous Air Pollutants in the Urban Metropolitan Area

A numerical simulation system is developed to predict the dispersion of hazardous air pollutants (HAPs) over a populated city due to accidental release. Ulsan, as one of the megacities in Korea, is chosen as an ideal testbed for the simulation, as it is located in complex terrain and hosts a national industrial complex on the outskirts of the city. The system is based on the California puff model (CALPUFF) for simulating a HAP’s dispersion, in which the three-dimensional atmospheric circulation derived from the observed weather station data is specified at a fine horizontal resolution of 200 m. A test scenario is developed for the accidental release of benzene during the daytime and nighttime, respectively, by a fictitious explosion of a storage container, and the injection amount is determined arbitrarily yet comparable to those in the past accidents. In attempting a quantitative assessment and zoning the level of potential risk over the impacted area, multiple simulations have been conducted each day with different hourly varying meteorological conditions in August. The dispersion characteristics of the air pollutant depend largely on the local wind patterns that vary substantially from day to day. Nevertheless, the composite analysis sufficiently identifies the impacted area by the HAP’s dispersion due to the local prevailing wind such as the land–sea breeze circulation. An immediate hazardous area is determined based on the vulnerability map constructed by zoning the level of risk determined by the spatial distribution of the HAPs’ concentration and the harmfulness standard to the human body.


Introduction
Modern industry uses various kinds of chemicals, which can be not only the source of air pollutants gradually degrading ambient air quality but also often cause hazardous accidents. The dense populations in urban cities tend to increase vulnerability, even with narrow pollutant dispersion from nearby accidents [1]. There were several notorious accidents, including the case in Bhopal, India, in 1984, which resulted in more than 200,000 casualty losses. Another was the natural gas leak in Porter Ranch, California, in 2015-2016, which continued over four months and led to the relocation of up to 5000 residents.
In South Korea, about 42 accidents occur on average per year based on the national statistics reported by the Korea Research Institute of Chemical Technology for 2003-2018 [2]. The accident that happened in Gumi in 2012 was the worst, where the hydrogen fluoride (HF) gas leaked from a chemical-producing factory, causing 23 casualty losses. The devastating effects included killing crops

Atmospheric Dispersion Model
The model used in this study is the California puff model (CALPUFF) [20], which has been used widely for solving various atmospheric dispersion problems in previous studies (e.g., [9,21,22]). In [23], showed that the model was able to simulate the dispersion of H 2 S differently according to the different meteorological conditions specified for flat land, land-sea, and mountain-valley geographies, and [24] modified the model to reflect the impacts of heterogeneous land surface conditions on the dispersion Atmosphere 2020, 11, 477 3 of 18 and local concentration of air pollutants. The models are based on the non-steady-state Gaussian puff model, which represents a dispersed plume consisting of several puffs (or lumps) of pollutants. Each puff is transported by the specified three-dimensional wind fields, and the size of the puff changes over time by the internal diffusion process with a Gaussian distribution. The concentration equations are formulated by [25] as below: where C is the ground-level concentration, Q is the pollutant mass in the puff, σ x and σ y are the standard deviations of the Gaussian distribution in the along-, cross-, and vertical-wind directions, respectively, and d a and d c are the distances from the puff center to the receptor in the along-and cross-wind directions. In Equation (1), g determines the vertical dispersion as where σ z is the standard deviation of the vertical Gaussian distribution, H e is the effective height above the ground off the puff center, and h is the mixed layer height. The three-dimensional wind fields are derived from the discrete meteorological data through statistical interpolation and extrapolation from the California Meteorological Model (CALMET), which is a diagnostic model for processing meteorological data required for CALPUFF. CALMET is available for processing various meteorological data from observation stations or model simulation outputs to the gridded three-dimensional data. More detailed information is provided in the user's guide for CALMET [26] and CALPUFF [25], respectively.

Experiment Settings
This study selects benzene as a representative chemical among the HAPs, which is the first-class carcinogenic chemical broadly used in our local study area. Benzene is a stable and non-polar chemical with low reactivity. The biggest reaction is the substitution of hydrogen, but it requires catalysts and a sufficient amount of heat. Therefore, this study disregards a secondary chemical reaction process under the normal outdoor atmospheric conditions and treats it as a pure dispersive gas by meteorological conditions. Table 1 shows the chemical properties of benzene specified in this study, following the Korean Ministry of Environment [27].  Figure 1 shows the modeling domain, in which the urban area with the industrial complex is surrounded by high mountain ranges and ocean. The surface elevation and the land cover type were specified based on the ASTER global digital elevation model version 2 (ASTGTM2) and Global Land Cover Characterization (GLCC), respectively. To produce the meteorological fields, CALMET processed the surface data from the local Automated Surface Observing System (ASOS) and unmanned automated weather stations (AWS), and the rawinsonde data in Pohang, which is the nearest available upper-level observation site. The source point in the study domain was chosen as one of the industrial businesses located at 35.504 N, 129.353 E. This region is identified as the largest benzene emission area within the industrial complex, based on the annual emission and waste amount total reported officially to the government. Although the worst scenario of an accident is the entire leakage of all contained gas in the storage, this information is unavailable in public in Korea. The potential storage amount was estimated based on the waste transport, emission, and the displacement for each business industry reported to the government database of the Pollutant Release and Transfer Register (PRTR) [28]. Based on the record, the total waste of benzene transported in Ulsan in 2013 was 416,204 kg. The annual emission of benzene from the targeted point source was also recorded as 26,097 kg yr −1 , which accounts for 39% of the total emission of 66,908 kg yr −1 in Ulsan. It is assumed that the waste amount is also proportional to the emission amount and this study applied 39% of the 416,204 kg yr −1 total in Ulsan, which is 162,319.56 kg yr −1 or 13,526.63 kg mon −1 . Given the monthly waste total, this study further assumes that the storage will be replenished once for a month, then this amount of waste can be assumed as the amount of monthly storage. Our scenario is 60% leakage of the monthly storage to the atmosphere within 12 h, which is estimated at 94 g s −1 , or 339 kg h −1 . For comparison, Table 2 shows the record of past accidents that occurred in Ulsan based on the Integrated Chemical Information System (ICIS) [29]. In Table 2, the emission rate is simply estimated based on the duration period and total amount of leakages from the information system, which shows lots of variation. By comparing these with the emission scenario in this study, the value of 339 kg h −1 is the one likely to occur. The source point in the study domain was chosen as one of the industrial businesses located at 35.504 N, 129.353 E. This region is identified as the largest benzene emission area within the industrial complex, based on the annual emission and waste amount total reported officially to the government. Although the worst scenario of an accident is the entire leakage of all contained gas in the storage, this information is unavailable in public in Korea. The potential storage amount was estimated based on the waste transport, emission, and the displacement for each business industry reported to the government database of the Pollutant Release and Transfer Register (PRTR) [28]. Based on the record, the total waste of benzene transported in Ulsan in 2013 was 416,204 kg. The annual emission of benzene from the targeted point source was also recorded as 26,097 kg yr −1 , which accounts for 39% of the total emission of 66,908 kg yr −1 in Ulsan. It is assumed that the waste amount is also proportional to the emission amount and this study applied 39% of the 416,204 kg yr −1 total in Ulsan, which is 162,319.56 kg yr −1 or 13,526.63 kg mon −1 . Given the monthly waste total, this study further assumes that the storage will be replenished once for a month, then this amount of waste can be assumed as the amount of monthly storage. Our scenario is 60% leakage of the monthly storage to the atmosphere within 12 h, which is estimated at 94 g s −1 , or 339 kg h −1 . For comparison, Table 2 shows the record of past accidents that occurred in Ulsan based on the Integrated Chemical Information System (ICIS) [29]. In Table 2, the emission rate is simply estimated based on the duration period and total amount of leakages from the information system, which shows lots of variation. By comparing these with the emission scenario in this study, the value of 339 kg h −1 is the one likely to occur.
The minimum detection range of benzene in the air is known as 0.64 µg m −3 [30]. The annual standard for benzene in Korea is 5 µg m −3 , and the US EPA defines the dangerous level of inhalation is 30 µg m −3 or higher according to the Integrated Risk Information System (IRIS) [31] and the Chemical Summary Report of the benzene [32] These values are used in the health risk assessment in Section 3.3. The numerical experiments were conducted for the domain for 1-31 August 2013 at a 200 m spatial resolution. The reason for choosing August was to consider the effect of the spread of HAPs due to the regional circulation specific to Ulsan. The temperature of Ulsan is highest this month during the year and the local daytime sea breeze wind develops at its peak phase. Therefore, the dispersion to the adjacent residential area is more likely to occur in summer. Although this study focuses on summer, it is noted that the modeling framework developed in this study can be generalized easily to other seasons such as winter. Although not shown here, much stronger northwesterly winds prevail in winter, which tends to slow down the sea breeze and weaken the transport of atmospheric pollutants toward the inland region. The experiments were carried out twice a day. Each experiment was run for 12 h, started from the identical time of the day at 14:00 local standard time (LST), with an identical emission scenario of benzene. Therefore, the experiments have 31 simulation outputs, all of which have different chemical distributions in time and space due to the difference in meteorological conditions only. Hourly outputs were then averaged for a month for each hour of the day to depict the time-averaged dispersion characteristics driven by the local observed meteorological conditions. Chaotic weather variability was averaged out and the diurnal change in local circulation driven by the land-sea contrast became evident.

Validations for Wind and the Test Cases
This study first validates the simulation of local wind by CALMET against the observations. Figure 2 compares the local wind variation, which is averaged each hour over August 2013. The observed local wind is westerly at noon and it turns into southeasterly at around 15:00 LST and to southerly until 20:00 LST, due to the development of the sea breeze in the study domain. By the weakening of the sea breeze and developing land breeze, the wind changes to northwesterly after 21:00 LST until 11:00 LST. In general, the observed and the simulated wind vectors are in good agreement in terms of direction and magnitude, which demonstrates that CALMET can represent the local observed wind variation realistically. When the daily mean wind components are removed and only the diurnal components are compared, the simulations also exhibit good agreement with the observations, although the overall magnitude is weaker than the observed. The CALMET diagnostic model adjusts the sparse wind observations onto the model grids at 200 m resolution, not only based on statistical interpolation using an inverse-distance method but also by considering the kinematic effects of the terrain, slope flows, blocking effects, and three-dimensional divergence minimization [22]. Therefore, the simulated wind is not exactly the same as the observed, although the difference is negligible. Note that the southerlies are relatively stronger than the northerlies, suggesting that the sea breeze is more pronounced than the land breeze in this region. August is characterized by strong diurnal wind variation in this coastal region, typically in summer. During winter, the observed winds are mostly northerlies with much suppressed diurnal wind variation (not shown). on statistical interpolation using an inverse-distance method but also by considering the kinematic effects of the terrain, slope flows, blocking effects, and three-dimensional divergence minimization [22]. Therefore, the simulated wind is not exactly the same as the observed, although the difference is negligible. Note that the southerlies are relatively stronger than the northerlies, suggesting that the sea breeze is more pronounced than the land breeze in this region. August is characterized by strong diurnal wind variation in this coastal region, typically in summer. During winter, the observed winds are mostly northerlies with much suppressed diurnal wind variation (not shown).  Figure 3 compares the sample cases of benzene dispersion in time during the daytime and the nighttime accidents. These cases are chosen as the typical cases in the summertime weather in Ulsan, characterized by the local land-sea breeze circulation. The chemical dispersion pattern depends largely on the specified wind condition. The initial direction of the plume and the level of surface concentration is quite different between the daytime emission-started and the nighttime emissionstarted cases. When the emission starts in the afternoon at 14:00 LST, the emission plume tends to spread to the northwest by the developing sea breeze in the one hour after. It tends to move further northwestward by the enhanced sea breeze in the three hours after. The plume also becomes wider in the cross-wind direction by the diffusive nature of Gaussian puffs. Six h later at 20:00 LST, the sea breeze slows down and the overall wind becomes calm. The surface concentration shows a circular pattern with a slight stretch in the zonal direction by background winds. At 12 h after at 02:00 LST, the wind becomes northerly, suggesting the development of the land breeze in the nighttime, and the plume is advected toward the ocean. The surface concentration in the downstream tends to increase a lot by the nighttime suppression of the atmospheric boundary layer. Compared with the southerly sea breeze, the land breeze is much weaker and this is consistent with the results shown in Figure 2.  Figure 3 compares the sample cases of benzene dispersion in time during the daytime and the nighttime accidents. These cases are chosen as the typical cases in the summertime weather in Ulsan, characterized by the local land-sea breeze circulation. The chemical dispersion pattern depends largely on the specified wind condition. The initial direction of the plume and the level of surface concentration is quite different between the daytime emission-started and the nighttime emission-started cases. When the emission starts in the afternoon at 14:00 LST, the emission plume tends to spread to the northwest by the developing sea breeze in the one hour after. It tends to move further northwestward by the enhanced sea breeze in the three hours after. The plume also becomes wider in the cross-wind direction by the diffusive nature of Gaussian puffs. Six h later at 20:00 LST, the sea breeze slows down and the overall wind becomes calm. The surface concentration shows a circular pattern with a slight stretch in the zonal direction by background winds. At 12 h after at 02:00 LST, the wind becomes northerly, suggesting the development of the land breeze in the nighttime, and the plume is advected toward the ocean. The surface concentration in the downstream tends to increase a lot by the nighttime suppression of the atmospheric boundary layer. Compared with the southerly sea breeze, the land breeze is much weaker and this is consistent with the results shown in Figure 2.
In the nighttime emission-started case, the initial plume moves eastward by weak background westerlies in the one hour after. The concentration tends to increase further in the downstream in the three hours after, with a more stretched plume in the along-and cross-wind directions. Six h later at 08:00 LST, the surface concentration tends to decrease, presumably due to the dilution by the development of the atmospheric boundary layer in the morning. The plume changes drastically in direction and it moves to the northwest by the strong sea breeze circulation in the afternoon at 14:00 LST, in the 12 h after emission.
Overall, CALMET represents the wind reversal by the local land-sea breeze circulation. Accordingly, the chemical dispersion pattern is realistic, reflecting the important role of local wind patterns. In the nighttime emission-started case, the initial plume moves eastward by weak background westerlies in the one hour after. The concentration tends to increase further in the downstream in the three hours after, with a more stretched plume in the along-and cross-wind directions. Six h later at 08:00 LST, the surface concentration tends to decrease, presumably due to the dilution by the development of the atmospheric boundary layer in the morning. The plume changes drastically in direction and it moves to the northwest by the strong sea breeze circulation in the afternoon at 14:00 LST, in the 12 h after emission.
Overall, CALMET represents the wind reversal by the local land-sea breeze circulation. Accordingly, the chemical dispersion pattern is realistic, reflecting the important role of local wind patterns.
This study further examines the mechanisms for a much higher benzene concentration in the nighttime compared with that in the daytime. Figure 4 shows the vertical cross-section of the benzene dispersion at the source point. The horizontal axis is aligned following the surface wind direction at the source point, whose direction is indicated separately for the daytime and the nighttime emissionstarted cases in Figure 1. As shown in Figure 3, the southeasterly sea breeze wind near the emission source transports the chemical to the northwest in the daytime emission-started case in Figure 4a. The dispersion pattern also shows a vertical expansion, thereby contributing to decreasing the concentration level near the surface. Note that the mixing height in the boundary layer for the daytime and nighttime is significantly different. Considering that the typical height of the atmospheric boundary layer is about 1-2 km [33], benzene tends to spread in the vertical direction due to the development of the atmospheric boundary layer during the daytime. In the meantime, surface heating during the daytime tends to induce a vertical upward motion, thereby driving a more enhanced vertical transport of the chemicals. The vertical motion is more accelerated at the increase in height due to less turbulence momentum mixing. On the other hand, strong subsidence develops on the opposite side, which is southeast to the emission source and toward the ocean. The wind pattern in Figure 4a shows the local sea breeze circulation developed in the daytime, with an upward motion over the land and downward motion over the ocean. It shows the sea breeze at the lower This study further examines the mechanisms for a much higher benzene concentration in the nighttime compared with that in the daytime. Figure 4 shows the vertical cross-section of the benzene dispersion at the source point. The horizontal axis is aligned following the surface wind direction at the source point, whose direction is indicated separately for the daytime and the nighttime emission-started cases in Figure 1. As shown in Figure 3, the southeasterly sea breeze wind near the emission source transports the chemical to the northwest in the daytime emission-started case in Figure 4a. The dispersion pattern also shows a vertical expansion, thereby contributing to decreasing the concentration level near the surface. Note that the mixing height in the boundary layer for the daytime and nighttime is significantly different. Considering that the typical height of the atmospheric boundary layer is about 1-2 km [33], benzene tends to spread in the vertical direction due to the development of the atmospheric boundary layer during the daytime. In the meantime, surface heating during the daytime tends to induce a vertical upward motion, thereby driving a more enhanced vertical transport of the chemicals. The vertical motion is more accelerated at the increase in height due to less turbulence momentum mixing. On the other hand, strong subsidence develops on the opposite side, which is southeast to the emission source and toward the ocean. The wind pattern in Figure 4a shows the local sea breeze circulation developed in the daytime, with an upward motion over the land and downward motion over the ocean. It shows the sea breeze at the lower levels up to a few hundred meters, although there is uncertainty in the simulation for the vertical depth of the sea breeze circulation.
The dispersion pattern in the nighttime is quite contrasting to the case in the daytime. Due to the westerlies near the emission source in the nighttime, as shown in Figure 3, benzene spreads to the east in the figure. The surface cooling and nocturnal stable boundary layer tend to suppress the vertical dispersion of benzene. Further, the vertical motion in the lower atmosphere maintains weak subsidence, which tends to further inhibit the vertical transport of benzene. As a result, the dispersion sticks to the ground and the benzene concentration tends to be larger near the ground-level in the nighttime compared with that in the daytime. This explains why the impacted area by the high level of the benzene concentration becomes larger in the nighttime, as shown in Figure 3.  The dispersion pattern in the nighttime is quite contrasting to the case in the daytime. Due to the westerlies near the emission source in the nighttime, as shown in Figure 3, benzene spreads to the east in the figure. The surface cooling and nocturnal stable boundary layer tend to suppress the vertical dispersion of benzene. Further, the vertical motion in the lower atmosphere maintains weak subsidence, which tends to further inhibit the vertical transport of benzene. As a result, the dispersion sticks to the ground and the benzene concentration tends to be larger near the ground-level in the nighttime compared with that in the daytime. This explains why the impacted area by the high level of the benzene concentration becomes larger in the nighttime, as shown in Figure 3.

Time-Averaged Dispersion Pattern
The time-averaged dispersion patterns obtained from the numerical experiments for the entire 31 days in August represent the time evolution of the vulnerable area at the chemical accident in the summertime. For the daytime emission-started case (Figure 5), the wind variation shows much resemblance to that in the sample case of 15 April (Figure 3), and it shows the signal of the land-sea breeze circulation driven by the differential heating between the land and the ocean in the diurnal timescale more clearly by eliminating random weather variability. The sea breeze develops from the south during the afternoon to the early evening and the off-shore wind develops in the nighttime from the west in the late evening and from the north after midnight. The model also represents the impact of surface friction exerted on the wind realistically by slowing down the wind speed over the elevated mountain area.  Figure 1 for the daytime and the nighttime emission-started cases, respectively. Wind (arrow) and the source point (red triangle) are indicated in each figure. Note that the vertical motion is exaggerated by multiplying the vertical wind component by 100.

Time-Averaged Dispersion Pattern
The time-averaged dispersion patterns obtained from the numerical experiments for the entire 31 days in August represent the time evolution of the vulnerable area at the chemical accident in the summertime. For the daytime emission-started case ( Figure 5), the wind variation shows much resemblance to that in the sample case of 15 April (Figure 3), and it shows the signal of the land-sea breeze circulation driven by the differential heating between the land and the ocean in the diurnal timescale more clearly by eliminating random weather variability. The sea breeze develops from the south during the afternoon to the early evening and the off-shore wind develops in the nighttime from the west in the late evening and from the north after midnight. The model also represents the impact of surface friction exerted on the wind realistically by slowing down the wind speed over the elevated mountain area. Figure 5 also shows the time evolution of the time-averaged benzene concentration during the daytime emission-started cases. Each figure uses the color scheme with just three concentration levels of 0.64, 5, and 30 µg m −3 , which are the minimum value of detection by the instrument, the annual standard of benzene concentration in Korea, and the dangerous level defined by US EPA, respectively. The averaged dispersion patterns show the hint of a few rapid dispersion cases affected by strong wind in specific dates, indicating that the one-month average is not sufficient enough to obtain generalized dispersion patterns. Nevertheless, the time evolution of the dispersion patterns still captures the basic features affected by the diurnal variation of the local circulation, particularly reflected in the higher concentration levels. Benzene dispersion shows a gradual expansion in the area over time of the continuous emission at the source point that started in the afternoon. Moreover, when the southerly winds prevail in the afternoon, benzene spreads to the north following the sea breeze. Then, it spreads southeastward by the weak land breeze from the late evening to the nighttime. Based on the simulations, the detectable range of benzene can be observed at the 5 km distance from the emission source in 1 h, 10 km in 3 h, and 15 km in 6 h. The benzene concentration level exceeds the annual Atmosphere 2020, 11, 477 9 of 18 standard value in Korea at the 5 km distance in 3 h and 10 km in 6 h, and 15 km in 12 h. The area under the dangerous benzene level stays within a 2 km radius from the emission source until 3 h after the emission. However, the area of dangerous level spreads quickly in the late evening and after midnight, reaching the 5 km distance after 6 h.
Atmosphere 2020, 10, x FOR PEER REVIEW 9 of 18   The wind variation during the nighttime emission-started cases ( Figure 6) also shows much resemblance to the sample case of 16 August (Figure 3). The land breeze develops in the nighttime until it changes to a sea breeze in the afternoon. According to this, benzene tends to spread southeastward until the morning, and then northwestward in the early afternoon. Compared with the daytime emission-started cases, benzene tends to spread more quickly in the initial few hours in the nighttime emission-started cases. The detectable range of benzene can be observed mostly within the radius of the 10 km region from the source in one hour, and even at the 15 km distance it is affected in a few cases. In three hours, it is detectable beyond the 10 km radius in most of the downstream region. The spreading region becomes narrow in direction and more oriented toward the ocean, reaching out to the 15 km radius in 6 h (08:00 LST), as the land breeze develops more strongly at this time near the source point. When the sea breeze develops after 12 h (14:00 LST), benzene is not detectable in the southeastern area and it is carried to the land area and is detectable beyond the 15 km radius. The benzene concentration exceeds the annual standard and even reaches the dangerous level at the 5 km distance immediately within one hour after the emission. The dangerous area expands further and gets close to the 10 km radius in the downstream in 3 h. Compared with the case at 05:00 LST, the areas exceeding the annual standard and the dangerous levels tend to reduce within 6 h (08:00 LST), and this seems to be related to the dilution of benzene in the atmospheric boundary layer in the morning. Both the areas exceeding the annual standard and the dangerous level are confined within the 5 km radius within 12 h after the emission (14:00 LST). The comparison between the daytime and the nighttime emission-started cases suggests that the dispersion can be more dangerous in the nighttime cases, especially for the initial few hours. This can be explained by the difference in the three-dimensional wind patterns. The air pollutant spreads more quickly and reaches further in the nighttime (c.f., Figures 5b and 6b). Note that the near-surface wind is stronger in the daytime due to the development of the sea breeze, which would otherwise transport the pollutant further from the emission source than in the nighttime. However, a strong vertical motion (Figure 4a) tends to increase the vertical transport of the pollutant and therefore decreases the near-ground-level concentration in the daytime. The daytime heating over the land can induce a vertical motion. In addition, the sea breeze tends to decelerate as it penetrates through the land, and The comparison between the daytime and the nighttime emission-started cases suggests that the dispersion can be more dangerous in the nighttime cases, especially for the initial few hours. This can be explained by the difference in the three-dimensional wind patterns. The air pollutant spreads more quickly and reaches further in the nighttime (c.f., Figures 5b and 6b). Note that the near-surface wind is stronger in the daytime due to the development of the sea breeze, which would otherwise transport the pollutant further from the emission source than in the nighttime. However, a strong vertical motion (Figure 4a) tends to increase the vertical transport of the pollutant and therefore decreases the near-ground-level concentration in the daytime. The daytime heating over the land can induce a vertical motion. In addition, the sea breeze tends to decelerate as it penetrates through the land, and this induces horizontal convergence and a vertical motion. The daytime heating and the development of the atmospheric boundary layer also contribute to the vertical dilution of the air pollutant. On the other hand, the upward motion is quite suppressed in the nighttime. Further, the stable atmospheric boundary layer reduces the vertical dilution and it tends to increase in the near-ground-level concentration. During the nighttime, the overall wind magnitude is weak and the land breeze is not as pronounced as the sea breeze in the daytime. This also makes the prediction of the dispersion direction difficult. There are a few cases when the dispersion direction is opposite to the prevailing wind direction, such as the case shown in Figure 6a, and this indicates that the nocturnal wind variation is not as regular as that in the daytime.

Risk Assessment
Previous time-averaged concentration patterns do not necessarily represent the vulnerable area by the chemical accident. When the wind is changing drastically during the accident, the time-averaged concentration tends to be low for the specific location, although that spot may be temporarily exposed to the dangerous level of 30 µg m −3 in the concentration. In this regard, the accumulated hit rate for each benzene concentration exceeds the dangerous level from the initial time to the designated hour and is examined for the daytime emission-started cases and the nighttime cases, whose results are shown in Figures 7 and 8, respectively.
Overall, these accumulated hit rate patterns show more smoothed patterns in space when compared with their counterparts of the time-averaged patterns (i.e., Figures 5 and 6). During the daytime emission-started cases (Figure 7), the vulnerable area tends to expand with time. In one hour, the hit rate near the 5 km radius from the emission source is less than 3%, and then it becomes 3% in the downstream region. In 6 h, the hit rate reaches as high as 10% at the 5 km radius, particularly in the downstream side. In 12 h, the hit rate tends to increase by more than 15% for those regions. In the meantime, the results indicate that the risk area more than 3% of the hit rate is mostly within the radius of 10 km from the emission source. Assuming that 3% is the threshold for taking action such as the enforcement of the evacuation of residents, the area within the 10 km radius from the industrial complexes should be carefully monitored for an immediate response to a chemical accident that may occur in the daytime.
Throughout the nighttime emission-started cases (Figure 8), the hit rate tends to increase in the downstream side, following the land breeze until the morning. The most contrasting feature to the daytime cases is the rapid increase of the hit rate for the initial 3 h. The hit rate becomes as high as 20% near the 5 km radius from the source in one hour, and then 20-30% in 6 h. In 12 h, the hit rate tends to increase more than 35% in the downstream side at the 5 km distance. In 12 h, the circulation is reversed and the hit rates tend to be steady or slightly decrease in the windward side. The result suggests that the risk area with more than 3% of the hit rate can expand up to the 15 km radius within 6 h after the accident. Therefore, a broader region should be monitored for an accident that may occur in the nighttime.
Previous time-averaged concentration patterns do not necessarily represent the vulnerable area by the chemical accident. When the wind is changing drastically during the accident, the timeaveraged concentration tends to be low for the specific location, although that spot may be temporarily exposed to the dangerous level of 30 μg m −3 in the concentration. In this regard, the accumulated hit rate for each benzene concentration exceeds the dangerous level from the initial time to the designated hour and is examined for the daytime emission-started cases and the nighttime cases, whose results are shown in Figures 7 and 8, respectively.  The risk assessment based on the accumulated hit rate in the above highlights that the vulnerable area expands more quickly in the nighttime emission-started cases in the downstream side, thereby increasing the potential health risk in the affected region. To examine this difference in a quantitative way, this study compares the accumulated hit rate averaged over land within the 5, 10, and 15 km radius from the emission source ( Figure 9). In one hour, most of the exposures are concentrated within the 5 km radius after the daytime emission starts. The exposures over the area within the 5 km radius are less frequent in the nighttime emission-started cases, with the probability of about half of the daytime emission-started cases, but it shows the exposures even at 10 and 15 km. This indicates that the nighttime dispersion spreads much farther than the daytime dispersion so that it could affect the residents in the farther area. The exposures increase steadily for the daytime emission-started cases at all distances. On the other hand, the probability of exposures seems to be saturated after 3 h in the nighttime emission-started cases, showing almost comparable values to the values after 6 h. Then, it tends to decrease slightly within 12 h, when the local circulation changes into the sea breeze in the afternoon. Overall, these accumulated hit rate patterns show more smoothed patterns in space when compared with their counterparts of the time-averaged patterns (i.e., Figures 5 and 6). During the daytime emission-started cases (Figure 7), the vulnerable area tends to expand with time. In one hour, the hit rate near the 5 km radius from the emission source is less than 3%, and then it becomes 3% in the downstream region. In 6 h, the hit rate reaches as high as 10% at the 5 km radius, particularly in the downstream side. In 12 h, the hit rate tends to increase by more than 15% for those regions. In the meantime, the results indicate that the risk area more than 3% of the hit rate is mostly within the radius of 10 km from the emission source. Assuming that 3% is the threshold for taking action such as the enforcement of the evacuation of residents, the area within the 10 km radius from the industrial complexes should be carefully monitored for an immediate response to a chemical accident that may occur in the daytime.
Throughout the nighttime emission-started cases (Figure 8), the hit rate tends to increase in the downstream side, following the land breeze until the morning. The most contrasting feature to the The results suggest that the probability of exposures can change much depending on the starting time of emission. It also depends on the distance from the emission source, which makes it difficult to decisively state which accident is more dangerous. Within the initial 6 h after the emission starts, the nighttime event shows more probability of exposures, especially in the distant area. However, the probability of exposures is reserved in 12 h, when the diurnal land-sea breeze circulation is reversed. Therefore, the duration of the accident is also important for assessing the risk quantitatively.
As shown in Figure 1, the urban downtown area of Ulsan is located within the 5 km radius northwest of the emission point. Within one hour of emission, the daytime emission cases are more dangerous because of the sea breeze blowing over the land area. The probability of risk is almost twice as large than the case in the nighttime emission cases. As the emission lasts longer, the probability of exposures becomes comparable between the daytime and the nighttime emission-started cases.
of the daytime emission-started cases, but it shows the exposures even at 10 and 15 km. This indicates that the nighttime dispersion spreads much farther than the daytime dispersion so that it could affect the residents in the farther area. The exposures increase steadily for the daytime emission-started cases at all distances. On the other hand, the probability of exposures seems to be saturated after 3 h in the nighttime emission-started cases, showing almost comparable values to the values after 6 h. Then, it tends to decrease slightly within 12 h, when the local circulation changes into the sea breeze in the afternoon. The results suggest that the probability of exposures can change much depending on the starting time of emission. It also depends on the distance from the emission source, which makes it difficult

Discussion
The CALPUFF simulation results in this study show that the dispersion of HAPs is affected significantly by horizontal transport driven by local circulation in the diurnal timescale and also by vertical diffusion by the diurnal variation of the atmospheric planetary layer. These dependencies are demonstrated well by sufficient dispersion tests driven by different weather conditions. Previous studies often assessed the impact assessment based on a few selected cases of numerical simulation, which is hard to draw a generalized conclusion. This study constructed the system and tested for the case of benzene dispersion, which can be easily applied to other HAPs.
It is worth comparing the results obtained from this study with the guidelines in the 2016 Emergency Response Guidebook distributed by the United States Department of Transportation (ERG2016) [34]. ERG2016 is the latest version intended for use by emergency responders during the initial stage of the leaking and atmospheric dispersion of hazardous materials due to various types of accidents. Large spill cases in their scenario involve the release of more than 208 L for liquids or 300 kg for solids, which approximately corresponds to our emission scenario (339 kg h −1 ). The guidebook provides the protective action distance (PAD) for each airborne material for the daytime and the nighttime, respectively, which is the predicted size of the downwind areas being affected by a cloud of toxic gas. The people in this area should be evacuated and/or sheltered inside buildings. The PAD in the nighttime is mostly more than doubled compared with the distances in the nighttime, due to the less vertical mixing of the toxic plume with the atmosphere. On average from our experiments, the dangerous level of benzene concentration (i.e., 30 µg m −3 ) can be detected in the downwind areas at a 4.0 km distance from the source point during the daytime, and 8.0 km at the nighttime, being consistent with the guidebook. ERG2016 recommends the initial isolation distance for fire and large spills of benzene as 800 m in all directions. Even though the guidebook indicates the isolation distance can be more than doubled in the case of strong nocturnal inversion, it may still not be sufficient in our tested cases for Ulsan in the summer. The dangerous level could be detected even at 6.2 km in the daytime, and 14.8 km in the nighttime from our tests on specific days. This suggests that the local meteorological condition is critical in the emergency response to chemical accidents. As this numerical simulation system based on the CALPUFF model can reflect the real-time meteorological condition, it is useful not only for assessing the potential risk and establishing long-term measures but also in responding to actual accidents and enforcing instant safety measures for residents.
It is noted that the simulation results and quantitatively assessed risk are also season-dependent, and the result is quite different when the simulations are conducted for the winter when the diurnal circulation becomes weak and the northwesterly winds prevail. Prevailing wind will disperse HAPs more quickly toward the ocean, which maintains a low probability of exposures in the inland area. Indeed, the diurnal variation of wind enhanced in the summertime due to a stronger land-sea heat contrast tends to increase the risk of exposure, particularly in the initial few hours. Then, the ventilation impact of a strong sea breeze tends to offset the increase of exposure probability after 3 h of emission.
Limitations of this study also exist. Although this study has validated the simulated wind fields with the observations, it is not sufficient to verify the detail of the circulation pattern on the local scale. In particular, the CALPUFF model used in this study has great difficulty in representing the mountain-valley circulation. To improve this, wind observations in more locations within the computation domain are needed. The weak land breeze in the nighttime could be partially related to the insufficient cooling over the mountainous land area in the nighttime.

Conclusions
This study developed a numerical simulation system to predict the dispersion of HAPs and assess the risk of exposures over a populated megacity after a chemical accident in the industrial complex nearby. The system is based on CALPUFF, which is a non-steady-state Gaussian puff model capable of representing the Gaussian diffusion and advection of pollutants by three-dimensional winds. A test scenario was the accidental release of benzene by the fictitious explosion of the storage container, whose amount was arbitrarily specified but based on the amounts in the past accidents. Multiple simulations were conducted for 1-31 August 2013 at the horizontal resolution of 200 m to obtain generalized dispersion patterns and the hazardous area. Each experiment was conducted for each day for 12 h, and started from 14:00 LST for the daytime emission-started case and 02:00 LST for the nighttime case. The patterns of the time-averaged concentration and the accumulated hit rate eliminated the day-to-day weather variability and showed the impacts from the diurnal variation of the local land-sea breeze circulations and the dilution in the atmospheric boundary layer.
The dispersion direction is much affected by the local circulation changes in the diurnal time scale. For the daytime emission-started case, the on-shore wind develops from the south from the afternoon to early evening, driven by differential heating between the land and ocean, and then the off-shore wind develops from the west and north in the late evening to midnight. Benzene shows the northward transport toward the residential area by sea breeze as well as a gradual dispersion in time. The dispersion direction is reversed in the nighttime when the sea breeze ceases and the land breeze develops in the nighttime. For the nighttime emission-started cases, land breeze tends to spread benzene southeastward from the initial time to the morning, and then the transport moves northwestward when the sea breeze develops in the afternoon. When the cases are compared between the daytime and the nighttime emission-started cases, the HAPs dispersion can be more dangerous in the nighttime, especially for the initial few hours due to the stable atmospheric boundary layer and the increase in the near-ground-level concentration. During the nighttime, the overall wind magnitude is weak and the land breeze is not as pronounced as the sea breeze in the daytime.
The probability of HAPs exposures to the harmful level determined by the US EPA standard changes much in time and space, depending on the starting time of emission. The exposures also depend on the distance from the emission source. Within the initial few hours after the accident, the nighttime event shows a higher probability of exposures. On the other hand, the probability is reserved in 12 h, when the circulation changes in the diurnal time scale. In this regard, how long the accident may last is another important factor to assess the risk quantitatively.
The potential risk of HAPs dispersion in the test area in Ulsan turns out to be serious as it is located within the 5 km radius from the center of the chemical complexes. In the summertime, when there exists a strong diurnal sea breeze, the toxic pollutants can impact the downtown area within an hour after the accident, and the probability is almost twice as large as the case in the nighttime emission. As the emission continues, the probability of exposures becomes comparable between the daytime and the nighttime emission-started cases.
There should be much to improve in the dispersion modeling systems in CALPUFF. In particular, more detailed wind patterns should be provided to resolve the microscale circulation such as the mountain-valley diurnal winds. This can be done by either substantially increasing the meteorological observation sites or by providing high-resolution meteorological reanalysis data to the CALMET system. There should be much room for further improvements to the developed system in this study, and the quantitative risk assessment made here could be revised.