key: cord-0817418-k6ovqpsn
authors: Cai, Hao; Li, Xianting; Chen, Zhilong; Wang, Mingyang
title: Rapid identification of multiple constantly-released contaminant sources in indoor environments with unknown release time
date: 2014-06-17
journal: Build Environ
DOI: 10.1016/j.buildenv.2014.06.006
sha: 8ffe2548d4867c792ea343f9f2536bdea1f29d57
doc_id: 817418
cord_uid: k6ovqpsn

The sudden release of airborne hazardous contaminants in an indoor environment can potentially lead to severe disasters, such as the spread of toxic gases, fire, and explosion. To prevent and mitigate these disasters it is critical to rapidly and accurately identify the characteristics of the contaminant sources. Although remarkable achievements have been made in identifying a single indoor contaminant source in recent years, the issues related to multiple contaminant sources are still challenging. This study presents a method for identifying the exact locations, emission rates, and release time of multiple indoor contaminant sources simultaneously released at constant rates, by considering sensor thresholds and measurement errors. The method uses a two-stage procedure for rapid source identification. Before the release of contaminants, only a limited number of time-consuming computational fluid dynamics (CFD) simulations need to be conducted. After the release of contaminants, the method can be executed in real-time. Through case studies in a three-dimensional office the method was numerically demonstrated and validated, and the results show that the method is effective and feasible. The effects of sensor threshold, measurement error and total sampling time on the source identification performance were analysed, and the limitations and applicability of the method were also discussed.

Indoor environments in modern society are persistently threatened by a huge variety of airborne hazardous contaminants, such as toxic, inflammable, or explosive chemicals, harmful microbes, and radioactive substances [1e4] . The sudden release of these contaminants, accidentally or intentionally, can potentially result in exposure to toxic gases, fire, explosion, epidemic, or radiation. All these incidents would lead to significant casualties and property losses if no effective response measures were taken. In recent years, several tragic incidents, such as the Tokyo subway sarin attack, severe acute respiratory syndrome outbreak and the Fukushima nuclear leak, have aroused the concern of indoor environmental safety. To prevent and mitigate disasters due to the sudden release of hazardous contaminants, it is critical to promptly and accurately identify the characteristics of the contaminant sources, such as locations, emission rates, and release time. Source identification is an important premise for developing effective response measures, such as source elimination, dilution ventilation, air purification, evacuation, and shelter-in-place [5, 6] .

With known locations and emission rates of contaminant sources, the transient distribution of a contaminant can be well predicted using a dispersion model. However, in most practical applications the information on contaminant sources is incomplete or unknown, and needs to be identified using the sensor measurements of the airflow and concentration fields. This type of identification constitutes an inverse problem, which has been widely studied for several decades in the fields of heat transfer [7e9], groundwater contamination [10, 11] , soil pollution [12, 13] , and atmospheric constituent transport [14e16] .

Due to increasing concern about indoor environmental safety, research on identifying indoor contaminant sources has developed rapidly in recent years. The methods available for identifying indoor contaminant sources can be roughly categorised as backward and forward methods, as shown in Table 1 . To date, most backward and forward methods have shown their applicability in identifying a single indoor contaminant source [17e40] . In contrast, very few studies have discussed the scenarios involving multiple indoor contaminant sources, which are frequent and common in real-world applications [18, 38] . In residential or working environments, chemical contaminants are ubiquitous and continuously emitted from building materials, furniture, and equipment. In emergency situations related to earthquakes, blasts or impacts, hazardous chemicals may simultaneously escape from many leakage points. During an epidemic outbreak, many infected people may be found in hospital or other buildings. In terrorist attacks, chemical or biological agents may be deliberately released from many locations. Therefore, research on identifying multiple indoor contaminant sources is of great practical significance.

For identifying multiple indoor contaminant sources, a major challenge to the backward methods is how to ensure the uniqueness of the identification results. When multiple contaminant sources are released, different combinations of sources may result in the same or close sensor measurements at a given time point. In other words, different causes may lead to the same outcome. Because the backward methods use the same outcome (sensor measurements at a given time point) as the initial condition, it is difficult to differentiate various possible causes. Besides, the backward methods would require a significant reduction in computing time for solving inverse computational fluid dynamics (CFD) models. This potentially limits their application in real-time source identification.

As shown in Table 1 , forward methods use sensor measurements during a given time period, rather than at a single moment, as inputs. In indoor environments, it is quite rare that different release scenarios will produce the same or close sensor measurements at all time points. Therefore forward methods are more promising for finding a unique solution. In addition, the two-stage procedure of the forward methods (see Table 1 ) also makes real-time source identification possible. A major constraint of the forward methods is that most of these methods only work well using a fast contaminant transport model, such as a multi-zone model, because a large number of simulations are usually required before a release incident to cover all possible scenarios. The multi-zone model generally assumes a uniform contaminant distribution in one zone [41, 42] . Therefore the forward methods using the multi-zone model cannot identify the exact locations of contaminant sources in each zone or room.

