A Digital Twin Framework for Environmental Sensing with sUAS

This paper proposes a digital twin (DT) framework for point source applications in environmental sensing (ES). The DT concept has become quite popular among process and manufacturing industries for improving performance and estimating remaining useful life (RUL). However, environmental behavior, such as in gas dispersion, is ever changing and hard to model in real-time. The DT framework is applied to the point source environmental monitoring problem, through the use of hybrid modeling and optimization techniques. A controlled release case study is overviewed to illustrate our proposed DT framework and several spatial interpolation techniques are explored for source estimation. Future research efforts and directions are discussed.


Introduction
Environmental sensing and monitoring includes a variety of applications such as measuring local air quality, soil remediation, oil spills, ground water quality-tracing, etc. These applications can include both fixed and mobile sensors. For example, sensing ecosystem flux can be done with fixed sensors in the soil, above ground sensors, with mobile sensors (e.g. ground vehicles or sUAS [36]) or remote sensing techniques (such as with satellite imagery or manned aircraft [47]). Typically fixed sensor approaches work well to provide precise ecosystem flux measurements (e.g. chambers, autochambers, and eddy covariance towers [13]). In terms of site level emissions, fixed open path lasers have been used to evaluate and compare with accepted Environmental Protection Agency (EPA) methods (e.g. dairy farms) [10]. Within the past few years, several works in the area of environmental sensing with sUAS have been done (see the comprehensive review by [18]). Typically the sensors consist of low cost chemical sensors (e.g. amperometric, metal-oxide, non-dispersive infrared, or photo-ionization detectors) or higher cost tunable diode laser absorption spectroscopy (traditional open/closed cell -TDLAS or backscatter based -sTDLAS). Other techniques like optical gas imaging, which uses thermal camera based sensing have been used with sUAS. For instance, the sTDLAS (e.g. Pergam Methane mini-G (SA3C50A)) has been used for determining spatial methane in permafrost [38]. In the oil and gas industry, the TDLAS (e.g. open path laser spectrometer -OPLS [22], or path integrated) has been used for detection [50,51] and for quantification [27], such as in quantifying the Alberta Methane Field Challenge (AMFC) site emissions [52]. The use of fixedwing sUAS have been used as well in both detection [15] and quantification [8].
Due to the large variety of approaches and applications of environmental sensing, the focus of this paper will be on ground-to-air emissions (as an extended versions of the conference paper [30]) with an application to gas plumes (e.g. fugitive methane). The motivation of studying fugitive methane (CH 4 ) gas, stems from the fact that methane has a large green house gas (GHG) warming potential relative to CO 2 and occurs from both natural and anthropogenic sources. Modeling gas plumes in general can be done using a variety of tools (e.g. deterministic and stochastic approaches) which provide information about how the plume may develop given a set of initial and boundary conditions. However, these visualization require a certain degree of model fidelity to capture realistic behavior (i.e. plume meandering and trapping). Computational fluid dynamics approaches can be applied but suffer from high computational cost. Stochastic methods have shown better computational efficiency but for a given initial condition the output is an ensemble (due to the inherent randomness) and presents difficulty in implementing continuous time. On the other hand, hybrid modeling (which encompasses a combination of stochastic and deterministic methods) can provide computational efficiency while maintaining highest possible fidelity. In practice, computational efficiency of these models would need to be suitable for deployment on a laptop. Thus, there exists a trade-off between model fidelity and run-time. This is important for applying concepts such as DT to further improve on ES work.
What is a DT? Some of the earliest work of using physical twins (or hardware twins) can be seen as early as National Aeronautics and Space Administration's (NASA) Apollo program [16]. In this case there existed at least two space vehicles or rovers at once. This allowed for the continuous testing of code, maneuvers, and conditions before sending them to the deployed version. Software can also be introduced to the hardware twin, as in the 'iron bird' which combined a physical interface of a plane cockpit with an aerodynamics simulation for training purposes. The simulation was used for a visual aide as well as to compute the meteorological conditions and resulting forces on the plane [16]. Some of the early reports of the current concept of digital twin (as digital equivalence of a physical product) showed up in a product lifecycle management course at the University of Michigan in 2003 taught by Michael Grieves. Since then a white paper was written in reference to virtual factory replication [28]. It is explained that the digital twin (DT) in this scenario contains three parts: (1) physical products in real space (2) virtual products in virtual space and (3) and connections that tie together the real and virtual products. The use cases for DT are: (1) conceptualization -visualization of the physical and virtual products, (2) comparison -compare physical product behavior with the virtual product to improve results and (3) collaboration -solutions found can be applied immediately to other factories (or systems). For ES, the environment can be full of complexity and having a DT can provide a deeper understanding between the sensors and physical parameters such as source rate, location, and atmospheric conditions. Additionally, DT's have the ability to update the system parameters, over time, as well as give insight on performance metrics. This paper proposes a DT framework for ES from a point source emission as an extended version of the conference paper in [30]. A controlled release case study is used to solve an optimization-based behavior matching problem of the DT framework, using the estimated concentration signal, c m = f (x, v|θ) (where x is the position of the sUAS, v represents the wind field and θ is the DT parameters), and the measured concentration signal, c m , with a single sUAS measurement system. The mass balance method is applied to experimentally quantify the source emission rate. This paper extends the work in [30] by: defining DT levels and capabilities; adding additional details and figures of the DT model, approach, sensitivity of parameters, and how the observation problem is effected by sensor location (i.e. measurement operator); exploring additional spatial and spatial-temporal interpolation methods to quantify the source and compare data between the DT simulation and physical experiment; proposing the use of Mittag-Leffler function into the spatial-temporal gas distribution modeling framework. The paper is organized as follows: Section 2 describes the DT framework, Section 3 overviews the controlled release experiment, setup and behavior matching optimization, Section 4 we conclude on the simulation results, and Section 5 provides a discussion on future work to further this field. This paper proposes a DT framework for ES with a gas plume source (see Fig. 1).

