Implementation of a Lagrangian Stochastic Particle Trajectory Model (LaStTraM) to Simulate Concentration and Flux Footprints Using the Microclimate Model ENVI-Met

the Abstract: The number of studies evaluating ﬂux or concentration footprints has grown considerably in recent years. These footprints are vital to understand surface–atmosphere ﬂux measurements, for example by eddy covariance. The newly developed backwards trajectory model LaStTraM (Lagrangian Stochastic Trajectory Model) is a post-processing tool, which uses simulation results of the holistic 3D microclimate model ENVI-met as input. The probability distribution of the particles is calculated using the Lagrangian Stochastic method. Combining LaStTraM with ENVI-met should allow us to simulate ﬂux and concentration footprints in complex urban environments. Applications and evaluations were conducted through a comparison with the commonly used 2D models Kormann Meixner and Flux Footprint Predictions in two different meteorological cases (stable, unstable) and in three different detector heights. LaStTraM is capable of reproducing the results of the commonly used 2D models with high accuracy. In addition to the comparison with common footprint models, studies with a simple heterogeneous and a realistic, more complex model domain are presented. All examples show plausible results, thus demonstrating LaStTraM’s potential for the reliable calculation of footprints in homogeneous and heterogenous areas.


Introduction
The number of study sites for the measurement of the surface-atmosphere exchange of energy and matter such as greenhouse gases (i.e., carbon dioxide, methane, nitrous oxide) has grown substantially during recent years [1][2][3]. Flux data are used to increase process understanding of the terrestrial carbon cycle, to constrain models, and to inform policy makers, e.g., in the context of climate change research and for modeling purposes. The stateof-the-art method which is used at measurement sites is the eddy covariance technique. For an accurate interpretation of the flux data, it is necessary to estimate the surface area that influences the signal measured at the detector location, which is the so-called flux footprint or flux source area [4][5][6]. The concentration or flux footprint of a quantity is a function of the characteristics of the Earth's surface, the state of the atmosphere as well as the height of the detector above ground level. At the surface, the distribution of sources and sinks (homogenous/heterogeneous), the surface roughness length, and the topography are important factors in shaping the footprint's location and dimensions. Atmospheric stability, wind speed and direction, as well as turbulent intensity are atmospheric factors influencing the footprint. For this purpose, several models have been published since the 1990s which can be distinguished into four categories: (1) analytical models [7][8][9][10][11][12]; (2) Lagrangian stochastic particle dispersion models (LS) [13,14]; (3) large-eddy simulations (LES) [15][16][17]; and (4) closure models [18][19][20].
Most of the available footprint models are valid only in case of a homogeneous source and in flat topography. This includes all models based on analytically derived mathematical expressions to calculate advection and diffusion of the quantity of interest. LES footprint models [16] and closure models [20] overcome this limitation, i.e., they are able to simulate fluxes in heterogeneous environments, e.g., high vegetation, forests, and urban areas. One characteristic of the first mentioned type of models is the high computational demand. This means that only computer clusters with hundreds of processors are able to calculate footprints in reasonable time. In contrast, Reynolds-Averaged Navier-Stokes (RANS) models have lower computational demand and can hence run on personal computers. However, as Hellsten et al. [16] state there is no backward-integration approach for any building-resolving LES-LS footprint predictions available yet. Backward-integration approaches represent a backward tracking of particles from the target (detector) to the source (emission source). Compared to forward tracking from source to target, where a multitude of particles is tracked without ever reaching the detector zone, backward approaches feature a significantly higher computational efficiency as the number of particles eventually reaching the detector zone is already determined at simulation start. That allows us to define smaller detector zones to be thus comparable with real measurement devices [14].
In the present study, the microclimate model ENVI-met [21], a RANS model, was coupled to a newly developed Lagrangian model called LaStTraM (Lagrangian Stochastic Trajectory Model) to simulate particle trajectories with a backward-integration approach based on the microclimate's three-dimensional output data in order to close this gap. Since ENVI-met features a sophisticated plant physiology stomata model (A-gs; [22]) to simulate transpiration and carbon uptake of plants, the influences of heterogeneous sources and source distributions (e.g., trees, shrubs, grass) on the footprint can be modelled in complex topographies (e.g., within dense and complex urban areas) [23,24]. LaStTraM works with all ENVI-met versions including the free LITE version.
To evaluate the capability of the newly developed post-processing Lagrangian model to simulate footprints in a homogeneous scenario, the model results are compared to a wellknown Lagrangian model from Kljun et al. [9,10] and an analytic model from Kormann and Meixner [11]. Furthermore, the influence of heterogeneous sources in more complex model areas on the concentration and the flux footprint is discussed.

