## Abstract

Since highway traffic has become one of the major emission sources of air pollution, air pollution prediction near roadway tunnel portals is a very important subject. Although many models have been suggested to predict pollutant concentrations near roadways, almost all models can be applied to only at-grade or cutoff straight highways. Therefore, a numerical model applicable to the site near roadway tunnels in complex terrain has been developed.

The first stage of this study is to make a database of air quality and meteorological conditions near roadway tunnel portals. The second stage is a screening of several wind field models. The third stage is an evaluation of the numerical schemes for the advection equation, mainly carried out based on the results of the rotating cone problem.

In this limited comparative study, the most accurate and high-speed computing scheme was the Taylor–Galerkin scheme. Next, a three-dimensional model based on this scheme was developed by operator splitting of locally one-dimensional calculations.

The final stage is a validation study of the proposed model. The composite model consists of a wind field model, a model for the jet stream from a tunnel portal, and a model for the diffusion and advection of pollutants. The calculated concentrations near a tunnel portal have been compared to air tracer experimental data for two actual tunnels: the Ninomiya and the Hitachi Tunnels. Good evaluation scores were obtained for the Ninomiya Tunnel. Since predictive performance for the Hitachi Tunnel was not sufficient, some additional refinements of the model may be necessary.

## Introduction

Highway traffic has become one of the major emission sources of air pollution, and especially automobile exhaust gas pollution near roadway tunnel portals can sometimes become a serious environmental problem. Therefore, air pollution prediction is a very important subject, and many models have been proposed to predict pollutant concentrations near roadways. Although these models can usually be applied to at-grade or cutoff straight highways, they are usually not well applied to tunnels that are located in an area of complex terrain. Well-known numerical models for complex terrain that have been developed in the last decade, for example, by Tesche et al. (1987) and Uliazs (1993), are mainly prepared for larger-scale areas than the objective area concerned here, namely, near a tunnel portal.

The first stage of this study is to carry out air tracer field experiments in order to prepare a database. The second stage is a screening of wind field models, and the development of a suitable wind field model. The third stage is an evaluation study on the numerical methods for the advection equation to incorporate the composite model and the development of the three-dimensional numerical model suitable to the air quality simulation near a tunnel portal. The final stage is the validation study of the proposed numerical model.

The paper describes the development of a three-dimensional model for air pollution simulation near roadway tunnel portals. The procedure and results with regard to the wind field model screening appear in Kimura et al. (1996) and only a brief summary is written in this paper. The matters concerning field experiments are outlined in section 4.

## Selection of numerical method for advection equation

The numerical method for the advection equation may be the most important factor in composite atmospheric models because an inappropriate descretization of the advection terms can sometimes have a devastating effect on a numerical model and should be selected carefully. There are many review papers concerning the numerical methods, such as Rood (1987), McRae et al. (1982), and Chock (1991), and it was decided to carry out an evaluation study to confirm their conclusions. The candidate numerical methods chosen for this evaluation study are four methods as follows: simple upwind differencing method of second order, Taylor–Galerkin method selected from the schemes of Galerkin-type finite element method, cubic spline method selected as a representative of semi-Lagrangian-type schemes, and particle method because this is sometimes used for Monte Carlo simulation. The spectral and pseudospectral methods were not chosen because of difficult correspondence with complicated boundary condition over complex terrain.

### Description of methods

#### Upwind differencing method

The upwind differencing is one of the most popular numerical models, although it is largely concerned with numerical diffusion. Its advantage is that it is very easy to understand and that it suits computer programs. We have chosen this method of the second-order accuracy for the underlying evaluation study. The upwind differencing of a space derivative is shown in Eq. (1):

where *C* is concentration, *x* is a spatial coordinate, and Δ*x* is a spatial increment.

#### Taylor–Galerkin method