The Proposed DT Framework
In order to discuss a framework there needs to be some agreed upon definition of DT. Over the years there have been several variants of the DT definition such as those shown in Table 1. The Mechatronics Embedded Systems and Automation Lab has since added a more general definition which is as follows, "A Digital Twin is the combination of multiple, individual, and detailed simulation models (continuous, discrete, hybrid), where its interconnection represents the dynamics of a complex system, which is updated periodically (windowed or real time) with the system information in order to reflect the system current status as well as predict its future behavior and possible faults". For ES, the model needs to be selected such that there is a trade off between fidelity and computational speed. Considering the quote from statistician George Box, "all models are wrong, but some are useful", we need the desirables of the DT model to be capable of: (D1) exhibiting relevant digital artifacts with suitable level of fidelity, (D2) short or long time-scale DT is a virtual, digital equivalent to a physical product Grieves [28] 2003 Up-to-date representation of an actual physical asset in operation Mathworks [5] 2019 Dynamic virtual representation of a physical object or system, usually across multiple stages of its lifecycle. It uses real-world data, simulation or machine learning models, combined with data analysis, to enable understanding, learning, and reasoning.

IBM [2] 2020
Software representations of assets and processes that are used to understand, predict, and optimize performance in order to achieve improved business outcomes.

GE Digital [6] 2019
DT is a perfect digital copy of the physical world: a digital twin. This twin would enable you to collaborate virtually, intake sensor data and simulate conditions quickly, understand what-if scenarios clearly, predict results more accurately, and output instructions to manipulate the physical world..

Deloitte [3] 2020
An integrated Multiphysics, multiscale, probabilistic simulation of an as-built vehicle or system that uses the best available physical models, sensor updates, fleet history, etc., to mirror the life of its corresponding flying twin

