2023 Volume 48 Issue 4 Pages 187-201
Fugacity models are used widely to predict the time-dependent behaviors of chemicals in environments containing several media (e.g., air, sediment, soil, and water). However, these fugacity models work on the assumption that the concentration of a chemical in each medium is uniform, so they cannot describe the spatial distribution of the chemical. We developed a new fugacity model, termed InPestCFD, incorporating computational fluid dynamics to describe both the time-dependent distribution and the spatial distribution of a chemical in a medium. InPestCFD was used to calculate the behavior of an insecticide released from an aerosol canister in a room. Indoor airflow and aerosol particle behavior were calculated via computational fluid dynamics and using a Lagrangian dispersion model. Transport of the insecticide among media (aerosol particles, air, ceiling, floor, and walls) was calculated using the fugacity model. The time-dependent distributions and spatial distributions of the insecticide in the media agreed well with real measurements.
When a chemical is released into the environment, the chemical will become distributed between the media (e.g., air, sediment, soil, and water) in the environment, and humans can be exposed to the chemical through contact with these media.1) When assessing the risks a chemical poses to humans, the level of exposure to the chemical is compared with the tolerable level of exposure to the chemical determined using toxicological data. If the exposure level is lower than the tolerable level, it is concluded that there is no concern for human health.2,3) It is therefore essential to understand the behavior of the chemical in the environment. This is usually achieved by determining the concentrations of the chemical in relevant environmental media. However, it is practically impossible to achieve a detailed understanding of the behavior of a chemical in the environment by only making such measurements because of the enormous resources required and the likelihood of disturbing the environment. Many models that predict the behaviors of chemicals in the environment have therefore been developed.4–6) The fugacity model, developed by Mackay,7,8) is widely used to predict the concentrations of chemicals in environmental media from the physicochemical properties of the chemicals.9–12)
A fugacity model13) is incorporated into the EPI Suite developed by the US Environmental Protection Agency and Syracuse Research Corporation to predict the distribution of a chemical between air, sediment, soil, and water after the chemical has been released into the environment. Fugacity models take biota into consideration, so can be used to assess bioaccumulation of chemicals discharged into the environment14) and accumulation of chemicals in humans through the food chain.15) Fugacity models have been used to assess spatial and temporal variations in chemical concentrations in indoor air and on ceilings, floors, and walls after insecticides have been released from aerosol canisters,16) electric vaporizers,17) and hand wand sprayers.18) The US Environmental Protection Agency has used a concept similar to a fugacity model to predict concentrations of pesticides in indoor media.19) Shin et al.20) used a fugacity model to predict the concentrations of various chemicals in carpets and wall materials. Li et al.12) used a fugacity model at the screening level to assess the risks to human health posed by a chemical volatilizing in a room by predicting the concentrations of the chemical in human organs after exposure through dermal contact, ingestion, and inhalation.
In a fugacity model, the concentration C (mol/m3) of a chemical in a medium is determined by multiplying the fugacity f (Pa) of the chemical by the fugacity capacity Z [mol/(m3 Pa)] of the medium. Fugacity is the tendency of a chemical to escape from one medium and enter another. The fugacity capacity is the ability of a medium to hold onto a chemical, which is usually determined from the physicochemical properties of the chemical (octanol–water partition coefficient, vapor pressure, and water solubility). There are four fugacity model levels, with the mass balance equations becoming more complex moving up a level.21,22) The level I model is used for a closed system assuming that all compartments are at equilibrium and there are no chemical inputs or losses. The level II model is used for an open system assuming that all compartments are at equilibrium. The level II model takes continual emissions and advection into consideration. The level III model is used for an open system at steady state without assuming that all of the compartments are at equilibrium. The level IV model is used for an open system not at steady state. The level IV model describes time-dependent changes in fugacity using differential equations. The level IV fugacity model allows transfer of a chemical among media and the time-dependent concentrations of the chemical in the media to be predicted. However, these models do not indicate areas in a medium with high or low concentrations of the chemical of interest.
People and various organisms move around the environment. It is important when performing a risk assessment for a chemical in the environment to determine the concentrations of the chemical of interest in different parts of each medium in the environment. ChemFate is a multi-media dynamic fate and transport model for predicting the environmental fate of organic chemicals by using a fugacity model.23) Iwami et al.24) described spatial and temporal variations in the concentrations of dioxin and dioxin related compounds on soil surfaces using a fugacity model to describe a system with dioxin and dioxin related compounds discharged from an incineration plant stack and deposited on the soil. A Gaussian puff model was used to describe diffusion in the atmosphere. Pivato et al.25) predicted the concentrations of pesticides in air using a Gaussian plume model and fugacity model to describe the behaviors of pesticides volatilized from a vineyard. These models contribute to performing a comprehensive exposure assessment by predicting concentrations of various environmental media. However, each model has the following problems. ChemFate assumes that the concentration of the chemical of interest in a medium is uniform. It is therefore not possible to predict the spatial distribution of a chemical in a medium. The puff model is a statistical model for describing diffusion of a chemical in concentric areas around a source of emissions, so does not include the effects of advection. The plume model can take the wind direction into account but cannot predict short-term temporal changes in chemical concentrations.
Computational fluid dynamics (CFD) is a simulation method in which equations for fluid flows are solved numerically.26) CFD is currently used in various fields, including aircraft design, aerospace engineering, and weather forecasting. Accurate boundary conditions need to be set to allow precise CFD calculations to be performed. There remain problems parameterizing physical processes such as interactions between airflow, heat, humidity, and radiation.27,28) The behaviors of chemicals in rooms have been analyzed by CFD. However, the analyses only considered decreases in chemical concentrations caused by ventilation29) or both ventilation and deposition onto surfaces.30) Volatilization of deposited chemicals from surfaces to the air or decomposition of chemicals in the air and other media in the room should be considered.
In this study, a model of chemical transport in a room (termed InPestCFD) was developed by including CFD in a level IV fugacity model. The CFD part described airflow (which is not considered in existing fugacity models) in the indoor air compartment. The fugacity (which is not used in existing CFD models) modeled chemical transport between different media. This is the first time that CFD has been incorporated into a fugacity model. Combining CFD and fugacity models allows calculations to be performed that would be difficult or impossible using only one of the models. A fugacity model cannot indicate which areas in a medium will contain high or low concentrations of the chemical of interest. A CFD model cannot accurately describe transport of the chemical of interest between different media without precise boundary conditions being defined. We combined a fugacity model and CFD model to describe both time-dependent behavior of the chemical of interest in various media and spatial distribution of the chemical in each medium. The new model InPestCFD was used to predict the behavior of an insecticide in a room after the insecticide had been released as an aerosol from a canister. The predictions of InPestCFD were compared with the results of an existing model and measurements of the behavior of an insecticide in a room. This model enables us to conduct more-detailed exposure assessments by predicting not only the temporal transfer of the chemical between the different media but also the spatial distribution of the chemical in each medium.
We used the results of two wind tunnel experiments to verify indoor airflow calculated using the CFD model.31) In the first experiment, the vertical wind velocity and turbulent kinetic energy profiles were measured 13 m from the air inlet along the center-line of the test section. Boundary layer flow was produced in a wind tunnel at Niigata Institute of Technology (height 1.8 m, width 1.8 m, length 16 m), which did not contain any obstacles. The x-axis was defined as the mainstream direction perpendicular to the y-axis (the spanwise direction) and the z-axis (the vertical direction). The measurement position was x=0 (Fig. 1a). Figure 1b shows the measured velocity and turbulent kinetic energy. In the second experiment, a rectangular acrylic parallelepiped obstacle (height 0.2 m, width 0.2 m, length 0.05 m) was mounted so that the center of the base was at x=0 and y=0 (Fig. 1c). The average wind speed was measured along the vertical cross section at x=−0.075, −0.025, 0, 0.025, 0.05, 0.1, 0.2, 0.3, 0.4, and 0.55 m along the centerline through the obstacle (Fig. 1d).
We used the results of experiments performed at the Sumitomo Chemical Co., Ltd. to verify the indoor behavior of an insecticide predicted using InPestCFD. The room was a typical size for Japan and had a wooden floor coated with polyurethane resin (Fig. 2). The walls and ceiling were covered with polyvinyl chloride wallpaper. The ventilation rate was controlled by a fan on the outlet wall and an air inlet with vertical slits on the opposite wall to give an exchange rate of 1.5 hr−1. The temperature was kept at 25°C by an air conditioner, which was in a large room surrounding the sprayed room. The air leaving through the outlet was directed out of the large room and not allowed to return to the sprayed room or large room. A space aerosol canister (300 mL) containing 0.45 g of d-tetramethrin (Neo-Pynamin Forte, C19H25NO4, (1,3,4,5,6,7-hexahydro-1,3-dioxo-2H-isoindol-2-yl)methyl (1R)-cis-trans-2,2-dimethyl-3-(2-methyl-1-propenyl)cyclopropanecarboxylate) as an active substance was used. Four canisters of the aerosol were sprayed slightly upward at the same time through four small windows at a rate of 0.82 mL/s for 2.5 sec through each window. Thus, the air enters through the inlet, circulates in the room, and exits through the outlet. The sprayed aerosol particles circulated in the flowing air, settled as a result of gravity, and became smaller because of evaporation of the solvent in the particles. The insecticide in the particles adhered to the room materials, re-volatilized to the air, and penetrated the room materials. The behavior of the insecticide in the room after spraying was investigated by performing three separate experiments.
The parameters and variables used are shown in Supplemental Table S1. The equations used were similar to the equations used by Matoba et al.16) and are shown in Supplemental methodology.
2.1. CFD modelAirflow in the room was described using a CFD program constructed using Eqs. 1–8. The equations governing the behavior of an incompressible fluid are shown in Eqs. 1 and 2.26)
 | (1) |
 | (2) |
