Modeling turbulent transport of aerosols inside rooms using eddy diffusivity.

One major approach to modeling dispersion of pollutants inside confined spaces describes the turbulent transport of material as the product of an eddy diffusivity and the local concentration gradient. This paper examines the applicability of this eddy diffusivity/gradient model by (1) describing the conditions under which this approach is an appropriate representation of turbulent transport, and (2) re-analysis of data provided in studies that have successfully applied gradient transport to describe tracer concentrations. We find that the solutions of the mass conservation equation based on gradient transport provide adequate descriptions of concentration measurements from two studies representative of two types of sources: instantaneous and continuous release of aerosols. We then provide the rationale for the empirical success of the gradient transport model. The solutions of the gradient transport model allow us to examine the relationship between the ventilation rate and the spatial and temporal behavior of the dose of material associated with aerosol releases in a room. We conclude with the associated implications on mitigation of exposure to aerosols such as airborne virus or bacteria.

Relatively simple dispersion models based on gradient transport can be used to estimate impact of aerosol emissions inside naturally ventilated rooms.

Ventilation rates based on tracer decay experiments provide useful estimates of transport eddy diffusivities.

Because dose increases with ventilation rate at certain distances from the aerosol source, models such those described here can provide useful guidance on designing ventilation rates.

INTRODUCTION

There is increasing evidence that the Novel Coronavirus (SARS‐CoV‐2) can spread infection through aerosols generated by an infected individual through coughing, sneezing and speaking loudly in an enclosed space shielded from outdoor air motion. Models for the transport and dispersion of such aerosols inside enclosed space have used three approaches to describe the spatial and temporal distributions of aerosol concentrations. The first divides the room into several compartments each of which is treated as a well‐mixed reactor. Turbulent transfer of material between these compartments is modeled using an exchange rate parameter. The model reduces to a set of coupled ordinary differential equations whose solution provides the concentrations at the locations of the compartments. The second approach describes the turbulent flux of material as the product of an eddy diffusivity and the local concentration gradient. The third approach, which is computationally intensive, simulates turbulent transport of aerosols by tracking the motion of aerosol “particles” in a turbulent velocity field. The velocity field in the second and third approaches are commonly generated with Computational Fluid Dynamics (CFD) models, which in most cases, represent turbulent transport of momentum and energy using eddy diffusivities , in the corresponding governing equations. These approaches yield a continuous concentration field within the room.

The objective of this paper is to determine the usefulness of the gradient transport model by examining (1) the assumptions that underlie this description, and (2) the performance of gradient transport models in describing observations provided in the literature. We focus on past studies in which transport by the mean flow plays a minor role, and in which transport by turbulent velocity fluctuations is formulated with the gradient transport model. , , , ,

This paper does not examine models for instantaneous concentrations as in Cheng et al., (2014). We focus on models that estimate ensemble averaged concentrations. We do not examine transport of aerosols that have appreciable settling velocities, or emissions associated with initial momentum, such as those during coughing and sneezing, which have been studied elsewhere. , Our focus is on dispersion of aerosols that are small enough to act as passive tracers of air motion. There is evidence that transport of viruses by such aerosols makes a significant contribution to Covid‐19 infections that occur indoors.

EDDY DIFFUSIVITY AND MEASUREMENTS OF INDOOR TURBULENCE

The representation of turbulent transport with the product of concentration gradient and eddy diffusivity is based on the analogy with molecular transport, in which diffusion is modeled with molecular diffusivity. At the molecular scale, the diffusivity is based on the valid approximation that at the length scales of interest, transport is mediated through collisions between molecules. The length scale governing these collisions, referred to as the mean free path, is several orders of magnitude smaller than the length scale of the concentration distribution. In this case, one can show that the transport of mass (as well as momentum and heat) is the product of the concentration gradient and the molecular diffusivity. The magnitude of this diffusivity is governed by the product of the mean free path and the root mean square of molecular velocities.

The gradient transport concept used to describe transport by molecules is transferred to transport in turbulent flows by assuming that the equivalent diffusivity is the product of a dominant length scale of turbulent “eddies” in the flow and a velocity scale. The length scale is associated with the dimensions of the physical process that generates turbulent motion and the velocity scale is proportional to the root mean square of the turbulent velocity fluctuations. So, we can express the eddy diffusivity K as K = αul, where u is the velocity scale and l is the length scale. This expression provides an order of magnitude estimate of the eddy diffusivity, but requires a value for the constant, α, for numerical predictions.