NASA[26] 2012
A digital twin is a multi-faceted dynamic set of smart digital models of a system or a subsystem along with all its constituents, which accurately represent the design of a product, production process or the performance of a product or production system in operation.
Dufour et al. [23] 2018 DT are precise, virtual copies of machines or systems driven by data collected from sensors in real time, these sophisticated computer models mirror almost every facet of a product, process or service.
Tao et al. [46] 2019 characteristics evolves along with the real system given the same conditions, and (D3) can be used to solve relevant problems for the system or improving methodology [30].
The DT concept can also be broken into four different stages or levels based on their current implementation: (L1) Pre-Digital Twin, (L2) Digital Twin, (L3) Adaptive Digital Twin, and (L4) Intelligent Digital Twin. Each stage has increasing levels of model sophistication. Aside from L1, all of the levels have a physical twin. As the DT levels go from L2 to L4 the data acquisition from the physical twin increases. Different system health and performance parameters can be extracted for predictive maintenance (from batch updates to real-time updates). Additionally, these stages can start to incorporate machine learning (ML) in operator preferences as well as for the system/environment [34].
The proposed DT framework for a point source first starts on L1 with a proposed model for representing the complex system. In the case of gas concentration, the system is governed by a partial differential equation such as the advection diffusion equation, where the C represents the concentration of the gas, D is the diffusion coefficient, and v = [u, v] T is the wind vector field. The source can be represented using Kronecker delta function δ s = δ(x − x s )δ(y − y s )δ(z − z s ), where x s = [x s , y s , z s ] T is the location of the point source. The wind field parameters can be solved generally by using the something like the incompressible Navier-Stokes equation, Here ν is the kinematic viscosity, w is the thermodynamic work, and g is the gravitational acceleration. Solving these systems consecutively can be expensive and thus, some assumptions need to be applied to simplify the models. One such approach to solve this problem was proposed by Farrell et al. [25] and is summarized here for completeness. The first assumption is to break the concentration dynamics into discrete filaments (or packets of molecules) based on different length scales (i.e. large scale advection (a), intermediate scale turbulent mixing and stirring (m), and local small scale diffusion (d)). The transport of the k-th concentration molecule can then be described by, However, the small scale diffusion term can be absorbed into the shape of the filament, such that the filament's i-th location can be described as, Now the only thing left to resolve is the wind vector field. One approach is to look at average wind vector such that v =v + v , where the overbar (·) represents the average and the prime ( ) denotes the deviation from the average. Substituting the decomposition of the wind vector into (2) and using the following assumptions the Navier-Stokes equation can be drastically simplified: 1. Coriolis forces, geostrophic winds, and molecular viscosity can be deemed small and neglected; 2. Measurements are conducted close to the ground where winds are relatively constant w.r.t. the altitude.
The wind field can then be represented by the following equations that resembles the viscous Burgers equation, This is realized by using the simplest K-closure method, which relates the diffusivity (K x,y ) to, u u = −K x ∂u ∂x and similarly for the v v , and v u . For the simulation aspect, the wind field meandering can be produced by perturbing the flow with colored noise using a second order transfer function approach, H (s) = Ga/(s 2 + bs + a) driven by white noise [25]. This colored noise is updated on the corners of the finite grid and adjacent boundary conditions are linearly interpolated (see Fig. 2). The meandering behavior (see Fig. 3) can be adjusted by scaling a, b, G and K x,y . The interior nodes can be solved implicitly [37] or explicitly by finite differences [40].
The large scale advective term, v a , can be calculated by using a bilinear interpolation between the nearby wind vectors points or estimated using the nearest neighboring point relative to the filament position, x i . The grid size used to calculate the the wind field can remain course (large Δx and Δy) as it reflects the larger length scales of the model. Additionally the large spatial separation between points helps to increase simulation runtime. To maintain efficiency we used the nearest neighboring point for this work. To calculate the intermediate scale term, v m , we use a random process satisfying, where ω ∈ R n and η is a white noise random process with spectral density σ 2 η . The matrices A , B , C , and D are appropriately sized. The resulting standard deviation grows as Gσ η √ t. In the simplest 2D case D = 0, (B ) T = C = [1, 1], and A = I , where I is the identity matrix.
To calculate the small scale diffusion, v d the use of a growth term, R(t), is needed. This growth term represents the standard deviation of the packet of molecules inside the discrete filament, where r i (t) = ||x − x i || is the Euclidean distance from the sensor and the i-th filament, R i0 = R i (0) and γ where Q is the filament source rate in molecules per filament. The total concentration measured, C(x, t), can be calculated given an effective sensor area, A e , where N represents the number of filaments inside A e = πr 2 e .
The effective area can be approximated as a circle with radius r e . The sensor can be modeled as a low pass filtered system with thresholding,ċ = f BW (C(x, t) − c). The sensor bandwidth is f BW and threshold is λ. The measured concentration, c m would then be, This current model can represent the gas plume as 2D flow. However, in most scenarios the gas needs to be represented as a 3D flow usually with building resolving capabilities. To start we will examine just rural case with no buildings. Furthermore solving the wind field for 3D case results in an increase in computation because the number of discrete points grows as 1/Δ 3 In [30], the model was extended to 3D by utilizing a detection probability based on experimental observations in [31,45]. The experiments showed that the detection probability varies as a function of altitude and wind speed. An example of this can be seen in Fig. 4, where the sigmoid function can be fit to the probability of detection given different wind conditions or atmospheric stability classes.  (2017) under controlled methane release of 5 standard cubic feet per hour (SCFH) at 30m downwind [31] By using this knowledge one can essentially extend the 2D model to 3D. Given that sensors do not always detect the target gas, we can introduce a conditional detection probability P d . Letting the P d be a function of downwind distance d, altitude z, and wind velocity |v a |. Making the assumption that the lateral or cross wind position does not affect the probability with respect to the simulation, the relation becomes, P d (δ d |z) = f (d, z, |v a |). Given the nature of the mass balance method, the downwind distance d does not readily change, and therefore can be approximated as Here M is the slope at the half probability point, and z bias is the altitude at which half probability occurs. From the sigmoid curves in [31] (see Fig. 4 we can alleviate this limitation. Although this only works if the source location and downwind distance are known ahead of time. Additionally, there is a relationship between the value of M, z bias , and the atmospheric stability classes that needs to be further investigated. The atmospheric stability classes (A-F), first introduced by Pasquill in 1961 and later reformulated by Gifford, are typically referred to as the Pasquill-Gifford (PG) curves [32]. The PG curves relate wind speed and solar irradiance to the dispersion coefficients that govern the spread of the plume. This is represented in the Gaussian plume model, or sometimes referred to as the infamous 'bell curve', where d l is the lateral (or crosswind) distance from the plume centerline, H = h + ΔH is the effective height of the plume, h is the stack height, Q p is the source rate in g/s and ΔH is the plume rise [7]. The plume rise term can be calculated using models such as Holland's formula which is dependant on temperature and wind speed. However, Holland's formula was designed for exhaust stacks and also takes into consideration the pressure, stack diameter, and temperature of the ambient air and stack. As a result if the stack velocity is small, such as with sub-surface leaks that diffuse through the ground, the plume rise term tends toward zero. Therefore other factors should be considered, such as lapse-rate, which takes into consideration how the temperature changes as a function of altitude. Furthermore, for small leaks, even surface temperature may come into effect by introducing vertical momentum through natural convective forces. Holland found that there is a small correlation between the plume rise and the temperature gradient near the ground. Briggs mentioned that the plume rise is only dependent on the temperature gradient of the air to which the plume is rising. It is also mentioned that the gradient near the ground is not a good representation of the gradient at higher altitudes [17]. Thus, this part of the DT framework will be reserved for future research. The simulation source rate,Q * , given in molecules per second, and is determined by the number of filaments that are discretely released from the point source location. The particular way the filaments are released depends on two parameters, the number of filaments per puff n f and the number of puffs per filament n p . The filament source rate can then be defined by Q =Q * /(n f n p ), which is used in (5). The different combinations of these parameters can result in qualitative changes in the measured signal from the sensor. This qualitative behavior can be matched to the measured signal of the physical system and sensor (i.e. behavior matching). The behavior matching can be undertaken given an observable data set that captures these qualitative changes within the model. The key point here is observable data and thus the question becomes, how to get it? The answer depends generally on two things: sensing in the right place and sensing at the right time. Given a set of observations s we want to infer the states or plume field, x i . This observation problem can be formulated as a state estimation problem using orthogonal decomposition methods based on singular value decomposition(SVD) [24]. Essentially the problem can be represented as a L 2 minimization, where represents the orthogonal modes and ν the associated coefficients. This approach can yield solutions that are often ill-posed and not necessarily unique. However, it can be shown that the estimation error depends on the approximation basis (othrogonal modes) as well as the measurement operator, H, The symbol + represents the Moore-Penrose pseudoinverse. The measurement operator depends entirely on the sensor location and not the orthogonal modes. Therefore, the choice of sensor location is an important task. For the single sUAS measurement system, this comes in the form of path planning and is subject to disturbances from changes in the stability of the wind. When applied to the mass balance method, path planning involves scanning back and forth at different altitudes inside a region of interest. The size of this region can incorporate twice the expected plume dispersion (deviation from the plume centerline) σ d l and σ z , such as in [31] to ensure a likely encapsulation of the plume (see Fig. 5). This approach aims to capture rich signal data for conducting the mass balance method. However, this process takes time to complete with a single sensor system, and thus can present spatial-temporal issues when reconstructing the mass flux.
To advance the pre-Digital Twin model (L1) to a Digital Twin (L2) we need to have the physical twin exist and have the virtual system model of the physical twin. This step requires a system identification like optimization of the pre-Digital Twin model, given a dataset with observable properties and tunable parameters to minimize on. Once this minimization is complete the Digital Twin will allow for testing different path planning strategies to improve mass balance method results. The tunable parameters (see Fig. 6) in this DT model consist of: the effective radius r e , the number of filaments per puff n f , the number of puffs per second n p , the intermediate scale mixing and stirring dispersion σ η , and a correction factor c f where the corrected simulation source rate is (Q * ) c = c fQ * . The value ofQ * is chosen to be the same as the physical source rate, such that Q = (Q * ) c /(n f n p ), and c f is determined through behavior matching optimization.

Case Study
The experimental data used to convert the L1 Digital Twin to L2 was based off a controlled release single source field experiment in the Merced Vernal Pools and Grassland Reserve just north of the University of California, Merced campus. The source setup included a pure methane bottle connected to a pressure regulator, followed by several meters of PFTE (Polytetrafluoroethylene) tubing, where a rotameter flowmeter gauge was placed for controlling the release rate. Additional PFTE tubing was used to connect the output of the flowmeter to a 5 Gallon bucket filled with coarse rocks (≈ 1 cm in diameter). The rocks served to slow the stack velocity down and let the wind carry the methane, simulating a surface point source. The source rate was initially set to 10 SCFH, whereas after laboratory testing the source rate was calibrated to be Q * = 12.9 SCFH. The landscape of the site was relatively flat with little topology and no surface obstructions (see Fig. 7a). A sUAS, namely the DJI M210 (see Fig. 7b), was equipped with the OPLS sensor, a highly sensitive spectrometer (detection rates of 10ppb s −1 ) and acquisition rates of 5Hz [22]. The atmospheric stability during the experiment was determined to be C and B (slightly to moderately unstable) using the PG curves [32]. The wind speed and direction was measured with a RMYoung 3D ultrasonic anemometer. The experiment consisted of 16 downwind flight paths of horizontal transects ascending from ground level till a clean air fetch was reached, referred to as a curtain (see Fig. 5).
Once the experimental data is collected, the most data information rich flights are to be selected to perform behavior matching with the L1 digital twin. In this Fig. 6 The timeseries behavior and sensitivity with tuning parameters θ given one pass through the plume Fig. 7 (a) Experiment location at Merced Vernal Pools and Grassland Reserve. The yellow star indicates the source, the green triangle is the RMYoung anemometer, and the flight area is shown red dotted line. (b) DJI Matrice 210 (M210) equipped with a sensitive methane sniffer developed by [22] and attached to the DJI Skyport. (c) Example mass balance measurement for one curtain flight experiment, some of the curtain flights (such as flights 13 and 14) resulted in no detected emissions and cannot be used to behavior match the model (see Table 2). After analysing the 16 flights using the mass balance method (see Fig. 7c), flight 1 showed stable wind conditions, consistent methane detection throughout the ascending transects, and good estimation of the true source rate. Thus, flight 1 was chosen for behavior matching the sensor signal of the L1 digital twin to the experimental sensor signal (see Fig. 8). The simulation source rateQ * was set to (Q * = 2.4 × 10 24 molecules/s). Using the PG curves given the current meteorological conditions (i.e. atmospheric stability) the initial condition of σ η was chosen to be 8 m. The wind speed and direction measured from the experiment was used as a direct input into the L1 digital twin, using the mean wind field approximation. Using the behavioral characteristics of the hyperparameters shown in Section 2 and in an attempt to reduce complexity, let n f = 1, such that n f n p has units as filaments per second. Heuristically tuning the hyperparameters, to give a best guess of θ = [c f , r e , n p , σ η ], the following initial conditions were chosen: c f = 0.2, r e = 5 m, n p = 40 puffs/sec, and σ η = 8 m. It can be also seen that there is a lag present in the system (see Fig. 8) which can be attributed to model inconsistencies as well as initial conditions of the plume. Since there is currently no implementation of plume rise (future research), the average altitude where concentration detection occurred was used to determine the z bias term and the vertical dispersion was set to σ z = 5 m. Fig. 8 The concentration time series in ppm during flight curtain 1, where n p = 40 puffs/sec and the symbols τ gi , τ pi , and τ lag represent the i-th gap length, i-th pulse length, and the time lag between the simulation and experimental data [30] If we represent the L1 digital twin given hyperparameters θ as a function at time t, with the sUAS position (i.e. sensor position) x and the wind speed and direction measurement using the ultrasonic anemometer (RMYoung) v. The output consists of a concentration measurementĉ m (t) and filament positions x i . In the physical experiment the locations of the The flow rates (experimental Q * and simulationQ * ) are in SCFH, the subscripts denote the spatial interpolation method, and the superscripts are the associated parameters, where K is Kriging, I DW is inverse distance weighting, and KDM is Kernel DM+V plume filaments are not observable given that methane gas is not visible to the human eye. The plume can potentially be observed using a large number of high quality sensors, of which, is expensive, difficult to set up, and only be inferred using spatial interpolation techniques. Thus, the only output available to optimize over is the concentration measurement, c m (t), in the least squares sense. In order to compensate for the stochastic behavior of the plume dynamics, the cost function was designed to average n T = 20 trials. Unfortunately, this leads to multiple possible solutions to the minimization. Therefore, additional constraints on the hyperparameter space, Ω, and regularization can be added to constrain the magnitude of the concentration measurement such that, where c max = max(c m ). The cost function weights were chosen to balance the effects of each term and were determined from trial and error. The solution to (15) can be found using a variety of optimization methods and techniques, such as Gradient Descent or Stochastic Gradient Descent. Faster convergence can be achieved using Hessian based approaches but due to the cost of the forward model calculation and considering complexity of implementation these methods were not considered. Other gradient free approaches are more desirable for this problem such as Genetic Algorithm, Extremum Seeking Control, or Pattern (Direct) Search can be applied. With simplicity in mind, the Pattern Search Optimization was used here to behavior match the digital twin to the experimental data. After 65 iterations and over 390 cost function evaluations (7880 forward model runs) the gradient of the cost function dJ /dt ≈ 0 (see Fig. 9). The model parameters were found to be: c f = 0.11, r e = 7.5 m, n p = 78 puffs/s, and σ η = 8 m.
To calculate and perform the mass balance method, sometimes called the box method, requires a way to measure the flux entering or leaving a control volume. Typically, upwind and downwind curtain flights are used to measure the flux. However, in the case of a controlled release the downwind curtain flight is sufficient to calculate the flux, granted the plume is encapsulated in the flight path (see Fig. 5). Using the timeseries vector c m = [c m (Δt), c m (2Δt), . . . , c m (n t Δt)] T with the corresponding measurement locations, {x} n t 1 , a concentration matrix, C m can be estimated on an n by m spatial grid through spatial interpolation methods. The flux calculation can be then given as the integration over curtain area with the average wind vector multiplied by the cosine angle of the normal vector to the curtain with the enhanced concentration, Here C b represents the background concentration measurement. An example of a flux plane calculation can be shown in Fig. 7c. Several kinds of techniques can be applied to the timeseries data to generate C m , such as Ordinary Kriging, Inverse Distance Weighting (IDW) or statistical gas distribution modeling (e.g. Kernel DM+V/ W). For the case of Kriging, a linear unbiased estimator, can be solved using a  . 11 Here the turbulence intensity, T i = σ u /ū, of the wind is plotted against the normalized flow rate, which trends towardQ * /Q = 1 as T i → 0. A boxplot of the source estimation for each flight (set of 10 simulations), given a spatial interpolation method, is shown for each T i condition observed from the experimental data. A smaller boxplot (shown within each boxplot) contains tails that represent the uncertainty. We can observe 1) several flights agree with experimental data 2) majority of flights underestimate the flow rate in short due to increased T i and 3) the remaining flights that disagree are likely not capturing vertical plume behavior variogram (or semi-variogram), Generally, the redundancy and closeness covariance matrix are calculated from fitting the experimental variogram to a function with respect to the spatial distance or lag, such as the exponential or Gaussian, and knowing the covariance of the measured points (referred to as the sill). This method requires second order stationarity in the spatial observations [49]. Unfortunately this is only approximately the case during stable atmospheric conditions. Methods such as IDW look at the inverse distance when determining There has been some efforts to improve this by looking at minimum error variance within IDW [12]. Wind information has been incorporated into Kernel DM+V/W approach in 3D [42]. Another adaptation to statistical gas distribution mapping approach was implemented in [29] by using the Gaussian plume kernel, which incorporates the wind information into the covariance function in determining the weights. The downside to this approach is that it assumes the source location. This is not ideal for determining the flux in this case. The spatial-temporal effects on modeling have been looked at as a separation between the spatial covariance and the exponential decay between temporal measurements [11,53], where N is a Gaussian weight function, (·) (k) represents the k-th grid point, {x} i are the sensor measurements at time t i and ϕ is given as, However, when dealing with complex systems such as in turbulent environments there may be several factors that determine what value of a ϕ is. If we look at the weight summation of exponentials, with elements increasing towards infinity, it can be described by fractional order dynamic behavior [35]. The Mittag-Leffler function [44] can be used to describe this behavior, such that when α = β = 1, E 1,1 (t) = e t . Given that β = 1 and letting α decrease we can see in Fig. 10 that the tail behavior decays much more slowly than the exponential function. By combining the approach in (19) with (21) we can apply a new weighting scheme with tuning parameter, α, The value of t * can be chosen for each grid point based on the time t i of the closest measured point that minimizes |{x} i − x (k) |. For the single sensor mapping problem we test the change in sensitivity for α = 1, 0.5, and 0.1. After optimization of the L2 digital twin, using the experimental position and wind measurements, each flight was simulated 10 times. The spatial interpolation was carried out using Kriging, IDW, and Kernel DM/ V to see a comparison in quantitative results across the different approaches. It can be observed that all the spatial interpolation methods tend to perform similarly. Also, for the spatial-temporal method using (22) does not show noticeable performance improvement. Future work can be done to address the optimal configurations of these methods and their hyperparameters. It can be observed that the source rate estimations fluctuate in value depending on current weather conditions and turbulence (see Table 2). Given small deviations in the wind, the sUAS curtain flight path can appropriately capture the spatial data through the mass balance plane. When the wind deviates slightly, the data we capture can sometimes be spread out over the flux plane (i.e. like a blurred image). As a result, this can lead to higher source rate estimations or to poor detection and characterization of the plume. Since the wind direction fluctuations correlate to different stability classes it makes sense that there should be some correlation with the turbulence intensity as well (see Fig. 11). Thus, mindfulness of the current atmospheric conditions are important when planning and conducting these flights. Furthermore, the wind behavior, to some extent, is always changing in magnitude and shifting in direction. It seems this behavior is very rarely calm and stable, which may put physical limitations on the sensitivity of accurate spatial sampling with single mobile sensor systems (at least in the case for air). Comparing the experimental data to that of the digital twin one can observe qualitative similarities in the distribution of sensor detections, as seen in Fig. 12. Additionally, larger less frequent spikes at lower altitudes were observed that are consistent with what was experimentally observed (see Fig. 13).