where um (m/s) is a velocity vector component (m=1, 2, and 3 indicate the x-, y-, and z-directions, respectively). ρ is the air density (kg/m3). ν is the molecular dynamic viscosity (m2/s). νt is the eddy viscosity (m2/s). 〈 〉 indicates the ensemble average. The summation convention is assumed for m and n. The standard k−ε turbulence model was described by using Eqs. 3–8.
 | (3) |
 | (4) |
 | (5) |
 | (6) |
 | (7) |
 | (8) |
where Cμ=0.09, Cε1=1.44, Cε2=1.92, and σε=1.3.33) kt is the turbulent kinetic energy (m2/s2). ε is the energy dissipation rate (m2/s3). Pk is a production term. Dk and Dε are diffusion terms for kt and ε, respectively.
2.2. Aerosol dynamicsThe velocity vector for an aerosol particle vp (m/s) consisted of an average velocity vector for indoor airflow at the position of the particle u̅ (m/s), a turbulent velocity vector vt (m/s), and a velocity vector derived from an external force (i.e., gravity) F (m/s), as shown in Eq. 9.
 | (9) |
where σu is the standard deviation of each indoor air velocity component. ζ is a random number with a Gaussian distribution (mean 0 and variance 1). σu was defined as 
because the turbulent kinetic energy was defined as 1/2 of the sum of the squares of the standard deviations of each component. vt(0) at time=0 was derived from the velocity at which the aerosol was released from the canister. vt(t) after t=0 was calculated by using the random walk model.34) a and b are the contributions of turbulent diffusion calculated by using the Markov chain model in Eq. 10.34)
 | (10) |