Although the Galerkin finite element method using a chapeau function is an attractive scheme for the discretization of a space derivative, the combination with time differencing makes the scheme unstable because of computational pseudonegative diffusion. The Petrov–Galerkin and Taylor–Galerkin methods were developed to overcome these defects. The Taylor–Galerkin method was proposed by Donea (1984) and Donea et al. (1987). Chock (1991) carried out an extensive evaluation study for the numerical methods and concluded that the Taylor–Galerkin method is one of the best choices for solving the advection equation.

The one-dimensional advection equation descretized by Taylor–Galerkin method is shown in Eq. (2):

where

where *t* is time, Δ*t* is the time increment, and *μ* is the Courant number. Moreover, the wind derivatives in Eq. (3) are approximated by first-order center differencing.

#### Quasi-Lagrangian cubic spline method

In the field of air quality simulation, the quasi-Lagrangian cubic spline method was introduced by Pepper et al. (1979). This method uses the following basic concept: the concentration of place *x*_{i} at time (*n* + 1)Δ*t* is identical to that of place *x*_{i−1} at time *n*Δ*t.* To estimate the concentration of place *x* (*x*_{i−1} < *x* < *x*_{i}), the following interpolation scheme was used:

where *D*_{i} is a spatial differential and can be obtained by solving the following equation:

where Δ*x*_{−} = *x*_{i} − *x*_{i−1} and Δ*x*_{+} = *x*_{i+1} − *x*_{i}.

#### Particle method

The dispersion phenomena can be simulated by the movement and distribution of a large number of particles, and the atmospheric turbulence is modeled by a pseudorandom number generated by a computer. This method is known as particle method or random-walk method.

In the underlying comparative study this method was selected as a candidate to solve the advection and diffusion equations. However, to evaluate only the advection term, the diffusion process was omitted and movement of a particle was expressed by using only the mean wind field. The initial number of particles was set proportionally to the initial concentration field.

### Evaluation of methods

The computational performance of the numerical methods for the advection equation can be measured by the rotating cone problem. In this test, the initial concentration field is represented by the cosine function, and its center is biased from the center of the circulating flow field. The two-dimensional advection equation is shown below:

where *C* is concentration, and *u* and *υ* are wind components for *x* and *y* directions, respectively. Chock and Dunker (1983) and Chock (1985, 1991) have carried out an extensive evaluation study for the numerical methods, using the rotating cosine hill test in a 33 × 33 grid. Rotating cone testing was carried out to reconfirm the evaluation results given by Chock. The experimental conditions shown in Figs. 1 and 2, including the initial concentration field and the wind field, were the same as that of Chock (1991).

Seibelt and Morariu (1991) used the deformational flow field to test their semi-Lagrangian numerical advection methods. Their purely deformational flow field is expressed as follows:

where *α* is a constant.

In this case, the deformation *F* cannot be zero, even though the divergence becomes zero:

In this flow field the calculated mass after *n* time increments can be analytically obtained. This analytical result suggests that a continuous growth of the total mass may occur in the ordinary Lagrangian model in which the first-order difference of the spatial derivative of wind is used.

To evaluate the numerical methods, this deformational flow field was also used. The initial concentration field was set for a rectangular-shaped block of 8 × 8 grid elements in the center of a 32 × 32 computational domain, which is the same as that of Siebert et al. (1991). The flow field for the former half computational cycle is shown in Fig. 3, and the flow for the latter half cycle is in reversed direction. Therefore, the rectangular block stretches along the *y* axis. After that, the flow becomes reversed, stretches along the *x* axis, and returns to the shape of initial conditions after one cycle of computation.

As for the rotating cone and deformational flow problems, the two-dimensional equation shown in Eq. (5) was numerically solved by operator splitting (McRae et al. 1982). In the rotating cone problem, the maximum Courant number *μ* was set at 0.4. This scheme was inappropriate (only for the upwind differencing) for *μ* = 0.4, and so the Courant number was set at 0.2. The number of particles used in the particle method was 500. The calculated concentration distributions of the rotating cone after two revolutions are shown in Fig. 4. The results of the deformational flow problem after five cycles are shown in Fig. 5. Good performance was obtained for cubic spline and Taylor–Galerkin methods. The upwind differencing method reveals the worst in this comparison because the concentration field after two revolutions reveals a quite different shape from the initial condition. In the particle method, a little bias of the center position and some spikes were observed. These errors may have occurred because the first-order differencing of time derivative was used, and the particle moved along the tangential line of airflow at the position of this particle. If a higher-order time difference is adopted, these error may be suppressed.

