Method of Source Identification following an Accidental Release at an Unknown Location Using a Lagrangian Atmospheric Dispersion Model

: A computationally efficient source inversion algorithm was developed and applied with the Lagrangian atmospheric dispersion model DIPCOT. In the process of source location estimation by minimizing a correlation-based cost function, the algorithm uses only the values of the time-integrated concentrations at the monitoring stations instead of all of the individual measurements in the full concentration-time series, resulting in a significant reduction in the number of integrations of the backward transport equations. Following the source location estimation the release start time, duration and emission rate are assessed. The developed algorithm was verified for the conditions of the ETEX-I (European Tracer Experiment — 1st release). Using time-integrated measurements from all available stations, the distance between the estimated and true source location was 108 km. The estimated start time of the release was only about 1 h different from the true value, within the possible accuracy of estimate of this parameter. The estimated release duration was 21 h (the true value was 12 h). The estimated release rate was 4.28 g/s (the true value was 7.95 g/s). The estimated released mass almost perfectly fitted the true released mass (323.6 vs. 343.4 kg). It thus could be concluded that the developed algorithm is suitable for further integration in real-time decision support systems.


Introduction
The identification and characterization of sources of radioactive substances that have been detected, without or with limited prior information about their origin, has attracted the attention of researchers for a long time due to regular incidents of this type [1][2][3]. In such situations, the first questions that arise are about the origin of the release and the quantity of the released substance. The application of atmospheric dispersion models (ADMs) for the analysis of the potential sources of the detected contaminants is one of the important tools (often called 'inverse modelling' or 'source inversion') that help in solving the source identification problem. Even in the case of the Chernobyl accident in 1986, during the first three days, before the accident was confirmed by the USSR, the Chernobyl origin of the radionuclides detected in Sweden was at first identified by means of modelling and trajectory analysis [4]. The advanced and mathematically rigorous approach of solving inverse atmospheric dispersion problems was developed by Marchuk [5] based on adjoint equations of atmospheric transport, while Pudykiewicz et al. [6] applied this method for the global atmospheric transport problem following an hypothetical nuclear test. Following those pioneering works, the scientific community has made substantial progress towards solving the problem of locating an unknown source releasing a hazardous atmospheric pollutant through developing source inversion methods in particular for cases of regional-scale dispersion [3,[7][8][9][10][11][12]. These approaches are implemented by combining the available measurements, atmospheric transport models, and data assimilation algorithms. Presently, the source inversion methods have reached a relatively high level of sophistication, and during the recent incidents involving radioactive substances, such as detection of Ru-106 over the territory of Eurasia in the fall of 2017, they were successfully applied by researchers from different countries [8][9][10][11][12].
Besides the applications of source inversion methods to real accidents-in which the source parameters are never known exactly-their proper evaluation requires comparisons of their results with data of controlled atmospheric transport field experiments on the regional scale, in which source parameters are known. One of the most famous such experiments is the European Tracer Experiment-ETEX [13,14]. A considerable amount of work evaluated the source inversion methods against data of ETEX [11,[15][16][17][18]. However, most of them concentrated on estimating the emission inventories and time parameters of the release, while only Tomas et al. [11] reported on the estimated source location obtained after application of their model to the ETEX experiment. The estimated source location reported in [11] was at a distance of about 156 km from the true location (distance evaluated in this paper from the coordinates presented in Table 1 of [11]). Therefore, it is clear that additional studies concentrating on evaluation of estimated source location are needed.
Besides the strong public interest for information about the unknown source of detected radionuclides, radiation protection authorities are highly interested in the operational applicability of computational methods for this purpose. Some progress has been achieved in an automated solution of the source inversion problem in the case of a known source location-assessing in this case the quantities of the released substances [19,20]. In case of an unknown source location, an operational version of a simplified source inversion method was developed by Wotawa et al. [7] for the Comprehensive Nuclear-Test-Ban Treaty Organization (CTBTO). However, the automation of source inversion methods and their integration in real-time computerized decision support systems remains a challenging task because of the ill-conditioned nature of the related mathematical problem, i.e., the possibility of non-uniqueness of the solution or its sensitivity with respect to input parameters. An additional complexity in applying source inversion methods in the framework of real-time decision support systems arises from the fact that significant computational resources are required to run backward transport simulations. In standard source inversion approaches, the number of backward integrations of the atmospheric transport equation is equal to the number of measurements that are used for solving the source inversion problem. Currently, the application of inverse modelling methods requires a high user expertise and, in practice, such methods are usually applied by the model developers.
In a previous work [20], Kovalets et al. integrated an algorithm of estimation of timedependent emission rates for the case of a known source location in the Lagrangian ADM DIPCOT (DIsPersion over COmplex Terrain) [21,22], functioning within the EC-funded nuclear emergency response system RODOS (Real-time On-line Decision Support system for nuclear emergency management) [23]. The goal of the present work is to develop a computationally efficient source inversion algorithm for the case of an unknown source location applied with the DIPCOT ADM suitable for further integration in real-time decision support systems, such as RODOS. The present study is also motivated by the described above need for more studies of source inversion methods, with emphasis on evaluation of estimated source location against data from field experiments on atmospheric transport. This paper is organized as follows: in Section 2, the source inversion method is described; in Section 3, the results of the evaluation of the proposed method against the 1st release of the European Tracer Experiment-ETEX-I [13,14]-are presented; and the last section presents the conclusions drawn from the presented research study.

Assessment of Source Location
In this work, we use the correlation-based cost function for solving the problem of the identification of a source location as in [10,24]. Assuming a constant release rate, this approach allows for separating the solution of the problem of identification of the source coordinates, start time, and duration of the release from the problem of emission rate estimation. Moreover, Tomas et al. [11] successfully applied a time-integrated correlationbased cost function to the problem of estimation of source location. It could be shown that the usage of such a cost function is equivalent to the assumption of a stationary continuous source. Thus, the problem of source location estimation was further separated from the problem of estimation of the release start time and duration. The assumption of continuous source is also used in this work at the stage of source location identification. The correlation-based cost function and the corresponding minimization problem used for estimation of source location are defined as follows: where the arithmetic averaged variables are denoted by overbars; , , and are the source coordinates to be estimated by the minimization procedure; is the vector of the observed concentrations; and is the vector of the respective simulated concentrations. The observed concentrations may be time-integrated before being used in the minimization procedure, as it is further discussed below.
To evaluate the concentration values for different source coordinates, we integrate the equations of the Lagrangian puff/particle model DIPCOT in time [21]. Let us consider for example the forward equation of change of the -coordinate of puff p: where is the Reynolds-averaged x-component of the wind velocity vector, evaluated at the position of puff ( , , ), and represents the turbulent velocity fluctuation component. The latter is calculated by the stochastic Langevin Equation (3). In (3) the coefficients and ("drift" and "diffusion" coefficient, respectively) are functions of space, velocity, and time and depend on atmospheric stability conditions, while ( ) are increments of a Wiener process with zero mean and variance dt, uncorrelated in time. The relationships that are used in DIPCOT for calculating the coefficient are deduced from the solution of the Fokker-Planck equation containing the Eulerian velocity probability density function. The form of the latter depends on the atmospheric stability conditions. In the vertical direction for unstable conditions, turbulence is assumed to be inhomogeneous and skewed, while for stable conditions, it is assumed inhomogeneous Gaussian. In the horizontal direction, turbulence is assumed to be homogeneous isotropic Gaussian for all stability conditions. The coefficient is calculated by requiring consistency of the Lagrangian velocity structure function estimated by Equation (3) with the one obtained from the Kolmogorov's similarity theory. The expressions for calculating the standard deviations of the turbulent velocity fluctuations and the Lagrangian time scales that enter the relationships for the coefficients and depend on the stability conditions, the mixing layer height, and the height above ground. For specific details, the reader is referred to [21]. The model-puffs are created in accordance with the distribution of the sources volumetric density ( , , , ). In the case of a single source located at point ( , , ), acting with release intensity during the time interval [ , + ∆ ] ( is the start time of release, and ∆ is the release duration), is defined by the following relationship: where (. ) are delta-functions, and [ , +∆ ] ( ) is the 'indicator' function: Equations (2) and (3) are solved from the start simulation time = 0 to the end of simulation period = . In case of a source inversion problem, T represents the time of the last measurement.
To avoid a solution of Equations (2)-(5) of the forward model for every possible source position, we integrate the backward equations of puffs movement that for the xcoordinate of puff p could be written as follows: where 'backward' time ′ is defined as ′ = − and ranges from ′ = 0 to ′ = and the velocity fluctuation is defined by Equation (3) but with substituted by ′ . Similar equations are written for y-and z-coordinates of puffs. The backward equations are integrated separately for each measurement, and puffs are created in accordance to the distribution of the corresponding measurements' probe function , defined as follows. If a particular measurement was taken at spatial point ( , , ) and collected during time interval [ − ∆ , ] ( is the end time of measurement interval, and ∆ is the measurement duration) then the probe function is given by the relationship: where the indicator function [ −∆ , ] ( ) is defined analogously to Equation (5).
The backward Equation (6) in which the first (advection) term on the r.h.s. changes sign compared with the corresponding term of the forward transport Equation (2) while the second (diffusion) term remains unchanged is constructed by analogy to the Eulerian adjoint advection-diffusion transport equation [6]. Concentrations * calculated from puff coordinates ( ′ , ′ , ′ ) following the solution of the backward equations are thus analogous of those calculated from the solution of the adjoint equations. Even though strictly speaking Lagrangian and Eulerian models of atmospheric transport are not equivalent [25], it is assumed that forward and adjoint concentrations, calculated from puff coordinates after the integration of Equations (2),(5)-(7) approximately satisfy the same Lagrange duality relationship, as for the case of solutions of forward and adjoint equations of Eulerian atmospheric transport model [6]: Substituting relationships (4) and (7) in Equation (8), we obtain the following: where index k is the position of the respective measurement in measurement vector . Thus, provided that the backward Equations (6) and (7) was solved, the model value to be compared with the corresponding measurement is easily evaluated for arbitrary source location by using the adjoint concentration field. The r.h.s. expression of Equation (9) is further denoted as source-receptor function (SRF). As it was noted above, at the stage of source location identification, it may be appropriate to use the assumption of continuous release to separate the problem of source location identification from the problem of release start time and duration estimation. With the assumption of continuous release, the relationship (9) for SRF is changed to the following: Relationship (10) does not contain the unknown start time and duration of the release and ∆ , and it is further used in the process of source location identification. Since the finite-duration release is substituted by the continuous release in the source location estimation, it can further be assumed that it is reasonable to substitute the measurements time series at a station by a single time-integrated measurement at that station. In that case, vectors and are vectors of positions for all detector stations and not time. Indeed, even if the use of the continuous release assumption formally separates the start time of the release from the minimization problem (1), the plume arrival times at the measurement stations still affect the values of the cost function. Thus, overfitting is likely to happen, i.e., release at the wrong location and time could potentially better fit the plume arrival time at a station than the release at the true location and time. The use of time-integrated measurements excludes the influence of the correctness of the plume arrival time on the value of cost function (1). It should be noted that, in the case of stationary meteorology, the use of time-integrated observations is fully justified mathematically: it could be shown that, in stationary meteorological conditions, time-integrated concentrations following a finite duration release are proportional to concentrations following a continuous release with a proportionality constant independent of sensor's location [26].
Assuming that the measurement intervals for a given station do not overlap and their respective number is , the simulated value corresponding to time-integrated measurements is evaluated by a substitution of probe function (7) with the sum of probe functions defined for all measurement intervals: Then, the source-receptor function (10)  where index k is now the number of observation station. The adjoint concentration * was obtained from the solution of the backward Equation (6) with puffs created according to the probe function (11). Thus, the use of time-integrated measurements leads to a reduction in the computational time by a factor of 1/ −1 ( could be different for different stations; therefore, the overbar denotes averaging of reduction factors over different stations). This reduction factor is approximately equal to the average number of measurements at a single station. As it is shown below, the use of time-integrated measurements does not reduce the accuracy of the algorithm in calculating the source location.
For implementing the solution of the backward Equation (6) in the ADM DIPCOT, it is sufficient to change the signs of horizontal velocity components in the input meteorological data before processing them with the meteorological pre-processor (MPP) [22]. As the MPP uses a divergence minimizing procedure, the vertical velocity components is automatically inverted. The time order of the calculated meteorological fields by the MPP is inverted in accordance with the definition of backward time ′ entering Equation (6). The solution of the backward transport Equations (6) and (11) corresponding to different measurement sensors can be performed by a single ADM run with multiple sources located at the positions of different sensors and releasing different species. The calculated adjoint concentration of a particular species represents the solution of the backward transport equation for the corresponding sensor.

Assessment of Release Start/End Times and Inventory
After estimation of the release location, the start-time and duration ∆ of the release are assessed by minimizing the same cost function (1), using concentration measurements with a maximum possible time resolution instead of time-integrated measurements. Model concentrations are calculated in the course of the minimizing procedure with the already estimated source location but with different start times and durations of release. To do so, a forward run is performed with a continuous source at the estimated location and with a reference release rate (for example, = 1). The release period is split into = ⁄ subintervals represented by different groups of puffs ('species'), where is the release duration corresponding to a single group of puffs. Concentrations created by different groups of puffs at a given point represent the contribution of the corresponding release subinterval to the total concentration. The range of the possible values for is (0, , ), where , < is considered the time of the earliest detection of the contaminant at a measurement station. The range of possible values for ∆ is (0, − ), where the upper limit corresponds to a continuous release, lasting until the end of the simulation period. In the minimization algorithm, the ranges of the possible values for and ∆ are discretized with the time step . In practical applications, it may be appropriate to set this time step equal to the minimum time resolution of measurements, as it is below (Section 3.1). By selecting appropriate groups of puffs it is then possible to calculate concentrations created by releases with different start times and durations.
The pseudo-code of the minimization algorithm for finding and ∆ is the following: It should be noted that, at the stage of assessment of and ∆ , the value of the correlation-based cost function is still independent of the particular value of release rate . The latter is estimated after and ∆ have been defined by multiplication of the reference release rate by the ratio of the average observed over average simulated concentrations with the already evaluated and ∆ :

Method Evaluation against ETEX-I
Evaluation of the proposed method was performed for the conditions of the 1st release of the European Tracer Experiment-ETEX-I [13,14]. In the next paragraphs, we describe the model setup for the conditions of the experiment and the source inversion scenarios and present the results of the calculations.

DIPCOT Setup for the Conditions of ETEX-I
ETEX-I was performed on the territory of several European countries (Figure 1  In a previous work [27], DIPCOT was successfully evaluated for the conditions of the ETEX experiment. In this study, the DIPCOT dispersion simulations were performed using ERA5 reanalysis meteorological data downloaded from the Climate Data Store (https://cds.climate.copernicus.eu/#!/home, last accessed on 22 January 2021). The data have been retrieved from the following datasets: "ERA5 hourly data on pressure levels from 1979 to present" and "ERA5 hourly data on single levels from 1979 to present". The spatial resolution of both datasets is 0.25° × 0.25°. The original meteorological data have been processed by the MPP to make the conversion towards the computational grid of the ADM and to calculate missing variables that are needed by the ADM. The MPP performs a divergence minimizing procedure to ensure mass conservation, during which the three wind velocity components are adjusted by taking into account the underlying topography. The Lagrangian-puff mode of operation of DIPCOT has been employed. A non-depositing passive tracer has been assumed. The time interval between puff releases was set to 10 s. The model was set to calculate concentrations at 10 min intervals, and from these results, time-integrated and time averaged concentrations were calculated. The solutions of the backward equations were performed with the same settings of DIPCOT as in the forward runs. For the backward runs, the meteorological fields have been processed by the MPP changing the signs of the horizontal velocity components and by reversing their time sequence as described in the previous section. In the backward DIPCOT runs, puffs were released from the locations of the measurement stations as described in Section 3.2. The computational domain of DIPCOT is Cartesian. The geographical coordinates of its lower-left corner for the ETEX-I were set at (3.05° W, 38.52° N) and that of its upper-right corner were set at (30.45° E, 60.63° N), resulting in a dimension of 2600 km in the westeast direction and 2400 km in the south-north direction in a Lambert Conformal Projection. This computational domain covered an extended area around the ETEX-I domain.
The adjoint concentrations were saved on a computational grid with a horizontal spatial resolution of 10 km at 1.5m above ground. Thus, in inversion runs, it was considered known that the release happened near the ground surface. Before using the adjoint concentrations in the source inversion algorithm, the validity of the Lagrange duality relationship (8) was checked by comparing concentrations calculated at locations of stations by adjoint concentrations and Equation (9) with the concentrations calculated in a forward run, assuming release from the true source location. The correlation coefficient between the concentrations calculated in the two different ways was 0.89. Thus, it was concluded that adjoint concentrations could be then used in the source inversion algorithm.
In the solution of the forward transport equation for estimation of the release start time and duration according to algorithm (13), the time interval between groups of puffs representing different times of release was set to = 3 h, equal to the time resolution of concentration measurements in ETEX. Concentrations in the forward runs were calculated and saved at the exact locations of measurement stations.

Results of Simulations
A basic aim of this study was to evaluate the feasibility to use a single time-integrated concentration measurement at each station instead of all of the measurements available in the concentration time series at the respective station. For this purpose, we performed two sets of source inversion (SI) calculations. In the first set (SI-1) we used all of the 3 h measurements from a reduced set of the stations that are available for ETEX-I. Thus, we selected 46 stations (appearing as black rhombuses in Figure 1) out of the full set of ETEX-I stations ( Figure 1) and performed backward dispersion simulations to calculate adjoint concentrations with probe functions defined by Equation (7). Two inverse computations were performed in this set to determine the source location: in the first run (SI-1a), the original 3 haveraged observed concentrations were used for source inversion (in total, 419 measurement values). In the second run (SI-1b), the sum of the observed concentrations at a given station multiplied by the averaging period of a single measurement = 3 h was used in source inversion: For the stations where concentrations were continuously monitored, i.e., when there were no time gaps in the measurements, the above value corresponds to the time-integrated concentration. The total number of time-integrated measurements used in SI-1b was thus equal to the number of the selected stations (46)-almost by one order of magnitude less than the number of measurements in SI-1a. In backward runs for the scenario SI-1b, the puffs were released according to the distribution of the probe function (11) corresponding to the sum (15). The estimated source locations in both cases SI-1a, b are shown in Table 1. As it is seen from the table, the difference between estimated source locations in the two runs is small: distances between estimated and true source locations in runs SI-1a,b are 99.3 and 97.0 km, respectively. Moreover, the distance between the two estimated source locations is 23 km, which is by a factor of more than 4 less than the abovementioned distances between estimated and true source locations. It should be also noted that, as presented in Table 1, the maximum correlation coefficients achieved in runs SI-1 were quite different: 0.41 in case of SI-1a and 0.76 in case of SI-1b. Thus, in the case of using time-integrated measurements, higher values of correlation coefficients are obtained, which is expected considering that slight discrepancies in plume arrival times lead to deterioration of correlation coefficient in the case of using time-dependent measurements. Table 1. True and estimated source locations in SI-1 in cases of using all observations from the selected stations (SI-1a) and using time-integrated observations from the selected stations (SI-1b).  Figure 2 shows spatial distributions of the correlation coefficient, normalized by its maximum value (Table 1), calculated in runs SI-1a, b. It should be noted that, despite the similar shapes of the obtained distributions, in the case of using time-integrated measurements (SI-1b), the area of possible source locations (defined for instance as the area covered by the 0.5 isoline of normalized correlation) is greater compared with the case of using time-dependent measurements (SI-1a). Nevertheless, because of the closeness of the estimated source locations in runs SI-1a and b noted above, it could be concluded that it is reasonable to use time-integrated measurements for the estimation of a source location in operational applications of a source identification method. A reduced set of observations for runs SI-1 was used for two main reasons: to reduce the required computational resources and, because in an emergency situation, the number of measurements is usually limited, so testing with a realistic number of measurements is necessary. The stations for the runs SI-1 were chosen based on two criteria: (1) the area covered by the selected stations was comparable with the whole experimental area and (2) time breaks in measurements at the respective stations were absent. However, because the accuracy of the obtained estimation of source location was almost the same in runs SI-1a,b as discussed above, it was concluded that calculation with the raw 3 h averaged concentrations for the whole ETEX-I dataset is not necessary.
Therefore, in the second run (SI-2), only time-integrated measurements from all ETEX-I stations that had given valid measurements (168 values) were used for estimation of all source parameters-location, release start time and duration, and release rate according to the method described in Section 2. Table 2 presents a summary of the obtained results from this run. The estimated source location is at distance of about 108 km from the true source, which is slightly worse compared with the set of runs SI-1. Although the number of measurements used in the SI-2 run was about 3.7 times greater than in the SI-2b run, the distance between the estimated and the true source location changed only by about 11%. Therefore, the source inversion method is not sensitive to the number of measurements used, at least when this number exceeds several dozens of values. At the same time, the results of source inversion in the SI-2 run are of comparable and even better accuracy than the results of source inversion reported for the same experiment by Tomas et al. [11]. It is also interesting that the distance between the estimated and the true source location is comparable with the resolution of the measurement network, estimated as the square root of the area covered by the measurements divided by the number of measurement stations and equal to ≈104 km.  Figure 3 shows the spatial distribution of the normalized correlation coefficient obtained in the process of source location estimation based on the run SI-2. Qualitatively this figure is similar to the distribution in runs SI-1. To proceed with the uncertainty analysis in run SI-2 the confidence limits of the estimated source location were calculated by application of the bootstrap method [28]. To do that, the MATLAB "bootci" function was used with the measurement vectors as input together with the pre-calculated adjoint concentrations and combined with the subroutine implementing the minimization process of the cost function (1). The number of iterations (i.e., runs of the minimization process with randomly excluded records from the measurements vector) used for obtaining the 95% confidence intervals of the estimated source coordinates was set to 10,000. The estimated confidence limits for source coordinates were −5.61 ≤ lon ≤ −0.24 deg. and 47.61 ≤ lat ≤ 51.87, as shown in Figure 3 by the square with the indicated coordinates of corners. The domain within the confidence limits is quite large; however, its area is by a factor of about 30 less than the area of the computational domain. Thus, the inverse method allows for reducing the search area of the unknown source by more than one order of magnitude.
The estimated start time of the release in run SI-2 (23 October, 15:00) is very good, keeping in mind that the difference with the true start time is only about 1 h, i.e., less than the averaging time of measurements-3 h, which is the best possible accuracy of estimate of the start time. The release duration obtained in run SI-2 (21 h) overestimates the duration of the true release by a factor of almost 2 (12 h). This increase in the release duration is compensated by a decrease in the estimated release rate (4.28 g/s as compared to 7.95 g/s) so that the estimated released mass fits almost perfectly the true released mass (323.6 vs. 343.4 kg).

Conclusions
In this work we developed a computationally efficient source inversion algorithm for the case of an unknown source location and applied it with the Lagrangian atmospheric dispersion model DIPCOT. The algorithm is based on a previously developed source inversion approach that uses a correlation-based cost function, but with the additional assumption of continuous release at the stage of source location estimation. The latter assumption allows for separating the processes of source location estimation from the assessment of release start time, duration, and released substance quantity. The backward equations of atmospheric transport are solved to calculate adjoint concentrations, which are subsequently used for the rapid calculation of the concentration values to be compared with measurements, for different locations of the unknown source. The computational efficiency of the algorithm of the source location estimation is further improved by using only the values of the time-integrated concentrations at the respective stations instead of the full concentrations time series. This approach, which was proposed for the first time in this work, allows for reducing the number of integrations of the backward transport equations by a factor approximately equal to the average number of measurements at a single station. Following source location estimation, the release start time and duration are evaluated by minimizing the same correlation-based cost function but using the assumptions of finite duration release and full concentration time series measured by stations. The emission rate is then readily obtained by normalizing the unit (reference) release rate with the ratio of the average observed over simulated concentrations.
The developed algorithm was evaluated for the conditions of the famous ETEX-I experiment on regional-scale atmospheric transport. It was shown for the first time that using the values of the time-integrated concentrations at the respective stations instead of the full concentrations time series did not increase the error in the estimated source location. Using time-integrated measurements from all of the ETEX stations, the distance between the estimated source location and true location was 108 km, which is comparable with the results of other published algorithms. The estimated start time of the release was only about 1 h different from the true start time-a difference that is less than the averaging time of measurements (3 h), which is the best possible accuracy of estimate of the start time. The release duration evaluated by the calculations (21 h) was overestimated by a factor of almost two compared with the true value (12 h). This increase in release duration was compensated by a decrease in the estimated release rate (4.28 g/s compared with 7.95 g/s) so that the estimated released mass almost perfectly fitted the true released mass (323.6 vs. 343.4 kg).
It can thus be concluded that the developed algorithm is suitable for solving source inversion problems of atmospheric transport of the temporal and spatial scales similar to the conditions of ETEX experiment-release duration up to 24 h, spatial scale of up to about 1000 km. To aid decision makers in dealing with such problems, the method could be integrated in real-time decision support systems for hazardous substances atmospheric releases. The application of the method to problems of larger spatial and temporal scales needs further testing, especially with respect to the possibility of using time-integrated measurements.
Author Contributions: Conceptualization, methodology, software, and validation, S.A. and I.V.K.; writing-original draft preparation, I.V.K.; writing-review and editing, S.A. All authors have read and agreed to the published version of the manuscript. Data Availability Statement: The study used ERA5 reanalysis meteorological data downloaded from the Climate Data Store (https://cds.climate.copernicus.eu/#!/home, last accessed on 22 January