Airborne Transmission Risk in Urban Buses: A Computational Fluid Dynamics Study

The Heat, Ventilation, Air Conditioning (HVAC) systems in urban buses promote high particle motion. In this paper the effect of centralized roof-top HVAC on virus transmission by droplets caused by sneezing and coughing was studied with the Lagrangian approach while the transmission through aerosols caused by breathing and talking was modeled with the Eulerian tracer method. It was found that the large droplets ( > 200 µ m) can travel more than 3 m without being affected by the airflow. On the other hand, the small droplets ( < 5 µ m) are easily dragged, dispersed


INTRODUCTION
The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has spread fastly around the world. Three transmission mechanisms have been identified: direct contact with contaminated surfaces (fomites), ballistic droplets, and aerosols. However, recent investigations suggested that fomites play a minimal role in SARS-CoV-2 community transmission (Harvey et al., 2021;Lewis, 2021). Regarding ballistic droplets, they are emitted by coughing and sneezing, and occasionally by talking and breathing. The large droplets cause local contamination by deposition around the emitter, and generally contain high viral load because they originate in the deep part of the respiratory system. The droplet size largely depends on the subjects, ranging from 1 µm to 1000 µm for sneezing Chao et al., 2009). However, more than 85% of droplets emitted during coughing, breathing, and talking are around 1 µm (Papineni and Rosenthal, 1997). For instance, coughing droplets range from 0.58 to 5.42 µm, and 82% of them are between 0.74 and 2.12 µm (Yang et al., 2007). These small droplets completely evaporate in air. The effect of the airflow over the droplets trajectory can be estimated by calculating the terminal velocity of droplets. For instance, droplets of 100, 200 and 500 µm have terminal velocities of 0.35 m s -1 , 1.4 m s -1 , and 9.2 m s -1 , respectively, but this is 4 × 10 -5 m s -1 for droplets of 1 µm. These small droplets are easily dragged by the airflow and need more than 10 hours to deposit. On the contrary the dynamics of large particles is mainly governed by the initial velocity, gravity and drag force. Busco et al. (2020) observed that droplets larger than 500 µm can travel more than 3 m until deposited on the floor. They also found that the motion of the head during sneezing affects the trajectory and maximum distance.
The importance of the third transmission mechanism (aerosols) was not clear until massive virus transmission events showed that virus is also transmitted by asymptomatic individuals through the exhalation of aerosols. Indoor environments with poor ventilation, and longer exposure time (Miller et al., 2021) were the common factors. In fact, the odds of transmission in Vuorinen et al. (2020) concluded that it should be necessary to inhale hundreds of aerosols to contract the disease. The WR model was also applied by Miller et al. (2021) using a zero-dimensional approach assuming homogeneous virus distribution in the domain (hereinafter called as "0D WR"). Based on real infection cases they proposed that the quantum emission rate (E) should be from 1 q h -1 to 2 q h -1 for breathing, from 1 q h -1 to 25 q h -1 for talking, and from 25 q h -1 to 125q h -1 for singing. Coarsely, these values should be taken as a reference. For instance, Ou et al. (2022) studied two real transmission cases where a person infected 8 passengers in a coach bus trip of 3 hours, and then infected 2 passengers more in a minibus trip of 1 hour. Ou et al. (2022) determined that the averaged ventilation rate in the coach and minibus were around 1.7 and 3.2 L s -1 per person (8-10 L s -1 per person recommended in offices and public spaces), respectively. For these two situations the 0D WR model estimated that the quantum emission rate should be between 35 q h -1 and 58 q h -1 .
The WR model was also applied by Dai and Zhao (2020), who identified that if people wear masks, then normal ventilation systems or natural ventilation was found to be sufficient to reduce the infection probability to 0.5% on buses, even at the top of the estimated range for E. However, if people did not wear face masks, then normal ventilation could only achieve 2.0% infection probability. Similarly, Krishnamurthy et al. (2020) applied the model and suggested that the probability of infection on a 2 hours bus travel running at 50% capacity with a single infected passenger should be around 19%. Finally, Shen et al. (2020) also applied the WR model. In relation to ground public transport, taxis were estimated to have the highest infection probability (3.2%), followed by coaches and school buses (2.9% and 2.2%, respectively), with the lowest infection probabilities in subway cabins and urban buses (both 0.6%).
The current paper aims to investigate the virus transmission mechanisms in urban buses. Sneezing and coughing are analyzed through Lagrangian methods in order to visualize the spread of droplets. On the other hand, the spread of aerosols due to breathing and talking is addressed through Eulerian methods combined with the WR model. Several locations for the infected passenger are considered, and the effect of the HVAC and the air renewal ratio is studied in depth by considering different operative conditions.

