A New Model Algorithm for Estimating the Inhalation Exposure Resulting from the Spraying of (Semi)-Volatile Binary Liquid Mixtures (SprayEva)

The spraying of liquid multicomponent mixtures is common in many professional and industrial settings. Typical examples are cleaning agents, additives, coatings, and biocidal products. In all of these examples, hazardous substances can be released in the form of aerosols or vapours. For occupational and consumer risk assessment in regulatory contexts, it is therefore important to know the exposure which results from the amount of chemicals in the surrounding air. In this research, a mechanistic mass balance model has been developed that covers the spraying of (semi)-volatile substances, taking into account combined exposure to spray mist, evaporation from droplets, and evaporation from surfaces as well as the nonideal behaviour of components in liquids and backpressure effects. For wall-spraying scenarios, an impaction module has been developed that quantifies the amount of overspray and the amount of material that lands on the wall. Mechanistically, the model is based on the assumption that continuous spraying can be approximated by a number of sequentially released spray pulses, each characterized by a certain droplet size, where the total aerosol exposure is obtained by summation over all release pulses. The corresponding system of differential equations is solved numerically using an extended Euler algorithm that is based on a discretisation of time and space. Since workers typically apply the product continuously, the treated area and the corresponding evaporating surface area grows over time. Time-dependent concentration gradients within the sprayed liquid films that may result from different volatilities of the components are therefore addressed by the proposed model. A worked example is presented to illustrate the calculated exposure for a scenario where aqueous solutions of H2O2 are sprayed onto surfaces as a biocidal product. The results reveal that exposure to H2O2 aerosol reaches relevant concentrations only during the spraying phase. Evaporation from sprayed surfaces takes place over much longer time periods, where backpressure effects caused by large emission sources can influence the shape of the concentration time curves significantly. The influence of the activity coefficients is not so pronounced. To test the plausibility of the developed model algorithm, a comparison of model estimates of SprayExpo, SprayEva, and ConsExpo with measured data is performed. Although the comparison is based on a limited number (N = 19) of measurement data, the results are nevertheless regarded as supportive and acceptable for the plausibility and predictive power of SprayEva.


Introduction
Room spraying as well as spraying onto rigid walls occurs in many industrial, professional, and consumer applications and can lead to high inhalation exposure to vapours and aerosols. At the same time, dermal exposure via the deposition of aerosol on the skin and via other pathways can also play a relevant role [1,2]. Hence, the components is quite flexible and reflects a wide variety of exposure situations and pattern related to the spraying of (semi)-volatile binary mixtures.

General Mass Balance Equations
In essence, the conceptual model for aerosol exposure described above can be expressed in mathematical terms by a set of mass balance differential equations. These mass balance equations are based on the assumption that continuous aerosol release can be subdivided into a number N l of pulses l released at discrete time steps and that overall exposure to aerosol droplets in a room with a volume V room results from summation over all release pulses with initial aerosol concentration A init,c and droplet size classes c, where each droplet size class is characterised by an initial droplet diameter d init,c . When the interval ∆t l between the pulses decreases, a quasi-continuous application can be achieved.
Taking into account the differential mass loss due to evaporation dm ev (t) dt , the mass loss due to settling dm Fsed (t) dt and the mass loss due to ventilation dm vent (t) dt as the relevant transport mechanisms, the differential change of the initial aerosol concentration A l,c (t l ) in the air of a single release pulse l can be expressed as: Given that (numerical) solutions to Equation (1) are available (which will be demonstrated in Sections 2.2 and 2.3), the total aerosol exposure at the time t is then obtained by summation over all release pulses N l and droplet size classes N c : Analogously, the differential change dC room,i dt in the vapour concentration of a component i can be expressed taking into account the vapour mass loss rate of component i due to ventilation Q vent ·C room,i (t) and the emission rates resulting from aerosol evaporation dm Aev,i (t) dt and eventually from evaporation from the floor  . The term C vent,i describes the vapour concentration of component i in the supply air and Q vent the ventilation rate.
From these mass balance considerations, it is clear that the physical processes behind the emission and mass loss rates are key to understanding the time dependency of aerosol and vapour exposure. Therefore, the corresponding processes will be described in more detail in Sections 2.4-2.9.