In a previous study we presented a theoretical model for quickly identifying the exact locations and emission rates of multiple constant contaminant sources indoors, using ideal sensors [38] . By using an analytical expression of indoor contaminant dispersion, the theoretical model only requires a limited number of CFD simulations. Because CFD simulations can provide detailed information of contaminant dispersion, the exact locations of contaminant sources can be identified by the theoretical model. The major constraint of this model is that it cannot work by using real sensors, mainly because the model cannot work when there is no response of sensors due to sensor thresholds. In real-world situations, the use of real sensors will make the problem of source identification much more difficult, primarily due to three issues: (1) no response of sensors due to sensor thresholds; (2) inevitable measurement errors; and (3) unknown release time of contaminant. This research aims to fill the gap by developing a more sophisticated method that can promptly identify the locations, emission rates, and release time of multiple indoor contaminant sources simultaneously released at constant rates, by considering the sensor threshold and measurement errors.

The source identification problem is specified by introducing the following assumptions. Assumption 1. The contaminant is dispersed as a passive gas in steady indoor airflow. The passive gas refers to the gaseous contaminant that is present in sufficiently low concentration that the effects of its dispersion on the indoor airflow field and air density can be neglected. Assumption 2. The number of potential sources is limited and their locations are known. The current work does not attempt to Optimization methods based on a limited number of CFD simulations Cai et al. [38e40] handle the situations where the number of the potential sources is unlimited and their locations are completely unknown. Assumption 3. The emission rates of contaminant sources are constant. In practice, a contaminant may be released continuously or instantaneously. The current study only aims to identify indoor contaminant sources that are continuously released at constant rates.

Assumption 4. All contaminant sources are released simultaneously. That is, the current work only deals with multiple contaminant sources in synchronous release.

The overall framework of the identification method is depicted in Fig. 1 . The method is composed of three essential components. Component I is an analytical expression of transient contaminant dispersion, which can predict the temporal-spatial distribution of a contaminant in real-time with known source information, such as source locations, emission rates, and release time. Component II is a source identification model that can identify the locations and emission rates of multiple contaminant sources with a known release time. Component III is an algorithm for identifying the release time of contaminant sources.

In our previous work, an analytical expression of transient contaminant dispersion was presented to predict the transient distribution of contaminants in real-time [43] . For the dispersion of a passive gas in a steady-state airflow field, the analytical expression can be expressed as:

where C p (t) is the contaminant concentration (mg/m 3 ) at an arbitrary indoor point p and moment t, C S,k is the contaminant concentration (mg/m 3 ) at the kth inlet, C 0 is the initial contaminant concentration (mg/m 3 ) indoors, S i is the emission rate (mg/s) of the ith contaminant source, Q is the airflow rate (m 3 /s), K is the number of indoor supply air inlets, N is the number of potential indoor contaminant sources, a Sk,p (t) is the transient accessibility of supply air (TASA) from the kth inlet to p at moment t, and a Ci,p is the transient accessibility of contaminant source (TACS) from the ith source to p at t. The TASA and TACS indices are dimensionless numbers that quantify the effects of the supply air and contaminant source on an indoor position at different moments, respectively. TASA is a function of flow characteristics, regardless of contaminant type and source. In a ventilated space, assume that only air from the kth inlet contains a hypothetical contaminant, and there is no source of the contaminant in the space. The TASA from the kth inlet to point p at moment t was defined as [43] :

TACS is a function of both the flow characteristics and source location, regardless of the emission rate and contaminant type. Assume that the initial concentration was zero, all inlet concentrations were zero, and only the ith contaminant source existed. The TACS from the ith source to p at t was defined as [43] : 

where C e,i is the average exhausted concentration (mg/m 3 ) under steady-state conditions. Using CFD simulations, the temporal-spatial distributions of TASA and TACS from each supply air inlet and contaminant source can be obtained in advance. With the TASA and TACS data, transient contaminant distributions under different combinations of supply air concentrations and source emission rates can be calculated in real-time through the simple algebraic calculations in Eq. (1). This feature of the analytical expression enables the real-time simulation of contaminant dispersion.

When an initial concentration and all supply inlet concentrations are zero, Eq. (1) can be reduced to:

Assuming the identified emission rate of the ith contaminant source is S * i , then