Model Description
In an offline coupling, all necessary atmospheric parameters for the calculations of LaStTraM are taken from simulation outputs provided by the microclimate model ENVI-met. The model ENVI-met is a prognostic, three-dimensional, high resolution microclimate model. With its physical fundamentals that are based on the principles of fluid mechanics, thermodynamics, and the laws of atmospheric physics, it is able to calculate three-dimensional wind fields, turbulence, air temperature and humidity, radiative fluxes, and building physics [25,26].
One of the key features of ENVI-met is the detailed modeling of vegetation. With its high spatial resolution (<5 m), ENVI-met allows the simulation of the individual photosynthesis rates, taking into account local solar radiation, air temperature and humidity, wind speed, and CO 2 concentration, among other parameters [23].
The spatial and temporal evolution of the wind field is calculated by applying the non-hydrostatic three-dimensional Navier-Stokes equation. Using the Boussinesq approximation, the density of air can be assumed to be constant [21,25,27]. Turbulence in ENVI-met is parametrized using a E − ε 1.5 order closure model. The E − ε model basically consists of two prognostic equations, one describing the production of turbulent kinetic energy (TKE) and the other its dissipation. In contrast to first order closure models, the E − ε model allows the simulation of advective processes in horizontally inhomogeneous environments without as much computation time as closure models of higher order [25,28]. The trajectory model horizontally and temporally interpolates all atmospheric parameter outputs from ENVI-met, making it quasi-independent of the ENVI-met model resolution. Both the interpolation as well as the advection of the particles are computed in parallel, making use of all available processor cores. The particles are released simultaneously at the start of a simulation. The interval in which simulation outputs are generated can be set by the user.
To simulate particle trajectories, LaStTraM operates with a backward-integration approach. The particles are "emitted" from the detector and are followed until predefined conditions are met, e.g., particles reach coordinates where plants or sources of different compounds, such as CO 2 , are located. The particles' starting position can either be set to the center of the detector or can be randomly generated inside the detector's dimensions. The number of particles to be simulated can be freely set by the user. All particles reaching certain areas are individually saved with various pieces of information such as their current three-dimensional position as well as their take off and touch down velocities.

Dynamic Interpolation of the 3D Model Domain
To reduce the dependencies of the simulation on the ENVI-met model resolution, all ENVI-met model output variables are dynamically interpolated for every particle and every timestep. The interpolation is carried out using inverse distance weighting (IDW) of all adjacent atmospheric ENVI-met grids to the particle in a 3D Moore neighborhood analysis, i.e., a maximum of 27 neighboring grids are taken into the interpolation [29]: where i, j, k are the ENVI-met grid coordinates of the particle,R int as the interpolated value, R(i, j, k) are the ENVI-met model value of the in the grid i, j, k and w(i, j, k) are the weight corresponding to the grid cell i, j, k calculated by: where d(i, j, k) are the Euclidean distance calculated using the 3D Pythagorean Theorem between the particles' actual position and the center of all adjacent atmospheric ENVI-met grid cells. Summing up all individual weights, the result must be equal to 1.
In case the particles' vertical position lies below the center of the lowest ENVI-met atmospheric grid cell, the interpolation of the wind vectors u, v, w is carried out using a logarithmic approach: where u k 0 , v k 0 , w k 0 are the wind speed at the lowest ENVI-met atmospheric grid cell, z 0 is the roughness length of the surface, and h k 0 ,c is the height of the center point of the lowest ENVI-met grid cell, k = 0.