Conclusion
In this work we develop a L2 DT model for environmental sensing application of a point source emission on the ground using sUAS. We give background details on spatial interpolation and modeling of sparsely measured concentration points, and also propose the use of the Mittag-Leffler function for capturing fractional order dynamics of the spatial-temporal gas distribution modeling. We compared the results from [30] with different spatial interpolation and modeling techniques. The results (shown in Fig. 11) indicate that the choice of spatial interpolation and modeling technique is less sensitive to method choice and more dependent to the DT model itself. The average mass balance calculation (excluding flights 13 and 14) of the source rate for the experiment using Kriging was determined to be 6.4±2.4 SCFH, and the average mass balance calculation of the source rate for the simulation results are 7.7±2.6 SCFH for Kriging, 8.1±3.4 SCFH for Kernel DM+V, 7.8±2.7 SCFH for IDW(p = 2), 7.8±2.7 SCFH for IDW(p = 1.5), 7.9±2.7 SCFH for IDW(p = 1), 8.1±3.3 SCFH for TD Kernel DM+V(α = 1), 8.1±3.4 SCFH for TD Kernel DM+V(α = 0.5), and 8.1±3.4 SCFH for TD Kernel DM+V(α = 0.1). Although the DT simulation results do not reflect the true source rate, they do reflect the experimental emission estimates, which is the goal of the DT. The qualitative behavior of the DT was also observed to be comparable to the mass balance calculations of the experiment as a function of turbulence intensity, with more accurate estimates when atmospheric conditions are stable.
This work also provided an extension to the conference paper in [30] by: (1) defining the DT and the associated levels and capabilities; (2) extending the details on the L2 DT model used in this study including the figures showcasing the flight path, probability of detection per transect pass, and sensitivity of model parameters in the timeseries behavior matching process; (3) giving an explanation of the behavior matching problem in the context of optimization and the measurement operator H ; and (4) extending the details on the case study setup and choice of flights for behavior matching.