Discretisation Approach
As the mass balance for a single release pulse can be described by a set of only two time-varying differential equations (Equations (1) and (3)) for each of the involved compounds, solving this system of ordinary differential equations (ODEs) analytically seems to be feasible at first sight. On the other hand, the calculation of the total exposure requires summation over all release pulses, which prevents analytical solutions to this problem. Hence, in order to obtain approximate time-dependent solutions to Equations (1) and (3), an extended Euler method was used that was developed in previous work [13].
In the case of wall spraying, a moving spray cone is assumed that sweeps a certain surface area W during the application. Hence, the treated surface area and the corresponding emission rate increase over time. As the sprayed substances may differ in their volatility, spatial concentration gradients within the sprayed liquid film on the surface can occur. In order to address this problem, we break down a certain area W into a sequence of N l small area elements ∆W that are numbered consecutively using index l. That is, the size of ∆W = W/N l determines the spatial resolution of the proposed approach. We assume that fractions of the product are applied instantaneously, spray pulse by spray pulse, where the applied amount reflects the work rate of the worker. However, in principle, this approach can also take arbitrary values for the applied amount, or even zero to reflect periods when the worker stops the application (breaks, postapplication phase).
The next step in the discretisation procedure is to define the temporal resolution that is needed to adequately describe the time-dependent aerosol and vapour concentrations resulting from each release pulse. This can be accomplished by dividing the overall exposure time t expo into N t time steps so that ∆t = t expo /N t . In brief, the time ∆t l between two release pulses l can be the same as ∆t. However, for reducing computational time, it may also be useful and justifiable to reduce the number of release pulses N l by allowing ∆t l to adopt integer multiples IM of ∆t, that is ∆t l = ∆t·IM and N t = N l ·IM. For clarification, in the results chapter, the value of IM will be stated explicitly for the example calculations.
In Figure 1, the discretisation procedure is exemplified for aerosol concentration A assuming IM = 1. The discretisation of the related variables of interest (droplet diameter, mole fractions, vapour concentration) and their substance-specific time-dependent derivates is not shown but can be expressed analogously. For each time step k and each release pulse l, the aerosol concentration at the beginning of time step k is represented by A l,c (t k ). For instance, A 2,1 (t 3 ) means the aerosol concentration of release pulse 2 in droplet size class 1 at the beginning of time step 3. The diagonal elements (that is k = l) of the triangular matrix in Figure 1 represent initial aerosol concentrations resulting from release pulse l at the time t k=l = ∆t (k − 1). The temporal iteration process starts with time step k = 1 and proceeds successively for each release pulse l to the maximum number of time steps N t .
In the case of wall spraying, a moving spray cone is assumed that sweeps a certain surface area W during the application. Hence, the treated surface area and the corresponding emission rate increase over time. As the sprayed substances may differ in their volatility, spatial concentration gradients within the sprayed liquid film on the surface can occur. In order to address this problem, we break down a certain area W into a sequence of small area elements ∆ that are numbered consecutively using index l. That is, the size of ∆ = / determines the spatial resolution of the proposed approach. We assume that fractions of the product are applied instantaneously, spray pulse by spray pulse, where the applied amount reflects the work rate of the worker. However, in principle, this approach can also take arbitrary values for the applied amount, or even zero to reflect periods when the worker stops the application (breaks, postapplication phase).
The next step in the discretisation procedure is to define the temporal resolution that is needed to adequately describe the time-dependent aerosol and vapour concentrations resulting from each release pulse. This can be accomplished by dividing the overall exposure time texpo into Nt time steps so that Δt = texpo/Nt. In brief, the time Δtl between two release pulses l can be the same as Δt. However, for reducing computational time, it may also be useful and justifiable to reduce the number of release pulses Nl by allowing Δtl to adopt integer multiples IM of Δt, that is Δtl = Δt·IM and Nt = Nl·IM. For clarification, in the results chapter, the value of IM will be stated explicitly for the example calculations.
In Figure 1, the discretisation procedure is exemplified for aerosol concentration A assuming IM = 1. The discretisation of the related variables of interest (droplet diameter, mole fractions, vapour concentration) and their substance-specific time-dependent derivates is not shown but can be expressed analogously. For each time step k and each release pulse l, the aerosol concentration at the beginning of time step k is represented by , ( ). For instance, , ( ) means the aerosol concentration of release pulse 2 in droplet size class 1 at the beginning of time step 3. The diagonal elements (that is k = l) of the triangular matrix in Figure 1 represent initial aerosol concentrations resulting from release pulse l at the time tk=l = Δt (k − 1). The temporal iteration process starts with time step k = 1 and proceeds successively for each release pulse l to the maximum number of time steps Nt. That is, each release pulse l of droplet size class c is assumed to instantaneously apply an initial amount to the room volume at time leading to an initial aerosol ( ) . Correspondingly, the initial release rate of the That is, each release pulse l of droplet size class c is assumed to instantaneously apply an initial amount to the room volume at time t k=l leading to an initial aerosol concentration ·∆t l . Correspondingly, the initial release rate of the spray is then For wall spraying, the term dm l,c (t k=l ) dt has to be replaced by the release rate of the overspray dm os,l,c (t k=l ) dt (see Section 2.8). Since all equations specified above refer to a droplet size class c, we now need to define the overall distribution that includes all the droplet size classes. Although there is no fundamental theoretical reason why droplet size data should approximate the lognormal distribution, it has been found to apply to most single-source aerosols [14]. Assuming a lognormal droplet size distribution with median mass diameter d m and geometric standard deviation GSD, a number of N c droplet size classes are defined, each in turn specified by a median mass diameter d m,c . The borders of each droplet size class c with median mass diameter d m,c are chosen in a way that it includes an equal proportion of 1/N c of the total mass release rate rr = dm l (t k=l ) dt of the spray. Then, the release rate for each droplet size class dt ·∆t l , describing the initial aerosol concentration of a spray pulse. An example of a cumulative mass distribution assuming lognormality is shown in Figure 2 with d m = 100 µm, GSD = 2 and N c = 7. The mass mediansc d m,c and boundaries of each droplet size class are depicted as green and black dots, respectively, each representing 1/N c of the total released mass. The symbol CDF denotes the cumulative mass distribution function.
define the overall distribution that includes all the droplet size classes. Although there is no fundamental theoretical reason why droplet size data should approximate the lognormal distribution, it has been found to apply to most single-source aerosols [14]. Assuming a lognormal droplet size distribution with median mass diameter and geometric standard deviation GSD, a number of droplet size classes are defined, each in turn specified by a median mass diameter , . , describing the initial aerosol concentration of a spray pulse. An example of a cumulative mass distribution assuming lognormality is shown in Figure 2 with = 100 µm, GSD = 2 and = 7. The mass medians , and boundaries of each droplet size class are depicted as green and black dots, respectively, each representing 1/ of the total released mass. The symbol CDF denotes the cumulative mass distribution function.

Numerical Approximation Using Euler's Method
To approximate the time-dependent solutions of the system of ordinary differential equations (ODE) described above and more specified in Sections 2.4-2.9, the explicit Euler method is used. Although the Euler method is straightforward and hence not very precise [15], it is regarded sufficient as in occupational exposure modelling, high numerical precision is not required due to other sources of uncertainty in the exposure assessment process. Starting with the initial value y(t0) = y0, the Euler method iterates approximate solutions using the equation ( ) = ( ) − ∆ • ( ) where = ( − 1) • ∆ with ∈ ℕ, = 1, . . . . The formulas of Euler's algorithm for a single ODE do not change when the method is used to solve a system of ODEs. Hence, the iteration formula applied to Equation (1) for a release pulse l and droplet size class c for instance gives

Numerical Approximation Using Euler's Method
To approximate the time-dependent solutions of the system of ordinary differential equations (ODE) described above and more specified in Sections 2.4-2.9, the explicit Euler method is used. Although the Euler method is straightforward and hence not very precise [15], it is regarded sufficient as in occupational exposure modelling, high numerical precision is not required due to other sources of uncertainty in the exposure assessment process. Starting with the initial value y(t 0 ) = y 0 , the Euler method iterates approximate solutions using the equation The formulas of Euler's algorithm for a single ODE do not change when the method is used to solve a system of ODEs. Hence, the iteration formula applied to Equation (1) for a release pulse l and droplet size class c for instance gives The iteration process over time for each release pulse l of droplet size class c at the time t k=l starts with the initial aerosol concentration A l,c (t k=l ). It should be noted that the initial value can change over time, so that the modelling of intermittent application patterns is possible. The total aerosol concentration is obtained by summation over all release pulses and droplet size classes: The substance-specific aerosol concentration can be calculated by multiplying the right side of Equation (5) with the mass fraction w i,l,c (t k ) of the respective component i: where and M i denotes the molecular weight of component i and x i,l,c the mole fraction of component i in the droplet size class c of the aerosol. In terms of exposure, it is important to note that not all sizes of droplets suspended in the air can be breathed in with 100% efficiency. Thus, any exposure assessment should focus on the particles that can actually enter the respiratory tract. The International Standards Organisation (ISO) [16] has agreed upon a so-called inhalable fraction. For wind velocities less than 4 m/s, the inhalable fraction is given by The inhalable fraction of the total aerosol concentration is then obtained by multiplying the right sides of Equations (5) and (6) with the right side of Equation (9).
The other variables of interest (x i , d, and C i ) can be obtained analogously by taking into account the mass losses of each of the three relevant mass transport processes (ventilation, evaporation, sedimentation). The corresponding mass balance equations needed for the iteration procedure will be described in the following.
It should be noted that the programming of the extended Euler algorithm of SprayEva is straightforward, essentially using three Do loops nested over k, l, c. In the Supplementary Materials, a piece of pseudocode based on Mathematica ® is given that is intended to illustrate how the iteration algorithm works for the room spraying of a binary mixture.