where C m p ðtÞ is the measurement of sensor (mg/m 3 ) at p and t; C * p ðtÞ is the calculated concentration (mg/m 3 ) at p and t by substituting S * i into Eq. (4); and e is the difference between C m p ðtÞ and C * p ðtÞ (mg/ m 3 ). Ideally, e is equal or close to zero.

Using Eq. (5), the emission rate of each potential source can be identified by solving a linear programming model, which can be expressed in a canonical form as follows:

where f(x) is called the objective function, Ax ¼ b and x ! 0 are called the constraint functions, x represents the vector of variables to be determined, c and b are vectors of known coefficients, A is a known matrix of coefficients, and (,) T is the matrix transpose. The linear programming model identifies the emission rate of each potential source by minimizing the objective function with constraint functions. The objective function is the sum of absolute errors between concentrations calculated by Eq. (4) and those detected by sensors. The constraint functions are derived from Eq. (5). The detailed derivation procedure of the linear programming model is given in Appendix section. After solving the linear programming model, the emission rate of each potential source can be determined, and then the number, locations, and emission rates of indoor sources can be identified, because it is assumed that the locations of all potential sources are known. A major constraint of this model is that it can only work with a known release time of contaminant sources. In practice, the release time is usually unknown and must be identified. The linear programming model excludes negative emission rates by using constraints x ! 0. Ideally, any potential source with an identified emission rate greater than zero can be thought as a source. However, the identification results will inevitably contain some errors in practice because of the errors of numerical simulations and sensor measurements. If the identified emission rate of a potential source is a tiny number, a quantitative criterion is required to help determine whether the potential source actually release any contaminant.

In a ventilated space, if only the ith contaminant source existed and all the inlets concentrations would be zero, then the average exhausted contaminant concentration under steady-state conditions is:

where S i is the emission rate of the ith potential source (mg/m 3 ). According to Eq. (7), a quantitative criterion is proposed as:

where S thd is the emission rate threshold (mg/s) under which the release of a source is suggested to be neglected; z is a factor far less than 1 (e.g., z ¼ 0.01); C thd is the concentration threshold of the contaminant (mg/m 3 ) above which a threat to people or property may exist during a specific time period. The emission rate threshold S thd is a function of C thd and Q, and is significantly scaled down by a small number z. If the identified emission rate of a potential source is less than S thd , it suggests that the release of the potential source can be neglected because it alone will not cause a severe threat to people or property. For most contaminants, the value of C thd can be determined by referring to public exposure guidelines, such as AEGLs (Acute Exposure Guideline Levels), ERPGs (Emergency Response Planning Guidelines), and TEELs (Temporary Emergency Exposure Limits) [44, 45] .

In some special cases, if the value of C thd cannot be easily determined, it somewhat makes sense to use a number small enough as the emission rate threshold. The reason is that if the identified emission rate is small enough, whether or not considering this source will not significantly affect the concentration field and corresponding response measures. In addition, the emission rate threshold can also be set as zero for conservative reasons.

The release time of a contaminant is denoted by t r , which refers to the time when the contaminant release event starts. After the contaminant release, the sensor that first detects the contaminant is called the first responding sensor (FRS), and the time when the contaminant is detected by the FRS is called the alarming time of the FRS and denoted by t d . The time interval between t r and t d is called the lag time of the FRS, and is expressed as t FRS ¼ t d À t r , as shown in the following figure. Fig. 2 In practice, t d is known, while t r and t FRS are to be determined. The presented algorithm begins by assuming t FRS ¼ 0 s, and continues by increasing t FRS at a time step Dt FRS . As shown in Fig. 1 , with an assumed release time and the measurements of all sensors over a period t after t d , the contaminant concentrations at the FRS can be predicted by sequentially executing Components II and I. An index called SLT is proposed to quantify the match degree between the predicted and measured concentrations at the FRS over period t.

The SLT index is defined as: Fig. 2 . Timeline of the contaminant release event.

where SLT j is the SLT index (m 3 /mg) when t FRS ¼ j  Dt FRS ; C 1,k is the kth measured concentration at the FRS (mg/m 3 ); C j 2;k is the kth predicted concentration at the FRS when t FRS ¼ j  Dt FRS (mg/m 3 );

and N FRS is the number of measured concentrations at the FRS. The denominator of Eq. (9) is essentially a widely used concept, namely the distance between two data sets. A smaller distance between the predicted and measured data will lead to a larger value of SLT, which further indicates a higher match degree between the two types of data. In practice, the errors in numerical simulations and sensor measurements are inevitable, and thus the distance between the predicted and measured data is always greater than zero, which ensures that the definition of SLT is always meaningful.

