Zonal modeling of air distribution impact on the long-range airborne transmission risk of SARS-CoV-2

A widely used analytical model to quantitatively assess airborne infection risk is the Wells-Riley model which is limited to complete air mixing in a single zone. However, this assumption tends not to be feasible (or reality) for many situations. This study aimed to extend the Wells-Riley model so that the infection risk can be calculated in spaces where complete mixing is not present. Some more advanced ventilation concepts create either two horizontally divided air zones in spaces as displacement ventilation or the space may be divided into two vertical zones by downward plane jet as in protective-zone ventilation systems. This is done by evaluating the time-dependent distribution of infectious quanta in each zone and by solving the coupled system of differential equations based on the zonal quanta concentrations. This model introduces a novel approach by estimating the interzonal mixing factor based on previous experimental data for three types of ventilation systems: incomplete mixing ventilation, displacement ventilation, and protective zone ventilation. The modeling approach is applied to a room with one infected and one susceptible person present. The results show that using the Wells-Riley model based on the assumption of completely air mixing may considerably overestimate or underestimate the long-range airborne infection risk in rooms where air distribution is different than complete mixing, such as displacement ventilation, protected zone ventilation, warm air supplied from the ceiling, etc. Therefore, in spaces with non-uniform air distribution, a zonal modeling approach should be preferred in analytical models compared to the conventional single-zone Wells-Riley models when assessing long-range airborne transmission risk of infectious respiratory diseases.


Introduction
Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) and the disease it causes (COVID-19) pose a threat to health security worldwide and has caused more than a million deaths (Feb. 2022) since the beginning of the pandemic . SARS-CoV-2 can infect a person through exposure to respiratory fluids carrying the infectious virus. This can happen either by being exposed to infectious droplets and particles (aerosols) in the surrounding air or touching surfaces contaminated by infectious respiratory fluid [2] . The evidence for the airborne transmission of SARS-CoV-2 contained in aerosols has grown as the pandemic progressed [3][4][5][6][7] . It is now widely accepted that airborne transmission of SARS-CoV-2 may be the leading cause of super spreading events that are recognized as the pandemic's primary drivers [8][9][10] . The preventive measures in indoor environments should include decreasing the risks because of airborne transmission, among others obtained by ventilation. To this end, validated and well-performing predictive mathematical models and risk assessment create fundamental tools for understanding and planning effective strategies to minimize risks associated with airborne transmission.
Wells-Riley is the classical and widely used model to quantitatively assess airborne infection risk (Wells [11] and Riley et al. [12] ). The Wells-Riley model presents a quick and straightforward evaluation method that implicitly calculates the airborne infection risk using the concept of quantum -one quantum is defined as the number of inhaled infectious droplet nuclei or the infectious dose required to infect 63.2% of susceptible persons in an enclosed space. The Wells-Riley model has extensively been used to evaluate the airborne infection risk of respiratory diseases, such as tuberculosis [13] , measles [14] , influenza [15] , H1N1 [16] , and more recently, SARS-CoV-2 [ 17 , 18 ]. However, the original Wells-Riley model is based on the well-mixed room air and steady-state generation of airborne pathogens, which may not be the case in real environments. Consequently, the airborne infection risk could be under-or overestimated [ 19 , 20 ]. Non-steady state solutions for the Wells-Riley model have been developed in the form of differential equations by Gammaitoni and Nucci [21] and alternative equations in a study by Rudnick and Milton [22] that introduced the concept of exhaled air volume fraction, which requires measuring CO 2 concentrations. These non-steady-state solutions still adopt the well-mixed assumption that is not the principle of advanced and novel ventilation systems such as displacement ventilation as they do not create a uniform indoor air distribution [ 23 , 24 ]. In buildings with such systems, the risk of airborne transmission examined by the conventional Wells-Riley model can be led to flawed interpretations [25] .
Various modifications have been proposed to overcome the well-mixing limitation to account for the spatial distribution of infection risk when using the Wells-Riley model [26][27][28][29][30] . One of the first solutions to this problem was presented by Ko et al. [26] . They used a theoretical model expansion of the Wells-Riley model, in which the enclosed space of an airliner was divided into multiple cabins (or zones). In this way, the degree of exposure to airborne infectious pathogens could be varied from cabin to cabin and better estimated. While this zone-model approach seems operational and reasonably simple, it does not consider the type of airflow distribution. Therefore, it is not easy to extrapolate this concept to other applications than the original use in the aircraft cabin. The Wells-Riley multi-zonal approach has been extensively applied to hospital environments in several studies by Noakes et al. [27][28][29][30] . While these studies consider the role of interzonal airflow distribution on the risk of airborne infection in multizonal departments such as adjacent hospital wards and corridors, the proposed models share the following limitations i) steady-state computation of quanta concentration in each zone is applied to the Wells-Riley models ii) the assumption of the interzonal-mixing factor value was selected based on CFD simulations for specific case-scenarios. In addition, modeling advanced airflow distribution methods such as displacement or protected zone ventilation was not within the scope of these studies.
Recent studies that consider the spatial distribution of infection risk calculated by the Wells-Riley model use numerical simulations (CFD) [31][32][33][34][35][36][37] or tracer gas measurements [38] . However, the CFD approach requires comprehensive information on the room configuration and ventilation design conditions and has limitations on a time-consuming simulation process. The tracer-gas measurements, on the other hand, require both an experimental setup and measurement equipment and have lower predictive capacity.
Therefore, the objective of the present study is to extend the Wells-Riley model so that it can be used for quick and straightforward infection-risk calculations in spaces where complete air mixing is not present. The models presented in this study are based on the zonal modeling approach presented in works by Noakes et al. [27][28][29][30] . However, the novelty of the multi-zonal approach presented in this study compared to the previous works [27][28][29][30] include: (i) transient-state computation of quanta concentration is presented for each zone (ii) the interzonal mixing factor is calculated based on the ventilation efficiency values derived from literature (iii) solutions for advanced airflow distribution methods such as displacement ventilation and protective zone ventilation are presented.