Eulerian Governing Equations
An Eulerian framework was used to model the airflow inside the bus. Simulations were performed using the reactingParcelFoam solver into the OpenFOAM-7.0 suite. This solver allows calculating the Lagrangian particle tracking as well as the Eulerian scalar transport of the gas exhaled. The continuity equation for the gas read as: where ρ is the air density, t is the time, and v is the gas velocity. The gas can be a mixture of several species. For the current problem evaporation is possible. Hence, air and vapor are considered. Each one is solved through a passive scalar transport equation: where Yi is the mass fraction of the specie i, and µeff is the effective viscosity, which is obtained from the κ-ω SST turbulence model. The momentum transport equation is: where p is the pressure, g is the gravitational acceleration, and Fs is a momentum source accounting for the drag force between the continuous and dispersed phases. The energy conservation equation of the continuous phase is written as: where h is the specific enthalpy, αeff the effective thermal diffusivity, Dp/Dt is the total derivative of p, and Fh is a heat source associated with evaporation. The effect of the body thermal plume was included in the model. A typical adult person at rest releases around 100 watts, but only 30% of this total power is directly transferred to the air through conduction/convection. The rest is transferred by radiation to the other surfaces and evaporation. Therefore, in the current model a power of 30 watts was homogeneously distributed over each mannequin.

Droplet Model: Lagrangian Equations
Due to their small size and low velocity, the droplets were assumed as spheres of variable size due to evaporation. Droplet fragmentation and coalescence were neglected, and the presence of non-volatile mass in droplets was considered in order to avoid complete evaporation. The particle trajectories were computed by Newton's second law. Only the drag (FD) and gravity (FG) forces were considered. The drag force was computed as: where dp is the particle diameter, and ρc is the air density. The drag force coefficient was computed through the stokes law: where Re is the Reynolds number. The evaporation mass rate Fm was calculated as: where M is the molecular weight, kc is the mass transfer coefficient, and Cs and C∞ are the molar concentrations of vapor on the droplet surface and in the bulk of the air, which are calculated from the ideal gas equation: where, psat and Tsat are the pressure and temperature of saturation, Rv is the constant of gasses, and xc is the concentration of vapor. Finally, the mass transfer coefficient kc was derived from the Sherwood number Sh for spherical droplets: ( ) Special Issue on Air Quality in a Changed World: Regional, Ambient, and where Dv is the diffusion coefficient of vapor in air, and Sc is the Schmith number, which is calculated from the Ranz-Marshall correlation.

Bus Model
Fig. 1 displays the bus studied, which is 11.1 m × 2.2 m × 3 m in length, width and height, respectively. The cabin has a total volume V of 47 m 3 with 40 seats and additional space for around 20 standing passengers. The driver is partially isolated by a vertical panel on the back. The bus has two floor levels, with the roof at 2.3 m in the front side, and at 1.9 m in the rear side.
The HVAC system is a centralized roof-top unit. This HVAC configuration is massively used in urban buses in Argentina. The opening of windows was not considered, so the air entered at 22°C through one inlet vent of 0.6 m × 1.6 m (around 1 m 2 ), and left the bus through two outlet vents of 0.6 m × 0.3 m (see Fig. 1). Many operating conditions were evaluated, considering air flow rates Qair of 0, 0.5, and 1 m 3 s -1 , and recirculation air ratios ϒrec of 100% (null fresh air), 75% (common in land transport media), and 0% (complete fresh air renewal). HVAC parameters were adopted from the commercial ThermoKing ATHENIA AM II 700, although the maximum ϒrec is restricted to 75% in the real device.
The RH was set to 40% and the bulk temperature was 27°C. In order to have suitable initial conditions, the problem was simulated for along 5 min before introducing the particles and tracers. For Lagrangian simulations the re-entry of particles trapped by the HVAC was not implemented. On the other hand, an ad hoc dynamic boundary condition was implemented for Eulerian simulations in order to consider the air recirculation with the corresponding renewal ratio.
Five mannequins representing the driver and passengers were included in the model. The computational domain was meshed using hybrid grids with local refining. A total of 886,544 cells were required for the whole bus, with cell sizes around 50 mm in the bulk and 10 mm close to the mannequins. The use of the face masks was only considered for the risk model.  Fig. 2(a) displays the particle injection angle used to represent the head motion. The velocity of the droplets and the exhaled air is also shown in Fig. 2(a). It was built by combining the dynamic pressure measurements by Busco et al. (2020) and the peak velocity reported by Bahl et al. (2020). The particles were injected using a cone injection model with an angle of 25° (Busco et al., 2020). As far as the authors know, the total liquid spread during sneezing has not been reported. Most CFD studies assume values around 5 mg (Busco et al., 2020;Löhner et al., 2021), which is similar to the reported by Xie et al. (2009) from coughing tests. However, we assumed a higher value of 50 mg with nonvolatile volume fraction of 6% (Aliabadi et al., 2010).