By increasing t FRS continuously, and running Components II and I iteratively, a curve of SLT j changing with t FRS can be plotted. The peak point of the curve indicates a best match between the two data sets. With the SLT curve, the release time t r can be identified by finding the peak point of the curve. For example, if the peak point corresponds to a lag time l  Dt 1 , then the identified release time was t r ¼ t d À l  Dt 1 . Furthermore, the source locations and emission rates can be determined, which were the outputs of Component II using t FRS ¼ l  Dt 1 as the assumed lag time of the FRS.

In practice, the presented identification method operates in a two-stage procedure, as illustrated in Fig. 3 . In the first stage, before contaminant release, CFD is used to simulate a limited number of benchmark scenarios of contaminant dispersion. In each scenario only one of the potential sources is released at a constant emission rate. Therefore the number of CFD simulations is equal to the number of potential sources. After the simulations, the TACS data from each potential source to each sensor are calculated and stored in a database.

In the second stage, after contaminant release, the source locations, emission rates and release time are identified using the measurements of all sensors over a specific period after the release time as the inputs. In this stage, Components I and II are operated iteratively to identify the release time of the contaminant. Nevertheless, this stage will be completed in a few seconds, because Components I and II are very fast. In the identification procedure, the analytical expression (Component II) plays an important role in reducing the computational load. For an indoor space with a limited number of potential contaminant sources, there are theoretically an infinite number of contaminant dispersion scenarios, by considering all the possible combinations of source parameters, such as source numbers, source locations, and emission rates. Because of the analytical expression, only a limited number of time-consuming CFD simulations are required for covering all the possible scenarios, which makes it possible to use CFD in the forward methods of source identification.

A three-dimensional office was used to validate the presented identification method (Fig.  4) . The office (X  Y  Z ¼ 9.6 m  3.2 m  5 m) was ventilated by a mixing ventilation (MV) system and a personalized ventilation (PV) system. The MV system was used to control the temperature, humidity and contaminant concentration of the whole room. The air supplied by the MV system entered the room from two supply air inlets (0.4 m  0.4 m) at 16 C and 0.4 m/s and exited from an exhaust air outlet (0.8 m  0.4 m). The air exchange rate was 3 ACH (air changes per hour) or 0.128 m 3 /s. The PV system was used to supply fresh air directly to the six people in the office through six table-mounted supply air terminal devices (ATDs). For each person sitting behind a table, there was an ATD that is close to his/her head. In the office, each computer, person, lamp, and the single window contributed 108, 75, 34, and 220 W heat to the room, respectively. The four walls, ceiling, and floor were set as adiabatic boundaries.

It is assumed that a toxic gas was intentionally released into the PV system in a terrorist attack. The toxic gas could potentially be delivered to the breath zone of the indoor people from one or more of the six ATDs. In this case, all the ATDs were potential contaminant sources, and denoted as CS1eCS6. Five sensors (SR1eSR5) were installed to identify the contaminated ATDs that were releasing the toxic gas (Fig. 4) . The positions of the potential sources and sensors are listed in Table 2 . A toxic gas spreading scenario in which the toxic gas was released from CS3 and CS5 at 50 mg/s was studied to test the presented identification method.

A commercial CFD program, Airpak, was used to simulate the indoor airflow field and contaminant dispersion. Airpak is customised from a general-purpose program, Fluent, for indoor environment simulations [46] . Airpak has been widely used and validated in indoor airflow and contaminant dispersion studies, as reported by Xu [47] .

An indoor zero-equation turbulence model was used to simulate the indoor turbulent flow [48] . An ideal gas law model was used to account for temperature-dependent density effects, and a surfaceto-surface radiation model was used to calculate radiative heat transfers [46] . For simplicity and illustration purpose, each contaminated ATD was modelled as a contaminant source, rather than a supply air inlet, in CFD simulations. After the test of mesh independence, the office room was divided into 121,038 hexahedral control volumes. The momentum equations were solved on nonuniform staggered grids using a SIMPLE algorithm [49] . A secondorder upwind scheme and linear under-relaxation iteration was used in the solution process.

In the presented identification method, Component II identifies the emission rate of each potential contaminant source by solving a linear programming model (Eq. (6)). Linear programming was first presented by Dantzig in 1947 to refer to the optimization problems in which both the objective function and the constraints are linear [50, 51] . After decades of development, linear programming has become mature and widely used. The most popular algorithm in linear programming is called the simplex algorithm [51] . In this study, the simplex algorithm provided by Matlab optimization toolbox 2.3 was used to solve the linear programming model.

The validation procedure includes four major steps.

Step I: prepare TACS data First the steady-state airflow field indoors was simulated using CFD. Next, with the same airflow field, six benchmark scenarios of contaminant dispersion were simulated using CFD. In each scenario, only one of the six potential sources ( Table 2 ) was released at a constant rate. Finally, using the concentration data simulated, the TACS data from each potential source to each sensor were calculated by Eq. (3).