When Taylor–Galerkin and spline methods were used, small ripples occurred behind the cosine hill, and some negative values appeared. The adoption of a numerical filter may be useful to suppress these ripples, and an additional study seems necessary. As for nearly nonreactive tracers, this defect does not seem to be fatal. More detailed considerations are made in the evaluation of the composite model.

For the quantitative evaluation, the mass conservation ratio *R*_{1} and the mass distribution ratio *R*_{2} were calculated as follows:

where *C*_{ij}(0) are initial concentrations and *C*_{ij}(*t*) are concentrations at the time of evaluation.

The evaluation scores for each numerical method are shown in Figs. 6 and 7, and the computational time is shown in Fig. 8. These scores suggest preference for the use of the Taylor–Galerkin method, which is comparable with Chock (1991).

### Example of Taylor–Galerkin model: Comparison with Gaussian model

The Taylor–Galerkin method turns out to be the most attractive method in this limited comparative study. Therefore, this method was selected for the basic numerical scheme, and the three-dimensional composite model was produced by operator splitting based on the locally one-dimensional Taylor–Galerkin method for the diffusion and advection equation.

This model was initially applied to the most simple case of the effluent diffusion from a tall stack where the wind and diffusion coefficients were constant. An analytical equation was derived that is identical to the Gaussian plume equation:

where *H* is source height, and *Q* is source strength; *K*_{y} and *K*_{z} are diffusion coefficients in *y* and *z* directions, respectively.

The calculated concentration fields based on the Taylor–Galerkin model and the Gaussian model are shown in Fig. 9. The results show sufficient consistency between these models.

## Model overview

### Model structure

The objective of the proposed composite model is to simulate automobile exhaust gas diffusion from a tunnel portal. The model consists of two major modules and three submodules. The first major module is a wind field module for calculating the three-dimensional wind components over a complex terrain; it contains a tunnel submodule to simulate the jet stream from a tunnel portal. The second major module is a diffusion module to calculate the concentration field of automobile exhaust gas from the tunnel portal; it contains two submodules:a meteorological preprocessing submodel to calculate the diffusion coefficients and an emission submodule that gives the emission rate for roadways.

The schematic diagram of this simulation model is shown in Fig. 10. Basic input data for this model concerns the terrain, meteorology, traffic, and tunnel configurations. The terrain data include the terrain elevation and land surface representation.

The NO_{x} emission intensity for a roadway is estimated by the emission submodule based on the traffic volumes for each automobile type and traffic conditions, and details for the estimation methods for NO_{x} emissions are shown in Okamoto et al. (1990).

### Wind field estimation

Only a brief summary of the comparative study for wind field models is presented here. The model names and evaluation scores are shown in Table 1. For this comparative study, three companies have cooperated. Each company was supplied with terrain data and wind data at one point (the most representative anemometer site shown in Table 2) for two runs of the Ninomiya experiment. They submitted the calculated three-dimensional wind components for all grid points. The calculated wind at 10 anemometer sites were extracted from the grid point data and compared with observed data. The companies did not participate in the process to evaluate the wind field model.

Table 1 shows that some fluid dynamic models do not reveal a good performance. The main reason may be that the data supply was not satisfactory (some fluid dynamic models required more detail terrain information about outer computational domain), rather than that the models themselves were not sufficient.

According to the comparative study on the wind field models (Kimura et al. 1996), a mass consistent (MASCON) model has been selected for the wind field model. Since we could not find a perfectly fitting fluid dynamic model, and the difference of the predictive performance between the candidate models was not so large, we adopted the model with the least computational requirements.