Mass Loss Due to Ventilation
In the well-mixed room (one-box) model, the sprayed material is assumed to be dispersed immediately and homogeneously through the room air [17][18][19]. At first sight, this assumption seems to be too simplistic given the fact that during spraying, the concentration of the sprayed material is highest near the source and forms a concentration gradient throughout the room. On the other hand, air turbulence is generated by the momentum flux of the sprayed aerosol and ventilation air as well as thermal convection. Thus, homogenisation inside a room can be achieved in minutes [20] as long as the room volumes are not too large and the room is well mixed due to either natural or technical ventilation. Koch et al. [12] concluded from a number of model calculations and comparison with measured data that the average aerosol concentration in smaller rooms (<150 m 3 ), can be described sufficiently well.
From a mathematical point of view the mass loss rate due to ventilation dm vent,l,c (t k ) dt for an aerosol release pulse l of droplet size class c at time t k can be expressed in quite simple terms [18]: Analogously, the differential change in the overall vapour concentration C room,i (t k ) of component i at time t k due to mass loss by ventilation is determined by the term Q vent ·C room,i (t k ). This simplification is justifiable when room volumes are not too large and when the room is well mixed due to either natural or technical ventilation [19,21].

Mass Loss Due to Droplet Settling
The vertical motion of the droplets and their deposition onto horizontal surfaces is determined by settling in the earth's gravity field. The velocity v sed at which droplets fall to the floor is given by Stokes's settling velocity: This formula applies to droplets with diameters up to 50 µm. For larger droplets, a correction is necessary [14]. The changing droplet diameter d(t) and density ρ dr (t) of the droplet due to evaporation make the settling velocity a time-dependent quantity decreasing with increasing residence time, where g denotes the gravitational constant and η the air viscosity. The mass density of the droplet ρ dr,l,c (t) can be calculated using the liquid molar volumes V i of the pure components and their mole fractions x i , where N i denotes the total number of components: Taking into account the surface area of the floor, F, the mass loss rate dm Fsed,l,c (t k ) dt due to sedimentation for an aerosol release pulse l of droplet size class c at time t k can then be expressed as: The substance-specific settling rates can be obtained by multiplying the right side of Equation (14) with the mass fraction w i (t) (see Equation (7)) of the respective component i: The summation of the sedimentation flow of component i over all release pulses l and droplet size classes c at the time t k then results in The left side of Equation (16) denotes the total mass flow of component i, depositing on and potentially evaporating from the floor (see Section 2.7).

Evaporation of Two-Component Droplets
The evaporation of the (semi)-volatiles from the liquid phase, i.e., the vapour-droplet partitioning, depends on their effective equilibrium vapour pressures, which are determined by the chemical compositions and thus the concentrations and activity coefficients of the constituents in the spray formulation. Commercial preparations are often complex mixtures of many substances with different physical properties. In addition, a wide variety of use scenarios is possible, making the evaporation of multicomponent droplets a complex process. Hence, for a basic understanding of multicomponent droplet evaporation, it is necessary to investigate droplet evaporation using simplified models. A number of models have been proposed for the evaporation of substances from aerosol droplets, each of which makes different assumptions and simplifications [22]. This study is based on the simplified rapid-mixing model (SRMM) for two-component droplet evaporation that was suggested by Wilms [22]. It is called a rapid-mixing model because it assumes a rapid mixing of the different liquids so that there is no concentration gradient within the droplet. The following differential equations are given for the evaporation rate, the size of a droplet, and its composition for binary mixtures. By extension to Wilms, vapour concentrations C room far away from the droplet are not assumed to be zero but can take finite values. This is regarded as an important modification as backpressure effects may influence the evaporation significantly. The evaporation mass flux of a single droplet then results for a component i in where d is the droplet diameter, D i the effective binary diffusion coefficient of the component i in the air, M i the molecular weight, p * i the saturation vapour pressure, x i the molar fraction, and γ i the activity coefficient of component i. The activity coefficient γ i , which depends on the molar fraction of component i, was introduced in Equation (17) to account for nonideal mixtures. Due to the consumed heat of evaporation, the effective temperature of the droplets is lower than the air temperature. Therefore, a reference temperature T r was used in this work (for Equation (17)) approximated using Hubbard's 1/3 rule [22,23], where T dr is the temperature at the droplet surface, estimated according to Hinds [14], who proposed an empirical equation for water droplets, and T room the room temperature.
Equation (17) only describes the mass flow of a single droplet and needs to be extended for the calculation of a spray pulse probably consisting of thousands of droplets. For this purpose, it is assumed that the droplets are spherical and noninteracting and that the number of droplets N dr in the room volume is related to the aerosol concentration A(t) and the density of the droplet ρ dr via the following equation: Assuming a two-component droplet (N i = 2) and combining Equations (13), (17), and (18), the following expression can be derived for the vapour emission of a release pulse l of component 1 at the time t k : The emission rate of the second component can then be determined by rearranging Equation (19) using the equation The summation of the evaporation flow of component i over all release pulses l and droplet size classes c at the time t k results in Since Equations (19) and (20) include d l,c (t) and x 1,l,c (t) as unknown functions, the set of equations needs to be extended with corresponding differential equations for . Following Wilms [22] (chapter 3.6) and using the simplified rapid-mixing model (SRMM) for two-component droplets, the following equations are obtained, where V i denotes the molar volume of component i and (23) Given that (numerical) solutions to Equation (23) are available, the mole fraction of the second component can be determined via the equation

