## 2 Mathematical Model  
### Governing Equations and Numerical Algorithm  
The airflow and pollutant diffusion at the tunnel portals are influenced by various factors. In order to explore the inherent law of the flow field, the following assumptions are put forward in this paper: (1) vehicles are regarded as external and accidental factors, with no influence on the airflow; (2) natural wind is regarded as the external boundary at a velocity of 0 m\(\cdot\)s\({}^{-1}\); (3) the airflow temperature at all inlet and outlet boundaries is 300 K, without any volume source; (4) as the pollutant components are scalar with similar distributional patterns, CO concentration distribution is used to represent the pollutant diffusion.  
The airflow velocity at the portals is much slower than the sound velocity, so it can be regarded as a steady turbulent flow of incompressible fluid [14]. This paper calculated the airflow at the inlet and outlet based on the Standard \(k\)-\(\epsilon\) turbulence model. The pollutant diffusion was described through the species transport equation. The buoyancy-driven effects caused by density changes due to varied components of mixed gas were not considered. The transport equations of the turbulent kinetic energy \(k\) and turbulent dissipation rate \(\varepsilon\) in the Standard \(k\)-\(\varepsilon\) turbulence model and species transport equation are shown in Equations (1)-(3).  
\[\frac{\partial}{\partial x_{i}}(\rho ku_{i})=\frac{\partial}{\partial x_{j}} \Bigg{[}\left(\mu+\frac{\mu_{t}}{\sigma_{k}}\right)\frac{\partial k}{\partial x _{j}}\Bigg{]}+G_{k}-\rho\varepsilon \tag{1}\]  
\[\frac{\partial}{\partial x_{i}}(\rho\varepsilon u_{i})=\frac{\partial}{ \partial x_{j}}\Bigg{[}\left(\mu+\frac{\mu_{t}}{\sigma_{k}}\right)\frac{ \partial\varepsilon}{\partial x_{j}}\Bigg{]}+C_{1\varepsilon}\frac{ \varepsilon}{k}G_{k}-C_{2\varepsilon}\rho\frac{\varepsilon^{2}}{k} \tag{2}\]  
\[\frac{\partial(\rho c_{s}u_{i})}{\partial x_{i}}=\frac{\partial}{\partial x_{ j}}\Bigg{[}D_{s}\frac{\partial(\rho c_{s})}{\partial x_{j}}\Bigg{]}+S_{s} \tag{3}\]  
In the above formula, both the values of \(i\) and \(j\) range from 1 to 3; \(x_{i}\) represents the coordinate components; \(\rho\), \(\mu\), \(\
u\), and \(u\) are the mass density of air, turbulent viscosity, coefficient of kinematic viscosity, and air velocity; \(c_{8}\) refers to the volume density of \(s\); \(D_{s}\) refers to the diffusion coefficient of \(s\); \(S_{s}\) refers to the quality of \(s\) generated from the chemical reaction in unit time and volume, namely, the productivity; \(C_{\mu}\), \(\sigma_{k}\), \(\sigma_{e}\), \(C_{1\varepsilon}\), \(C_{2\varepsilon}\), and \(C_{3\varepsilon}\) are model coefficients that follow the principle of asymptotic proportionality [16], and their values are 0.09, 0.24, 0.15, 1.44, 1.92, and 0.2.  
Equations (1)-(3) and the Reynolds-averaged Navier-Stokes equation [17; 18] constituted the governing equations for the pollutant diffusion of closely spaced twin tunnels. The finite volume method [19] was used to discretize the governing equations. For the discretization of derivatives, the skewness-corrected linear upwind interpolation scheme (second-order accuracy) was applied to the convective terms, and second-order linear interpolation was used with a nonorthogonal correction scheme for the diffusive terms [20]. The solution of the velocity field and pressure field was conducted through a fully coupled block algorithm. Due to the simultaneous solution of momentum and continuity equations, implicit block coupling of pressure and velocity variables leads to faster convergence with respect to classical, loosely coupled, segregated algorithms of the SIMPLE family of algorithms [21].