Approach
The model presented in this study is loosely based on the zonal models developed by Nicas [39] and Sandberg [40] , in which the enclosed spaces are divided into zones with uniformed mixed air. In this study, the zonal modeling approach is used to demonstrate the impact of different airflow distribution methods on the infection risk in an enclosed space occupied by one infected and one susceptible person. The development of this model is carried out in the following manner ( Fig. 1 ): the enclosed space is first divided into vertical or horizontal zones depending on the position of airflow supply and exhaust according to the considered airflow distribution methods. The next step is to express the interzonal airflow mixing factor between the two zones β as a function of the ventilation efficiency ε. This allows solving the set of differential equations representing the time-dependent distribution of infectious quanta in each zone by using the interzonal mixing factor β expressed as a function of the ventilation efficiency ε. Finally, based on the calculated zonal quanta concentrations from the previous step, the zonal infection risk may be calculated according to the Wells-Riley equation.

The Wells-Riley model extension assumptions
For this study, we assumed two-zone airflow distribution methods with horizontal and vertical zoning in the enclosed space and the airflow distribution methods classified as incomplete mixing ventilation, displacement ventilation, and protected zone ventilation ( Fig. 2 ).   For each of the airflow distribution methods represented by the different ventilation systems, the following assumptions were made: (i) The only emission source of SARS-CoV-2 is from the infected individual within the room who emits SARS-CoV-2 quanta at a constant rate. (ii) There are four removal mechanisms of the infectious quanta: deposition by gravitational settling, virus inactivation by biological decay, respiratory absorption, and ventilation with recirculation; resuspension rate R is neglected in the model( R ≈ 0).
(iii) Indoor airborne transmission is associated with aerosol droplets ≤ 5 μm in the dehydrated state [41] .
(iv) Only long-range ( > 1.5 m) airborne transmission is considered [42] (v) The volumetric flow rates of air into the room from outside and out of the room to the outside are assumed to be equal and constant for the duration of the analysis, so the systems are considered balanced concerning the inward and outward flows. (vi) Infectious respiratory airborne droplets become evenly distributed throughout each zone considered after reaching steady-state conditions, so for each zone complete mixing is assumed. (vii) No virus-laden airborne particles enter the supplied air (from the outside) ( n sup = 0). (viii) There is no prior source of infectious aerosol in space.
(ix) The viral content of a saliva droplet produced by an infected person is proportional to its initial volume [43] . (x) There is no simulated sunlight indoors (ultraviolet B solar irradiance ≈ 0 W/m 2 ).