Sneezing parameters
Regarding the droplet size, Han et al. (2013) reported that unimodal and bimodal distributions are equally possible. Therefore, both of them were considered. From these distributions and the total liquid mass adopted, we calculated the number of droplets versus particle size, which is  shown in Figs. 2(b) and 2(c). For the unimodal and bimodal distributions, the total number of particles were 878 and 11,500, respectively. Chao et al. (2009) measured the droplet size distribution during coughing, while Xie et al. (2009) measured the average total liquid mass emitted, which was 4.25 mg and 7.7 mg for female and male, respectively. Based on these data, we calculated the number of particles versus particle size shown in Fig. 3(a), which represents a total of 1,660 particles. Regarding the velocity, the measurements from Gupta et al. (2009) tests were averaged (see Fig. 3(b)). As noted, the velocity profile is similar to the sneezing one.

Breathing and talking parameters
The use of Eulerian approaches for modeling the concentration of particles is based on the small size of aerosols (Vuorinen et al., 2020;Ou et al., 2022). The average tidal airflows are around 10 L min -1 for breathing and 13 L min -1 for talking (Gupta et al., 2010). Therefore, the same airflow (Q = 0.6 m³ h -1 ) was assumed for both respiratory activities. Considering that the risk model is based on the virus concentration around susceptible subjects, the inhaled virus can be approached without considering the real inhalation/exhalation cycle. Hence, only exhalation with a constant gas flow rate and unitary concentration scalar tracer was considered.

Risk Estimation Model
The WR model estimates the probability of contagious R as an exponential function of the cumulative inhaled virus n during a period of time D: where n depends on the virus concentration C, which obeys a differential equation taking into account the volume of the domain V, the virus removal rate λ, and the effective quantum emission rate Eeff, including the effect of facial mask with an efficiency ξex to trap the virus during exhalation (Eeff = E (1 -ξex)).