The MASCON wind field model was developed according to Dickerson (1978) and Sherman (1978). Although the wind field model developed here is three-dimensional and similar to Sherman’s MATHEW model, the terrain-following coordinate system was employed. The first step is to obtain a guess field by using the weighted-average interpolation scheme from several meteorological points, in which the weight is proportional to the inverse of the square distance from the grid point to each measuring point, and the vertical extrapolation is made by using the power law of the wind profile (Turner 1994). There were 10 wind measurement points for each experimental site, and wind measuring was carried out at height of 6–10 m above terrain.

The second step is to modify the guess field so as to satisfy the continuity equation. The coordinate system (*x, y, z*) was converted to the terrain adjusted coordinate (*ξ, η, ζ*), and the continuity equation in this system is shown in Eq. (12):

where *h*(*x, y*) is the terrain elevation at the place (*x, y*), and

In the variational wind field model, the airflow through the obstacles varies depending on the ratio of the coefficients *α*_{1}/*α*_{2} [the same notation was used as in Sherman (1978).] When the atmosphere is stable, the airflow does not climb over a ridge but tends to follow the same altitude. This situation is simulated by using a larger *α*_{2} value. When *α*_{1} is equal to *α*_{2}, this flow is equivalent to the potential flow. Ueyama et al. (1984) studied the relation between the calculated wind field and the ratio *α*_{1}/*α*_{2} in their variational model, and they concluded that the most appropriate value for the ratio *α*_{1}/*α*_{2} in moderately complex terrain is 0.7 for stable and unstable conditions. Therefore, this value was chosen for all cases of the simulations at this stage. A detailed description of these procedures is given in Kimura et al. (1996).

For the lowest layer on the roadway, the vehicle-induced wind plays an important roll compared to the ambient wind. The model describing this effect is explained together with the model for jet stream from a tunnel portal. The coordinate system and grid spacing for this MASCON model is equivalent to that of tracer concentration calculation, and details are given in section 4.

### Jet stream from a tunnel portal

The velocity of the jet steam at a tunnel portal can been calculated with a traffic piston equation:

where

*ζ*_{e}: tunnel entrance loss coefficient,*λ*: tunnel wall friction loss coefficient,*L*: tunnel length (m),*D*: tunnel diameter (m),*A*_{r}: tunnel cross sectional area (m^{2}),*V*_{t}: traffic speed in tunnel (m s^{−1}),*U*_{0}: velocity of jet stream at the portal (m s^{−1}),*n*: number of vehicles in tunnel, and*A*_{m}: equivalent resistance area of the vehicles (m^{2}).

The velocity of a jet stream from a tunnel portal was calculated as follows. The first step is to calculate the jet stream at the portal. The second step is to calculate the decrease of the jet velocity (*x* axis: direction along the roadway) from the portal. Next, the distribution of the jet velocity in the *y*-axis and *z*-axis directions is calculated. This distribution is assumed to be Gaussian and shown in Eq. (14) as referred by Ueyama (1985):

where

*U*: flow velocity at location (*x, y, z*) (m s^{−1}),*U*_{0}: flow velocity at portal (m s^{−1}),*κ*: distance decrement coefficient,*σ*_{y}: standard deviation of wind in horizontal direction (m), and*σ*_{z}: standard deviation of wind in vertical direction (m).

The distance decrement coefficient *κ* varies by ambient wind direction, wind velocity, and traffic conditions; the values for this coefficient were determined based on the experimental results (Ueyama 1985). The functional forms of the *σ*_{y} and *σ*_{x} were estimated by the results of a scale model experiment published in Takaso et al. (1978).

The decrease of the jet stream was empirically modeled based on the Mikkabi Tunnel experimental data. The most typical cases are shown in Fig. 11. The exponential decay was observed up to the distance of five times the portal diameter and approached the asymptotic values that may be the function of the ambient wind speed, direction, and traffic conditions. The asymptotic values are large for nearly calm conditions, and relatively small for perpendicular and strong wind conditions. Therefore, these asymptotic values are assumed to be 3 m s^{−1} for the former case and 1 m s^{−1} for the latter case. Farther away from the portal, the jet stream from the portal decreased and merged with the vehicle-induced wind in the open field, and these effects cannot be clearly separated.

Based on these data, the portal jet Eq. (14) can be used up to the distance where the calculated velocity of the jet stream becomes equal to these asymptotic values. When the vehicle’s size is larger and driving speed is faster, the traffic-induced wind speed on the road is stronger. The traffic-induced wind is equivalent to this asymptotic value.

The jet stream velocity for the lower two grid points just above the roadway was fixed and did not vary under MASCON calculation. However, as for the other grid points, the jet stream velocity was combined with the guess ambient wind field calculated by the inverse-square-distance weighted averaging method, and its field was modified to satisfy the continuity equation in MASCON calculation.

### Diffusion coefficients

The diffusion coefficient *K*_{y} is assumed to be a function of Pasquill’s atmospheric stability and reveals the same value throughout the covered area, and *K*_{z} is an important parameter that has been set based on the boundary layer meteorological model.

The coefficient *K*_{y} has been calculated in Eq. (15), by using the diffusion parameter *σ*^{′}_{y} of the Pasquill–Gifford (PG) diagram. Since the objected area for calculation was within 100 to about 200 m from the emission source, the value of *σ*^{′}_{y} was taken from the PG chart of 100-m downwind distance,

where the assumed value of wind speed *u* in Eq. (15) is 3 m s^{−1}.

The estimation methods for the turbulent diffusion coefficient *K*_{z} in the boundary layer have been proposed by many researchers and one of the most simple and reliable methods may be the one described in Shir and Shieh (1974). The computational process can be summarized as follows: the stability parameter *s* proposed by Shir and Shieh is calculated by incoming solar radiation and wind speed, or with sky cover and wind speed. The stability length *L* is computed with stability parameter *s* and surface roughness parameter *z*_{0}, where *z*_{0} for each surface grid can be obtained by the terrain database. The nondimensional wind shear Φ_{m} is computed with stability length *L* and height *z.* The mixing height *H* is computed with the integrated incoming solar radiation. The mixing length *l* is computed with the von Kármán constant *k* and the mixing height *H.* The friction velocity *u** is computed from the von Kármán constant, the natural wind speed *u,* and the integrated value Ψ_{m} of the nondimensional wind shear Φ_{m}. Calculating Φ_{n} based on stability length *L* and height *z,* the diffusivity coefficient *K* is computed from friction velocity *u** and mixing length *l.*

These procedures were developed according to the reference of Shir and Shieh (1974). The turbulent diffusion on the roadway is heavily influenced by the vehicles. Therefore, the vehicle-induced diffusion coefficient *K*_{z0} was estimated based on the data published in Eskridge et al. (1979). This value *K*_{z0} is set 2.0 m^{2} s^{−1} for the lowest two levels of grid cells within the roadway.

### Numerical scheme

The basic mathematics of the diffusion module consist of the three-dimensional advection and diffusion equation on the terrain adjusted coordinate system. According to the review on numerical methods by McRae et al. (1982), the operator splitting method was used to solve the three-dimensional equation. The concentration at time level *n* + 1 can be expressed by

where *Ax* and *Ay* are the horizontal advection and diffusion operators, respectively: *Az* is the vertical advection, diffusion, source, and sink operator.

The operators *Ax, Ay,* and *Az* were all formulated by the locally one-dimensional Taylor–Galerkin numerical method. The boundary conditions were as follows:

top of the boundary: ∂

*C*/∂*z*= 0,terrain surface:

*K*_{z}∂*C*/∂*z*= 0,outflow side boundary: ∂

*C*/∂*x*= 0,inflow side boundary:

*C*=*C*_{background}, andsurface of the structure inside:

*K*_{n}∂*C*/∂*n*= 0, (*n*= perpendicular to the surface).

The *x* axis was made consistent with the direction of the roadway at the portal. The horizontal grid size is variable: 5 m for center 20 × 20 grid points and 30 m for the most outer grids. The vertical grid size is small near the surface (Δ*z*_{1} = 3 m) and large for higher altitude, and the top elevation of the highest grid cell is 200 m from the reference heights (origin of the coordinate system). The maximum number of the horizontal (*x, y* direction) and vertical (*z* direction) grid cells is 51, 51, and 18, respectively.

The time increment Δ*t* is set so as to maintain the numerical scheme stable, that is, 0.25–0.5 s depending on the Courant number. Since the area considered covers several hundred square meters, the stable concentration may be obtained after the appropriate time integration from the initial concentration field, usually 360 time increments (i.e., 180 s) under 1–3 m s^{−1} wind conditions. This stable concentration is considered as an hour-average concentration in this simulation.

## Database for the validation

The basic idea of this simulation model is formulated based on the experimental data at the Mikkabi Tunnel of the Tomei Expressway. This experiment was conducted in 1980 and 1981, and its scale may be the largest of its kind as an extensive field experimental project for the air pollution near tunnel portals (Ueyama 1985).