Mixing ventilation
A schematic representation of the theoretical model assessing the impact of the airflow distribution method on the quanta emission rate of SARS-CoV-2 for infection risk assessment in case of the mixing ventilation is shown in Fig. 3 .
The differential Eq. (1) represents the mass balance single-zone model for a completely mixed mechanically ventilated room: Since tracer gas has been confirmed as a suitable surrogate for the exhaled droplet nuclei < 5 μm when studying the airborne transmission in the built environment [ 44 , 45 ], we will use the Eq. (2) [42] to express the ventilation efficiency for removing droplet nuclei represented by a tracer gas concentration C from an indoor space under steady-state conditions, i.e. from the moment the quanta concentration in the occupied zone n ( t ) has reached the steady-state level: For a completely mixed condition, the ventilation efficiency is ε = 1. Under this condition, the quanta level in the exhaust equals the quanta level in the room ( n exh ( t ) = n ( t )) and by assuming that viruses in the supplied air are absent ( n sup = 0 the differential Eq. (1) can be reduced to To solve Eq. (4) in the form of a first-order differential equation dn (t) dt + n (t) · a = b, the equation may be rewritten as follows: The unique solution of quanta concentration n(t) in an indoor environment with complete mixing ventilation at time t is: where n 0 is the initial quanta concentration ( quanta To perform calculations with Eq. (5) and to predict indoor concentrations of quanta at time t , appropriate expressions for the source term S , deposition rate D , inactivation rate k , and absorption rate ζ and must be first be known.
To determine the probability of infection ( P , %) as a function of the exposure time ( t ) of susceptible people, the quanta concentration must be integrated over time through the Wells-Riley equation as follows:

Incomplete mixing ventilation
The contaminant distribution modeling approach in indoor spaces with incomplete mixing is based on two-zone models, in which the space under consideration is divided horizontally into two perfectly mixed zones, Fig. 3 . This phenomenon occurs due supplying heated air from a ceiling-mounted supply inlet [ 40 , 46 ]. Besides the inter-zonal mixing factor, the main unknown variable that needs to be determined in this modeling approach is the volume of the two zones. Previous studies have shown relatively uniform concentrations and airflow distributions established in the volume that may be identified as the occupied zone, i.e. up to around 1.8 m height from the floor [ 40 , 46 , 47 ]. But no general means of identifying the exact volumes of the two zones have been proposed. Therefore, the model for incomplete mixing ventilation in this study is also based on the assumption that the volume of an incompletely mixed room will be divided into the lower or occupied zone having a height of H = 1.8 m and the upper zone which consists of the rest of the room volume. Fig. 4 shows a schematic presentation of a two-zone exposure model for the room with an incomplete mixing loosely based on the two-zone models for the incomplete mixing used earlier in the literature [ 30 , 36 ]. In the exposure model, the space under consideration is divided horizontally into two perfectly mixed zones with uniform quanta concentrations in each zone n i ( t ) and n j ( t ): the occupied zone i reaching h occup . = 1.8 m above the floor and the rest is the unoccupied zone n j ( t ). The volume of the occupied zone can be expressed as The airflow mixing between the two zones is presented by the inter-zonal airflow rates Q i − j and Q j − i , where i-j denotes the zone i -to-zone j movement and j-i denotes the reverse air movement. By the principle of continuity of airflow in each zone, the inter-zonal airflow rates for this scenario must be where β is a dimensionless measure of the degree of mixing (mixing factor) between zones i and j in both directions. The quanta balance for the lower occupied zone i can be expressed in the following form: and after regrouping: The quanta flow rate balance for the upper zone j: And after regrouping Eq. (9) : The novelty in the zone exposure model is that the value of the mixing factor β is calculated as a function of the steady-state ventilation effectiveness, i.e. when As the outdoor air is contaminant free n sup ≈ 0 and no recirculation can be considered the ventilation efficiency can be expressed as ε = n exh (t) −n sup (t) . By dividing both sides of Eq. (10) with n i it becomes possible to derive an expression for the mixing factor β as a function of the ventilation effectiveness ε in the following manner: We used measurements (n = 25) from four studies in mixing ventilated rooms with supply located in the ceiling and exhaust in the sidewall [ 40 , 46-48 ] to determine a linear regression model expressing the relationship between the ventilation efficiency ε and the supply airflow, the room height and the and the difference between the supply and exhaust temperatures expressed as . This correlation is based on A fairly good fit ( R 2 = 0.76 was obtained ( Fig. 5 ).

Displacement ventilation
In a space with displacement ventilation, cool air is supplied at floor level so warm contaminated air is rising due to buoyancy forces and extracted at ceiling level. Theoretically, due to the thermal stratification, two zones will be formed in an indoor space: a lower, cooler, and cleaner zone and an upper warmer and contaminated zone [49] . The stratification point is located at the point where the airflow rate to the lower zone equals the airflow out from the lower zone [50] . As the heat load inside the room also creates rising convective flow, the stratification height y n (m) will depend on the amount of airflow supplied and the convective flow from the thermal load inside the considered space and can be determined from the airflow balance for zone i : The convective airflow rising from I number of human adults at height h = y n can be calculated according to the following model derived by Dokka [51] : The surface area of an average adult A per can be set at 1.8 m 2 [52] . The height of an average adult h per seated and standing can be respectively set at 1.2 m and 1.75 m. At normal room temperature conditions with normal indoor clothing and activity the convective part of the sensible heat output P per of a person can be estimated at 50 % [54] .
The neutral height y n has to be iterated in Eq. (13) until the sum of the convective airflows from the thermal load equals the supply airflow rate, i.e., h = y n Q con v = Q sup . As the stratification height across the inter-zonal cross-section will vary, we will define a share factor that defines the amount of the quanta that is emitted in the lower zone · S and the amount of the quanta emitted from the infected person in the upper zone (1 − ) · S . When y st > breathing/exhalation zone of the infected person then = 0, and vice versa: when y n < breathing/exhalation zone of the infected manikin then = 1. The quanta flow rate balance for the lower zone i: The quanta flow rate balance for the upper zone j: (17) or when expressed in the form of a first-order differential equation: Similarly, as for incomplete mixing ventilation, an expression for the displacement ventilation mixing factor as a function of the contaminant removal effectiveness ε was derived ( n sup ≈ 0): We used data (n = 30) from seven studies [53][54][55][56][57][58][59] to create a polynomial linear regression model expressing the relationship between the reported normalized concentrations A fairly strong fit ( R 2 = 0.71) was obtained ( Fig. 7 )

Protected zone ventilation
The protected zone ventilation [60] separates an indoor space into two well-mixed subzones i and j of equal volume 2 ) by using a downward plane jet as shown in Fig. 6 . We will assume that due to the uniformity of the airflow distribution from the jet plane, each zone will be supplied by an equal amount of the total airflow supply rate, i.e. Q i,sup = Q j,sup = Q sup 2 = Q 2 and that an equal amount of air will be extracted from each zone Q i,exh = Q j,exh = Q sup 2 = Q 2 . As a consequence, the inter-zonal airflow rates will be equal, i.e.
where β is the mixing factor between zones i and j . A simplified two-zone model concept of contaminant distribution for protected zone ventilation transformed into a two-zone exposure model for assessing the long-range airborne transmission risks in indoor environments is shown in Fig. 8 . The quanta flow rate balance for the polluted zone i: and after regrouping becomes using expression (21) becomes: The quanta flow rate balance for the protected zone j: For protected zone ventilation the mixing factor can be expressed as a function of the ventilation efficiency, as follows: We used an experimental study by Aganovic & Cao [61] with n = 9 data and the air change rates Q V to express the relationship between the ventilation efficiency. A linear regression model was fairly strong ( R 2 = 0.71), Fig. 9 .

Solutions for quanta concentrations for the two-zonal ventilation systems
The differential equations for the change in the quanta concentrations in zones i and j, i.e. the pair of Eqs. (8) and (10) for incomplete/imperfect mixing ventilation, Eqs. (16) and (18) for displacement ventilation, and Eqs. (22) and (24) for protected occupied zone ventilation can be written in the following forms: Where the constant coefficients A 1 , A 2 , B 1 , B 2 , C 1 and C 2 for the ventilation systems are presented in Table 1: Fig. 9. A linear regression fit model for expressing the ventilation efficiency ε for protected zone ventilation.

Table 1
Constant coefficients for solving the coupled system of differential Eqs.

Type of ventilation system
Incomplete mixing ventilation Fig. 10. Schematic representation of a simple indoor air mass-balance model in a three-zone model.
The unique solutions to this set of first-order differential Eqs. (27) and (28) : n j ( t ) = K 1 · e r 1 ·t + K 2 · e r 2 ·t + Where The coefficients K 2 and K 1 can be calculated using initial conditions n i (0) = 0 and n j (0) = 0:

Three-zonal theoretical ventilation model
While the two-zone models described in previous chapters are characterized by one single supply airflow inlet, we here introduce a three-zone ventilation model capturing the airborne risk infection dynamics in an occupied room with three inlets and one outlet. The modeling approach presented in this chapter is based on a recent experimental study on the airborne transmission risk of SARS-Cov-2 using tracer gas measurements in an enclosed space with the same layout of airflow supply inlets and exhaust outlets [62] . The three-zone model based on the study is presented in Fig. 10 .
Given that the supply airflow rates to each zone are equal Q i,sup = Q j,sup = Q k,sup , and each zone has an equal volume V i = V j = V k = V the interzonal flow rate between zones i and j can be determined by Q j − i = Q i − j + Q i,sup = β · Q + Q and the interzonal airflow rate between zones j an k can be determined by Both interzonal mixing factors β and γ were adjusted to fit the results based on the reported experimental measurements from the study [62] .
The quanta flow rate balance for zone i: and after regrouping: The quanta flow rate balance for zone j: The quanta flow rate balance for zone k: The differential equations for quanta concentrations in zones i, j, and k can be written in the following terms: The equations can be written in the following form: The unique solution to this equation given initial conditions n i ( t ) = n j ( t ) = n k ( t ) = 0 can be expressed as Then the general solution is n c ( t ) = K 1 · v 1 · e λ 1 t + K 2 · v 2 · e λ 2 t + K 3 · v 3 · e λ 3 t The particular solution is found as : (43) 2.9. The source and removal terms in the extended Wells-Riley model 2. 9

.1. The source term S
The pollutant source emission rate S is defined as the quanta emission rate of SARS-CoV-2 generated by infected persons and can be defined by [18] : The input data for c v was based on recent data [63] for the mean viral of the delta variant, thus c v = 6 . 91 · log 10 RNA ml was used as input in Eq. (4) . But this value can also be determined for any other variant and virus type for which the viral load is known and the calculations are made.
c i -conversion factor is defined as the ratio between one infectious quantum and the infectious dose expressed in viral RNA copies (quanta/RNA). Schijven et al. [64] developed a dose-response relationship for SARS-CoV-2 based on a SARS-CoV-2 isolate (hCoV-19/Netherlands/Zuid-Holland_10 0 03/2020). It was estimated that 3.13 · 10 9 RNA replicates per mL match 5.62 · 10 7 tissue culture infectious doses that infect 50% of the cells (TCID 50 ) per ml. Based on these data the fraction of RNA viral copies that corresponds to 1 plaque-forming unit (PFU), given that 1 PFU = 0.69 TCID 50 [65] conversion ratio is given by: 3 . 13 · 10 9 = 0 . 125 Given that Haas [66] recommended using the dose-response data for human coronavirus 229E (HCoV 229E) as representative for SARS-CoV-2 and that for each HCoV 229E eac hPFU has a probability of 0.055 to cause illness enables us to calculate that on average 1 0 . 125 · 0 . 055 = 1440 viral copies per infectious dose of SARS-CoV-2. Therefore, we estimated c i = 1 quanta 1440 RNA = 6 . 94 · 10 −4 quanta RNA . IR -( m 3 h ) -inhalation rate, i.e., the product of breathing ( N br ) and tidal volume ( V br ) -are both functions of the activity level of the infected subject. The inhalation rates for resting and standing averaged between males and females are equal to 0.49 and 0.54 m 3 h , respectively [67] . N i -droplet number concentration in the i th bin, particles cm 3 V i -the mean volume of a single droplet (mL) in the i th bin.
where D max and D min denote the bin's lower and upper diameter values, according to Nicas [68] . i -size bin of the droplet distribution.
The size distribution of droplet aerosols ≤ 2 μm generated during talking is determined experimentally by the works of Morawska et al. [69] and for respiratory droplets ≥ 2 μm by Chao et al. [70] . Both studies measured the size distribution of droplets for talking/voice counting at a distance of 10 mm from the participant's mouth opening. Therefore, the measured concentration of droplets represents the original size of the droplets at the mouth opening or the mass equivalent diameter of the particle D eq (m) at the temperature and RH in the respiratory tract (37 °C and RH = 99.5%). The total volume of droplets was calculated by multiplying the droplet number distribution by the mean volume corresponding to each diameter in the size distribution [71] .

Virus inactivation rate/biological decay constant k
To characterize the impact of relative humidity on the inactivation rate k , experimental data [72] on the survival time of SARS-CoV-2 in aerosols were aggregated for mean values of k (min −1 ) at two indoor relative humidity values (RH = 20% and RH = 53%) k RH = 20% = 0.0053 min −1 and at k RH = 53% = 0.0101 min −1 .

The deposition rate D
The deposition rate D of a virus-laden droplet can be expressed as follows: H person -the average height of the infected person(s), m The gravitational settling velocity of the droplets v s and m s can be determined from the following: The relationship between the RH and equilibrium droplet diameter ( D eq ) can be derived from the separate solute volume additivity (SS-VA) model for multi-component particles by Mikhailov et al. [73] . For the sake of brevity, the equations are not repeated here. All the equations can be found in a recent publication on the relationship between indoor RH and infection risk using the Wells-Riley model [71] .

Respiratory tract absorption rate, ζ
Respiratory tract absorption rate , ζ ( 1 h ) is a function of droplet diameter and tidal volume size [74] and can be calculated according to the following equation: IR is the inhalation rate of the exposed subject (which was assumed to be the inhalation rate for resting and standing averaged) at 0.52 m 3 h , and DE (%) is the total deposition number that can be approximated by following analytical expression as a function of the equivalent droplet diameter: Where b = 5.788, s = 2.574, d 0 = 1.2 and d 1 = 4.307 are constant values for an average tidal volume of 500 ml which is considered in this study [75] .

The infection risk calculation
To calculate the infection risk for the exposed person for complete, incomplete, and displacement ventilation we used n ( t ) = n i ( t ) while for protected occupied zone ventilation, we used n ( t ) = n j ( t ) in the Eq. (6) . IR is the inhalation rate of the exposed person (which was assumed to be the inhalation rate for resting and standing averaged at 0.52 m 3 h , and T is the total exposure time (h).
The following scenarios are considered: • 0.5, 2.0, 4.0, and 8.0 ACH for incomplete mixing ventilation at four different temperature differences between supply and exhaust air ( T = 0 K , 2 K , 5 K and 10 K ) • 0.5 and 2.0 ACH for displacement ventilation ( Section 3.2 ) at four different temperature differences between supply and exhaust air ( T = 0 K , 2 K , 5 K and 10 K ) and different standing and sitting positions for the infected and susceptible person • 4.0 and 6.0 ACH protective zone ventilation • 2.0 ACH for a large three-zone ventilation model.
The modeling approach is applied to a small room with one infected and one susceptible person present. The results are shown for a 40 m 2 x 3 m sized room with two occupants, one infected (RNA = 10 6 . 91 RNA ml ) and one susceptible person that is distanced by at least 1 m.

Incomplete mixing ventilation
The results in Fig. 10 compare the incomplete mixing zonal models for three different temperature differences between the supply and exhaust air ( T = 2 K,5 K, and 10 K) with the results of the complete mixing zonal model, i.e. the temperatures of the exhaust and supply air are equal ( T = 0 K). Fig. 11 shows that the temperature difference has a notable impact on the infection risk when the air is heated compared to the isothermal air supply. Although the differences between infection risks are least noticeable for a high ventilation rate of 8.0 ACH, the relative difference can still have an important effect when comparing the results for T = 0 K and T = 10 K. Increasing supply temperature to T = 10 K higher than exhaust air relatively increases infection risk up to more than 15 % for low ventilation rates (0.5 ACH) and up to 10 % for higher ventilation rates (6 ACH) after 90 min compared to complete mixing ( T = 10 K).

Displacement ventilation
For the case of displacement ventilation, the most important effect on the infection risk will be determined by the height of the neutral zone and the position of the breathing zone of both the infected and the susceptible person. Fig. 12 illustrates the infection risk when both persons are standing in a room with displacement ventilation. An increased ACH will significantly increase the neutral height position at higher supply airflow rates, resulting in lower infection risks. The impact of the temperature difference of supply and exhaust air on is more pronounced at lower ventilation rates (0.5 ACH).   At 2.0 ACH, when the neutral zone is above the breathing zones of both persons, the impact of the temperature difference is very low, while at higher ventilation rates of 6 ACH the impact of T on infection risk is almost insignificant. The importance of considering the specific type of displacement ventilation instead of complete mixing ventilation is best illustrated in Fig. 13 when the breathing zone of the standing infected person is located in the upper contaminated zone while the breathing zone of the sitting susceptible person is located in the lower clean zone. In this case, the temperature difference has a very strong impact on the infection risk, reducing it by up to more than two times for high-temperature differences between supply and exhaust air ( T = 10 K). However as the airflow rate increases, the height of the neutral zone increases. Most interestingly, when there is a temperature difference between the supply and exhaust air, the impact of increasing the ventilation rate by a factor of four (0.5 ACH → 2.0 ACH) will not considerably lower the infection risk.
In fact, at a higher temperature difference ( T = 10 K) the infection risk might be lower at lower ventilation rates (0.5 ACH) compared to the calculated infection risks at 2.0 ACH. Also, as in the previous case, when both the breathing zones of the infected and susceptible person are in the same zone that happens for higher ventilation rates, the impact of the temperature difference between the supply and exhaust air is almost negligible.

Protected zone ventilation with exhaust located in the ceiling
The case of protected zone ventilation in Fig. 14 shows that the infection risk in the protected zone with the susceptible person will be reduced compared to complete mixing conditions and that the tendency of reducing the infection spread is increasing with increased ventilation rate (0.5 ACH → 2.0 ACH → 6 ACH).
The accuracy levels were expressed in form of 95 % confidence intervals as calculated for the derived mixing factors for the two-zone models' Eq. (11) for incomplete mixing ventilation, Eq. (19) for displacement ventilation model, and Eq. (25) for protective zone ventilation. For all two-zone ventilation systems, the results showed that the range values of the calculated 95 % CI for the modeled infection risk were either entirely above the infection risk for complete mixing ventilation (as shown in Fig. 11 for incomplete mixing ventilation, Fig. 12 for displacement ventilation when the susceptible person is standing and Fig. 14 for the polluted zone for protected zone ventilation) or the 95 % CI range values were entirely below the complete mixing infection risk (as shown in Fig. 13 for displacement ventilation when the susceptible person is sitting and Fig. 14 for the protected zone for protected zone ventilation). This implies that the wide range of 95 % CI values was still predicting the infection risk consistently for the predicted regression Eqs. (11) , (19) , and (25) shown by the colored data bars in each of Figs. 11-14 .   Table 2 shows the airflow distribution method compared to completely mixed single-zone airflow distribution. The relative difference (%) of the infection risks calculated is expressed as P two−zone − P compl etel y mixing singl e zone P compl etel y mixing singl e zone after 360 min of exposure for the considered scenario. T is the temperature difference between the supply and exhaust/room temperatures. The single-zone Wells-Riley model based on a completely mixing assumption relatively overestimates the infection risk by at least 34 % for the case when the supply air is heated (i.e. incomplete mixing ventilation) regardless of the ventilation rate. The overestimation may also be higher by at least 13 % compared to displacement ventilation when the susceptible person is standing and at low ventilation rates (0.5 ACH). With higher ventilation rates the relative overestimation drops below 6 % at 2.0 ACH and becomes almost negligible at 6.0 ACH ( < 0.1 %). The underestimation of the infection risk is high for low ventilation rate (0.5 ACH) displacement ventilation when the person is standing, ranging between 37.7 and 59.7 % depending on the supply of air and temperature difference. Once again as in the case of when the susceptible person is sitting, the relative overestimation is below 6 % at 2.0 ACH and becomes almost negligible at 6.0 ACH ( < 0.1 %). For protected zone ventilation the Wells-Riley model based on completely mixing conditions may underestimate the infection risk for at least 10 % at 0.5 and 2.0 ACH, and up to more than 30 % at higher ventilation rates (6.0 ACH).

Three zonal ventilation model
The three zonal ventilation models are shown in Fig. 15 . It depicts the infection risk for three different positions of the infected person in each of the three zones. It is shown that the infection risk in the zone with both supply and exhaust (zone k) is equal for all three scenarios, regardless of the infected person's position and it is equal to the infection risk that would exist for complete mixing conditions. However, the infection risk in the other two zones with only supply inlets (i and j) are strongly influenced by the position of the infected person. Only when the infected person is in the zone with both supply and exhaust, will this type of ventilation system reduce the infection risk in the other two zones compared to complete mixing conditions. For the other two cases shown on the left and middle part of Fig. 14 , the infection risk will be higher in the or best case scenario equal to the infection risk for complete mixing ventilation conditions. Thus, when calculating the infection risk in larger spaces with a non-uniform distribution of the supply and exhaust outlets in the ceiling, using the Wells-Riley model based on the assumption of complete mixing may considerably over-or underestimate the infection risk.

Model limitations
There are several limitations of our model: 1 Our model applies only to the long-distance airborne transmission of SARS-CoV-2. Short-range is not dealt with. 2 The impact of convective flows within the space caused by thermal sources such as human thermal plumes on the flow field has been ignored. As a result, our method may not be suitable for highly occupied spaces and needs to be extended for such situations.
3 It is assumed that occupants remain where they are most of the time i.e. the occupant activities are not considered. Moving occupants can increase mixing between zones on the one hand and also increase the infection risk generating directional transportation of exhaled virus-laden aerosols. 4 We have transformed a three-dimensional problem into a two-dimensional problem. 5 The linear regression relationships developed for correlating the ventilation efficiency and indoor environment parameters (12), (20), and (26) may be subject to inaccuracy as they are compiled from studies with different experimental setups. 6 The zones have been firmly defined and mixing within a single zone may not be achieved completely. Other room sizes have not been taken into account. 7 The assumption that expelled aerosol droplets are evenly distributed in the air of the room implies that there is an immediate dilution of the expelled virus concentration. In reality, dilution will not occur instantaneously; it highly depends on the movement of the air in the room. Even in a well-mixed room, an exposed person standing directly in front of the infected person may inhale a much larger dose of airborne particles than an exposed person physically distanced at least 1 m apart. 8 The resuspension effect has been excluded as a removal mechanism in our version of the Wells-Riley model as previous studies on the airborne transmission of respiratory diseases have shown that disease transmissions could depend on the resuspension of floor dust [76] . 9 The values calculated with this model could vary significantly as a function of the activity levels of both the infected subject and the viral load in the sputum of the infected subject. 10 The infected individual was assumed to talk constantly, which may present an unrealistic overproduction of the number of respiratory droplets expelled. 11 The probability of infection calculated according to Eq. (6) is strongly dependent on the dose-response relationship c i . To allow for a more absolute comparison between infection risks from different studies, other assessment methods, which are independent of the dose-response data [ 74 , 77 ] should also be considered in future studies. 12 The validation of any Wells-Riley model for calculating airborne infection risk would require detailed input parameters from an actual outbreak and this is complex for several reasons. This method requires considerable input information from the observed outbreaks, including the type of ventilation system, ventilation rates, space volume, and exposure time of infected persons but also additional building system details such as the amount of recirculated air, filter efficiency, etc.. As there is no source providing this information the airborne risk calculation cannot be sure. Future outbreak reports must include this information for validated retrospective infection risk assessments.

Conclusions
The findings of the zonal analytical modeling of the infection risk in rooms occupied by an infected and susceptible person and the exposure time up to 6 h can be summarized in the following key points.
• Incomplete mixing • The results of our model showed that increasing the temperature difference between supply and exhaust air T increases the infection risk for incomplete mixing conditions compared to completely mixing conditions ( T = 0 K). • The single-zone Wells-Riley model may relatively underestimate the infection risk by a substantial amount compared to the two-zone incomplete mixing model. • Displacement ventilation • The relative difference to complete mixing conditions is mostly caused by the position of the neutral plane that depends on the heat load, amount of supplied air, and temperature difference between supply and exhaust air. • The single-zone Wells Riley model may considerably underestimate the infection risk when the susceptible person is standing at low ventilation rates (0.5 ACH). On contrary, when the susceptible person is sitting the conventional complete mixing model may overestimate the infection risk to a great extent. Regardless of whether the susceptible person is sitting or standing, the relative underestimation is low at 2.0 ACH and becomes almost negligible at 6.0 ACH ( < 0.1 %). • Protective zone ventilation • Protective zone ventilation decreases the infection risk in the protected zone with the susceptible person while it increases the infection risk in the polluted zone compared to completely mixing conditions. • For protected zone ventilation the Wells-Riley model based on completely mixing conditions may overestimate the infection risk to considerable extents at all ventilation rates of 0.5 to 6.0 ACH. • Three-zone model • A case study with a three-zone model shows the importance of the zonal distribution of ceiling supply inlets and exhaust outlets compared to completely mixing conditions. In this case, it is also important to consider the relative position of the infected person.
In conclusion, this study shows that using the Wells-Riley model based on the assumption of completely mixing air may overestimate the long-range airborne infection risk compared to some high-efficiency ventilation air distribution systems such as displacement ventilation and protected zone ventilation, but also underestimate the infection risk in a room heated with warm air supplied from the ceiling. Therefore, when assessing long-range airborne transmission risk of infectious respiratory diseases a zonal modeling approach should be preferred in analytical models compared to the conventional single-zone Wells-Riley models in indoor spaces where the room air and supply air is not completely mixed.

Data availability
Data will be made available on request.