The Lagrangian Stochastic Model
Based on Weil et al. 2004 [30], the particles advection is calculated by: P pos(x,y,z),t+∆t = P pos(x,y,z),t + P V,t ·∆t (4) where P pos(x,y,z),t+∆t is the new position of the particle after the time interval, ∆t, P pos(x,y,z),t is the original position of the particle, and P V,t is the velocity of the particle. The velocity consists of two parts: a resolved, deterministic part (P v res ) that is directly simulated based on the interpolated 3D wind vector fields of ENVI-met and a stochastic, random part (P V LS ): The resolved velocity part P v res is calculated by interpolating the u, v, w wind vectors for every particle position (see Equation (3)). The stochastic term of Equation (5), P V LS , is derived from the work of Weil et al., 2004 [30], where the considerations of Thomson, 1987 [31] are adapted. With these considerations and in accordance to Maronga et al., 2015 [32] and Steinfeld et al., 2008 [33], the difference in the particle's stochastic velocity is calculated by: where C 0 is a universal constant of the structure function (here set to 3; as in [14]), ε and e s as the inverse distance interpolated (see Equations (1) and (2)) dissipation and TKE, respectively, dξ as a random displacement vector whose components are independent from each other and in time [32][33][34]. The random number is calculated using a normal (Gaussian) distribution where the interval is limited to ±5.0σ. The parameter f s defines the subgrid-scale TKE fraction of the total TKE [30,33].
For the conditions at the boundary surfaces of buildings and digital elevation models (DEM), the simulation can be operated in two modes: an iterative solution mode, where in case a particle's new position was located inside a surface (buildings, DEM) the advection (Equation (8)) would be calculated with new random terms (( f s C 0 ε) 1 2 dξ) until a new position outside a solid structure is found or the maximum of 300 attempts is reached. In case no solution is found, the particle's position will be set to the previous time step. In the second mode, the termination mode, all particles that enter a solid surface (buildings, DEM) are no longer tracked and removed from the simulation. Regardless of the boundary surfaces mode the model operates with, the particles that would enter ground surface will be set to a minimum height of the user's choice (default value 0.0001 m). All particles leaving the lateral or top boundary of the model domain are removed.
While the resolved, deterministic velocity part of particles could be calculated in quite large time steps depending on grid size and velocity (e.g., CFL-Condition), the calculation of particles' Lagrangian stochastic velocity demand much smaller timescales [30].
To ensure small time steps, an upper limit of at least dt LS /40 is defined [32].

Footprint Calculation
A range of information is stored for every particle that reaches its target. This information includes the starting position within the detector, the final position within an emitter, travelled distance, elapsed time, mean velocity, as well as initial and final wind vector velocity. Since all data are stored for every particle reaching its target, the resulting file might be quite large (e.g.,~5 GB for around 20,000,000 particles reaching a source area).
To determine the relative contribution of a source to a measured concentration or flux at a predefined detector area, calculations can be performed that determine which surfaces/objects in a model area contribute to a measured concentration and flux at a given volume in space using the data information of the particles. The calculations of the concentration footprint C and the flux footprint F are carried out in accordance to Vesala et al. [5] and Flesch [35]: where x, y, z are the grid coordinates, Q is the source's strength, N is the total number of particles, n is the number of particles reaching the detector zone.
represents the sum of all particles reaching the target zone, w Det,i the initial w component at a particle's release at starting point within the detector zone and w Emit,i is the w component of wind velocity at touch down at source area [5,36].