When applied to transport of tracers in the atmospheric boundary layer, K derived from its definition, Flux/Gradient for momentum or heat is often applied to the dispersion of inert tracers. The validity of this approach is established posteriori if it yields useful predictions. For example, Nieuwstadt and Ulden (1978) showed that the eddy diffusivity for heat can be used in the mass conservation equation to estimate concentrations of a tracer (SO2) downwind of a point release in the near surface boundary layer (Prairie Grass Study, Barad, 1958). This success of the eddy diffusivity approach is rationalized by saying that eddy diffusivity works when the length scale of turbulence is much smaller than the scale governing the gradient of the tracer concentration. Or equivalently, eddy diffusivity can be used to model turbulent transport, if the turbulence length scale, l, is much smaller than the distance between the source and the receptor of interest. Although this type of argument cannot be readily defended from theoretical considerations (See Tennekes and Lumley ), the eddy diffusivity concept, with empirical modifications, yields useful results in some applications (Wyngaard ).

Here, we examine the usefulness of the gradient transport model in describing concentrations of aerosols released in enclosed spaces. We can relate the eddy diffusivity, K, to the growth of the dimensions of a puff of material, σ, released instantaneously from a point. The mass conservation equation that governs the resulting spherically symmetric concentration field is where r is the distance from the source and the eddy diffusivity, K, can depend on the time from release, t. The solution of this equation is Where Q is the instantaneous release at t=0 and C(r,t)→0 at r→∞. By substituting Equation 2 into Equation 1, we find that . This implies that if K is taken to be a property of the flow field, the spread of the instantaneous release has to follow .

There are relatively few observations of flows within enclosed spaces that can be used to infer the length and velocity scales that govern turbulent transport. Baldwin and Maynard, (1998) and Kovanen et al., present results of measurements of velocities made in a variety of indoor environments with different types and rates of ventilation. They find that the over 50% of velocities are below 10 cm/s and turbulent intensities are about 30%. Most of the measurements made by King et al., in a ventilated room are less than 10 cm/s. They found no relationship between the distributions of these velocities and characteristics of the indoor space or ventilation rate or type. 85% of the velocity measurements made by Baldwin and Maynard (1998) were less than 30 cm/s, and 50% were below 10 cm/s. Kovanen et al., conclude that the “airflow in the occupied zone fluctuates randomly”.

Measurements made by Wasiolek et al. in two mechanically ventilated rooms indicate that mean and turbulent velocities vary over a wide range from a fraction of a cm/s to a fraction of a m/s, and the indoor turbulent intensities range from 40 to 60%. Spectral analysis of the data indicated that the associated length scale of the turbulent motion is of the order of meters.

These studies provide information on the length and velocity scales that govern turbulent dispersion in enclosed spaces but are not accompanied by tracer experiments that would allow calculation of the associated eddy diffusivities. We must turn to others that infer eddy diffusivities by fitting measured tracer concentrations to the solution of the mass conservation equation formulated in terms of eddy diffusivities.

PERFORMANCE OF MODELS THAT USE EDDY DIFFUSIVITY

We examine two studies here. Cooper and Horowitz used an eddy diffusivity model to describe measurements of doses associated with an instantaneous release, while Cheng et al. used an eddy diffusivity model to interpret measurements associated with a continuous release of CO. We chose these two studies for several reasons: (1) they do not involve transport by mean flow, which allows us to focus on turbulent transport modeled with eddy diffusivity, (2) they refer to distinct types of tracer releases, one instantaneous and the other continuous, and (3) the papers provide sufficient information to allow extraction of the measurements used to evaluate the models presented in this paper.

These two studies measured concentrations integrated over time intervals, referred to as doses in the literature. Dose rather than concentration is related to adverse health effects, such as infections from airborne viruses or bacteria (Riley et al). Note that these studies, like other similar ones, did not measure the velocities that governed the tracer concentrations.

Instantaneous release: Cooper and Horowitz