Evaporation from the Floor
The vertical motion of droplets in the earth's gravity field leads to a mass flow of (semi)-volatile components potentially depositing on the floor as a liquid film. At the same time, evaporation from the liquid film occurs leading to a mass loss of the film by vapour emission. Hence, both processes have to be taken into account for deriving the differential mass change While the sedimentation flow is governed by Equation (16), the evaporation of a component i released at time t k from the liquid film can be predicted using the model suggested by Gmehling and Weidlich [24,25]: In this model, it is assumed that evaporation is driven by the difference between the partial vapour pressure of an individual component i and its vapour pressure (backpressure) in the room air p room,i = C room,i ·R·T/M i , which is reflected in Equation (25) by the term C room,i (t k ). The evaporation rate of substance i, , is proportional to this pressure difference and depends on the surface area F of the mixture and the mass transfer coefficient ß i . The activity coefficient γ i as a function of x F,i (t k ) is used to take into account nonideal behaviour of the components in the liquid as defined by Raoult and Henry. The mass transfer coefficient ß i is a function of the diffusivity D i of the substance in the air and of the air flow v air over the product surface [24]: Parameter notation and corresponding symbols and units are listed in Table 1. Previous work has described this modelling approach in more detail [13].
Combining Equations (24) and (25) then gives the differential mass change where the mole fraction of component i is defined as Since Equations (25) and (27) are coupled by Equation (28), the evaporation rate dm Fev,i (t k ) dt and the differential mass change of component i can be calculated simultaneously. In the case that evaporation is faster than sedimentation, potentially leading to numerical results where m F,i (t k ) ≤ 0, the evaporation flow is assumed to be determined by the sedimentation flow, that is, . In other words, the droplets are assumed to evaporate instantaneously when reaching the floor.

Impaction Module
For spraying on surfaces, only the nonimpacting fraction of the droplet spectrum (overspray) is relevant for aerosol exposure. At the same time, the fraction of the spray deposited on surfaces can evaporate and can contribute significantly to the overall vapour exposure. The importance of spray penetration and impaction is well recognized and has been extensively studied experimentally and theoretically [26]. A rigorous theory of spray dynamics would be very complex, making a self-consistent modelling of these processes a major challenge. However, an in-depth understanding of all involved processes is not always essential as long as modelling arrives at results within the required accuracy boundaries. On many occasions, it is therefore justifiable to develop simplified models suitable for practical applications. Based on the work of Sazhin et al. [27], Medrano et al. [28], and Su and Yao [29], the following simplifying assumptions are made in this study for estimating the proportion of the overspray PO and the (complementary) proportion (1 − PO) of the spray deposited on surfaces. - The initial droplet velocity v 0 is equal to the velocity of the liquid at the injector tip. -A full-cone turbulent round jet ensuing from a circular nozzle is assumed. -Air entrainment beyond the vicinity of the nozzle is considered. - The droplet size distribution is constant during the transport process. - The approach is restricted to medium volatile mixtures that should not exceed the vapour pressure of water.
According to Su and Yao [29], an inertial impaction parameter K for impacting sprays is defined in terms of the momentum conservation to reveal the relationship between spray transportation and aerodynamic conditions. The impaction parameter is defined in Equation (29) as the quotient of the stopping distance λ(x) based on the flow mean velocity v(x) to the spray diameter D(x) at the cross section of the target plate. The symbol ρ dr denotes the droplet density, η the dynamic viscosity of the air, and v(x), the flow mean velocity of the spray as a function of the impacting distance s. The droplet diameter d is assumed to be constant within the flying time and thus can be set equal to the initial droplet diameter d init,c , belonging to droplet size class c. These are simplifying assumptions, but given that the flying time of larger droplets that may reach the surface is much shorter than their evaporation time, they are probably reasonable [1]. In addition, smaller droplets predominantly belong to the overspray, where they evaporate, contributing not to the impacted mass but to vapour exposure.
In order to obtain the flow mean velocity v(s) needed in Equation (29) Sazhin et al. [27] and Medrano et al. [28] propose a two-phase-fluid model for a full-cone turbulent round jet with only the apex angle θ of the cone being a disposable parameter (see Figure 3).
is defined in terms of the momentum conservation to reveal the relationship between spray transportation and aerodynamic conditions. The impaction parameter is defined in Equation (29) as the quotient of the stopping distance ( ) based on the flow mean velocity _ ( ) to the spray diameter ( ) at the cross section of the target plate. The symbol denotes the droplet density, the dynamic viscosity of the air, and _ ( ), the flow mean velocity of the spray as a function of the impacting distance s. The droplet diameter d is assumed to be constant within the flying time and thus can be set equal to the initial droplet diameter dinit,c, belonging to droplet size class c. These are simplifying assumptions, but given that the flying time of larger droplets that may reach the surface is much shorter than their evaporation time, they are probably reasonable [1]. In addition, smaller droplets predominantly belong to the overspray, where they evaporate, contributing not to the impacted mass but to vapour exposure.
In order to obtain the flow mean velocity _ ( ) needed in Equation (29) Sazhin et al. [27] and Medrano et al. [28] propose a two-phase-fluid model for a full-cone turbulent round jet with only the apex angle of the cone being a disposable parameter (see Figure  3). Figure 3. Diagram of the mass entrainment process (adopted from [28]). Figure 3. Diagram of the mass entrainment process (adopted from [28]).
As a result, the dynamics of both droplets and entrained air can be described in terms of a two-phase flow with a zero relative velocity v(s) between air and droplets depending on the distance from the nozzle, s. Here, v 0 denotes the initial droplet velocity at the nozzle. v(s) = 2·v 0 When the result in Equation (30) is substituted into Equation (29), the impaction parameter K can now be calculated from readily available quantities.
The results obtained by Su and Yao [29] indicate that the transfer efficiency ξ(K) of impacting sprays onto the plate is a function of the impaction parameter. The predicted transfer efficiency is plotted against the impaction parameter K in Figure 1. It can be seen that almost all the points fall on the same curve. When K is greater than 0.3, the transfer efficiency is almost 1. In other words, if K exceeds the critical value of 0.3, according to [29], the droplets will be deposited; otherwise, they will be released into the air as overspray. The droplet diameter that corresponds to the critical value of K is called the critical diameter d crit .
Since a wall can be regarded as an infinite plate, the transfer efficiency ξ(K c ) of the impacting spray can be defined as the ratio of the rate of droplets actually deposited onto the wall dm Wd,l,c (t k ) dt to the overall rate is then obtained as: The substance-specific mass flow rate dm Wd,i,l,c (t k=l ) dt can be obtained by multiplying the right sides of Equations (31) and (32) with the mass fraction w i (t) (see Equation (7)) of the respective component i. Summation over all droplet size classes results in the overall impacting flow of component i at time t k .