Step II: simulate a toxic gas spreading scenario for testing the identification method A scenario in which CS3 and CS5 were simultaneously released was simulated using CFD. The simulated concentrations at the positions of five sensors (SR1eSR5) were sampled directly. These concentration data are called "exact measurements".

Step III: identify the contaminant sources To mimic the sensor measurements in a real-world experiment, the "exact measurements" were modified into "perturbed measurements" using a data pre-processing procedure. First, the exact measurements below a specific value, the sensor threshold, are eliminated. Next, the "exact measurements" were added with normally distributed random errors using a perturbance scheme, which is described in Section 4.4. The "perturbed measurements" were further taken as the inputs of the presented identification method to identify the source locations, emission rates and release time.

Step IV: evaluate the source identification performance

The source identification performance can be evaluated by comparing the identified results with those predefined in Step II. A series of indices, the scale of relative errors (SRE), were introduced to quantify the identification performance. The definitions of SRE indices are provided in Section 4.4.

The perturbance scheme used to generate random errors was defined as:

where C pert (t) and C exct (t) are perturbed and exact measurements (mg/m 3 ) at time t, respectively (mg/m 3 ); ε r is a random variable that follows a normal distribution with mean ¼ zero and standard

x is a fraction to adjust the level of measurement error (0 x 1); STD is the standard normal deviation.

Eq. (10) indicates that all concentrations were measured by sensors with the same relative error; i.e. a larger concentration may have a greater measurement error. To reflect the level of measurement errors more clearly, two parameters, d 1 and d 2 , were defined as:

where d 1 and d 2 are the maximum and average absolute relative error between perturbed and exact measurements, respectively; C pert and C exct,i are the ith perturbed and exact measurements, respectively; and M is the number of exact measurements. The SRE indices for performance evaluation were defined as follows [38] :

where S * i and S i are the identified and actual emission rates (mg/s) of the ith source, respectively; and N is the number of potential sources. SRE1 i indicates the absolute relative error between the identified and actual emission rate of the ith source. A lower value of SRE1 i indicates the more accurate identification of the ith source.

SRE2 indicates the average level of SRE1 i (1 i N). A lower value of SRE2 indicates the more accurate identification of all potential sources with respect to the overall average.

SRE3 indicates the maximum level of SRE1 i (1 i N) . A lower value of SRE3 indicates the more accurate identification of all potential sources with respect to the maximum relative error.

The steady-state airflow field indoors was first simulated, which shows a typical distribution pattern of mixing ventilation, as shown in Fig. 5 . With this airflow field, six benchmark scenarios of contaminant dispersion were further simulated over 180 s with a time step 0.5 s. In each scenario, only one of the six potential sources ( Table 2 ) was released at a constant rate. By using the simulated concentrations and Eq. (3), the TACS from each potential source to each sensor was calculated. Fig. 6 presents an example of TACS variation, from potential source CS3 to each sensor. The TACS curves from CS3 to SR3 and SR4 varied in magnitude and shape, while the curves from CS3 to SR1, SR2, and SR5 overlapped one another and had much lower values. This indicates that SR3 and SR4 can more easily detect CS3 than the rest sensors.

The toxic gas spreading scenario for testing the presented identification method was simulated over 180 s with a time step of 0.5 s. In this scenario, CS3 and CS5 were simultaneously released at a constant rate of 50 mg/s Fig. 7 shows the variation in concentration distribution in this scenario. Using the simulated concentrations, the lag times of each sensor with different sensor thresholds were obtained directly, and are summarised in Table 3 . For this scenario the FRS was sensor SR3 at sensor thresholds of 1.0, 5.0, and 10.0 mg/m 3 , and its lag times were 11, 14, and 15 s respectively.

The exact measurements at the three sensor threshold levels, 1.0, 5.0, and 10.0 mg/m 3 , were used to test the presented identification method. In each test, the total sampling time of the FRS was 150 s, the assumed lag time of the FRS was increased from 0 s at a step of 1 s, and the sampling interval of each sensor was 5 s. The measurements of all sensors were taken as the inputs to the identification method. Table 4 summarises the number of measurements from each sensor at different sensor thresholds. It is shown that the total number of measurements from all sensors decreased with the increase in sensor threshold. Fig. 8 presents the SLT curves calculated using exact measurements at different sensor thresholds. It is shown that each curve has a sharp single peak. By finding the horizontal coordinate value of each peak point, the lag time of the FRS at each sensor threshold was identified and listed in Table 5 . It is shown that the identified lag time of the FRS was identical to the actual ones for each sensor threshold.