In the study conducted by Cooper and Horowitz , a photographic developer component filtered for particles less than 40 µm in diameter was dropped from a height of 1 m to create a puff release. The material that was suspended was collected at distances of 1, 3, and 5 m from the release on filters connected to pumps, which collected samples over 30 min, in the time intervals, 0–1 min, 1–5 min, 5–15 min, and 15–30 min The mass of material collected on the filters combined with the pump flow rate yielded the doses of the material at the three distances.

The experiment, which yielded 15 dose samples, was conducted in a room with dimensions of . Six air inlets and three exhaust fans on the ceiling resulted in a relatively high ventilation rate, estimated to be 15 air changes per hour. Cooper and Horowitz (1986) indicate that “no measurable directional air currents were found in the areas of the experiments”. This suggests that transport of material from the source was dominated by random turbulent motion, which in principle can be modeled using the gradient transport model described next.

Consider a puff of material of mass, Q, released instantaneously at a point, and assume that it spreads by diffusion through the surroundings. The concentration when the puff spreads to a distance r from the source is proportional to Q/r. The dose delivered to a receptor at r is then this concentration multiplied by the time scale for the puff material to diffuse past the receptor. This time scale is r/K, where K is the eddy diffusivity, so that the dose becomes Q/(Kr); the transient dose delivered at r over a time t is then Q/(Kr)ϕ(tk/r) where ϕ→0 at small t and ϕ→1 at large t. This scaling analysis allows us to understand the relatively simple exact expression for the dose delivered at between t and t after release of the material, obtained by integrating the expression for the instantaneous concentration to obtain where the error function,. Equation 4 tells us that the dose falls off as 1/r; and it is delivered over a time scale of r/K.

For the sake of completeness, we present the expression for the dose delivered at r when the instantaneous release is advected past a receptor at a mean velocity of U (Seinfeld and Pandis, 1998): The time scale for advection across a room is L/U, where is L is the scale of the room. This scale is 50 times smaller than the diffusion time scale, , for K=10−2m2/s, U=10 cm/s and L=5m. So, in the presence of mean flow from source to receptor, the full dose of Q/(4πKr) is delivered quickly to the person downwind of the source.

Cooper and Horowitz integrated Equation 3 numerically to describe their measurements of dose of airborne powder measured at 1, 3, and 5 meters. Note that the expression for the dose, Equation 4, does not contain the ventilation rate explicitly. But its impact is contained in the eddy diffusivity, K, whose velocity and length scales are induced by the ventilation.

We fitted Equation 4 to the dose measured at 1 m sampled over 30 min using the eddy diffusivity, K, as the unknown parameter. The fitted value of K=0.025 m2/s yielded Q=1.34g, which represents the material that was suspended when the powder was dropped from 1 m; the suspended material is a fraction of the powder mass that was dropped (the paper suggests that about 1% of the dropped material was suspended). We accounted for “reflection” from the room's floor by placing an image source at the source height of 1 m below the floor.

Figure 1 shows the comparison between model estimates and measured doses using the single value of K=0.025 m2/s as a fitting parameter. The bias is computed as where the overbar represents an average. The fit between model estimates and corresponding measurements indicates that, in this particular case, the eddy diffusivity model provides a good description of the data. If we assume that the length scale is no larger than the distance to the closest receptor at 1 m, the corresponding velocity scale is no smaller than

FIGURE 1

Equation 1 fitted to measurements by Cooper and Horowitz using K=0.025m2/s. Solid circles in right panel correspond to measurements, and open circles to model estimates

Drivas et al. attempted to improve upon the model by including ventilation explicitly in the governing differential equation as a removal term. This cannot be justified because ventilation occurs at the boundaries of the enclosed space. They also included terms to account for “reflection” at the boundaries of the room, which is an approximate representation of the impermeable boundaries that are not affected by ventilation. Our inclusion of reflection from the room's floor through an image source yields comparable results.

Continuous release: Cheng et al.

Cheng et al. made measurements of CO concentrations averaged over half hour (essentially a dose over this period) associated with 30 min releases of CO in two rooms with natural ventilation. The point source of CO was located at the approximate center of each room, and the associated concentrations were measured with 37 (room 1) or 30 (room 2) real‐time monitors at different radial distances and angles. In each room, 16 monitors were placed at distances below 2 m: four each at 0.25, 0.5, 1, and 2 m from the source. In room 1, two monitors were placed at 3 and 4 m and one at 5 m from the source. In room 2, two monitors were placed at 2.8 m from the source and four monitors at the corners of the room at 3.56 m from the source. Measurements logged at 15 s intervals were used to derive half hour averaged concentrations, which were then combined at the different monitors to describe the average concentrations as a function of radial distance from the source.