Evaporation from Sprayed Walls
The mass balance equations for substance evaporation from a sprayed rigid surface can be formulated along the lines of evaporation from the floor. However, in contrast to the evaporation from liquid films deposited on the floor, spraying onto walls leads to an increase in the applied area that results in an increasing evaporation rate. As suggested in previous work [13], a logical approach to this problem is to break down the application area W into a sequence of N l small area elements, where each area element of size ∆W = W/ N l has its individual emission rate. The overall emission rate is then the sum of the individual emission rates.
For the differential mass balance of a single area element l, two mass flows must be considered for each component i, the impacting mass flow dm Wd,i,l (t k ) dt and the evaporation rate dm Wev,i,l (t k ) dt . The differential mass change dm W,i,l (t k ) dt of a component i in the liquid film on the wall at the time t k is then the difference of these flows: where analogously to Equation (25), the evaporation rate for an area element l can be expressed as If the evaporation rate is higher than the impaction flow, potentially leading to numerical results where m W,i,l (t k ) ≤ 0, the evaporation flow is assumed to be determined by the impaction flow, that is, . In other words, the droplets are assumed to evaporate instantaneously when reaching the wall.
Using Euler's method to estimate m W,i,l (t k ), x W,i,l (t k ), and C room,i (t k ), the total evaporation rate can then be calculated by summation over all area elements l: It should be noted that Equation (36) together with Equations (21) and (25), which are related to evaporation from aerosol droplets and from the floor, constitute the overall evaporation rate needed as indicated in Equation (3). Since evaporation from droplets is faster than evaporation from flat surfaces, the temporal resolution used in the Euler algorithm can be much lower when calculating the evaporation rate of flat surfaces. This fact was considered in the developed software by switching to longer discretisation times ∆t post if no aerosol generation had taken place (post-spraying phase) and exposure was dominated by evaporation from flat surfaces (wall, floor). In this way, the computing time could be significantly reduced.