ORIGINAL RESEARCH
Then, the total quanta inhaled n(D) is: where ξin the efficiency of the mask to retain particles during inhalation.
In this work the mask efficiencies for exhalation and inhalation were ξex = 0.5 and ξin = 0.3, respectively. These low values correspond to the minimum efficiency measured with homemade face masks (Fischer et al., 2020;Asadi et al., 2020;Li et al., 2021).
The virus removal rate λ is the sum of four terms (λ = λven + λfil + λdep + λdec) corresponding to ventilation (λven), filtering (λfil), deposition over surfaces (λdep), and inactivation in time (λdec). λven is the number of times per hour the air is completely renewed with fresh air, which is calculated as: Clearly, the best condition is for the higher airflow rate and null air recirculation (Qair = 1 m 3 s -1 and ϒrec = 0). In this situation λven is 63 h -1 .
Urban buses are often equipped with 7-8 MERV filters, which have only around 60% of efficiency to trap particles between 3 and 10 µm. Therefore, a low λfil of 1.0 h -1 is adopted. Finally, values of 0.3 h -1 and 0.62 h -1 are recommended for virus deposition (λdep) and inactivation (λdec) (Buonanno et al., 2020).
The WR model can be used to estimate the minimum ventilation rate to keep an acceptable infection risk. Coarsely, the quantum emission rate has to be known (Ou et al., 2022). It is clear that E is the most important parameter, depending on the particular characteristics of the emitter (e.g., viral load), and the respiratory activity (talking aloud, singing, breathing, etc.). Buonanno et al. (2020) suggested that E should be between 1 and 5 q h -1 for breathing, but based on recent real bus trip cases, Ou et al. (2022) calculated values between 35 and 58 q h -1 . Fig. 4 shows the risk Fig. 4. Wells-Riley 0D model estimations for different quantum emission rates (E) and operation conditions. estimation from the WR-0D model is drawn for different E and operation conditions. Unfortunately, the risk increases exponentially with E (compare the continuous blue line and the green dashed line). That is, to use a unit value for E will not allow scaling the results to other values, and a reasonable value should be adopted. The risk also increases exponentially with time, especially for null air renewal cases (Qair = 0). The use of face masks, even for low trap efficiencies, clearly reduces the risk. For this reason, in the current study we preferred to adopt the value of 58 q h -1 calculated for Ou et al. (2022), but due to the strong relationship between the risk (R) and E, from now on we will report the relative risk, calculated as the ratio between the calculated risk and the higher risk obtained for the whole bus with worse operative condition (Qair = 0 without air renewal). Thus, the results could be used to quantify the improvement on the traveling conditions obtained for the different operative conditions. . The colors of particles indicate the particle size, which ranges from 0 to 1000 µm. Large particles (> 200 µm) quickly cover the three seat rows in front of emitter. Many particles also deposit in the lateral walls and the roof. Simulations considering motionless airflow (Qair = 0) were also carried out, leading to conclude that the HVAC has negligible influence for the unimodal distribution.

Sneezing and Coughing
On the contrary, for the bimodal distribution the airflow caused by the HVAC has a significant effect on particle motion. Fig. 5 on the right displays the location of particles at 0.5, 7, and 21 s. In this case the colorbar is limited to 200 µm and most particles are below 40 µm. The small droplets are easily dragged by the airflow, and some of them are trapped by the outlet vents. At 7 s it is clear how the particles are pushed to the lateral wall before ascending to the roof. This particular motion is a direct consequence of the location of the inlet and outlet vents, which promotes large vortex patterns and high mixing in the central zone. Fig. 6 shows the final location of particles after sneezing of the four passengers and the driver (the colors of particles refer to the emitters). The result corresponds to the bimodal distribution, which is the worst in terms of particle dispersion. As noted, most of the particles emitted by the driver impacts over the windscreen and the board. In this zone the air moves slowly and the contaminated region is limited to the front of the driver. On the other hand, the effect of HVAC, spreading the small droplets, is quite significant for passengers 1 and 2. For passenger 3, the proximity with one of the outlet vents allows capturing the small particles, thus reducing the contaminated area. Finally, the particles emitted by passenger 4 are not captured and spread through the whole rear zone. The HVAC has poor efficiency to capture particles larger than 200 µm, and relatively low efficiency to capture particles less than 50 µm.
The scenario is quite different during coughing because the particles quickly reach a small size less than 10 m. This is caused by the small initial size added to the fast evaporation. Fig. 7 shows the particle position in time after coughing of passenger 2. Similar to the sneezing bimodal case, the particles initially move to the next seat dragged by the HVAC airflow, ascend to the roof, and are captured by the outlet vents. Due to the low inertia, only a few of them escape from the outlet vents and reach the front and rear sides of the bus. The outlet vents are able to capture the small particles showing the importance of using HEPA filters or increasing the air renewal ratio. Fig. 8 quantifies the number of particles that remain suspended in air, deposit in surfaces and escape through the outlet vents after sneezing and coughing of the passenger 2. The number of particles smaller than 20 µm is also included. As noted, the two particle size distributions for sneezing lead to completely different results. For the unimodal case there is not any particle less than 20 µm, and all particles are deposited after 2 s. On the other hand, for the bimodal case around 2,000 particles remain in the air after 10 s. At the end of simulation (100 s) around 85% of particles deposited over surfaces, and the rest escaped from the cabin. As a criterion of the authors, the bimodal distribution is the most dangerous in terms of particle spreading and contamination of surfaces. Regarding coughing the situation is quite different because most particles remain in the air after 10 s, but around 50% of them are removed from the cabin after 50 s. As a consequence of the small size, around 12% of particles are still suspended in the air after 60 s.