Simulating Footprints in Homogeneous Model Area-Comparison with Other Models
In order to evaluate the newly developed model, it is compared against two established and well-known footprint models: the Flux Footprint Predictions model (FFP) by Kljun et al. [14] and the analytic Kormann-Meixner model (K-M; Kormann and Meixner [11]). To prove the functionality of the newly developed model LaSTraM, the presented evaluation is held similar to the analysis of these models. It is hence based on two common meteorological cases (stable and unstable) and features the same corresponding detector heights (3 m, 20 m, and 40 m).
The analyzed model area is similar to the one used in Kljun et al. [14]. It consists of 400 × 400 × 40 grids with a horizontal resolution of 5 m and a vertical resolution of 2 m in a non-equidistant gridding, where the lowest grid cell is split into five subcells. The total extent of the area is thus 2000 m × 2000 m × 70 m. The model area features no building structures or other obstacles, the soil was digitized as "loamy soil", and the vegetation was set to grass of 0.2 m height throughout the whole model area, thus creating a homogeneous case.
The ENVI-met simulation was set to start at 6 am on June 23rd, the simulation ran for 24 simulation hours. The meteorological boundary conditions can be seen in Figure 1, the boundary conditions for wind speed and direction were constant and set to 2 m s −1 in 10 m height and 90 • over the whole simulation time. The post-processing analysis in LaStTraM was run using the model outputs of the ENVI-met simulations. Each simulation featured 10 8 particles. It took around 5.5 h to simulate the trajectories of the released 100 million particles in the 2 km × 2 km large area on a personal computer with an Intel Core i7-6700 CPU @ 3.40GHz (Table 1).    The footprint is sensitive to atmospheric stability, surface roughness, measurement height, and wind direction [5,14]. Since the atmospheric stability is a prognostic parameter of the microclimate model ENVI-met, the output files have to be analyzed to differentiate between unstable, stable, and neutral conditions. Analysis of the different output files showed an unstable atmosphere at 15:00 h and a stable atmosphere in the early morning hours around 3:00 h (see Figure 2). In order to compare the LaStTraM results with Kljun et al., 2002 [14] and Kormann and Meixner, 2001 [11], the output files of these hours have been used for the further analysis.  The footprint is sensitive to atmospheric stability, surface roughness, measurement height, and wind direction [5,14]. Since the atmospheric stability is a prognostic parameter of the microclimate model ENVI-met, the output files have to be analyzed to differentiate between unstable, stable, and neutral conditions. Analysis of the different output files showed an unstable atmosphere at 15:00 h and a stable atmosphere in the early morning hours around 3:00 h (see Figure 2). In order to compare the LaStTraM results with Kljun et al., 2002 [14] and Kormann and Meixner, 2001 [11], the output files of these hours have been used for the further analysis.  Apart from the friction velocity (u * ), all input parameters for the FFP and K-M model could be directly extracted from the output data of ENVI-met. The friction velocity (u * ) was then calculated by: where u is the wind speed, κ the von Karman constant (here 0.4), h is the height above ground, and z 0 as the roughness length. Table 2 shows all input parameters used to run FFP and K-M models.

Simulating Footprints in an Inhomogeneous Model Area
By coupling LaStTraM with ENVI-met's three-dimensional model outputs, the model is able to simulate footprints in model areas featuring obstacles such as buildings or vegetation. To evaluate the influence of obstacles on the footprint, a simple inhomogeneous and one complex realistic scenario were analyzed.
In the simple inhomogeneous scenario, a symmetrical obstacle (building) with a height of 20 m and a tree is added to the homogeneous model area. The detector is then placed in 3 m above ground at different distances to the obstacle (125 m and 820 m upwind distance from the obstacle) (see Figure 3).

Simulating Footprints in Realistic Model Area
In a realistic scenario, the performance of LaStTraM is evaluated in an urban model area featuring buildings, different soil, and surface materials, as well as vegetation (see Figure 4). While this scenario, due to the lack of empirical data, cannot be seen as proving  The simple inhomogeneous model run featured identical meteorological boundary conditions as the homogeneous model run. To reduce the computational cost, the model was only run for the unstable case at 15:00 h.

Simulating Footprints in Realistic Model Area
In a realistic scenario, the performance of LaStTraM is evaluated in an urban model area featuring buildings, different soil, and surface materials, as well as vegetation (see Figure 4). While this scenario, due to the lack of empirical data, cannot be seen as proving that the model is operating accurately in a complex environment, the scenario is designed to be a showcase visualizing the capabilities of the model.

