Stand-off detection of airborne chemical compounds has proven to be a useful method that is gaining popularity following technical progress. There are obvious advantages compared to in situ measurements when it comes to the security aspect and the ability to measure at locations otherwise hard to reach. However, an inherent limitation in many of the stand-off detection techniques lies in the fact that the measured signal from a chemical depends nonlinearly on the distance to the detector. Furthermore, the measured signal describes the summation of the responses from all chemicals spatially distributed in the line of sight of the instrument. In other words, the three dimensional extension of the chemical plume is converted into a two-dimensional image. Not only is important geometric information per se lost in this process, but the measured signal strength itself depends on the unknown plume distribution which implies that the interpretation of the observation data suffers from significant uncertainty. In this paper we investigate different and novel approaches to reconstruct the original three-dimensional distribution and concentration of the plume by implementation of atmospheric dispersion models and numerical retrieval methods. In particular our method does not require a priori assumptions on the three-dimensional distribution of the plume. We also strongly advocate the use of proper constraints to avoid unphysical solutions being derived (or post-process ‘adjustments’ to correct unphysical solutions). By applying such a reconstruction method, both improved and additional information is obtained from the original observation data, providing important intelligence to the analysts and decision makers.

For a week in April 2010 the ash emitted from the volcanic eruption at Eyjafjallajökull, Iceland, caused major disruptions in the European air traffic. There were in total 108.000 canceled flights affecting 10.5 million passengers and causing a $1.7 billion loss of revenue for the airline industry alone [

As the examples above show, deliberate or undeliberate, injection of gases or particles into the atmosphere may cause downwind hazards in terms of pollution, damage, discomfort, irrevocable injuries and in worst case death. Fortunately, many of these risks may be mitigated or reduced by taking precautionary safety measures such as donning gas masks, evacuating people or canceling flights as appropriate to the situation. Now, wearing gas masks or canceling flights when no danger is imminent carries a high cost too, hence there is a need for the ability to predict where these hazardous substances are transported following a release, and to be able to assess whether they still constitute a hazard. Predicting downwind concentrations and arrival times are done using atmospheric dispersion models, and translating these downwind concentrations to potential hazard is done using toxicological models. The dispersion model will require three categories of inputs: weather, topography, and a source. This implies that even though it is actually the plume that is of interest, we need to find these three inputs. When this is accomplished we are able to predict the plume. Now, the weather and the topography are most likely known with an acceptable certainty. The source, on the other hand, is often considerably more challenging to estimate. Properties such as the substance, location, time-dependent release rate, and initial plume dynamics need to be quantified which requires information provided by measurements or a priori information.

Measuring the flow and injection height of a volcano, or accessing conflict zones like Al-Mishraq and Khan Shaykhun are either difficult, dangerous or both. In lieu of measurements at the source and downwind point sensors, data in the form of stand-off detections may be used to inform us on the characteristics of the plume and, indirectly, the source. Stand-off detections may come from satellites (e.g. MetOp-A/B, Suomi, Aura, the French-German Space Lidar Mission MERLIN [_{3}, SO_{2}, NO_{2}, BrO, HCHO), from airplanes or helicopters carrying similar equipment, e.g. the DLR Falcon 20 or the Aerospace Corporation [

A typical example when stand-off measurements provide column densities, i.e. area concentrations, of the toxic plume. The area concentrations may be considered an image of the plume produced by the sensor of the stand-off detector. We are interested in reconstructing the plume that gave rise to the sensor response. Reconstruction of the plume is an ill-defined operation which is inherent to the type of problem that we are studying. The objective when reconstructing the plume is to find the plume that best explains the combined sensor response.

Many stand-off detectors provide us with an area concentration field, often in low spatial and temporal resolution. We want to increase the resolution and uncover the missing spatial dimension, in other words we want to use the low-resolution area concentration to reconstruct the high-resolution concentration fields of the toxic cloud. To do so, dispersion simulation is introduced and the source term is therefore required. Fusing sensor data to recover the best estimate of the source term is a maturing field of atmospheric dispersion modeling and is an inverse modeling problem. Most analyses found in the literature address the problem of estimating the source term based on in situ point measurements using either a probabilistic approach (Bayesian) or optimization techniques. There are results based on stand-off detections, in particular satellite observations of volcanic emissions, e.g. [

To summarize, we present a data fusion method that improves on inadequate information regarding concentration fields. The method is fed with empirical area concentration data in low resolution with no information regarding the time-dependent release rate or the initial plume dynamics (e.g. buoyancy effects). As an output the method provides the lacking source information as well as a high resolution concentration field. The key problem is to retrieve the source term which is the main focus in this work. Given the source term provided by the data fusion method, detailed analysis of the plume in the past and prediction of the plume in the future in an ongoing event may be obtained by a straightforward dispersion simulation.

The main idea described in this work is to fuse simulated concentration data with stand-off observations to acquire more information about the dispersing cloud than is provided by the observation data alone. The simulated data is provided as an output of an atmospheric dispersion model. The atmospheric dispersion model is based on the predominant physical phenomena that govern the dispersion of a chemical plume with the aim of emulating the reality. As mentioned above, the source term includes the following properties: substance, location, time-dependent release rate, and initial plume dynamics. The first two are known in many cases and for the remainder of this study it will be assumed that this is the case. Some key information is typically missing when studying a real case, in particular the time-dependent release rate, i.e. released mass per time unit, of the source and the coupled initial plume dynamics, i.e. the vertical distribution of the plume close to the source. However, since these properties scale in a well described manner with respect to the final concentration fields, they lend themselves to be tuned post-simulation to match the simulated plume with the observations of the real plume. This tuning constitutes the data fusion and is the core of the method presented here. While saving the details just for now, we establish that the main ingredients required to perform data fusion are:

Design a source term for simulations

Perform dispersion simulations

Formulate commensurable data sets

Solve for source term parameters

The reason why the third issue is raised is technical but important: to be able to conduct the data fusion the simulated data set and the measured data set have to be commensurable. To achieve this we have to bear in mind that the observation data can often be burdened with strong and unavoidable limitations. These limitations may typically lie in the temporal and spatial resolution of the measured data set. In contrast, the simulated data can be extracted in virtually any desired resolution in both the temporal and spatial domains which implies that it should be chosen to match the observation data. It is therefore the observation data that constitutes the constraints of the system and determines the construction of the commensurable data sets.

Let us now turn to the details, i.e. the purpose, options, implications, and possible problems related to each of the steps outlined above.

To be able to simulate the dispersion process at hand we need to design a suitable source term. In a typical synthetic dispersion simulation all properties of the source are defined as input to the simulation model. This is not possible in the type of cases studied here since the release rate and the vertical distribution of the released substance are unknown a priori. The problem with lacking source information is addressed by designing a source term that includes a set of parameters whose ‘true’ values will be obtained as the result of the data fusion process. This means that dispersion simulation will initially be conducted using guessed, or normalized, values for some of the parameters that may be tuned (altered) post-simulation. The tuning of these parameters is handled in the source parameter solving algorithm which is discussed later. The properties that need to be defined in the source design are: the known substance, the known source location, the known time interval of the release, a ‘temporal distribution method’, and a ‘vertical distribution method’. In the context of this work, we assume that the first property, the substance, is known. The location of the source is indeed often well known. One exception may be forest wildfires were the source location and extension actually changes with time which is a special case we will not discuss in detail here. The time interval of an entire release event refers to the time points when the release started and ended. In an analysis of an ongoing event, the time interval ends at the present time or at a time that relates to the last available observation. As with the location, this property is also often known with fair certainty. Even so, this is actually not pivotal information when the data fusion method is applied since the release rate is obtained as a result for all time points, hence the time interval can ideally be inferred with good accuracy from the results. However, all knowledge that guides the source term parameter solver as constraints will speed up the algorithm, improve the results and decrease the uncertainties. It is therefore recommended to provide the method with this information if available. We are then left with two important and unknown properties: the temporal distribution method and the vertical distribution method.

The term temporal distribution method refers to how the source is designed to manage the dynamics of the release rate, i.e. how the release rate changes over time. In reality the release rate is a continuous function with an autocorrelation over time. The autocorrelation may vary from case to case and over time but the source term becomes unrealistic as the autocorrelation approaches zero since that implies that the release rate may change freely from one second to the next. An estimate of a lower limit of the autocorrelation may be acquired from previous experience of similar events. There is no need to set an upper limit of the autocorrelation. There are two principally different methods that may be applied to the source design. The first method is to let the release rate vary continuously using a sum of chosen subfunctions, _{i}_{i}

Left, a small set of five subfunctions is shown as well as the release rate they add up to as equation 1 shows. Normalized Gaussian subfunctions have been chosen in this case. Right, illustration of the result from an arbitrary analytic plume rise function that provides the release rate density as a function of the current release rate, meteorology and height. The altitude of the initial plume is in general linked to the release rate, i.e. a high release rate implies a high altitude. The vertical integrals of the curves equal the release rates, i.e. the amount of released substance per time unit.

The second method applies discretization of the time into

_{0} to equal the start time of the release). Equation (2) is merely a trivial stepwise constant release function. Even though this method obviously misrepresents the continuousness of nature, it may be argued that the present uncertainties renders this approach ‘good enough’. In both methods it is strongly recommended that the autocorrelation over time is taken into consideration. This may be implicitly incorporated in the case of the continuous approach by the choice of subfunctions since they may eliminate all abrupt changes in release rate. In the method with a discretized timeline, the autocorrelation may be applied by the use of a regularization term which is discussed later.

As in the case with the temporal distribution of the release rate, the vertical distribution of the released substance is a result from the data fusion process but the method of representing the distribution must be defined in advance. It is well known that hot sources, such as large fires and volcanic eruptions, cause significant buoyancy which initiates a plume rise of the released chemical [

The analytic function usually incorporates a set of meteorological conditions, {Met}, that may change during the event but are always considered known, or deductible from the given meteorological input data. The standard analytic function for plume rise is the Briggs function [

The second method includes a definition of a limited number, _{k}

The set _{k}

In case of a continuous temporal method, _{k}_{k}

Here _{i}_{i,k}_{i,k}

The four different combinations of continuous and discretized time and heights. Yellow and blue colors indicate high and low release rate densities, respectively. Discretization in the vertical distribution is here made in five layers while the temporal discretization holds 30 different time intervals. These panels illustrate how the source term may be designed differently. The total masses in the fields are equal in all cases, i.e. the total release remains unchanged and independent on the choice of source term design.

The main principal difference between these two methods is the lack of direct dependence between total release rate and plume height in the discretized method, rather it is assumed that the dependence will emerge automatically from the source parameter solving algorithm. In other words, it is assumed that the actual plume rise will, due to the vertical variations in the meteorological conditions, influence the transport of the plume which will be observed. The source parameter solving algorithm will therefore force the vertical distribution in the source term to capture the correct initial plume rise for the simulated plume to mimic the observed data. To avoid physically unsound solutions, it is recommended that constraints and a regularization term are applied in the source parameter solving algorithm to smooth out too large variations in the relative release rate densities between the vertical layers of the source.

It might be worth to notice that the choices of temporal and vertical distribution methods show great resemblance. Either a continuous function or a discretization of the domains is applied. Even though all choices are valid, care must be taken in the choice of subfunctions, discretization, regularizations and constraints.

The four different combinations of continuous and discrete temporal and vertical distribution methods are illustrated in

With the source designed, regardless of method, a set of parameters (_{i}_{k}_{i,k}

Atmospheric dispersion simulation is the process of simulating the fate of released substances as they propagate in the atmosphere by the conveying wind field and disperse by turbulence. There are many different dispersion models that each handles meteorological information but also phenomena such as chemical degradation and dry and wet deposition differently. What they all have in common though is that they use the designed source term to deliver concentration fields at desired time points and at desired spatial resolutions. To simplify the nomenclature ahead we will assume that the dispersion model uses model particles (Lagrangian model). These model particles can represent both physical particles, aerosols and gases.

Here it is assumed that the dispersion process is linear which means that the evolution of the released chemicals is independent of the amount of released mass. Each model particle moves and evolves independently from the other model particles. Note that the linear assumption is not to be confused with the fact that the initial vertical distribution of the released mass do depend on the release rate. This relation lies within the source design and not the atmospheric dispersion model itself. The dependence of the vertical distribution upon the released mass will influence how much chemicals that are released at different heights, but not their movements in the air. This is a crucial property in the data fusion process since it allows for a direct recalculation of the concentration fields as a result of tuning the release rate and vertical distribution. That is, the position of each model particle at any given point in time will not change due to change in the parameters of the source.

Any set of source parameters give rise to a certain spatiotemporal release rate density as depicted in

Measurement data of an area concentration field is normally provided as data points on a two-dimensional grid or as scattered data points. In the latter case, the data points may be regularly or irregularly distributed over the geometric region of interest, i.e. the data set provides pointwise area concentration values distributed arbitrarily over the domain. Scatter data can be interpolated to a two-dimensional grid using methods such as bi-linear or bi-cubic interpolation [^{20} molecules per square meter. Commensurable area concentration fields can readily be rendered from simulated data by integrating the field column-wise as has been done by the detector. This implies that whenever a simulation provides a concentration field, the corresponding area concentration field is always available as well.

For the data sets to be truly commensurable they must be spatially represented in an equal manner as described above. Additionally, the sets must also contain data corresponding to area concentrations valid at the same time points. That is, the simulated area concentration fields must be collected at the same point in time as the measurements were conducted. This is often a trivial matter, but in cases where the time scale of the data measurements is on the same, or higher, order as the time scale of the dynamics of the plume, this issue must be addressed. When this occurs, it is possible that the plume has changed position and distribution from the beginning to the end of the measurement of one single area concentration field. Since area concentration measurements typically are conducted by sweeping measurements of parts of the area, each such part is measured at different time points. The solution is therefore to collect the simulated area concentration field at diverse time points that correspond to the time points of the measurements, i.e. the different parts of the simulated area concentration field are not to be collected at the same time point, rather their time points should mimic the sweeping mechanism of the detector. However, when the plume dynamics are significantly slower than the measurement procedure, this detailed data collection is unnecessary.

If the observed area concentration data is provided as scatter data points, a grid is constructed onto which the data is inserted. A two-dimensional (Cartesian) grid leaves few degrees of freedom. Any spatial translation of the grid will in general have little influence of the result and is ignored in this discussion. A more important setting is the grid resolution, i.e. the size of the grid cells. There are a couple of issues to consider while choosing the grid resolution. The plume must be resolved on the grid which sets a lower limit on the grid cell resolution. Though this limitation might be quite obvious, the other limit of the resolution scale might be less trivial to spot. In the source parameter solving algorithm proposed in this paper, the agreement between observed and simulated data is evaluated on the basis of comparison between individual grid cell area concentrations with no consideration to adjacent grid cells. This means that uncertainties in the prevailing winds may easily translate the simulated field one or more grid cells if the grid resolution is high. If this is the case, the source parameter solving algorithm will struggle to guide the solution in the right direction since there is a spatial mismatch in the field comparison. A high grid resolution may also require a large number of model particles depending on which method that is used to compile the simulated area concentration field. Common methods include straightforward box-counting [

In this final step of the data fusion process the observed and simulated area concentration data sets are compared with the aim of identifying the optimal values of the source parameters. That is, the values of the source parameters that provide simulation fields that best reassemble the observed fields are found. This set of parameters will then define the source term, which, as input to the dispersion model, will provide the full spatio-temporal development of the complete release event. Note that the final concentration field in high spatiotemporal resolution is directly available by asserting particle weight according to the optimal source parameters as discussed in Section 2.3. The reconstruction of the observed plume will then be complete.

The mechanism to find the optimal solution is principally straightforward. Any set of parameters will provide area concentration fields for all time points where observed data is available. Observed and simulated data sets are conformed on a suitable grid, rendering them commensurable, and the degree of disagreement is quantified using a scalar statistical measure

where _{obs,m}_{sim,m}

As a final note, we will discuss the timesaving trick referred to in the section Source design. When the source term parameters are changed it implies new weights of the model particles which translates into new area concentration fields. It may be time consuming to recalculate the area concentration fields over and over again as the search for the global minimum is under way. However, there is a shortcut in case of a source design utilizing discretization in both the temporal and vertical domains. In this specific setup there are many model particles with identical weights. One might capitalize on this fact by calculating one separate area concentration field for each area segment in the lower right panel of

As an alternative to using an optimization algorithm we may view the source estimation problem as a probabilistic one and cast the problem as such, that is, instead of searching for an optimum set of source parameters we look to assign a probability to each combination of parameters. Guided by these probabilities the source parameters may be chosen. Making this choice could be done in several ways, e.g. by picking the maximum likelihood estimate or the expected value. By far the most common method of assigning probabilities to parameters in source estimation problems is using Bayes method [

Bayes theorem thus relates the posterior distribution to the conditional probability, the ‘likelihood distribution’,

Non-detailed references to regularization terms and constraints have been given at several occasions above. A regularization term [

Please note that this case study is described in full detail in the article ‘The 2016 Al-Mishraq sulfur plant fire: Source and health risk area estimation’ [

October 20 2016, retreating Daesh troops ignite the extensive sulfur stock piles of the Al-Mishraq plant in northern Iraq. The fire spreads rapidly and a release of sulfur dioxide at the size of a minor volcanic eruption is soon a reality. The Swedish Defence Research Agency launches an investigation into the matter by employing satellite data from optical spectrometers measuring backscattered solar irradiance in the ultraviolet part of the spectrum to monitor, among other substances, airborne SO_{2}. As the satellites orbit the planet, they provide rapid measurements of SO_{2} in Dobson units, see

Measurement data from one instrument (MetOp-A) for October 24. Left, each band on the world map represents one polar orbit and the entire earth is covered in 24 h by two satellites working in pair. The plume over the Middle East is visible in the global map but even more clearly in the zoomed view over the Middle East to the right. Some data points are removed due to cloud coverage. Note that these values need to be scaled according to the actual height of the plume.

It is clear from photographs by Al-Sudani/Reuters [

The resulting source term presented in two different ways. Left, the release rate for the three different vertical layers. Right, this presentation uses the time intervals and the vertical heights of the virtual sources and indicates the release rate per volume. The two panels show the same information where the right panel includes information on the virtual source volumes. The right panel also functions as a depiction of the source design utilized in this case study. 25 time slots and three vertical layers result in 75 different virtual sources and an equal number of source parameters to be found. This figure is a reproduction from [

The optimization algorithm is built using the scalar statistical measure that is defined as the sum of the mean square error between the observed and simulated fields for all 12 time points, compare Eq. (7). This measure now depends on the 75 source parameters, i.e. the release rate densities for the virtual sources. Since the source term is defined using full discretization, individual plumes for the entire scenario can be produced for each virtual source as discussed above. Alteration of the source parameters are directly implemented as linear scaling of the plumes for each virtual source which allows for a relatively effective source parameters solving algorithm. A complication arises as the observed fields depends on information of the plume height, which is unknown a priori. This plume height is therefore estimated from the simulated plume given the current source parameters. This implies that both the simulated and the observed area concentrations become dependent on the source term parameters. Constraints on the daily difference in total airborne mass between the observed and simulated plume are applied together with a requirement of non-negative release rates. Moreover, a regularization term is implemented to moderate any unrealistic fluctuations of the release rate that may originate as artifacts from the low frequency of observation data. The regularization term and other case specific details are available in Appendix A coupled to the paper on this case study [

A variety of stand-off measurement instruments and systems provide area concentration data that represents column-wise integration of the concentration field along the line of sight. This is a problematic limitation for two reasons. First, the data is lacking information regarding the spatial distribution of the substance measured in the direction of the observation. Secondly, the data signal is often dependent on this unknown spatial distribution and this property gives rise to a significant and troublesome uncertainty for the measured area concentration. A general description of a method to address this issue is presented in this paper. The method allows for different choices in the approach both with regard to the source description and the source parameter solving algorithm. There is, however, a distinct path to follow regardless of the specific choices made by the data analyst. Implementation of such a method will provide extensive additional information of the physical reality hidden in the measurement data.

Four combinations for the source term representation have been discussed. It can be argued that the completely continuous or completely discrete setups are to be recommended. The case with a completely continuous source term, i.e. continuous in both the temporal and the spatial dimensions, is intrinsically the closest representation to a physical event. As a drawback this approach requires a well-working analytic plume rise model and a carefully thought-out set of release rate subfunctions. As a contrast, the completely discretized source term is a somewhat crude representation of the reality. However, compensation is to be found in the fact that this type of source term is easy to set up, fast to execute and can arguably be a sufficiently sophisticated design given the influences of uncertainties and errors from other factors in the analysis. A general note is that it might be computational heavy to locate a global minimum in a parameter space of high dimension, which is worth keeping in mind when considering the design of the source term.

It is important to consider, and probably include, constraints and regularization terms to avoid unphysical solutions. These solutions may arise due to sparse measurement data and biases in the simulated distribution due to errors in the meteorological predictions and/or the dispersion model. A generic constraint is to disqualify all negative release rates. It is worth noticing that this limitation should be included as a constraint in the source parameter solving algorithm. It might be tempting to take a shortcut and instead set all negative values to zero as post-processing. However, there is no guarantee that the actual global minimum will be identified and preserved using this shortcut.

In this paper we have put the source design method on a mathematical footing discussing limitations and possibilities of discretizing vs. keeping aspects of the source continuous, which stand in contrast to much of the literature where the discretized setup tends to be the starting point. In contrast to the literature covering satellite based stand-off detections to estimate source terms, the method presented here requires no prior assumptions on thickness/width of the plume – a property that highly affects how to interpret the measurements, instead this method iteratively determines these properties as part of the process of matching observed and simulated data. Finally we impose our constraints on the problem up-front instead of making post-processing ‘adjustments’. A key limitation of the method is that it subsumes the location of the source or the sources. Another limitation is that the method does not include a description of how to choose the best regularization term. Also, we have not performed a thorough study of under which conditions the method converges to the solution, but the Al-Mishraq case study suggests the method is fairly robust.

We note that atmospheric dispersion simulations target several shortcomings that are present using stand-off detection data. The highly resolved fields provided by the simulations help the analyst to interpret measured column data, increase the temporal resolution, adds an additional dimension to the concentration field and may be used in an ongoing event to forecast the future development of the plume and thereby assist authorities to establish health risk areas for the coming days.