Breathing and Talking
Simulations were performed to track the gas concentration from the different passengers and the driver during a typical trip of 20 min (mean trip time in half million population cities). Runs included three airflow rate conditions (Qair = 0, 0.5 and 1 m 3 s -1 ) and four recirculation ratios (ϒrec = 0, 50, 75, 100%). Fig. (9) shows the plumes of air emitted by the driver after 20 min with null  airflow (Qair = 0). The different plumes correspond to concentration thresholds of 0.005, 0.002, 0.0005 and 0.0001. The scale bar limit is fixed to 0.01 in all cases. The buoyant effect caused by the heat released by the manikins explain the upwards motion of the tracer. As noted, after 20 min the concentration in any place of the bus is higher than 0.0001 (0.01%), but it is around ten times more in the front side, where the driver is. A similar behavior is observed for all manikins. Fig. 10 displays the plumes of exhaled air concentration (threshold of 0.0001) for the passengers and the driver after 20 min. Note that the scale bar limit is fixed to 0.001. The pictures on the left correspond to null HVAC airflow (Qair = 0), which is clearly the worst condition because the concentration is higher than 0.001 in most of the bus for any emitter. On the other hand, the pictures on the right correspond to the best conditions with the maximum HVAC airflow (Qair = 1 m 3 s -1 ) and full fresh air (ϒrec = 0%). In these cases, the maximum concentration is around Fig. 9. Plumes of exhaled air concentration coming from the driver after 20 min without HVAC (Qair = 0). 0.0001, and mostly around the emitters. The pictures on the center correspond to the most common operation condition with Qair = 1 m 3 s -1 and air recirculation of ϒrec = 75%. In this case, the recirculation causes a fast mixing of the gas, although the fresh air helps to reduce the average concentration. As noted, the gas is spread through the whole domain. The maximum gas concentration is around 0.001 close to the emitter, but reduces up to 0.0003 in the rest locations. Clearly, the particular location of the HVAC inlet and outlet vents restrain the efficiency of the HVAC to remove the gas from the rear and front sides. Even for this optimal operation condition, the gas emitted by the driver and the passenger 4 hold around them.
The risk estimations based on CFD results were carried out using the following factor between gas concentration and virus concentration: Fig. 11 shows the virus concentration and transmission risk estimated with the WR model for the driver with two operation conditions (Qair = 0 and Qair = 1 m 3 s -1 with ϒrec = 0%). Although the risk estimations are strongly dependent on the quantum emission rate E (E = 58 q h -1 ), the comparison between the two extreme operation conditions, one with zero air renewal and the other with full air renewal, is of great interest. The graphs include the estimations using the tracer concentrations from CFD in each zone (see Fig. 1), the whole bus, and also the analytical WR 0D model estimation. As noted, the behavior is completely different for the case with and without HVAC. For the first case the virus concentration grows almost linearly with time in all zones, but specially in zone 1 (front of the bus). On the contrary, with HVAC the growth is only observed the first 5 min, and then the concentration, even in zone 1, remains almost constant. The 0D WR model predictions are always below the CFD predictions for the whole bus. This is because the 0D WR model takes into account the virus decay and deposition factors, which are not considered in CFD simulations. As noted, the local virus concentration in zone 1 is around two times more than the average. As expected from Eq. (10), a linear growth of the virus concentration (null air renewal) leads to an exponential increment on the risk. On the other hand, a constant virus concentration (full air renewal) leads to a linear increment of the risk. The maximum risk for the worst condition (Q = 0) after 20 min was around 3.2%, but this reduced up to 0.8% by renewing the air. Although these risk estimations seem to be low, the maximum risk quickly grows up to 9% without using face masks. Moreover, assuming a linear increment of the virus concentration (which is conservative), the risk estimation after 60 min should increase more than 4 times. Despite the relative reliability of the absolute values because of the uncertainty on E and the susceptibility of subjects to contract disease, the results are useful for comparison. Therefore, from now on the risk estimations will be reported as a relative risk of the maximum risk calculated for the worst case. Fig. 12 shows the relative risk in each zone and the whole bus. This is calculated as the ratio between the estimated risk and the maximum risk obtained for the whole bus for the worst condition (null air renewal condition). As expected, the higher values are for the null air flow rate. As noted, the local relative risk caused by the driver in zone 1 is nearly 2 times more than the average across the whole bus. In all cases the driver causes the highest risks in zones 1 and 2 Fig. 11. Virus concentration and transmission risk caused by the driver without air renewal (Qair = 0) and with full air renewal (Qair = 1 m 3 s -1 and ϒrec = 0%).
while passenger 4 is responsible for the highest risk in zone 5. For null airflow conditions the risk is a bit lower in the center of the bus, but the reduction is quite notorious if the HVAC is on with high air renewal ratios. As expected, the large reduction is found for ϒrec = 0% (full fresh air), but a significant reduction is also obtained even with ϒrec = 75%. The improvement is clearly noted mostly by analyzing the average relative risk in the whole bus (Avg), which is around 1 for null airflow, but reduces below 0.1 and 0.05 for Qair = 0.5 and Qair = 1.0 m³ s -1 with ϒrec = 0%, respectively. Fig. 13 shows the average and maximum relative transmission risk after 20 min due to each emitter. Results are for airflow rates of Qair = 0, 0.5, and 1 m 3 s -1 with ϒrec = 0 and ϒrec = 75%. The graphics on the left display the average values in the whole bus. As noted, for Qair = 0 (graph on the top left side) the average values are around 0.9, but these reduce below 0.1 and 0.05 for Qair = 0.5 and 1 m 3 s -1 with full air renewal (ϒrec = 0%). Despite full air renewal is not common in HVAC systems, the improvement is also very good even for partial air renewal of 25% (ϒrec = 75%), reducing the average relative risk below 0.4 and 0.3 for Qair of 0.5, and 1 m 3 s -1 , respectively (graph on the bottom left side). As noted, the improvement is not linear with Qair. As expected, the reduction in the maximum relative risk in the rear and front sides is not as good as that found for the average risk due to the particular distribution of the inlet and outlet vents. This is shown in the graphics on the right side. As expected, in most cases the maximums are given in the zone where the emitter is. Nevertheless, they remain less than 1.0 in most locations. Clearly, the analytical model cannot estimate the maximum local values (the WR 0D estimations are always averaged values), but it is noticeable the good agreement between this simple model and the average CFD results. Table 1 summarizes the local transmission risk estimated after 20 min with Qair = 0.5 m 3 s -1 and ϒrec = 75%. As shown, high local risks are still observed especially in the rear zone, but the average risk is reduced more than 50% in the whole bus.