Simulating Footprints in Realistic Model Area
In a realistic scenario, the performance of LaStTraM is evaluated in an urban model area featuring buildings, different soil, and surface materials, as well as vegetation (see Figure 4). While this scenario, due to the lack of empirical data, cannot be seen as proving that the model is operating accurately in a complex environment, the scenario is designed to be a showcase visualizing the capabilities of the model.  The meteorological boundary conditions were taken from an Energy Plus Weather File for Frankfurt, Germany. To resemble a hot summer day with clear sky conditions, July 7th was selected as the simulation date. The wind speed was set to 1.5 m/s in 10 m height coming from 230 • . The simulation was run for 24 h. Figure 5 shows the meteorological conditions for the whole simulation period provided by the EPW-file. The meteorological boundary conditions were taken from an Energy Plus Weather File for Frankfurt, Germany. To resemble a hot summer day with clear sky conditions, July 7th was selected as the simulation date. The wind speed was set to 1.5 m/s in 10 m height coming from 230°. The simulation was run for 24 h. Figure 5 shows the meteorological conditions for the whole simulation period provided by the EPW-file.

Results and Discussion
To give an overview on the capabilities of LaStTraM, model results of the homoge-

Results and Discussion
To give an overview on the capabilities of LaStTraM, model results of the homogenous, simple inhomogeneous, and complex model areas are presented. To evaluate the performance of LaStTraM, the results of the homogeneous model area are compared against the FFP model and the K-M-model. Since no other model was available to compare the simple inhomogeneous and complex model area outputs, their results are only visualized, and their plausibility is discussed.

Evaluation of Homogeneous Model Area Results
The evaluation of the ENVI-met based footprint simulation LaStTraM with two wellknown models, an analytical model and a model that is based on a Lagrangian model [9,11], demonstrates that the modeled footprints are comparable for detector heights of 3 m, 20 m, and 40 m in stable and unstable conditions despite the structural differences of the models (see Figure 6). In the unstable cases, the maximum footprint contribution ( f max ) of LaStTraM tends to be higher than the other two models while the length of the footprint (x max ) tends to be slightly shorter in LaStTraM. This phenomenon featuring both an underestimation of f max and an overestimation of the footprint length x max of the FFP and K-M model was also attested by Heidbach et al., 2017 [37]. Similarly, Van de Boer et al., 2013 [38] concluded that the FFP model and two other models rather overestimated the footprint length. Regarding the distance between the detector and f max (x max ), LaStTraM shows very similar results compared to the FFP and K-M model with slightly shorter distances in the unstable cases and distances between the FFP and K-M models in the stable cases (see Table 3). The footprints generated by LaStTraM show a behavior that is well-known from other models: the footprint length increases with atmospheric stability while f max decreases. Increasing the detectors' height leads to an increase in the footprints' width and x max but decreases f max .
Looking at the differences in x max between the three models, it becomes clear that LaStTraM predicts similar results to FFP (see Table 3). This stands in contrast to f max , where LaStTraM shows more similar values to K-M.
Aside from the two-dimensional crosswind integrated visualization, LaStTraM saves three-dimensional output files that can be visualized using ENVI-met's visualization tool LEONARDO. Figure 7 shows two-dimensional representations of the three-dimensional outputs of the homogeneous model area.
The footprints generated by LaStTraM show a behavior that is well-known from other models: the footprint length increases with atmospheric stability while decreases. Increasing the detectors' height leads to an increase in the footprints' width and but decreases . Looking at the differences in between the three models, it becomes clear that LaStTraM predicts similar results to FFP (see Table 3). This stands in contrast to , where LaStTraM shows more similar values to K-M. Aside from the two-dimensional crosswind integrated visualization, LaStTraM saves three-dimensional output files that can be visualized using ENVI-met's visualization tool LEONARDO. Figure 7 shows two-dimensional representations of the three-dimensional outputs of the homogeneous model area. The two-dimensional distribution of the footprint shows that the maximum of the three-dimensional footprint is much closer to the detector location in the unstable case than in the stable case (see Figure 8). Furthermore, the along-wind distance is increased in the unstable case compared to the stable case. In general, the footprint of the unstable The two-dimensional distribution of the footprint shows that the maximum of the three-dimensional footprint is much closer to the detector location in the unstable case than in the stable case (see Figure 7). Furthermore, the along-wind distance is increased in the unstable case compared to the stable case. In general, the footprint of the unstable case is broader due to more turbulence and thus higher atmospheric mixing during daytime. The two-dimensional distribution of the footprint shows that the maximum of the three-dimensional footprint is much closer to the detector location in the unstable case than in the stable case (see Figure 8). Furthermore, the along-wind distance is increased in the unstable case compared to the stable case. In general, the footprint of the unstable case is broader due to more turbulence and thus higher atmospheric mixing during daytime.
The evaluation of the new footprint model in comparison with the established models FFP and K-M shows that LaStTraM is capable of reproducing their results. The calculated footprints resemble the results of both models and even improve the imprecisions mentioned by Heidbach et al., 2017 andVan de Boer et al. (2013) [37,38], i.e., LaStTraM most likely does not underestimate and overestimate the footprint length.