However, to validate the proposed model, an independent database that was not used for the development of the model is necessary. Therefore, we carried out two field experimental programs to prepare a sufficient database. These programs consist of an air tracer experiment, meteorological observations, and measurement of the traffic conditions at experimental sites as shown in Fig. 12. The first program was carried out at Ninomiya Tunnel (*L* = 445 m) on Odawara-Atsugi Road, which has a traffic volume of about 30 000 vehicles per day. The second site is the Hitachi Tunnel (*L* = 2439 m) on the Joban Expressway, with 24 000 vehicles per day. Both tunnels have two tubes and two traffic lanes per tube. The Hitachi Tunnel has a jet fan to promote ventilation and is operated only for the case that the level of visibility or air quality in the tunnel decrease to a predetermined condition.

The air tracer experiments presented a most important database for this validation study. The tracer used was SF_{6} and released from the inside of the tunnel, about 150 and 300 m from the portal. Duration of the tracer release was about 160 and 170 h for the Ninomiya and Hitachi experimental sites, respectively. About 30 receptor points were connected to two multichannel FPD (flame photometric detector) sulfur analyzers with a precut filter for SO_{2} gas. At these receptor points, air sampling and analysis were continuously carried out throughout the duration of tracer release time. These receptor points were mainly located along the highway. At the other receptor points, bag air samplers were operated during each hour. The sampled air was collected in 12-L polyvinylidene chloride bags and analyzed by electron capture gas chromatography. This sampling was conducted for 21 and 18 times for each program, respectively, and meteorological and traffic conditions during these sampling are shown in Table 2.

About 10 meteorological sites were set around the tunnel portal, and the elevation height of the wind sensor was 6–10 m. Six anemometers with a one-direction wind component were also set along the highway, just outside the left lane and inside the tunnel to measure the velocity of the jet stream from the portal. The vertical temperature difference, insulation, and net radiation were also measured. Video tape recording for the highway lanes was carried out to measure the traffic volume and traffic conditions.

For the Ninomiya Tunnel experiment, the wind data were obtained by nine anemometers. The meteorological sites are shown in Fig. 12. However, for the Hitachi Tunnel experiment, eight anemometers are located within the area shown in Fig. 12 and two other anemometers are outside this map; one is the top of the mountain and the other is located on a farm at about 500 m distance from the portal.

## Evaluation of model performance

### Calculated velocity of the jet stream at the portal

The strength of the jet stream at the portal was estimated by Eq. (13), in which the empirical parameters were determined based on the Mikkabi Tunnel experiment. Therefore, the validation for this tunnel portal model is necessary. The calculated velocity *U*_{0} at the portal was compared to the observed jet stream measured by the one-directional anemometer located just inside the portal. The scatter diagrams for the Ninomiya and Hitachi Tunnels are shown in Fig. 13. The result for the Hitachi Tunnel is not much different from that of Ninomiya. This figure suggests that this tunnel portal model can present an appropriate jet stream velocity at the portal, and this calculated value can be used in the air quality simulation.