The ventilation rates in the two rooms were measured by first releasing SF6 in the room until the concentration was relatively uniform (the paper does not provide the criterion used to define uniformity) across the room. The release was then stopped and the concentration, monitored at two locations across the room, decreased as the tracer escaped the room. The time scale, τ, governing the decrease was obtained by fitting the equation, , where is the time after the tracer release was stopped. The inverse of τ is referred to as the air exchange rate, ACH=F/V, where F the rate of airflow out of the room, and V is the volume of the room.

Cheng et al. fitted the solution of the three‐dimensional eddy diffusion equation to the measurements using the eddy diffusivity as the unknown parameter. Following Drivas et al. , they assumed that the ventilation rate affected the entire plume as a bulk, which we believe cannot be readily justified, especially when the eddy diffusivity depends on the ventilation rate. They also accounted for reflection at the boundaries of the room. Their model provided an adequate description of the measurements. The agreement between measurements and model results was better in room 1 than in room 2; this could be related to the larger volume of room 1 (158 m3) compared to that of room 2 (59 m3) where boundary effects were more important.

Here we examine the data from Cheng et al. using the solution of the diffusion equation in a spherically symmetric environment: Where r is the distance from the continuous source, q, at r=0. The solution is for constant K and t is time after release; erfc is (1‐erf), where erf is the error function.

The concentration averaged over T (the dose divided by T) is, The dose at r over a time T after release is then . The solution does not include reflection at the boundary. We account for reflection by placing an image source at source height below the room's floor.

Figure 2 shows the comparison between model estimates (Equation 8) and corresponding measurements for room 1. The eddy diffusivity is determined by fitting the measured concentration at r=0.25 m to the model estimate of the steady‐state mean concentration, q/(4πKr).

FIGURE 2

Comparison of results from Equation 8 with measurements made by Cheng et al.

The measured average concentration is normalized with the time averaged concentration corresponding to a well‐mixed room, as in Cheng et al., given by Where q is the emission rate, and V is the volume of the room.

The model explains over 95% of the variance of the measurements, and the bias varies from close to zero to 0.26, with a mean of 0.08. As expected, the fitted eddy diffusivities increase with ventilation rate. We will examine the relationship in more detail in the next section.

The model (Figure 3) does not describe the measurements in room 2 as well it does in room 1 presumably because the mechanisms that invalidate the assumption of homogeneity and isotropy of turbulence built into the model has a bigger impact in the smaller room 2. Also, the higher ventilation rates of Experiments 5 and 6, where model bias is large, might indicate the presence of organized flow that is not accounted for in the turbulent diffusion model. “Reflection” at the boundaries is not likely to be significant because the inferred eddy diffusivities, shown in Table 1, imply that the time scale for diffusion to be affected by the wall at 3.9 m is over an hour.

FIGURE 3

Comparison of results from Equation 8 with measurements made by Cheng et al.25 (2011)

TABLE 1

Fitted Eddy diffusivities and corresponding ventilation rates for Room 2

Expt	K, m2/s	ACH, h−1
1	3.75E−03	0.19
2	4.71E−03	0.37
3	2.66E−03	0.41
4	2.60E−03	0.59
5	7.21E−03	2.08
6	1.43E−02	5.4

LENGTH AND VELOCITY SCALES ASSOCIATED WITH INDOOR EDDY DIFFUSIVITIES

The eddy diffusivities, presented in Tables 1 and 2, range from 2.3×10−3m2/s for a ventilation rate of about 0.2 ACH h −1 to 1.4×10−2m2/s for a ventilation rate of 5.4 ACH h −1. If we assume that the turbulent diffusion is the primary transport mechanism at the low ventilation rates of room 1, the inferred eddy diffusivity of 5×10−3m2/s translates into a turbulence velocity scale of 10 cm/s if we assume that the length scale is about 1/4th of the distance of 0.25 m to the monitor closest to the release location.

TABLE 2