With the identified lag times, the emission rate of each potential source was further determined and listed in Table 5 . Two indices, SRE2 and SRE3, were provided in Table 5 to help evaluate the identification performance. The results indicate that the identified emission rates were very accurate, even for the worst identification results corresponding to the sensor threshold 10.0 mg/m 3 . A possible explanation for the high accuracy is the use of exact measurements. The identification accuracy of emission rates also presented a declining trend with the increase in sensor threshold. This may be due to the fact that fewer measurements were used in source identification at a higher sensor threshold (see Table 4 ).

The identification method was further tested with two levels of measurement error, which were generated using the perturbance scheme described in Section 4.4. As an example, Fig. 9 shows the exact, highly perturbed, and lowly perturbed measurements of the FRS (SR3). According to Eqs. (10)e (12) , the level of measurement error was determined by x, and characterised by d 1 and d 2 . For the highly perturbed measurements, x, d 1 and d 2 were 0.1, 27.1% and 7.4%, respectively; whereas for lowly perturbed measurements these parameters were 0.01, 3.3% and 1.0%, respectively. Fig. 10 presents the SLT curves using highly perturbed measurements at different sensor thresholds. The SLT curves in Fig. 10 were much flatter than those in Fig. 8 . Moreover, multiple peak points appeared in the curve corresponding to the sensor threshold of 10.0 mg/m 3 . A comparison between Figs. 10 and 8 indicates that the measurement errors would decrease the accuracy of identifying the release time of the contaminant. Table 6 summarises the identification results using highly perturbed measurements. Compared with the identification results using exact measurements (Table 5) , when highly perturbed measurements were used the identified lag times of the FRS were generally less accurate, and the accuracy of identified emission rates significantly decreased. These results indicate that the measurement errors can significantly degrade the identification performance.

The identification method was further tested with lowly perturbed measurements. The SLT curves are plotted in Fig. 11 and the identification results summarised in Table 7 . By decreasing the level of measurement error, the SLT curves became steeper and more symmetrical, and each curve only had a single peak. The identified lag times of the FRS were generally more accurate, and the accuracy of the identified emission rates increased significantly. These results indicate that the identification performance can be improved by reducing the level of measurement errors.

In contrast to the results using exact measurements (Table 5) , the results using perturbed measurements (Tables 6 and 7) did not present the trend that the identification accuracy will decline with an increase in sensor threshold. In fact, the results corresponding to a sensor threshold of 10.0 mg/m 3 in Tables 6 and 7 were unexpectedly good. The difference between the results may be explained from two aspects. Firstly, the measurement errors were generated in a random manner, which resulted in the randomness of the identification results. Secondly, in the early stage of contaminant dispersion, the concentrations are normally low and changing rapidly, which can lead to large errors in the simulation results. The high sensor threshold may be helpful for eliminating these simulation errors.

The identification method was also tested by changing the total sampling times of the FRS. Table 8 summarises the total sampling time of each sensor. It is shown that one or many sensors were unable to provide concentration data when the total sampling time of the FRS was less than 150 s. Table 9 summarises the identification results using lowly perturbed measurements. The threshold and sampling interval of each a The toxic gas was actually released from CS3 and CS5 at 50 mg/s. First responding sensor (FRS) used a total sampling time of 150 s. Sampling interval of each sensor was 5 s. b SRE2, and SRE3 are performance indexes of identification results, referring to Eqs. (14) and (15). (14) and (15). sensor were 5 mg/m 3 and 5 s, respectively. The FRS was SR3, and its actual lag time was 14 s (see Table 3 ). When the total sampling time of the FRS was less than 50 s, the identification method failed to obtain accurate results. The failure of the model may be explained by two reasons. Firstly, except the measurements from the FRS, there are no measurements from the other sensors, or the measurements from the other sensors are not enough. Secondly, the measurements from the other sensors are probably important to the source identification. When the total sampling time of the FRS was equal to or greater than 50 s, all identification results, except that corresponding to 90 s, were very accurate and varied in a narrow range. As shown in Table 8 , with the increase in total sampling time of the FRS, more sensors provided concentration data for source identification. These results indicate that extending the total sampling time of the FRS or using measurements from more sensors can contribute to improving the identification performance. In practice, the source identification method can be executed repeatedly by increasing the total sampling time of the FRS. A minimum or appropriate total sampling time of the FRS can be determined when the identification results tend to be stable. The results corresponding to 90 s were much worse than expected. This may be because the measurement errors were randomly generated.

All calculations reported were conducted on a personal computer (CPU: Intel(R) Core(TM)2 T7200 @ 2.00 GHz). As shown in Tables 5e7 and 9, the computation times for each identification procedure in this study were only a few seconds. Therefore, the presented identification method is fast enough to be integrated into a real-time identification system.