CONCLUSIONS
The disease transmission in a typical urban bus was simulated through ballistic and aerosol transmission mechanisms considering the effect of the centralized roof-top Heat, Ventilation and Air Condition (HVAC) system. For ballistic cases (sneezing and coughing) the distance traveled by particles as well as the contaminated area around the emitters were determined. For sneezing cases two particle size distributions were considered: the first one characterized by large droplets and the other one by small droplets. It was found that large droplets (> 200 µm) can travel more than 3 m and deposit in less than 2 s without being captured by the HVAC outlet vents. On the other hand, the small droplets, characteristics of the bimodal distribution, were strongly affected by the HVAC airflow. The evaporation reduced the droplet size and around 15% of them became aerosols and remained suspended in the air for a long time. The bimodal distribution showed to be the worst in terms of airborne transmission.
Regarding coughing, the small droplets were easily captured by the HVAC when the emitter was close to the outlet vents. However, the studied HVAC shows to be very inefficient to renew the air in the rear and front sides of the bus.
The aerosols generated during breathing and talking were also studied, and the transmission  risk was calculated for many HVAC conditions. For null and low air renewal cases the aerosols quickly spread across the whole bus, reducing the local concentration. The transmission risk was estimated through the Wells-Riley model. Due to the uncertainty on the virus emission rate, a literature value was assumed (quantum emission rate of 58 q h -1 ) and the risk was reported as a relative value with respect to the maximum average one. This allows us to use the results to quantify the improvement on the indoor conditions for many HVAC operation conditions. It was found that a reduction of more than 50% is obtained with relative low HVAC airflow of 0.5 m³ s -1 with low air renewal ratios of 25%.