Fitted Eddy diffusivities and corresponding ventilation rates for Room 1

Expt	K, m2/s	ACH, h−1
1	2.31E−03	0.17
2	4.07E−03	0.51
3	4.82E−03	0.57
4	4.74E−03	0.78
5	7.21E−03	1.25

Cheng et al. (2011) show that . Figure 4 presents the data in Tables 1 and 2 to illustrate this relationship. This is consistent with the idea that ventilation induces turbulent diffusion across the room. As described earlier, the ventilation rate is computed from the SF6 decay rate, which suggests that the rate is essentially the rate at which SF6 diffuses away from monitor. Then, the inverse of the ventilation rate, ACH, is the time scale for diffusion across the room , where K is the effective eddy diffusivity, and L is a length scale relevant to the room, taken to be , where V is the volume of the room.

FIGURE 4

Relationship between Ventilation Rate and Fitted Eddy Diffusivity in Experiments described in Cheng et al. . The parallel lines next to the one‐to‐one line delineate a factor of two interval

Assuming that , we can estimate the magnitudes of the associated length and velocity scales by applying the turbulent kinetic energy (TKE, eg, Stull, ) equation, as suggested by Karlsson et al. . The TKE equation is given by Where u 2/2 is the TKE, Mech and Buoy are the mechanical and buoyant TKE production terms, PT is the pressure transport term, and Diss is the TKE dissipation rate.

Assume that pressure transport and buoyancy production are small compared to shear production of turbulence, and that all the mechanical energy, , that enters the room is converted into turbulent kinetic energy. Then the balance that governs the turbulent kinetic energy, in a room is given by. where , and are the variances in the three coordinate directions.

In Equation 11, u is velocity of the air entering the room, F is the associated volume flow rate and V is the volume of the room. The first term on the right is the dissipation term and second term is the outflow of turbulent kinetic energy. Now τ=V/F, and taking , Equation 11 yields the following expressions for u and l, We estimate C in Equation 11 using turbulent properties of the shear dominated surface atmospheric boundary layer in which the turbulence velocity scale is the surface friction velocity, u *, and the length scale is the height above ground, z. Taking u=2u*, K=ku *z, (Stull, ) we estimate c=0.036 where k=0.4 is the von‐Karman constant

Using Equation 12, we estimate values of u and l to examine their consistency with the application of eddy diffusivity to model turbulent transport inside a room. Figure 5 shows the variation of the length and velocity scales with ACH for L=4m, and two values of the velocity of air entering a room. We see that the inferred values of eddy diffusivities for ACH are compatible with turbulence length scales of less than 1 m, a result that provides tentative support for modeling turbulent transport using an eddy diffusivity and a concentration gradient.

FIGURE 5

Variation of turbulence velocity and length scales with ACH in a room. Blue line refers to length scale, and red line to velocity scale

The observed relationship suggests that the ventilation rate inferred from the decay rate of a tracer released in a room is simply an inference of the rate at which the tracer diffuses away from the point of release once it is stopped. Thus, an increase in ventilation rate does not always lead to a reduction in dose as can be seen from Equation 4 re‐expressed as the dose from an instantaneous release at distance r from the source over an exposure time, T, from release:

Increasing the ventilation rate increases eddy diffusivity which increases dispersion and thus decreases concentration: the term outside the parenthesis on the right‐hand side of the equation. But at the same time, it increases the fraction of the maximum dose delivered to the receptor: the error function complement, erfc, term on the right‐hand side. Figure 6 illustrates the effect of increasing ACH on dose normalized by the material released over time T in a room with a dimension of L=4 m. The dose is multiplied by a nominal breathing rate of 10−4m3/s to simulate the fraction of the total released mass inhaled by a person. We see that at both r=2 m and 4 m, the effects of the two opposing aspects of the ventilation rate on dose. At r=4 m from the source, the dose increases first with ACH up to before the effects of increasing dispersion become evident. At large enough ACH, which depends on room size and exposure time, the dose becomes inversely proportional to ACH.

FIGURE 6

The variation of dose from an instantaneous with ventilation rate. The room dimension L=4m, and exposure times T=0.5 and 1 hr. The breathing rate is taken to be 10‐4m3/s

CONCLUSIONS AND IMPLICATIONS FOR MITIGATION