Example Calculation
Due to its flexibility, the proposed model algorithm potentially covers a range of exposure situations including the spraying of mixtures. One area where spraying is often used as an application technique is the disinfection of rooms and surfaces with biocidal products. In order to demonstrate the possibilities and limitations of the SprayEva algorithm, an example calculation is presented in the following with hydrogen peroxide as an active substance.
According to the Regulation (EU) No 528/2012 product assessment report [30], there are two product types amongst others, PT 2, 3, and PT 3 involves the disinfection of animal housing by spraying aqueous solutions of 7.4% (w/w) of hydrogen peroxide. This type of disinfection is typically performed by spraying aqueous hydrogen peroxide solutions with different devices generating larger droplets, leaving the product on the treated surfaces and not rinsing it off.
The worker typically uses a high-performance spraying device with a handheld spray lance. According to Koch et al. [12], for spraying onto a surface, flat fan or hollow cone nozzles with a flow rate of 0.6 L/min and a pressure of 3 bar are used, and flat cone nozzles tend to generate coarse droplets. Analyses performed so far with other models and the existing results of measurements at workplaces have demonstrated that the droplet spectrum of the spraying method used has an impact on the exposure concentration [12]. In the example calculations, a mass median diameter of d m = 250 µm and a hollow cone nozzle was assumed. The number of droplet size classes characterising the over spray was set to N c = 5. The geometric standard deviation of the droplet distribution was fixed to 1.8 and was not considered by Koch et al. [12] to be a parameter that is subject to much variation.
The example calculation aimed to characterize a situation where the room volume is located at the lower end of possible sizes of animal housings (see Table 1). The product concentration w pr,i , the surface coverage SC, and the ventilation rate Q vent were assumed to follow the data given in the ECHA assessment report. Since settled droplets probably do not cover the whole floor in real workplaces if only one wall is sprayed, the area of the floor relevant for evaporation was assumed smaller than the total floor area F. This was accounted for by multiplying F by an arbitrary factor of 3 4 , resulting in correspondingly reduced evaporation rates from the floor. Other process-and scenario-related parameters such as the room height h room , release rate rr, sprayed area S, application time t ap and exposure time t expo , opening angle of spray cone 2 θ, distance nozzle-wall s, and nozzle diameter D 0 were estimated using the data given in [12]. Using Equation (29) together with Figure 4, the critical diameter could be estimated as d crit = 202 µm.
The substance-specific input parameters, including the activity coefficients γ, needed for the model calculation were adopted from previous work [13] and are listed in Table 1. The reference temperature T r for droplet evaporation was estimated using Hubbard's rule [23] in combination with the empirical equation proposed by Hinds [14] for calculating the temperature at the droplet surface T dr . Regarding surface evaporation, calculations are based on the room temperature T room . The initial vapour concentrations C init of water and hydrogen peroxide were assumed to be zero.
In addition to these physically based parameters, the quantities determining the temporal and spatial resolution of the iteration must also be specified. The example calculations were carried out using temporal resolutions of ∆t = 0.2 s for the spraying phase and ∆t post = 2 s for the post-spraying phase, resulting in a total calculation time of about 8 min. The multiple integer IM of ∆t, determining the spatial resolution, was set to IM = 20, which proved to be a good compromise between accuracy and computation time.
floor relevant for evaporation was assumed smaller than the total floor area F. This was accounted for by multiplying F by an arbitrary factor of ¾, resulting in correspondingly reduced evaporation rates from the floor. Other process-and scenario-related parameters such as the room height hroom, release rate rr, sprayed area S, application time tap and exposure time texpo, opening angle of spray cone 2 θ, distance nozzle-wall s, and nozzle diameter D0 were estimated using the data given in [12]. Using Equation (29) together with Figure 4, the critical diameter could be estimated as dcrit = 202 µm. The substance-specific input parameters, including the activity coefficients γ, needed for the model calculation were adopted from previous work [13] and are listed in Table 1. The reference temperature Tr for droplet evaporation was estimated using Hubbard's rule [23] in combination with the empirical equation proposed by Hinds [14] for calculating the temperature at the droplet surface Tdr. Regarding surface evaporation, calculations are based on the room temperature Troom. The initial vapour concentrations Cinit of water and hydrogen peroxide were assumed to be zero.
In addition to these physically based parameters, the quantities determining the temporal and spatial resolution of the iteration must also be specified. The example calculations were carried out using temporal resolutions of Δt = 0.2 s for the spraying phase and Δtpost = 2 s for the post-spraying phase, resulting in a total calculation time of about 8 min. The multiple integer IM of Δt, determining the spatial resolution, was set to IM = 20, which proved to be a good compromise between accuracy and computation time.
The proposed extended Euler algorithm was programmed using the Mathematica ® environment. Typical calculation times are in the range of seconds to minutes as long as the spraying times are in the range of seconds to few minutes. However, spraying times > 5 min can be time consuming in the current code version, resulting in calculation times even in the range of hours. In contrast, longer post-spraying phases are less demanding in terms of computing speed because Δtpost can be set much longer than Δt. This is because evaporation from flat surfaces is slower than the evaporation from droplets. The proposed extended Euler algorithm was programmed using the Mathematica ® environment. Typical calculation times are in the range of seconds to minutes as long as the spraying times are in the range of seconds to few minutes. However, spraying times > 5 min can be time consuming in the current code version, resulting in calculation times even in the range of hours. In contrast, longer post-spraying phases are less demanding in terms of computing speed because ∆t post can be set much longer than ∆t. This is because evaporation from flat surfaces is slower than the evaporation from droplets.
Using the input from Tables 1 and 2 in Table 3, the time-dependent modelled aerosol (total and inhalable) and vapour concentrations as well as the sedimentation and emission rates were determined separately for H 2 O and H 2 O 2 . For clarification, the diagrams also include the individual time-dependent curves of the release pulses and the corresponding time-averaged curves. Using the input from Tables 1 and 2 in Table 3, the time-dependent modelled aerosol (total and inhalable) and vapour concentrations as well as the sedimentation and emission rates were determined separately for H2O and H2O2. For clarification, the diagrams also include the individual time-dependent curves of the release pulses and the corresponding time-averaged curves.        Lege nd Starting with the aerosol concentration, the graphs show that the time courses for the spraying and post-spraying phases can be quite different. In particular, it becomes clear that exposure to H2O2 aerosol reaches relevant concentrations only during the spraying phase. Table 3, panel 1 also reveals an increase in the aerosol concentration over the whole spraying phase, which can be explained by an increase in the backpressure (see panel 6), Starting with the aerosol concentration, the graphs show that the time courses for the spraying and post-spraying phases can be quite different. In particular, it becomes clear that exposure to H2O2 aerosol reaches relevant concentrations only during the spraying phase. Table 3, panel 1 also reveals an increase in the aerosol concentration over the whole spraying phase, which can be explained by an increase in the backpressure (see panel 6), Starting with the aerosol concentration, the graphs show that the time courses for the spraying and post-spraying phases can be quite different. In particular, it becomes clear that exposure to H2O2 aerosol reaches relevant concentrations only during the spraying phase. Table 3, panel 1 also reveals an increase in the aerosol concentration over the whole spraying phase, which can be explained by an increase in the backpressure (see panel 6), Starting with the aerosol concentration, the graphs show that the time courses for the spraying and post-spraying phases can be quite different. In particular, it becomes clear that exposure to H2O2 aerosol reaches relevant concentrations only during the spraying phase. Table 3, panel 1 also reveals an increase in the aerosol concentration over the whole spraying phase, which can be explained by an increase in the backpressure (see panel 6), Starting with the aerosol concentration, the graphs show that the time courses for the spraying and post-spraying phases can be quite different. In particular, it becomes clear that exposure to H2O2 aerosol reaches relevant concentrations only during the spraying phase. Table 3, panel 1 also reveals an increase in the aerosol concentration over the whole spraying phase, which can be explained by an increase in the backpressure (see panel 6), Starting with the aerosol concentration, the graphs show that the time courses for the spraying and post-spraying phases can be quite different. In particular, it becomes clear that exposure to H2O2 aerosol reaches relevant concentrations only during the spraying phase. Table 3, panel 1 also reveals an increase in the aerosol concentration over the whole spraying phase, which can be explained by an increase in the backpressure (see panel 6), Starting with the aerosol concentration, the graphs show that the time courses for the spraying and post-spraying phases can be quite different. In particular, it becomes clear that exposure to H 2 O 2 aerosol reaches relevant concentrations only during the spraying phase. Table 3, panel 1 also reveals an increase in the aerosol concentration over the whole spraying phase, which can be explained by an increase in the backpressure (see panel 6), leading to slower evaporation rates (panel 3) and longer droplet lifetimes respectively. This is also reflected in panels 2 and 3 showing the settling and evaporation rates of H 2 O and H 2 O 2 as components of aerosol droplets.
Although the aerosol droplets are more or less restricted to the spraying phase, evaporation from sprayed surfaces takes place over much longer time periods. The temporal pattern of evaporation is again determined by the volatility of the components and by the resulting back pressure. The influence of the activity coefficients is not so pronounced. Briefly, there is a time shift to be expected in the release of H 2 O 2 because H 2 O 2 is less volatile than water, which is confirmed in Table 3, panels 4 and 5. It should be highlighted that this delay is enhanced if there are large surfaces involved, leading to a significant backpressure that delays evaporation. This effect is also confirmed in panel 6, which shows lower vapour concentrations of H 2 O 2 at the beginning of the exposure period.
The delayed occurrence of the concentration peaks of H 2 O 2 suggests that the worker should leave the room at least after the application phase to minimize exposure as according to the Regulation (EU) No 528/2012 assessment report of a biocidal product (family) [30] hydrogen peroxide is assigned the hazard statement H332, 'harmful if inhaled'.