The identification method presented was developed based on four assumptions. Recognition of the underlying assumptions is thus important to understand the limitations and applicability of the method. In the following, the underlying assumptions are analysed, and the limitations and applicability of the method are discussed.

The first assumption is that the contaminant can be treated as a passive gas and the indoor airflow field is in steady-state. For most indoor environments equipped with mechanical ventilation, the airflow field can reach a steady-state much faster than the contaminant dispersion, and thus the airflow field can be regarded as steady during the dispersion of contaminant. Normally, indoor airflow is turbulent and the contaminant concentration is low enough that the contaminant dispersion only has a trivial effect on the airflow field and air density [43] . This assumption is, therefore, reasonable for most indoor environments. For quiescent environments, the presented method may also be applicable. The reason is that a static environment can also be thought with fixed airflow mode, in which pure molecular diffusion transports contaminant. It is worth noting that, the presented method is not applicable to the situations where the dispersion of contaminant is significantly different from that of passive gas, such as the releases of particles, aerosols, and a massive leakage of dense or light gas. In addition, the method cannot be applied to the indoor environments where the airflow is always changing, such as the naturally ventilated environments.

This work also assumes that there are a limited number of potential sources and their locations are known. This assumption excludes the application of the method in situations where the number of the potential sources is unlimited and their locations are completely unknown. Nevertheless, a variety of contaminant dispersion scenarios that satisfy this assumption can still be found in practice, such as the emission of volatile organic compounds from furniture, leakage of toxic or inflammable gas from suspected leak points, viruses spreading from infected people at specific locations, toxic agents released from supply air inlets by terrorists.

The third assumption is that the contaminant sources are continuously released at constant rates. According to this Fig. 11 . Variation in SLT with assumed lag time of the FRS using lowly perturbed measurements at different sensor thresholds. (14) and (15). )  SR1  SR2  SR4  SR5   30  0  0  0  0  40  0  0  0  1  50  0  0  1  11  60  0  0  11  21  90  0  0  41  51  120  0  8  71  81  150  28  48  101  111 a Sensors with no measurement were highlighted in grey.

assumption, the presented method is applicable to constantlyreleased contaminant sources. In addition, the method may also be used for the situations where the emission rates of the contaminant sources are slowly changing, provided that the change over time is trivial and negligible during the short period of source identification. This assumption also indicates that the presented method is not applicable to the situations where the emission rates of the contaminant sources are instantaneous or rapidly changing. Finally, to simplify the source identification problem, it is assumed that all contaminant sources are released at the same time. This assumption implies that the presented method can only be used to identify multiple contaminant sources in synchronous release. For identifying multiple asynchronous releases, further research should be conducted in the future.

This study presents a method to quickly identify the exact locations, emission rates, and release time of multiple indoor contaminant sources simultaneously released at constant rates, by considering sensor thresholds and measurement errors. The method was numerically illustrated and validated through case studies, and the results show that the method is effective and feasible. The case studies indicate that the performance of the method was closely related to the sensor threshold, measurement error, and total sampling time. When the levels of the sensor threshold and measurement error were given, the total sampling time of the FRS should exceed a specific value to ensure the success of source identification. Besides, the accuracy of source identification can be improved by decreasing the levels of the sensor threshold and measurement error, extending the total sampling time of the sensors, and using measurements from more sensors at different positions. The case studies also show that the method is promising for integration into an online identification system, because of its capability for rapid identification.

The proposed method may be useful for preventing and mitigating indoor environmental disasters due to the sudden release of hazardous gas, such as industrial accidents, terrorist attacks, and epidemic outbreaks. Further validations of the method in a laboratory chamber and real indoor environments are ongoing and will be reported in the near future. Some major limitations of the method are also worth attention. Firstly, the method cannot be applied to situations where the dispersion of contaminant is significantly different from that of passive gas, such as the releases of particles, aerosols, and a massive leakage of dense or light gas. Secondly, the method is not applicable for identifying instantaneous releases or continuous releases with fast changing rates. In addition, the method can only be used to identify multiple contaminant sources that are simultaneously released. Therefore, further research is required to develop more sophisticated techniques that can be applied to more complex situations of contaminant dispersion. 