where tLxi (s) is the Lagrangian integral scale for the turbulent velocity of each component. The z-direction component of the velocity vector F3 (m/s) derived from the external force was represented by gravity and viscous forces, as in Eqs. S1 and S2 (see Supplemental methodology 1).35)
Each aerosol particle becomes smaller over time because of evaporation of the main solvent in the aerosol. The time-dependent diameter of a particle dp (m) was determined using Eq. S3 (see Supplemental methodology 1).35) Pd is the vapor pressure of the solvent (Pa). P∞ is the vapor pressure of the solvent (Pa) far from the particle. We assumed that P∞ was zero.16)
2.3. Fugacity modelThe behavior of the active substance sprayed into the room was described using a program based on the fugacity model using Eqs. 11–15. We used a level IV fugacity model to describe time-dependent fugacity. The equations were based on equations used by Matoba et al.16) However, each medium (air, ceiling, floor, and walls) was divided into grids to allow numerical values obtained by performing CFD calculations to be used.
The fugacity of the active substance in an aerosol particle l in cell (i,j,k) fp,l (Pa) changed with time and was calculated using Eq. 11 (see Supplemental methodology 2).
 | (11) |
where fair,i,j,k (Pa) is the fugacity of the substance in the air of each cell. Dp,l,air,i,j,k [mol/(s Pa)] is a parameter for transfer of the substance in each cell between the particle and air. Dp,l,air,i,j,k [mol/(s Pa)] was given by using Eq. S7 using kair,i,j,k (m/s), which is the velocity of air in cell (i,j,k). Both fair,i,j,k (Pa) and kair,i,j,k (m/s) were determined for each cell by performing CFD calculations. In contrast, Matoba et al. defined fair (Pa) as the fugacity of the substance in the air in the entire room and assumed that kair (m/s) was the sum of the ventilation rate and gravitational settling velocity.16) The first term on the right-hand side of Eq. 11 is the fugacity increase caused by the particle size decrease. The second term represents transfer of the active substance between the particle and air. The third term represents degradation of the substance in the particle.
The fugacity of the substance in an air cell fair,i,j,k (Pa) will change over time and is described by Eq. 12 (see Supplemental methodology 3).
 | (12) |
where Sct is the turbulent Schmidt number (0.7).36) S (mol) is the amount of substance transferred between the air and particles, and/or between the air and room materials. In Eq. 13, S is the amount of substance transferred between the air and particles.
 | (13) |
where N is the number of particles in an air cell not adjacent to any room materials (floor, ceiling, or walls). If a cell is adjacent to a room material, the amount of substance transferred (mol) between the air and room material is incorporated into S. For example, air cell (1,j,k) is adjacent to the inlet wall surface; thus, S is calculated with Eq. 14.
 | (14) |
where Dm,air,1,j,k is a parameter for transfer of the substance between the air and room materials [mol/(s Pa)]. Dm,air,1,j,k is given by Eq. S11.
The fugacity of the substance in a room material cell will change over time, as described by Eq. 15 (see Supplemental methodology 4) for room material cell (0,j,k) and air cell (1,j,k) adjacent to the room material.
 | (15) |