### Calculated concentration field

Air quality simulation near tunnel portals has been carried out by using the meteorological and traffic condition data for Ninomiya and Hitachi experiments. In this simulation, the wind field was made by a MASCON variational model using the complete set of available wind data: 9 meteorological stations for Ninomiya and 10 stations’ data for Hitachi.

Comparison of the observed and calculated surface concentrations are shown in Figs. 14a–d. These figures show the typical cases for Ninomiya and Hitachi data. In run R1-1 the wind direction was parallel to the tunnel portal jet, and the contour lines of high concentration were stretched along the roadway. This tendency was reproduced by this simulation. In run R1-10, the wind direction was perpendicular to the roadway. The contour lines for calculated and observed SF_{6} concentrations were stretched to the south side of the roadway.

Comparisons for R2-1 and R2-4 of the Hitachi experiment are shown in Figs. 14c,d. For the cases of the Hitachi experiment, the calculated concentration distributions were limited to a narrower area near the portal. The reason for this discrepancy may be the underestimation of the wind near the surface, which is caused by the wind field model.

The portal of the Hitachi Tunnel is surrounded by a steep ridge, and the wind direction at the top of the ridge sometimes differed with that of the base of the valley. Although wind measuring points were located on both the ridge and valley bases, this situation was not sufficiently simulated by this MASCON model, and it may be a reason why evaluation scores for Hitachi are not as good as those for Ninomiya.

### Statistical scores

At the receptor points, the calculated concentrations were compared to the observed values. As for the calculated concentrations, no significant negative values (below −0.01 ppb) were seen, even without the use of numerical filters. Scatter diagrams for hourly normalized SF_{6} concentrations are shown in Fig. 15. The simulation model results reveal an overprediction for the Ninomiya and Hitachi data. The correlation coefficients of 1-h Ninomiya and Hitachi data were 0.904 and 0.485, respectively. At the Ninomiya Tunnel, a better predictive performance was obtained than at the Hitachi Tunnel. The most important reason for this difference may be that for the Ninomiya experiment, the MASCON variational wind field model works well because the terrain at the Ninomiya Tunnel is less steep than at the Hitachi Tunnel.

There are many plots of underestimation in the upper-left corner of the scatter diagram of the Hitachi Tunnel. These plots mainly correspond to the site vicinity to the portal. One of the reasons may be an underestimation of the extent of the width of the jet stream for the Hitachi Tunnel. As for the other sites of the Hitachi experiment, the tendency of overestimation can be seen. This reason may be the deficiency of the wind field model in steeper terrain. The treatment of meandering during the averaging time of concentration may also be an important reason. To predict the 1-h-average concentration, the fluctuation of wind direction during the 1-h period should be considered.

## Conclusions

An air quality simulation model applicable to the roadway tunnel portals located in the complex terrain was developed. This model is a kind of three-dimensional numerical model and uses the operator splitting method based on the Taylor–Galerkin descretization. This numerical method was chosen by the limited evaluation study for numerical methods. The three-dimensional numerical model based on this Taylor–Galerkin method presents good performance in computing speed and accuracy and is easy to use.

Some empirical parameters in each submodule were determined by data archives. Therefore, the validation study is very important in increasing the reliability of this simulation model. The predictive performance was validated based on the Ninomiya and Hitachi Tunnels’ experimental data prepared separately for validation study. The results indicate that the proposed model can provide the appropriate concentration distributions. However, since evaluation scores for the Hitachi experiment are not sufficient, more extensive evaluation and refinement of the model may be necessary.

## REFERENCES

_{6}tracer gas. J. Appl. Meteor., 18, 401–412.

_{x}emissions from a roadway. Atmos. Environ., 24, 1535–1544.

_{2}distributions in the St. Louis metropolitan area. J. Appl. Meteor., 13, 185–204.

## Footnotes

*Corresponding author address:* Shin’ichi Okamoto, Tokyo University of Information Sciences, 1200-2 Yatocho, Wakaba-ku, Chiba 265-8501 Japan.