Comparison with Measured Data
Although this paper predominantly outlines the theoretical concept of the new approach, answering the question "How closely does the model approximate to real workplace exposures?" is essential as well. Since confirming observations increases the likelihood that the model is not fundamentally flawed, some model estimates were compared with measured data as part of this research. It proved to be difficult to find quality data for scenarios against which to compare SprayEva exposure predictions. However, Koch et al. [12] carried out a comparison of SprayExpo and ConsExpo 4.1 estimates with measurement data from a limited number of simulated spray experiments. The experiments were performed in a storage room (length: 11.2 m, width: 4.5 m, height: 4.0 m) under controlled conditions using different spraying devices and nozzles (flat fan nozzle, hollow cone nozzle, cold fogger, thermal fogger) characterized by different MMDs and release rates. The geometric standard deviation of the droplet distribution was fixed to 1.8 and was not considered by Koch [12] to be a parameter that is subject to large variation.
The applications were performed according to the same time scheme: • 5 min measurements prior to spraying; • 7 min spraying (4 × 1 min spraying with 1 min interruption in between); • 30 min follow-up measurements.
The spray droplets were released directly into the room without the interference of limiting surfaces, such that they could be simulated exactly with the program part 'room spray' in SprayExpo, ConsExpo and SprayEva.
In all these experiments, an aqueous formulation consisting of 0.2% dysprosium acetate and 1% sodium chloride was applied. While water is considered the volatile component (see Table 4), the inorganic dysprosium acetate and sodium chloride were assumed nonvolatile due to their extremely low vapour pressure. The data recorded by Respicon [31] and by an unspecified hygrometer were used for the analyses of the time aerosol concentration of the nonvolatile portion and the relative humidity (RH), respectively. After the termination of the experiment, the filters of the Respicon were analysed for the amount of trapped dysprosium, and this value was then extrapolated to calculate the total mass concentration of the nonevaporating ingredients. The measured relative humidity was used to calculate the vapour concentration of water in the room air. For comparison with model estimates, the time-weighted average over the spraying and post-spraying phases was calculated.
The input parameters such as release rate, temperature, initial humidity, and MMD as well as the individual measurements used for comparison with estimates from SprayEva are summarised in Table 4, where the inhalable aerosol and vapour concentrations are expressed in terms of the 37 min time-weighted average. It should be noted that the values for room temperature and initial humidity were provided via personal communication with Koch. Figure 5 shows the scatter plot of the Respicon measurements for the nonvolatile portion versus the model estimates of SprayExpo, ConsExpo 4.1, and SprayEva. The diagonal represents the 1:1 line, i.e., where the tool estimate is identical to the corresponding measurement value.
Looking at the deviations between the model values and the Respicon measurements, significant differences between the three models can be observed. For SprayExpo, underestimations are found only at low concentrations below 0.5 mg/m 3 . In contrast, SprayEva tends to overestimate exposure below 0.1 mg/m 3 . This partial lack of agreement is also confirmed numerically by the absolute and relative biases calculated according to Hornung [32], who defined bias as the mean difference between predicted estimates and measured exposure on a logarithmic scale. While the overall bias SprayExpo = −0.13 and the relative bias SprayExpo = −12.6% of SprayExpo were slightly negative, a positive bias SprayEva = 0.18 and relative bias SprayExpo = 19.6% was found for SprayEva. At the same time, Pearson correlation coefficients between the measurement results and tool estimates were r SprayExpo = 0.998 and r SrayEva = 0.992, demonstrating excellent correlation. Although the correlation coefficient of ConsExpo was excellent as well (r ConsExpo = 0.994), for some experiments, the model estimates deviated considerably from the measurement results. It seems that overestimation increases with decreasing concentrations. This visually apparent trend toward overestimation was also confirmed numerically by the overall bias (bias ConsExpo = 1.31) and the relative bias (rel.bias ConsExpo = 269.6%).  Looking at the deviations between the model values and the Respicon measurements, significant differences between the three models can be observed. For SprayExpo, underestimations are found only at low concentrations below 0.5 mg/m 3 . In contrast, SprayEva tends to overestimate exposure below 0.1 mg/m 3 . This partial lack of agreement is also confirmed numerically by the absolute and relative biases calculated according to Hornung [32], who defined bias as the mean difference between predicted estimates and measured exposure on a logarithmic scale. While the overall biasSprayExpo = −0.13 and the relative biasSprayExpo = −12.6% of SprayExpo were slightly negative, a positive biasSprayEva = 0.18 and relative biasSprayExpo = 19.6% was found for SprayEva. At the same time, Pearson correlation coefficients between the measurement results and tool estimates were rSprayExpo = 0.998 and rSrayEva = 0.992, demonstrating excellent correlation. Although the correlation coefficient of ConsExpo was excellent as well (rConsExpo = 0.994), for some experiments, the model estimates deviated considerably from the measurement results. It seems that overestimation increases with decreasing concentrations. This visually apparent trend toward overestimation was also confirmed numerically by the overall bias (biasConsExpo = 1.31) and the relative bias (rel.biasConsExpo = 269.6%).
While the experimental and measurement conditions were appropriate for the Respicon aerosol measurements, this was only partially the case for the water vapour concentration measurements. Firstly, the accuracy of the hygrometer used was comparatively low (+/−5%), and secondly, the atmospheric pressure was not measured; that is needed for the calculation of the mass concentrations. As a substitute, the mean atmospheric pressure of 101325 Pa was therefore used for the calculations. That is, weather-related changes in air pressure are not taken into account in the calculations. Figure 6 shows the scatter plot of the water vapour measurements versus the model estimates for SprayEva only; vapour exposure was out of the scope of the SprayExpo and ConsExpo spray models. While the experimental and measurement conditions were appropriate for the Respicon aerosol measurements, this was only partially the case for the water vapour concentration measurements. Firstly, the accuracy of the hygrometer used was comparatively low (+/−5%), and secondly, the atmospheric pressure was not measured; that is needed for the calculation of the mass concentrations. As a substitute, the mean atmospheric pressure of 101325 Pa was therefore used for the calculations. That is, weather-related changes in air pressure are not taken into account in the calculations. Figure 6 shows the scatter plot of the water vapour measurements versus the model estimates for SprayEva only; vapour exposure was out of the scope of the SprayExpo and ConsExpo spray models. Looking at the deviations between the model values and the measurement values, SprayEva tends to overestimate exposure (biasSprayEva = 0.14, rel.biasSprayEva = 15.5%). In addition, the correlation between the tool predictions of exposure and the measurement data is rSprayEva = 0.904, good but not excellent. This partial lack of agreement may be at least partly explained by the measurement limitations described above. Secondly, Looking at the deviations between the model values and the measurement values, SprayEva tends to overestimate exposure (bias SprayEva = 0.14, rel.bias SprayEva = 15.5%). In addition, the correlation between the tool predictions of exposure and the measurement data is r SprayEva = 0.904, good but not excellent. This partial lack of agreement may be at least partly explained by the measurement limitations described above. Secondly, SprayEva does not take into account the condensation of water vapour on surfaces, which may occur in reality. Taking these limitations into account, the results are nevertheless regarded as supportive and acceptable (due to the trend of conservative estimates) for the plausibility and predictive power of SprayEva.
Finally, it has to be noted that the simulation results presented in Table 4 were calculated using temporal resolutions ∆t in the range between 0.01 and 0.05 s for the spraying phase and 2 s for the post-spraying phase. The integer multiple IM for ∆t was set to 20, that is, the spatial resolution and the calculation time were reduced by a factor of 20, which proved to be a good compromise between accuracy and computation time.