where Ndep,1,j,k is the total number of particles deposited on room material cell (0,j,k) during period ∆t (s).
3. Computation and boundary conditions3.1. Wind tunnel experiment simulationFirst, we calculated airflow in a wind tunnel experiment to validate the CFD model. Airflow around an obstacle in a wind tunnel experiment was calculated using a staggered grid for a 2 m×2 m×1.2 m calculation domain. The domain was divided into an 80×80×48 grid (i.e., 307,200 sections) with a minimum mesh size ∆x=∆y=∆z=0.025 m (the minimum spacing between measurement points). The TVD method26) was used to solve the advection term in Eq. 1. The velocity and pressure fields were solved using the SMAC method.37) The code was written in Fortran 90 using the Microsoft Visual Studio 2017 user interface and compiled using the Intel® onlap Base Toolkit. Calculations were performed using the Windows Server 2019 standard 64-bit operating system using a two processor Intel® Xeon Gold 6234 @3.30 GHz system with 128 GB installed memory (RAM).
The velocity and turbulent kinetic energy measured in the wind tunnel experiment shown in Fig. 1b were applied to the inflow boundary condition of the model. A zero-pressure gradient was applied to the outflow boundary. The inflow boundary profile found in the wind tunnel experiment was described by the power law with an exponent of 0.25. Slip wall boundary conditions were applied to the side wall and top boundaries. No-slip conditions with a wall function were applied to the floor boundary.38)
3.2. Simulating the indoor behavior of the insecticideFlow in the room shown in Fig. 2 was calculated using the same program as for the wind tunnel experiment. The fugacity model was added to the program code to predict the behavior of the insecticide for 24 hr after the insecticide had been sprayed. The mesh size ∆x (=∆y=∆z) was 0.05 m. The room was divided into 186,624 sections with 72, 54, and 48 lines in the x-, y-, and z-directions, respectively. Calculations were performed until airflow appeared to be at steady state. Two types of tests were performed to determine whether indoor airflow became steady. In one test, the maximum and average differences in wind velocity during period ∆t in all the cells in the domain were assessed. In the other test, the time taken after the insecticide had been sprayed until the number of particles in the room halved was determined. The results are shown in Results section 1.2.
Huge computational resources were required to precisely describe the airflow and the substance behavior because very fine time steps ∆t were required to calculate the airflow and the behaviors of very many aerosol particles. Three measures were implemented to decrease the computational time. First, indoor airflow calculations were performed separately from the particle and fugacity calculations. Steady airflow was used to calculate the aerosol particle and substance behaviors in the room. Second, the number of aerosol particles was up to 108. Therefore, 1,000 or 10,000 particles were represented by one hypothetical particle. It appeared that this representation did not affect the particle behavior attributed to gravitational settling and airflow for the hypothetical particles because the diameter and mass of individual hypothetical particles were the same as for actual sprayed particles. In addition, dynamic memory was allocated to variable numbers for the particle parameters (diameter, fugacity, position in the air, and velocity). Third, a very short time step was required when the fugacity changed dramatically, e.g., when the particle diameter decreased or the substance began to permeate the room materials. At other times it was not necessary to set such a short time step to calculate the fugacity. A variable time step ∆t was therefore used and changed according to the difference in fugacity before and after the previous ∆t.
The input parameters used in InPestCFD are shown in Supplemental Table S1. The inlet airflow velocity uin (m/s) in the inlet cell was calculated from the ventilation rate ACH (h−1), room volume Vroom (m3), and inlet size used in the experiment (Eq. 16). Inlet size was set to 0.3 m×0.2 m because the mesh size ∆x was 0.05 m in this calculation.
 | (16) |