We have conducted an examination of the data available in the literature to examine the applicability of gradient transport/eddy diffusivity model to turbulent transport of emissions within a room. We have fitted analytical models based on gradient transport to measurements of tracer concentrations made in two studies, one with an instantaneous release and another with a continuous release. The models provide good descriptions of the data, which suggest that the gradient transport model of turbulent transport is applicable to the indoor environment when the ventilation rates are below 5 ACH h −1; higher rates are more likely to induce organized mean flow. The analytical models show that the dose delivered by emissions from a point source decreases as from the source, where K is the eddy diffusivity, and r is the radial distance from the source.

The preceding dispersion models tell us about the variation of the mean doses as a function of distance from the source but cannot be used to estimate a “safe” distance, if there is one, until we have an idea of the magnitude of the dose that can result in infection. However, these equations provide information on reducing risk of infection. Avoiding environments with potential high viral release, Q, or an emission rate, q, is the most obvious way of reducing risk: avoid places, such as bars, where people are likely to speak loudly. A high ventilation rate decreases risk in two ways. First, it removes air contaminated with emitted material from a room. It also increases turbulent velocities, which disperse emitted puffs and thus decrease the concentrations of material that they transport.

Equation (4), (8) 8 tell us that in an environment in which the mean flow is not organized, such as in a room that is naturally ventilated, the mean doses vary as1/r, where r is the distance from the source. Although this reduction in dose is relatively rapid spatially, the mixing occurs in all directions relative to the source. The emitted material is mixed throughout the room over a diffusion time scale of the order of . Taking and K=0.01m2/s, which is large for a naturally ventilated room, this time scale works out to be 40 min. So, staying a room with a source over a period of an hour is likely to bring everyone in the room in contact with the emissions. So, in a naturally ventilated room, time spent in the room is likely to be as important as distance from the source as a measure of risk of infection. The outbreak of Covid‐19 infections during a 2.5 h choir practice in a poorly ventilated church in Skagit County, Washington in March, 2020 supports this idea.

When the ventilation creates organized mean flow, Equation 5 tells us that the dose also decreases as 1/r, but the maximum dose is delivered quickly to the person in the direction of the induced large‐scale flow. So, in such a situation both distance and direction govern risk of infection. Time spent is less important if the “receptor” is located downwind of the source as first pointed out by Lu et al., who traced the infections that occurred during January to February 2020, among 10 people belonging to 3 families to their presence in a restaurant. The air conditioner in the restaurant induced strong flow in which the source of infection was seated at a table that was upwind of the tables that were affected by the infections; other diners and the serving staff were not infected.

Equation 5 also tells us that the dose received from an emitted puff downwind of a source at r decreases with increases in the mean velocity, U, and the eddy diffusivity K=ul. In an indoor environment, both U and u are of the order of cm/s, while outdoors, it is at least 10 times larger, u is about 10 cm/s, and U is several m/s. Thus, the risk of infection is much smaller outdoors than it is indoors. But even outdoors, it is important to avoid being directly downwind of the source of infection where the full dose of Q/(4πKr) is delivered over a time σ/U, where σ is the size of the puff.

One way of reducing risk of infection is to wear masks, which might dilute the aerosol concentration at the source by increasing the size of the emitted puff, which also has the effect of allowing the smallest scales of turbulence to distort and disperse the puff. The mask also prevents the emission of large droplets from inoculating a person close by. There is some evidence to indicate that a mask might also reduce the risk to the wearer (Gandhi and Rutherford 2020).

The results from the models presented here suggest that ventilation throughout the room is not always beneficial and might make matters worse by spreading emissions. On the other hand, local exhaust ventilation might reduce the risk. This type of ventilation, which has been suggested to reduce the risk of disease transmission , can be designed to remove pollutants close to the source. In a workplace or a bar, this would involve bringing in fresh air at each table close to the ground while exhausting air at roof level to take advantage of the convection induced by the warmer temperatures close to the ground. It is clear from the results presented here that mitigation is best achieved through ventilation systems that reduce horizontal transport of emissions and rely on vertical transport of air.

CONFLICT OF INTEREST STATEMENT

None reported.

PEER REVIEW

The peer review history for this article is available at https://publons.com/publon/10.1111/ina.12901.