Evaluation of Simple Inhomogeneous Model Area Results
The results of the simple inhomogeneous case indicate that the tree contributes significantly differently to the footprint compared to the grass vegetation at the same distance from the detector (Figure 8). In the simulation where the detector is positioned close to the building, less particles reach the tree compared to grass cells at the same distance. This is caused by the flow induced by the building that leads to a significant updraft downwind of the tree (see Figure 9). In comparison to the grass at the same location, less particles are thus traced back to the tree from the detector close to the building.

Evaluation of Simple Inhomogeneous Model Area Results
The results of the simple inhomogeneous case indicate that the tree contributes significantly differently to the footprint compared to the grass vegetation at the same distance from the detector (Figure 8). In the simulation where the detector is positioned close to the building, less particles reach the tree compared to grass cells at the same distance. This is caused by the flow induced by the building that leads to a significant updraft downwind of the tree (see Figure 9). In comparison to the grass at the same location, less particles are thus traced back to the tree from the detector close to the building. In the simulation, where the detector is placed further downstream of the building the same effect-uplifting of flow and subsequent down-drafting due to the buildingleads to more particles reaching the tree compared to the grass vegetation below the tree. Furthermore, it can be seen that despite the simple inhomogeneous model area and the Lagrangian Stochastic approach, where the particles velocity is calculated including a stochastic, random part ( in Equation (7)), the vast number of simulated particles (100 million) ensures the high symmetry of the footprint (see Figure 8). In the simulation, where the detector is placed further downstream of the building the same effect-uplifting of flow and subsequent down-drafting due to the buildingleads to more particles reaching the tree compared to the grass vegetation below the tree. Furthermore, it can be seen that despite the simple inhomogeneous model area and the Lagrangian Stochastic approach, where the particles velocity is calculated including a stochastic, random part (P V LS in Equation (7)), the vast number of simulated particles (100 million) ensures the high symmetry of the footprint (see Figure 8).
Both footprints show greater horizontal extents than in the homogenous results. This is likely caused by the much higher turbulence induced by the complex structures in the model area. Nevertheless, with increasing distance, the number of particles from emitter cells decreases, just as in the homogeneous case.

Evaluation of Realistic Model Area Results
The results for the realistic model area confirm the general functionality of the model in complex model domains. Despite having no measurement data for this scenario, the distribution of particles shows plausible results: the vegetation closer to the detector is reached by more particles than those further away and the horizontal extent increases with increased distance from the detector (see Figure 10). Additionally, vegetation with greater vertical extent, i.e., larger trees, is more likely to be reached by the particles released from the detector. This is especially the case for the trees directly located upwind of the detector as the flow passing through them is lifted and directed upwards towards the detector's location.