The turbulent kinetic energy kt (m2/s2) in the inlet cell was calculated using two methods. One method gave kt=0.005 m2/s2, calculated from the airflow in a clean room using the standard k-ε model.39) The other method gave kt=0.435uin2, calculated for slit flow using a direct numerical simulation40) in which uin (m/s) was the inlet velocity because vertical slits were used in the experiment room inlet.
No-slip conditions were used for all of the wall surfaces. The 1/7 power law for the wall function was used for the velocities of cells adjacent to the walls.39) The diffusion coefficient Dm (m2/s) of the substance in the room materials was calculated using the empirical equation Dm=10−9D. The diffusion coefficient D (m2/s) for the substance in air16) was determined using an equation published by Wilke and Lee41) and the physicochemical properties of d-tetramethrin (taken from the Hazardous Substance Data Bank42)). The oil content ratio rm for the room material was 0.04 for the floor coated with polyurethane resin (based on the urethane resin content) and 0.3 for the ceiling and walls covered with polyvinyl chloride wallpaper (based on the polyvinyl chloride content). The aerosol particles were divided into three groups, particles with diameters of 10 µm (10% of the particles), 30 µm (60% of the particles), and 50 µm (30% of the particles), as described in a previous publication by one of the co-authors of the present study.43) The total amount of the substance was calculated by multiplying the product application rate AR (mL/s) by the product application period tappl (s) and the concentration of the substance in the product Cproduct (g/mL).
The wind velocities for the system with a rectangular obstacle mounted on the wind tunnel floor were calculated using the CFD program. Sampling was performed at 10 points in front of and behind the obstacle. The results were compared with the measurements made during the wind tunnel experiment. The data for four points (x=−0.075, 0.05, 0.1, and 0.55 m) are shown in Fig. 3. The calculated and measured values for all of the points agreed well (data not shown). At x=−0.075 upstream of the obstacle, the approaching airflow was strongly deformed by the rectangular obstacle and backward flow occurred near the floor surface. Downstream of the obstacle, negative values were found at x=0.05 and 0.1 because of circulating flow behind the obstacle. However, the value at x=0.55 was not negative. The model reproduced these phenomena well, and the calculated and measured velocities for each point and height agreed well.
The airflows found 300 and 3600 sec after the calculations started are shown in Fig. 4. At 300 sec, air entered at the inlet, circulated throughout the room, and exited through the outlet. However, this was not the case at 3600 sec, and stirred airflow was found. The airflow did not change markedly after 3600 sec (data not shown). This indicated that indoor airflow was not at steady state at 300 sec. The maximum velocity differences for time step ∆t at 300, 1200, 1800, and 3600 sec were 1.6×10−5, 1.1×10−5, 8.3×10−6, and 1.0×10−5 m/s, respectively. The maximum velocity difference was slightly larger at 300 sec than later but changed relatively little after 1200 sec because the maximum velocity difference was derived from the velocity near the air outlet. In contrast, the average velocity differences for ∆t at 300, 1200, 1800, and 3600 sec were 8.8×10−8, 4.5×10−8, 4.2×10−8 and 2.7×10−8 m/s, respectively, i.e., the average velocity difference decreased over time (steady state was approached over time). When aerosol particles were sprayed into the air at 300, 1200, 1800, and 3600 sec, the time required for the number of particles with initial diameters of 10 µm (which would not have fallen to the ground within 1 hr) to halve were 23, 19, 18, and 18 min, respectively. The number of particles in the air in the room decreased over time because some particles exited the room through the outlet. Air circulation at 300 sec caused stagnation of particles in the room. The time taken for half of the particles to be lost was therefore longer when the aerosol was added at 300 sec than when the aerosol was added later. The particle half-life decreased over time and was 18 min when the aerosol was added at 2,700 sec (data not shown). The results indicated that airflow at 3,600 sec was effectively at steady state. Therefore, the parameters for airflow at 3,600 sec were used in subsequent calculations to describe insecticide behavior in the room.
Matoba et al.43) developed a model (later called InPest) based only on a fugacity model. Calculations were performed using InPest for the same experiment as described above. The time-dependent fugacity of the active substance in the media within 24 hr determined using InPest and InPestCFD are shown in Fig. 5. The fugacity values for the substance in small, medium, and large aerosol particles calculated using both InPest and InPestCFD increased with time because of evaporation of the solvent in the particles increasing the concentration of the substance in the particles. Both InPest and InPestCFD indicated that small particles were still present after 24 hr. However, InPest indicated that all of the medium particles had fallen to the floor after 5,050 sec (1.4 hr) and that the large particles had disappeared after 71 sec (1.2 min) (Fig. 5a). InPestCFD indicated that the medium and large particles remained present until 18,360 sec (5.1 hr) and 7,960 sec (2.2 hr), respectively (Fig. 5b). InPest only takes gravitational settling of particles to the floor into consideration, but InPestCFD also takes into consideration particles riding on the flowing air and diffusion using the random walk model. Some particles were therefore predicted to stay suspended in the air longer by InPestCFD than by InPest.
The fugacities of the substance in the air and room materials were described using one curve (Fig. 5a) because InPest assumes that the concentration of the substance is uniform in each medium. However, fugacity was calculated for each cell in InPestCFD, so multiple fugacity values were calculated for each medium (Fig. 5b). The fugacity in air calculated using InPest increased gently after spraying as the substance was transferred from the particles to the air. The fugacity in air increased sharply after 35 sec as large particles settled to the floor and the substance on the room materials re-volatilized to the air (Fig. 5a). In contrast, InPestCFD predicted that the fugacity in air increased gently rather than sharply until 600 sec (0.17 hr) as the particles gradually adhered to the room materials because of gravity, airflow, and diffusion (Fig. 5b). The fugacity in air predicted by InPest and InPestCFD then decreased because of ventilation and degradation .
The fugacities of the substance in the ceiling, floor, and wall materials predicted by InPest increased abruptly at 35 sec (Fig. 5a) because InPest calculated the amount of large particles falling to the floor and distributed that to the ceiling, floor, and walls according to the proportions determined from room material measurements. However, InPestCFD predicted that the fugacity for the floor increased gradually because gravity, airflow, and diffusion of the large particles were all taken into consideration (Fig. 5b). The fugacities for the ceiling and walls increased sharply within 10 sec of spraying because InPestCFD contained the assumption that the sprayed particles scattered and adhered to the ceiling and walls. The fugacities of the room materials predicted by InPest and InPestCFD decreased because of re-volatilization to the air, permeation into the materials, and degradation.
2.1.2. Concentrations in airThe indoor behavior of the insecticide calculated using InPestCFD was verified using the results of experiments performed by the Sumitomo Chemical Co., Ltd. The calculated and measured values for the center of the room for 24 hr after spraying are shown in Fig. 6. The measured concentration at each time was determined by dividing the amount of the active substance collected by the sampling tube by the volume of air sampled in 20 min. The calculated concentration was determined by dividing the sum of the amounts of the active substance in air and particles that remained for 20 min in the calculation cell at the center of the room by the cell volume. The experimental and InPestCFD results indicated that the maximum airborne substance concentration was reached immediately after spraying and then the concentration decreased rapidly because of particle settling. Airborne particles were continually discharged in air through the outlet. However, the substance concentration did not markedly decrease from 6 hr after spraying because of re-volatilization of the substance, mainly from the floor. InPestCFD reproduced this well.
The aerial concentrations were predicted using two different values for turbulent kinetic energy kt (m2/s2) in the cell adjacent to the air inlet. Airflow was calculated using kt=0.005 m2/s2 without a slit in the air inlet and using kt=0.435uin2=0.011 m2/s2 with a slit in the air inlet. These kt values were a factor of two different, but the aerial concentrations of the substance were not different (data not shown). We therefore used kt=0.011 m2/s2 to take the air inlet shape into consideration. Calculations were also performed using different numbers of aerosol particles (1,000 or 10,000 particles represented by one particle). However, the substance concentrations in the air were not markedly different (data not shown). We therefore treated 10,000 particles as one hypothetical particle in the model to avoid unnecessarily increasing the computation requirements.
2.1.3. Amounts of insecticide in indoor mediaThe time-dependent amounts of insecticide determined using InPest and InPestCFD are shown in Fig. 7. The time-dependent amounts of the substance determined using InPestCFD and InPest generally agreed, but there were some differences. For example, InPest takes only gravitational settling into consideration. All of the large particles (diameter 50 µm) rapidly disappeared from 70 sec after spraying. In contrast, InPestCFD takes both gravitational settling and airflow into consideration. Large particles remained suspended in the flowing air when the particles become smaller because of solvent volatilization. The substance in large particles therefore gradually disappeared after 70 sec (Fig. 7a). This agreed with experimental results indicating that gravitational settling becomes slower as the particle size decreases.39)
The amount of substance in the air is shown in Fig. 7d. InPest assumes that all aerosol particles are sprayed into the upper part of the room. It will take a little time for the particles to fall to the floor. Re-volatilization of the substance from the room materials occurs after the large aerosol particles have fallen to the floor. In contrast, InPestCFD takes motion caused by the turbulent velocities of the particles into consideration, meaning that particles can adhere to the ceiling and walls even immediately after spraying. The amount of insecticide in the air caused by re-volatilization therefore increased earlier in InPestCFD than InPest (Fig. 7d′).
InPest cannot express how the aerosol particles physically adhere to a ceiling and wall after the aerosol is sprayed slightly upward. InPest therefore calculates the amount of particles that fall to the floor and then distribute that into the ceiling, floor, and walls according to the proportion determined from measurements of the room materials. In contrast, InPestCFD can describe movements of particles caused by gravity, the aerosol release velocity, and airflow. The model therefore reproduced temporal variations in the amounts of the substance on the ceiling, floor, and walls without artificial manipulation. As shown in Fig. 7f, there was a large difference between the amounts of the substance on the walls calculated using InPestCFD (solid line) and InPest (dotted line). This difference is discussed in Results section 2.2.2, but it appeared that the amounts of the substance in the small areas around the inlet and outlet were not determined when the real measurements were made. When the amounts of the substance in the calculation cells adjacent to the inlet and outlet were subtracted from the amounts of the substance on the whole walls, the amounts on the walls calculated using InPestCFD (dashed line in Fig. 7f) and InPest were similar.
2.2. Spatial behavior2.2.1. Concentrations in airThe predicted and measured spatial distributions of the concentrations of the substance in air 6 hr after spraying are shown in Fig. 8. The measured concentrations were 0.05–0.1 µg/m3 and tended to be higher near the floor on the downwind side (B, D, I, and J in Fig. 8b) than elsewhere. This was because the substance was transported from the inlet (supplying fresh air) to the outlet and re-volatilized from the floor. This was described well by InPestCFD, and the calculated and experimental results generally agreed except for points E, F and H. The differences at points E, F, and H may have been caused by (1) insufficient reproducibility of the deposited residue on the floor, (2) uncertainty in the airflow or turbulence in the CFD model, or (3) experimental uncertainty. We concluded that the difference at point F was caused by (1). We assumed that little residue was deposited on the floor at point F (Fig. 10a) and so the rate of re-volatilization from the floor was low. However, for points E and H, the calculated values were higher than the measured values. We considered that the differences for points E and H were caused by (2), although the CFD part of InPestCFD contained viscosity and diffusion terms as well as turbulent kinetic energy terms validated in several studies.33,39,41) We assumed that the substance in air and re-volatilized from the floor was transported smoothly in the flowing air from the inlet along the floor to point E, as shown in the substance distribution along the plane (Fig. 8c2). However, side airflow was relatively weak, so the substance was retained at point H (Fig. 8c3).
The predicted (using InPestCFD) and measured spatial distributions of the substance on the ceiling, floor, and walls 6 hr after spraying were compared. The measured distributions are shown in Fig. 9. The measured concentrations were derived from the total residue on each cotton sheet divided by the size of the sheet (0.3 m×0.6 m). The calculated distributions of the substance on the floor and the wall containing the outlet are shown in Figs. 10a and 11a, respectively. The calculated concentrations were converted into the same form as the measured concentrations (Figs. 10b and 11b). To compare the calculated values, the concentrations on the floor and the outlet wall shown in Fig. 9 were extended symmetrically to give Figs. 10c and 11c, respectively. The residue concentrations divided by the sizes of the cotton sheets placed on the room materials are summarized in Table 1.
| Floor | Ceiling | Inlet wall | Side wall | Outlet wall | ||
|---|---|---|---|---|---|---|
| Measured value | Mean±S.D. (Min–max) | 247±158 (34–540) | 29±23 (2.7–66) | 12±4.6 (4.6–20) | 13±5.4 (4.5–27) | 26±22 (5.2–77) |
| Calculated value | Mean±S.D. (Min–max) | 294±291 (6.3–1058) | 34±21 (9.6–82) | 16±14a (1.4–65)a | 16±19 (1.8–44) | 233±613b (4.6–1922)b |
a 14±10 µg/m2 and (1.4–45 µg/m2) except for the cells around the inlet. b 19±30 µg/m2 and (4.6–113 µg/m2) except for the cells around the outlet.
The calculated concentration on the floor was higher in the center of the floor than near the walls, which generally agreed with the measured concentrations (Fig. 10). The measured concentrations on the floor were 34–540 µg/m2 and the calculated concentrations on the floor were 6.3–1058 µg/m2. The mean measured residue concentration on the floor was 247 µg/m2, which was similar to the calculated residue concentration on the floor (294 µg/m2), as shown in Table 1. The measured concentrations on the side walls (not containing the inlet or outlet) were 4.5–27 µg/m2 (mean 13 µg/m2) and the calculated concentrations were 1.8–44 µg/m2 (mean 16 µg/m2). The measured concentrations on the ceiling were 2.7–66 µg/m2 (mean 29 µg/m2), and the calculated concentrations were 9.6–82 µg/m2 (mean 34 µg/m2). These concentrations indicated that the calculated and measured concentrations roughly agreed.
The calculated and measured concentrations on the wall containing the outlet are shown in Fig. 11. The calculated concentrations around the outlet were extremely large (Figs. 11a and 11b) but the measured concentrations were not (Fig. 11c). It appeared that the cotton sheet could not be spread all the way to the outlet edge, so the amount of the substance near the outlet was not determined in the experiment. The calculated concentrations in the cells adjacent to the outlet were 1678 and 1922 µg/m2 (Fig. 11b). However, subtracting the concentrations in the cells adjacent to the outlet gave concentrations of 113 and 73 µg/m2, respectively, which were similar to the measured concentration of 77 µg/m2 (Fig. 11c). The mean calculated residue concentration on the outlet wall (233 µg/m2) was quite a lot higher than the measured residue concentration (26 µg/m2), as shown in Table 1, but the calculated mean concentration except for the cells around the outlet was only 19 µg/m2. The calculated distribution data for the wall containing the inlet are not shown, but the calculated concentrations for the cells below the inlet were higher (24 and 65 µg/m2) than the measured concentration (9.2 µg/m2 in Fig. 9). This would have been because air became caught around the inlet and circulating flow occurred in the simulation. Similar to the outlet, the experiment did not detect high substance concentrations around the inlet because the cotton sheet had limited coverage over the inlet edge. Subtracting the substance concentrations in the cells adjacent to the air inlet gave concentrations of 8.4 and 45 µg/m2. The mean calculated residue concentration on the inlet wall (16 µg/m2) was decreased to 14 µg/m2 by excluding the cells around the inlet, and the measured mean was 12 µg/m2. No significant differences were therefore found between the calculated and measured spatial distributions on the ceiling, floor, and walls.
2.2.3. Residues in the roomInPestCFD was able to calculate the amounts of the substance residues on the ceiling and walls as well as on the floor from particle movements caused by the aerosol release velocity, gravity, and indoor airflow. The ratios for the substance concentrations on the room materials 6 hr after spraying are shown in Fig. 12. The predicted concentrations on the inlet and outlet walls were calculated after subtracting the amounts of the substance on the cells adjacent to the inlet and outlet, respectively, as described in the 2.2.2. About a quarter of the substance (22% for the measurements, 30% for the predictions) remained on the floor because of gravitational settling of the aerosol particles. Several percent of the substance (3% for the measurements, 4% for the predictions) remained on the ceiling (because the aerosol was sprayed slightly upward). A small proportion of the substance (2% for the measurements, 4% for the predictions) was on the outlet wall because the particles were carried toward the outlet by the flowing air. The concentrations were lower on the side walls than the outlet wall. However, the amount of substance remaining on the side wall (2% for the measurements, 3% for the predictions) was similar to the amount of substance remaining on the outlet wall because the area was larger for the side walls (two walls each 3.6 m×2.4 m) than the outlet wall (2.7 m×2.4 m). The amount of the substance on the inlet wall was relatively low (0.7% for the measurements, 0.9% for the predictions) because of fresh air being drawn into the room through the inlet. These results were reproduced well by InPestCFD.
When an insecticide is sprayed into indoor air, humans are exposed to the active substance through the dermal and/or oral routes. The exposure level can be estimated from the substance concentration on the floor, which is usually calculated assuming that all of the substance is deposited uniformly on the floor.44–46) Therefore, no model has previously been developed to predict the substance concentrations on the ceiling and walls of such a room. The substance concentration distributions on floors have been determined in a few studies,47) but considerable effort is required to make such measurements. However, exposure through inhalation can be estimated from the substance concentration in the air, which is usually measured by drawing indoor air through a collection tube. This method disturbs the airflow near the sampling point. It is therefore not possible to determine the substance concentration distribution in the air in detail because many sampling points would be needed and there would be strong total airflow disturbance. InPestCFD can predict the spatial distributions of a substance in the air and on the ceiling, floor, and walls over time. This makes it possible to assess differences in exposure through the dermal and oral routes caused by an infant crawling on different areas of the floor and then licking their hands. It also makes it possible to assess differences in exposure through inhalation caused by a person standing or sitting in different areas in the room. InPestCFD therefore makes it possible to perform various personal exposure assessments to take into account behavioral patterns.
In the experiment on the distribution of d-tetramethrin on room materials, the room was structurally symmetrical in the plane containing the air inlet and outlet. The cloud of sprayed aerosol particles was virtually uniform immediately after spraying through the four windows. Cotton sheets were therefore placed uniformly on the room materials. The analysis results determined assuming that the room was bilaterally symmetrical are summarized in Fig. 9. However, assessing indoor air movement using the CFD method led us to conclude that the airflows were different at the different windows. InPestCFD clearly indicated that the residues on the room materials were not bilaterally symmetrical (Figs. 10a, 10b, 11a, and 11b). The model indicated that even the airborne concentration was asymmetrical (Fig. 8c1 and 8c3). The model makes it possible to describe the distribution of the substance in the air and on the room materials after the insecticide is sprayed at any particular position in a room.
Distribution of the substance in the room is strongly affected by indoor airflow and the aerosol canister’s method of spraying. For example, the indoor airflow is changed by the position and/or size of the air inlet and outlet, and/or the ventilation rate of the room. In this study, four aerosol canisters were sprayed from four directions at the same time. Thus, the concentrations in each medium were relatively uniform. If an aerosol canister was sprayed from one direction in a manner performed usually at home, some locations in the air and on the floor with higher or lower concentrations of the substance would be more clearly evident. When such spraying was repeated many times for a long term, a specific location in the room would become extremely high. Because the conventional fugacity model assumes that the substance is uniform in each medium, it is not possible to express these phenomena. In contrast, because InPestCFD can describes the distribution of the substance affected by indoor airflow and spraying method of the aerosol canister, it is possible to predict these phenomena and contribute a more detailed exposure assessment for both short and long terms.
In InPest, aerosol particles are roughly divided into three groups: 10% with diameters of 10 µm, 60% with diameters of 30 µm, and 30% with diameters of 50 µm. These values were based on both the fly knockdown effect and insecticidal effect being optimal when the average aerosol particle size is about 30 µm48) and on information about the aerosol particle sizes 30 cm from the canister nozzle. However, when the particle sizes were measured in a rigorous way, 10% of the particles were found to have diameters of 6 µm, 60% to have diameters of 20 µm, and 30% to have diameters of 45 µm. The particles were smaller than described above. When InPestCFD was used with these sizes, the predicted average amount of the substance on the floor (228 µm/m2) was slightly lower than the measured average amount of the substance on the floor (250 µm/m2). The calculated amounts on the ceiling and walls (40 and 22 µm/m2, respectively) were slightly higher than the measured amounts on the ceiling and walls (29 and 16 µm/m2, respectively). This was because the smaller particles would have made the aerosol particles less likely to fall to the floor and more readily remain suspended in the air. As mentioned above, how the particles are described at time 0 affects how they are predicted to be distributed in the room. However, it is not easy to precisely classify and describe particles in the air for two reasons. First, we assumed that P∞ (the vapor pressure of the solvent away from the particle) was 0 in Eq. S3 (see Supplemental methodology 1). However, the aerosol particles in a real aerosol cloud are very close to each other, so P∞ will not be 0. This means that the actual rate at which the particles become smaller will probably be lower than the predicted rate. Second, Zhang et al.49) found that the particle sizes in air increase as the distance from the nozzle increases because the aerosol particles will collide and aggregate. We only considered decreases in aerosol particle diameters and did not consider increases in the particle diameters. There is currently a limit to the P∞ value input and no suitable mathematical formula for collision and aggregation is included in InPestCFD. We will revise InPestCFD to address new findings when they become available.
We developed a new fugacity model incorporating CFD for the first time. The new model was used to simulate the behavior of an insecticide released from an aerosol canister in a room. The model is applicable for use in estimating the behaviors of other chemicals in a room because it uses the physicochemical properties of the chemical of interest as input parameters. The model is also applicable for use in predicting the behaviors of pesticides in greenhouses by using the fugacity capacities of foliage and soil instead of room materials. We believe that the model could potentially be used to predict the behaviors of various chemicals outdoors by taking variations in wind speed, wind direction, and temperature caused by climatic conditions into consideration.
We developed a new model incorporating CFD into a fugacity model for the first time. InPestCFD was used to reproduce the spatiotemporal behavior of an insecticide sprayed into a room. The predicted and measured behaviors generally agreed. InPestCFD allows more detailed exposure assessments than were previously possible to be performed, including assessments taking into account the position in the room in which a person stands or sits. We will next use InPestCFD to improve our understanding of the behaviors of insecticides under various conditions with different spraying positions and the room conditions (e.g., air inlet positions and ventilation rate). This will allow more useful assessments of the risks posed to humans by chemical substances than can currently be performed. We will attempt to increase the reliability of the model by using it to predict the behaviors of other chemicals outdoors as well as indoors. We will use the model to estimate exposure under various conditions to allow better risk assessments to be performed and therefore to improve the safety of humans and other biota when chemical substances are used.
The online version of this article contains supplementary material, which is available at https://www.jstage.jst.go.jp/browse/jpestics/.