where a i;j ¼ a j Ci =Q , b j ¼ C m j , x i ¼ S * i , i2 [1,N] , and j2 [1,M] , and C m j is the jth measurement from a specific sensor at a specific moment. a j Ci is the TACS from the ith source to a specific sensor at the specific moment when C m j was measured. In Eq. (16) , if the identified emission rate of each potential source x i approaches or is equal to its actual value, P M i¼1 je i j will be close to zero or a minimum. Therefore, x i or S * i can be identified using a nonlinear programming model, as follows: a The toxic gas was actually released from CS3 and CS5 at 50 mg/s. Threshold and sampling interval of each sensor were 5 mg/m 3 and 5 s, respectively. The actual lag time of FRS was 14 s, and lowly perturbed measurements were used. b SRE2, and SRE3 are performance indexes of identification results, referring to Eqs. (14) and (15) .

To reduce the solution complexity, the nonlinear programming model was further transformed into a linear programming model, as follows:

; (18) where x ¼ ðx 1 ; x 2 ; /; x N ; x 0 Nþ1 ; /x 0 NþM ; x 00 Nþ1 ; /x 00 NþM Þ T , The objective function of the linear programming model is the sum of absolute errors between concentrations calculated by Eq. (4) and those measured by sensors. By minimising the objective function, the emission rate of each potential source can be solved.

CFD analysis of dense gas dispersion in indoor environment for risk assessment and risk mitigation

Numerical simulation of combustible gas diffusion in the confined room with obstacles

Indoor air quality and health

Assessment of external gamma exposure and radon levels in a dwelling constructed with phosphogypsum plates

Decision analysis of emergency ventilation and evacuation strategies against suddenly released contaminant indoors by considering the uncertainty of source locations

Planning for protective action decision making: evacuate or shelter-in-place

Inverse heat transfer problems

A three-dimensional inverse heat conduction problem in estimating surface heat flux by conjugate gradient method

The method of fundamental solutions for the inverse heat source problem

State of the art report on mathematical methods for groundwater pollution source identification

Groundwater pollution source identification from limited monitoring well data: part 1dtheory and feasibility

Identifying sources of soil inorganic pollutants on a regional scale using a multivariate statistical approach: role of pollutant migration and soil physicochemical properties

Source identification of heavy metals in pastureland by multivariate analysis in NW Spain

Identification of pollution sources in urban areas using reverse simulation with reversed time marching method

Back-calculation of the strength and location of hazardous materials releases using the pattern search method

Utilizing state estimation to determine the source location for a contaminant

Identification of contaminant sources in enclosed environments by inverse CFD modeling

History recovery and source identification of multiple gaseous contaminants releasing with thermal effects in an indoor environment

History source identification of airborne pollutant dispersions in a slot ventilated building enclosure

Identification of contaminant sources in enclosed spaces by a single sensor

An inverse method based on CFD to quantify the temporal release rate of a continuously released pollutant source

Inverse modeling methods for indoor airborne pollutant tracking: literature review and fundamentals

Location identification for indoor instantaneous point contaminant source by probability-based inverse computational fluid dynamics modeling

Prompt tracking of indoor airborne contaminant source location with probability-based inverse multi-zone modeling

Experimental verification of tracking algorithm for dynamically-releasing single indoor contaminant

Principles and applications of probability-based inverse modeling method for finding indoor airborne contaminant sources

Identifying index (source) patient location of SARS transmission in a hospital ward

Inversely tracking indoor airborne particles to locate their release sources

Rapidly locating and characterizing pollutant releases in buildings

Systems approach to evaluating sensor characteristics for real-time monitoring of high-risk indoor contaminant releases

Influence of indoor transport and mixing time scales on the performance of sensor systems for characterizing contaminant releases

Towards improved characterization of high-risk releases using heterogeneous indoor sensor systems

Gaussian process emulator approach for rapid contaminant characterization with an integrated multizone-CFD model

Application of neural networks trained with multizone models for fast detection of contaminant source position in buildings

Real-time identification of indoor pollutant source positions based on neural network locator of contaminant sources and optimized sensor networks

Contaminant source identification within a building: toward design of immune buildings

A method to identify the point source of indoor gaseous contaminant based on limited on-site steady concentration measurements

Fast identification of multiple indoor constant contaminant sources by ideal sensors: a theoretical model and numerical validation

Rapid identification of single constant contaminant source by considering characteristics of real sensors

An optimization method of sensor layout to improve source identification accuracy in the indoor environment

An analysis of combined CFD and multizone IAQ model assembly issues

Multizone airflow modeling in buildings: history and theory

An analytical expression for transient distribution of passive contaminant under steady flow field

Perspectives on chemical hazard characterization and analysis process at DOE

The air toxics health effects database (ATHED)

Airpak 2.1 user's guide

Effectiveness of hybrid air conditioning system in a residential house

A zero-equation turbulence model for indoor airflow simulation

Numerical heat transfer and fluid flow

Reminiscences about the origins of linear programming

Engineering optimization: methods and applications

This study was supported by the National Natural Science