Future Efforts
Future research into making the digital twin (DT) smart is desirable as well as getting the digital twin to L3 (allows for real-time updates between the physical system and DT as well as utilize machine learning operators) and L4 (DT model utilizes adaptive techniques such as reinforcement learning between the physical system and DT for understanding system performance and features) [34]. For the DT to be smart it needs to have 5 characteristics such as Cognizant (aware of capabilities and limitations), Taskable (handle high-level, often vague, instructions), Reflective (learn from experience to improve performance), Ethical (adhere to a system of societal and legal norms), Knowledge Rich (reason over a diverse body of knowledge) [1]. Efforts in cloud computing can be leveraged with the L3 and L4 DT's to make the system smarter with tools such as Microsoft Azure and Amazon Web Services.
Just like the controlled release scenario in an open environment [30,45] or in a real world scenario with active well sites and obstructions [50][51][52] there is a need to understand benchmarks for a given environment. For example, a full scale experiment in Oklahoma City such as Joint Urban 2003 test [9], where a sulfur hexafluoride (SF 6 ) gas was released from a known source with several sensors located throughout the city. The resulting experiment was coupled with high fidelity modeling, such as FEM3MP [19], to explore source estimation techniques [21]. However, the data sets are not readily available for others to benchmark. Other test scenarios such as the Mock Urban Setting Test (MUST) have been used in similar explorations [33]. In both of these examples the sensors are fixed and there is little to no public availability of the data sets (i.e. only shared among participants). This presents a need for a benchmark data set or framework (such as the DT) to provide a scenario based environment to test source term estimation techniques etc. In addition, open source data archiving of the data set enables reproducibility within the research field.
Single mobile sensor systems (e.g. sUAS with OPLS) have inherent spatial-temporal issues that limit gas distribution mapping usage to certain (stable) atmospheric stability classes, otherwise the accuracy can degrade and uncertainty grow. Utilizing a grid of sensors will ultimately improve the field reconstruction efforts, as it can better separate the spatial-temporal effects and provide opportunities for sensor location optimization [48]. A group of mobile sensors can be implemented and controlled using techniques such as Central Voronoi Tessellations [20]. Furthermore, the desired sampling locations can potentially take advantage of adaptive sparse sampling methods [39] based on information metrics such as information entropy. Data driven approaches, given rich sampling locations can then be potentially applied (e.g. Dynamic Mode Decomposition [43], Compressive Sensing [14]). Even machine learning techniques may then be available such as shallow neural networks with sparse sampling [24] or physics informed neural networks (PINN) [41]. PINN's have been used in similar problems such as semantic inpainting [54]. Utilizing the power of machine learning with physical constraints combined with adaptive spatial sampling techniques could lead to more accurate reconstruction of flux planes, perhaps providing better source estimation and localization results overall.
Acknowledgements The authors would like to thank Jairo Viola for his helpful discussions on digital twins. The authors also thank the reviewers for helpful comments in improving this manuscript.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.