Discussion
Due to the complex relationships between substance properties, process parameters, work practices, and ventilation, currently there is a lack of mathematical models to provide useful predictions of worker exposure for spraying operations, especially for the spraying of (semi)-volatile substances. That is, exposure modelling requires strong simplifications taking into account only the most influential input parameters. Since this inevitably leads to model uncertainties, there is always a trade-off between the accuracy of a model, the number and availability of the required input parameters, and the scope of application of the model. Because the availability of the input parameters is limited in the risk assessment procedure under the Biocide and REACH regulations, the main goal when developing SprayEva was to keep the model as simple as possible without losing too much predictive power and scope. For instance, we decided not to go beyond the well-mixed room model because published studies [18,[33][34][35][36] and our own evaluation exercise revealed reasonable agreement with the measured data as long as room sizes were small to medium. However, exposure estimates assuming well-mixed room conditions may be questionable for larger rooms. In this case, more sophisticated models such as two-box models should be used to avoid the underestimation of exposure. It should be noted that the iterative approach of SprayEva can be readily extended to more compartments, taking into account technical control measures [37] such as local exhaust ventilation or containment. Finally, this study was restricted to binary mixtures because they lead to much simpler equations, highlighting the dominant parameters for droplet evaporation. At the same time, they allow for the much quicker calculation of the droplet evaporation since there are fewer equations to be solved. However, it should be noted that according to Wilms [22], mixtures with more than two components can addressed with SRMM as well. This can be a topic for future work.
This study primarily outlines the theoretical concept of the model without claiming to provide a complete empirical, measurement-based evaluation of the room-and wallspraying scenarios covered by SprayEva. Unfortunately, suitable measurement data on vapour concentrations during wall-spraying scenarios are not available. Therefore, the data used for comparison with model estimates are restricted to the room spraying of an aqueous solution of nonvolatile components. Such exposure situations are within the scope of the proposed approach but are not typical for SprayEva, which mainly targets (semi)volatile mixtures. While the comparisons of the aerosol concentrations of the nonvolatile components with exposure estimates showed good agreement, experimental evidence for accurate model estimates of SprayEva for (semi)-volatile components is still subject to uncertainty. Although some measurement data for water vapour are available, the informative value of the data is limited due to the nonoptimal measurement conditions. In particular, the lack of measuring accuracy of the hygrometer used should be mentioned here. Nevertheless, the correlations of the measured data with the model estimates are good and are only moderately positively biased, which might be explained by water vapour condensation on surfaces, which is not addressed by SprayEva. Thus, the model tends to overestimate vapour exposure, which is desirable from a safety perspective.
The example calculations of the spraying of H 2 O and H 2 O 2 mixtures reveal that exposure to aerosol droplets is more or less restricted to the spraying phase. The evaporation from sprayed surfaces takes place over much longer time periods, as illustrated by panels 5 and 6 in Table 3, where the temporal pattern of evaporation is determined by the volatility of the components. In brief, there is a delay to be expected because H 2 O 2 is less volatile than water, which causes water to evaporate first. The example calculations also show that evaporation from large surfaces creates a backpressure that in turn affects the time course of air concentrations by significantly delaying the release of substances. This suggests that the worker should leave the room after the spraying phase to minimize exposure. As already addressed in previous work [13], the nonideal behaviour of components should be considered in any case. Although nonideal behaviour is not very pronounced in the case of H 2 O 2 /H 2 O mixtures, it can play an important role with mixtures whose components differ greatly in polarity. Therefore, when using SprayEva, exposure assessors should be aware of deviations from ideal behaviour because the delayed or premature release of hazardous substances can determine the risks workers are facing.
Although the proposed mass balance equations in combination with the iterative approach provides a very flexible framework for modelling of exposure scenarios, there are still limitations the assessor has to consider when using SprayEva.
The strong simplification regarding the dispersion of the spray cloud by instantaneous diffusion differs heavily from the actual physical dispersion behaviour, in particular in high or very large rooms. In these settings, a concentrated particle cloud is initially created close to the source, spreading only very gradually through the whole room. Therefore, special conditions such as overhead spraying are not further differentiated in the well mixed room model. In addition, the impaction module is at present restricted to wall spraying and does not include spraying onto floor surfaces, where sedimentation and impaction contribute to the overall mass on the floor, making the mass balance difficult.
Furthermore, the influence of chemical reactions on the concentration course is not yet taken into account, which can play a significant role in practice (e.g., reaction of H 2 O 2 with biofilms). Although these effects can in principle be considered in the mass balance to some extent, the corresponding kinetics are regarded as complex and might be a topic for future research. This also applies to more sophisticated dispersion models, such as two-box models or models that can take into account additional technical control measures, such as local exhaust ventilation or containments. For instance, Ganser and Hewett [37,38] demonstrated that one-and two-box models can be extended to situations where various forms of a local exhaust ventilation are involved.
Finally, it has to be noted that calculations with the current code in the Mathematica ® [39] environment can be time consuming since the computing time increases quadratically with the number of iteration steps. For the post-spraying phase, computing speed is less demanding because the evaporation process from flat surfaces is slower than for droplets, allowing for longer discretisation intervals. Computer capacity and code optimization is therefore a crucial factor for reducing the computing time. First attempts in this direction using an R-environment showed that calculation time could be reduced by more than an order of magnitude.

Conclusions
An important finding from the example calculations and from the evaluation exercise is that the proposed iterative algorithm is very flexible, covers a wide range of scenarios including binary (semi)-volatile mixtures, and allows variable application patterns in roomand wall-spraying scenarios.
The example results suggest that situations where high evaporation rates from large surfaces and hence backpressure occur can influence the time course of airborne concentrations by significantly delaying substance release.
This also applies for mixtures where liquid-phase nonidealities play a role. Although nonideal behaviour is often less pronounced in practice, mixtures whose components differ greatly in polarity should be paid special attention, as deviations from ideal can lead to the delayed or early release of components. Both deviations from ideal and backpressure effects can significantly change the time course of exposure, leading to exposure after the original application is finished and thus possibly requiring the adjustment of risk management measures.
Finally, it has to be noted that the relationships between substance properties, process parameters, work practices, and ventilation are complex and require strong simplifications in model building. Therefore, exposure estimates are inevitably uncertain to some extent. Nevertheless, the small-scale evaluation exercise showed good to excellent correlations between measured data and model estimates and only a small to moderate tendency towards overestimation. Thus, the results are regarded as supportive for the predictive power of SprayEva and acceptable from a safety perspective. Promising applications of the new model are, for example, spray applications of disinfection and cleaning products as they often contain volatile and (semi)-volatile active substances. Data Availability Statement: Data is contained within the article. The data presented in this study are available in [12].

Conflicts of Interest:
The authors declare no conflict of interest. The authors of this manuscript work for the unit responsible for occupational exposure assessment in the German Competent Authority for REACH and biocidal products regulation.