Conclusions
Based on the three-dimensional model outputs of the microclimate model ENVI-met, the post-processing tool LaStTraM simulates probability distributions of particles using the Lagrangian Stochastic method. Combining LaStTraM with ENVI-met should allow us to simulate flux and concentration footprints in complex urban environments.
A footprint comparison for homogenous model areas between LaStTraM and two established models shows slight deviations for the maximum footprint contribution ( )

Conclusions
Based on the three-dimensional model outputs of the microclimate model ENVI-met, the post-processing tool LaStTraM simulates probability distributions of particles using the Lagrangian Stochastic method. Combining LaStTraM with ENVI-met should allow us to simulate flux and concentration footprints in complex urban environments.
A footprint comparison for homogenous model areas between LaStTraM and two established models shows slight deviations for the maximum footprint contribution ( f max ) and the footprint length (x max ). These deviations, however, might be a result of an underestimation of f max and overestimation x max by both the FFP and K-M model as demonstrated by earlier studies [37,38]. The results of the inhomogeneous model area and the complex model area shows reasonable results where the effects of obstructions alter the particle distribution significantly. While these scenarios provided promising results, a comparison with empirical data would be needed in order to verify the accuracy of the model in complex urban areas.
The great potential of LaStTraM lies in its ensemble with ENVI-met's averaged Navier Stokes model outputs, which should allow the study of footprints in heterogeneous and hilly environments. Our attempt to simulate a footprint in heterogeneous conditions demonstrated the capabilities of the model, i.e., incorporating three-dimensional emitters as well as the influence of obstacles on particle trajectories.
By incorporating ENVI-met's physiological aspects of vegetation, i.e., the cell based calculation of stomatal conductance and thus transpiration, the simulation of the footprint using LaStTraM can further be improved. Footprints of different compounds can also be studied separately, depending on the variable of interest. For example, the dependency of various compounds on the flux footprint can be studied by defining individual areas/volumes that act as sources of these compounds. All results can then be visualized with the built-in ENVI -met program LEONARDO.
It was demonstrated that 100 million particles within a 2 × 2 km 2 model area can be modelled within a few hours on a PC. The order of magnitude of particle numbers, the model area size as well as the grid resolution are comparable to numbers reported in LES footprint studies before [16], with the essential difference that the footprint model is run in backwards mode, which means that every modelled particle (that hits a source in the model domain) influences the footprint. In contrast to a forward in time trajectory mode, the backwards mode has the advantage that the detector does not have to be unrealistically large in order to gather a sufficiently large number of particles, as demonstrated before [16].
One drawback to eddy covariance users is that the ENVI-met model cannot be run routinely for each 30 min flux average, which it has in common with all RANS and LES models. The possibility of simulating whole years with the microclimate model ENVI-met could be to split the simulations into partially overlapping timeframes, e.g., 12 simulations, one for every month of the year, including a couple of days of overlap. These could be run on 12 PCs in parallel and the overlapping results could be stitched together in the end to create a year-long simulation [39].
Regardless of the timeframe, LaStTrAM is capable of reproducing typical atmospheric conditions (unstable, stable) that will help to elucidate the origin of measured fluxes in heterogeneous environments. The LaStTraM footprint model is therefore uniquely positioned for eddy covariance studies in complex urban terrains, where the premises of analytic models are generally not met. It also can be used to analyze flux footprint problems that frequently occur in rural locations when eddy covariance towers are positioned close to forest edges or when other steep changes in roughness occur.
The present work is far from exploiting the whole potential of LaStTraM, as larger and especially more complex areas could be modelled in future studies. Furthermore, it has the potential to model the footprint of reactive gases instead of only inert particles as shown in this study. Vesala et al. (2008) highlight the potential of ensemble-averaged closure models especially for sites with complex topography and canopy heterogeneity. This includes urban sites, which are of increasing interest in surface-atmosphere exchange studies. By combining LaStTraM with ENVI-met, model areas featuring buildings, sources, vegetation, and uneven terrain can be simulated. This makes it a viable tool to calculate and analyze three-dimensional footprints in complex environments.