Numerical Framework for the Computation of Urban Flux Footprints Employing Large-eddy Simulation and Lagrangian Stochastic Modeling

Conventional footprint models cannot account for the heterogeneity of the urban landscape imposing a pronounced uncertainty on the spatial interpretation of eddy-covariance (EC) flux measurements in urban studies. This work introduces a computational methodology that enables the generation of detailed footprints in arbitrarily complex urban flux measurements sites. The methodology is based on conducting high-resolution large-eddy simulation (LES) and Lagrangian stochastic (LS) particle analysis on a model that features a detailed topographic description of a real urban environment. The approach utilizes 5 an arbitrarily sized target volume set around the sensor in the LES domain, to collect a dataset of LS particles which are seeded from the potential source-area of the measurement and captured at the sensor site. The urban footprint is generated from this dataset through a piecewise post-processing procedure, which divides the footprint evaluation into multiple independent processes that each yield an intermediate result that are ultimately selectively combined to produce the final footprint. The strategy reduces the computational cost of the LES-LS simulation and incorporates techniques to account for the complications 10 that arise when the EC sensor is mounted on a building instead of a conventional flux tower. The presented computational framework also introduces a result assessment strategy which utilizes the obtained urban footprint together with a detailed land cover type dataset to estimate the potential error that may arise if analytically derived footprint models were employed instead. The methodology is demonstrated with a case study that concentrates on generating the footprint for a building-mounted EC measurement station in downtown Helsinki, Finland, under neutrally stratified atmospheric boundary layer. 15


Introduction
Micrometeorological measurements in densely built city environments pose an antipodal problem: they are essential in establishing the fundamental basis for the study of urban microclimates, but these measurements are endowed with pronounced uncertainties, which mainly originate from the topographic and elemental complexity of the urban landscape.The resulting noncompliance between the theory and practice in urban micrometeorological measurements undermines the study of how our cities interact with the surrounding atmosphere.At the very heart of this discord lies the problem concerning the determination of effective source areas, or footprints, of urban flux or concentration measurements.
The footprint is a concept used to describe the surface area that contains the sources and sinks which contribute to the measured quantity obtained by a sensor (Pasquill, 1972).In another words, it is such a sensor's "field of view" whose identification is essential in interpreting the obtained flux or concentration values in their correct spatial extent (Schmid, 2002).Mathematically, the footprint is a transfer function f , which relates the value of a measurement (of flux or concentration) η at location x M = (x M , y M , z M ) to the spatial distribution of Q from a volumetric domain of interest: Published by Copernicus Publications on behalf of the European Geosciences Union.
where f has dimensions of inverse of integration units (m −3 ).In the subsequent presentation the vertical dimension of domain is collapsed and thereby f has dimensions of inverse area (m −2 ).The footprint can also be interpreted as a spatial weighting function that expresses the probability with which a fluid element that coincides with an element of Q contributes to the measurement at x M (Pasquill and Smith, 1983).In accordance with Sogachev et al. (2005), this study does not adhere to the strict interpretation where the footprint is only a function of turbulent diffusion and sourcesensor location, but allows the possibility that, for instance, variations in source-area topography can influence the result.
In this context, topography refers to an elevation model of the land and buildings together.Consequently, the footprint should provide the critical link between the point measurement and the geographical distribution of sources, yielding a complete characterization of η with regard to its contents.
In an effort to achieve this, analytical closed-form solutions have been derived for the footprint functions -see Schmid (2002) for a comprehensive review -but only under the assumptions that (1) steady-state conditions prevail during the analyzed period, (2) turbulent fluctuations in the atmospheric boundary layer (ABL) are horizontally homogeneous, and (3) there is no vertical advection.These assumptions allow the governing equations to be reduced to a time-averaged balance between advection and turbulent diffusion which admits, with appropriate parametrization of the turbulent flow field, a closed-form expression for the footprint function.
The underlying assumptions are often acceptable in measurement sites where the sensors are mounted on towers that have been appropriately placed above homogeneous forested landscapes and well above the surface roughness sublayer height where the effects of the individual roughness elements disappear.However, due to practical regulations constraining measurement campaigns in densely populated cities, sufficiently tall flux towers cannot be erected above the skyline of central urban areas.It is often inevitable that if the urban microclimate is to be studied experimentally, the measurements must be obtained near the border of the roughness sublayer by sensors that are mounted either on low-rise towers or on top of tall buildings.In these suboptimal conditions, assumption (2) becomes strictly invalid and assumption (3) highly questionable because urban boundary layer (UBL) flows are typically characterized by developing and strongly heterogeneous flow conditions, particularly at lower elevations where individual buildings influence the turbulence.
Considering that the analytical footprint models effectively provide ellipse-shaped probability distributions for the source contributions without any regard to topographic heterogeneities, it becomes clear that the use of such source-area models becomes highly suspect in real urban conditions.This is an unacceptable state of affairs in the urban micrometeorology research and immediately calls for targeted efforts to alleviate the uncertainties associated with the invaluable urban flux-measurement data.Although, the first efforts by Vesala et al. (2008), utilizing the method by Sogachev et al. (2002), already explored topography-sensitive urban footprints, the applicability of the documented approach has not reached the scale and accuracy requirement of the urban footprint problems considered herein.
As a response, this works introduces a new numerical methodology to construct detailed topography-sensitive footprints for complex urban flux measurement sites by the means of pre-and postprocessing developments and a largeeddy simulation (LES) solver suite that features an embedded Lagrangian stochastic (LS) particle model.This coupled model will be referred to with the acronym LES-LS.The proposed methodology is designed to be first and foremost a postprocessing procedure, which exploits the current stateof-the-art LES-LS modeling framework in an urban setting with a minimal investment in the initial setup.
The principal objective is to provide a reliable computational framework, founded on a high-resolution LES-LS analysis, to generate the most accurate footprint estimates feasible without the need to conduct tracer gas experiments, which are nearly impossible to arrange in residential areas.These computationally generated footprints open up the possibility to study the appropriate placement of new measurement stations and to assess the magnitude of the potential misinterpretation which may arise from the application of closed-form footprint models to urban flux or concentration measurements.The proposed framework is also supplemented by a convenient technique to approximate this error with the assistance of a land cover classification dataset.
The methodology is demonstrated with a numerical case study, which is staged in Helsinki, the coastal capital city of Finland, and focuses on the eddy-covariance (EC) measurement site mounted on the roof of Hotel Torni (Nordbo et al., 2015;Kurppa et al., 2015), which is the tallest accessible building in the downtown region.The building height is 57.7 m and the EC sensor is situated 2.3 m above it corresponding to 74 m height above the sea level.Thus, the effective measurement height (a.g.l) is z M = 60 m -d = 45.1 m, where d = 14.9 m is the displacement height of the site according to Nordbo et al. (2013).The mean building height of the surrounding area is 24 m.The site belongs to SMEAR III (Station for Measuring Ecosystem-Atmosphere Relations, Järvi et al., 2009) and is also part of the urban network of atmospheric measurement sites (Wood et al., 2013).Its potential source area closely resembles a typical European city arrangement that features perimeter blocks with inner courtyards.
This study employs the PArallelised LES Model PALM (Maronga et al., 2015;Raasch and Schröter, 2001), which has been previously applied to footprint studies by Steinfeld et al. (2008) and very recently by Hellsten et al. (2015), who constructed footprints for an idealized city environment as a precursor study to this work.The presented contribution places special emphasis on the issue of composing footprints for flux measurement sites that are surrounded by arbitrarily heterogeneous topography and may be compromised by the fact that they are mounted on top of actual buildings instead of conventional radio-mast-like towers.Such a complex urban setting requires a new mechanism for constructing footprints, which is accompanied by a requirement that the associated LES-LS simulation is capable of resolving the relevant turbulent structures ranging from the street-canyonscale phenomena within the roughness sublayer to the larger ABL structures, while also accounting for the interaction between them (Anderson, 2016).

Numerical modeling framework
The PALM model utilized in this study is an open-source numerical solver for atmospheric and oceanic flow simulations.The software has been carefully designed to run efficiently on massively parallel supercomputer architectures and it is therefore exceptionally well suited for high-resolution UBL simulations considered herein.The LES model employs finite-difference discretization on staggered Cartesian grid and utilizes an explicit Runge-Kutta time-stepping scheme to solve the evolution of velocity vector u = (u, v, w), modified perturbation pressure π * , potential temperature θ , and specific humidity q v fields from the conservation equations for momentum, mass, energy, and moisture, respectively.The conservation equations are implemented in an incompressible, Boussinesq-approximated, non-hydrostatic, and spatially filtered form, which indicates that the conservation of mass is imposed by the solution to a Poisson equation for π * .The filtering refers to the separation of scales in LES where the turbulent scales containing the majority of energy are resolved by the grid while the diffusive effect of the unresolved subgrid-scale (SGS) turbulence is accounted for by a SGS turbulence model.To achieve closure in the final system of equations, PALM implements the 1.5-order SGS turbulence model by Deardorff (1980), modified according to Moeng and Wyngaard (1988) and Saiki et al. (2000).The model involves an additional prognostic equation for SGS turbulent kinetic energy (SGS-TKE) e.
The embedded Lagrangian particle model in PALM implements the time-accurate evolution of discrete particles (either with or without mass) through a technique that conforms to the LES approach: the trajectories are integrated in time such that the transporting velocity field is decomposed into deterministic (i.e., resolved) and stochastic (i.e., subgrid-scale) contributions.The deterministic velocity components are directly obtained from the LES solution, while the random components are evaluated according to Weil et al. (2004).Although LS modeling approaches that are less computationally expensive exist (Glazunov et al., 2016), warranting further investigation on their applicability to urban problems, the presented high-resolution urban flow problem is assumed to require the highest level of description also from the LS model; the interaction between the atmospheric wind and the cascade of multistoried buildings and street canyons gives rise to strongly anisotropic turbulence structures, which are not reliably amendable to parametrization.
While the LES-LS simulations are carried out in large supercomputing facilities, the preprocessing of the urban topography model and the postprocessing of the final footprint from raw data is performed on a personal workstation utilizing freely available numerical scripting and data visualization technologies.See the paragraph on code availability at the end of this paper.

Urban topography model
The urban topography model, used in describing the bottom wall boundary of the LES domain, is prepared from a detailed 2 m resolution laser-scanned dataset of the Helsinki area (Nordbo et al., 2015).The data are conveniently available in raster map format and, in addition to the height distribution h(x, y), also include a distribution of land cover types LC(x, y) ∈ {1, 2, . .., N LC } where N LC is the number of land cover classes in the dataset.Both raster maps are shown in Fig. 1.Access to similar surface data source is a critical prerequisite for the presented methodology.
The horizontal domain for the LES analysis extends L x = 4096 m in the mean wind direction and L y = 2048 m in the crosswind direction and is spatially oriented such that x axis is coincident with the geostrophic wind direction of the case study.The EC measurement site at Hotel Torni is pivotally located in the LES domain to facilitate the determination of its footprint.However, the extracted raster map has to be first purposefully preprocessed to attain a form that complies with the LES analysis-specific requirements.The following manipulations were applied to obtain the final topography model depicted in Fig. 2.
1.The first half of the topography model (where x < L x /2) is flattened for the purpose of generating physically realistic ABL conditions at the inlet through turbulence recycling technique (see below).
2. The lateral sides were made identical for cyclic boundary condition treatment by applying a zero-height margin that smoothly blends toward the values in the interior.
3. Immediately upstream of the outlet boundary, a margin with sloping terrain height is applied to force the highly www.geosci-model-dev.net/10/4187/2017/Geosci.Model Dev., 10, 4187-4205, 2017 turbulent flow (caused by the buildings near the end of the domain) to slightly accelerate before reaching the outlet boundary where reversed flow causes numerical difficulties.

Physical setup for the LES model
The meteorological conditions for the simulation are adopted from 9 September in 2012 when near-neutral ABL conditions were recorded with the EC measurements made on top of the Torni building.Lidar measurements (Wood et al., 2013) from the chosen time frame yielded |u g | = 10 m s −1 for the geostrophic wind in a southwesterly direction (α = 218 • ) and δ ≈ 300 m for the boundary layer height.The Coriolis force (corresponding to latitude 60 • N) is included to account for the turning of the flow within the boundary layer.
The meteorological conditions are conveyed to the simulation by means of a precomputed ABL solution over flat surface, which in this context represents the surface of the Baltic Sea bordering Helsinki from the south.The boundary conditions for the velocity solution in this precursor simulation were set such that a fixed value is applied at the top and a noslip condition at the bottom boundary of the domain while setting all the lateral boundaries as periodic.
For the precursor simulation the solver was run with an option that explicitly conserves the initial mass flow rate across the system, which was specified by initializing the velocity field with a constant value u t=0 = 0.95 u g .This initializa-tion value was determined by trial and error with the objective that the precursor solution would ultimately yield the desired geostrophic wind value at z = δ for the horizontally averaged velocity field ū pre .The boundary layer growth was controlled by initializing the potential temperature field with a vertical profile θ 0 (z) that features a strong inversion layer at 300 < z < 350 m.This θ 0 (z) profile is defined according to the following lapse rates: The precursor LES solution was computed on a grid that has the same resolution and vertical dimension as the principal urban LES grid, but its lateral dimensions are smaller by an integer division.Table 1 summarizes the respective grid characteristics.The study features a spatial resolution of 1 m, which is unprecedented at this scale.Giometto et al. (2016) found the same resolution to be sufficient to capture the relevant turbulence physics within a real urban roughness sublayer.However, the effect of grid resolution on the final result is not investigated in this work.The influence of the structural details of the urban surface (balconies, chimneys, ventilation ducts, stationary cars, small-scale vegetation, etc.) not included in the urban topography model are taken into account by specifying a uniform roughness length z 0 = 0.05 m on the bottom boundary surfaces (Letzel et al., 2012).The precursor simulation generates a highly resolved ABL solution that will be utilized, first, in a recursive manner to initialize the entire urban LES flow field with turbulence and, second, to aid construction of appropriate inlet boundary conditions though a technique labeled turbulence recycling, which is based on the method by Lund et al. (1998) with modifications by Kataoka and Mizuno (2002).The implementation of this boundary condition in PALM is presented in Maronga et al. (2015), but to aid discussion the description is also covered here with modified notation.
Denoting prognostic field variables by ψ = ψ(x, t) where ψ ∈ {u, v, w, θ, e}, the precursor solution is used to extract horizontally averaged vertical profiles ψ pre (z) for the turbulence recycling boundary condition.These stationary profiles are utilized at the inlet boundary in the urban simulation to conserve the global state of the mean flow, but in a manner that also incorporates physically sound turbulent fluctuations that occur in an ABL flow.This is achieved by specifying a recycling plane, that is, a yz plane at a windwise coordinate x rc , placed sufficiently far downstream from the inlet to prevent feedback of disturbances between the two planes.The fluctuations are obtained from the recycling plane through the following technique: where the spatial mean (in the crosswind direction) ψ y = ψ y (z, t) at the recycling plane is computed as a time dependent vertical profile that carries a dependence on N y .Finally, utilizing the precursor generated mean profiles, the turbulence recycling inlet boundary condition becomes In association with the turbulence recycling, the top boundary condition in the main simulation is specified as a slipwall.
In this study, the recycling plane is situated, as shown in Fig. 2, in accordance with the precursor domain length such that (x rc − x in ) = 1024 m ≈ 3.4 δ and the same distance is allocated from the recycling plane to the edge of the urban topography to ensure that disturbances originating from the urban terrain are not conveyed back to the inlet.The chosen turbulent inlet arrangement generated no observable feedback effect on the incoming turbulence field.

LS particle model setup for the footprint evaluation
The embedded LS particle model is employed such that, after the initial transients in the LES solution have subdued (after approximately 5 min of simulation), the release of particles is activated within the region outlined in Fig. 2. The release area extends 3030 m (≈ 41z M ) in the upwind direction and 780 m (≈ 10.5z M ) in both lateral directions from the Hotel Torni's EC site.The release area has been trimmed according to preliminary trial simulations to reduce the number of redundant particles in the domain.Denoting the Lagrangian coordinate vector of the lth particle by X l (t) = X l (t), Y l (t), Z l (t) , the release locations X l o = X l t=0 are uniformly distributed 2 m apart in the x and y directions while the vertical coordinate is set Z o = 1 m above the topography: The release height of one grid spacing at 1 m resolution is inferred to be a justifiably close to the surface to represent both the traffic emissions as well as the surface atmosphere exchanges.It also lowers the risk of accumulating a large number of particles within the first grid cell where the velocity values are dictated by the logarithmic wall function and the vertical advection of particles solely by the stochastic model due to Weil et al. (2004).Thus, the underlying assumption is that, at 1 m resolution, the release height of 1 m above solid surfaces does not significantly influence the footprint distribution, which is evaluated at 2 m horizontal resolution.
The raw particle data for constructing footprints through LES-LS modeling in an arbitrarily heterogeneous environment are obtained by setting a target volume around the specified sensor location x M and recording which particles hit this target.Although this approach appears natural and straightforward at first sight, a closer scrutiny reveals a number of problematic issues which arise with this setup, particularly when the flux sensor is mounted on a building (or close to one) instead of a tower.Purely from the perspective of particle data acquisition in the LES-LS simulation, setting a larger target volume would directly alleviate the computational effort required to gather a large enough dataset of particle hits, but this would clearly violate the formal premise that the footprint should be evaluated for the coordinate x M of the sensor.However, it turns out that the discrete setting of the LES-LS approach calls into question the relevance of seeking an urban footprint for a precise point near the surface of a solid structure.
Consider the problem of strictly concentrating on the exact location x M of the sensor.This effort becomes immediately futile as the spatial resolution with which the buildings are described in the topography model (which contains information on elevation changes only) cannot account for structural details that, in reality, influence the flow conditions at the precise location of the sensor.The same reasoning also extends to the LES flow analysis where the computational cost would become prohibitively expensive if the resolution would be set according to the ∼ 10 −1 m scale of structural detail of building facades and rooftops in the hypothetical situation that such datasets were available.Therefore, it is important that the methodology for evaluating footprints in urban environments comes with a prerequisite that the resolution demands of the LES-LS model are purely dictated by the turbulent structures within the urban canopy and not the fine details of the sensor site.On these grounds, the method to collect particle data in the LES-LS simulation is based on setting a finite target volume around the sensor location x M without strictly dictating the appropriate size.This is done understanding the fact that the flow around the sensor mounting building strongly interacts with the flow, resulting in strong gradients in the mean velocity field in the vicinity of the sensor.This is bound to further complicate the subsequent postprocessing of the flux footprint because the eddycovariance approach necessitates that the effect of the mean flow should be eliminated through the process of coordinate rotation (Aubinet et al., 2012), which is presented in the context of this study in Sect.2.3.1.Clearly, the discrete LES-LS approach in an arbitrarily complex urban environment is endowed with pronounced uncertainties.For this reason, the postprocessing procedure has to encompass a capability to conduct spatial sensitivity analysis on the intermediate footprint results and, according to its outcome, selectively exploit the particle dataset in the final processing of the result.
Adopting this strategy reduces the level of rigor required at the setup stage of the LES-LS analysis and simplifies the guidelines for the particle acquisition: the target volume should be centered at x M and its dimensions chosen to represent the sensor site proportionately (vagueness intended) to the dimensions of the building geometry.In all cases, it is important to acknowledge that, as a rule of thumb, more than 10 7 particle hits need to be recorded at the target volume during the course of the LES-LS simulation to gather a large enough dataset for flexible postprocessing.In general, it pays off to specify an oversized target and gather a large dataset, accepting that it contains certain percentage of particle hits whose contribution will be discarded.In this study, using L T ≈ 4 m to represent the horizontal length scale of Hotel Torni's apex structure on which the sensor is mounted, the target for monitoring particle hits is specified as a box of volume The box is centered at the apex (which closely coincides with the actual sensor location x M ) such that it extends 2.5L T in the crosswind and upward directions, 1L T in both streamwise directions and 0.5L T downward, entering partly into the building structure.And, to reiterate, these dimensions were chosen under the guiding principle that the target box reasonably represents the sensor site and enables particle hits to be gathered at higher rate.Figure 3 provides an illustration of the size and placement of the target box in relation to the surrounding urban topography.The monitoring is performed at 0.5 s intervals, which corresponds to approximately eight LES time steps.This allows the same particle to be recorded multiple times at different locations within the target box.This feature is intentional and desirable because of the chosen postprocessing strategy.

LES-LS analysis
The precursor simulation is run for 1.5 h physical time to develop the desired ABL profile.The initialization of the primary LES-LS computation with this precursor solution expectedly results in short-lived unphysical fluctuations around the urban topography, but after 3 min of simulation these overshoots have been advected away from the domain.The release of LS particles is initiated after 5 min of simulation, and from there on particles are released simultaneously in puffs at 10 s intervals such that two particles are seeded from each location at every instance.This translates into releasing approximately 2.36×10 6 particles every interval.The release schedule was determined by trial and error to best utilize the computational capacity of the supercomputer.Each particle is assigned a maximum lifetime T l max = 1200 s, which is long enough to guarantee that even the particles that are advected by the slowest ∼ 0.2 u g velocity scales manage to travel over 2 km during this time frame.The total number of particles in the whole domain converged to approximately 68×10 6 .Particles reaching any of the lateral boundaries or the top boundary are "absorbed", that is, deleted and deallocated from the computer's memory while the wall boundary below functions as an ideally smooth reflective surface for the particles.The simulation was run for 3 h physical time during which ca. 19×10 6 particle hits were recorded at the target volume.The computation cost of this simulation is comparable to running 3-4 urban flow simulations with the objective of studying turbulence.In absolute terms, the simulation took ca. 10 days on the Cray XC40 supercomputer "Sisu" (CSC -IT Center for Science, Finland) with 2048 CPUs which amounts to ca. 5.3×10 5 CPU hours.The LS model constituted merely 20 % of the total CPU time of the LES-LS simulation, which is an appreciably moderate value considering the high number of particles handled by the solver.

Piecewise postprocessing methodology for constructing the footprint
During the LES-LS simulation, the sampling of particle hits at the target volume V T entailed recording each particle's (identified as l) coordinate of origin X l o , incident velocity T at the target, and the associated sample location X l T (indicating where the particle hit the target), ultimately giving rise to a large dataset where N rp refers to the total number of released particles.
According to the issues discussed in Sect.2.2.3, the postprocessing of S is now required to account for the spatial uncertainty and facilitate a sensitivity study on the obtained result.This is achieved by introducing a piecewise processing strategy where the principle idea is that the original dataset S is split into smaller subsets according to a Cartesian discretization of the target volume V T .See an example illustration in Fig. 4. Thus, the target is divided into subvolumes V i,j,k , satisfying V T = n x k n y j n z i V i,j,k where i, j , and k are the Cartesian indices of the subvolumes.The number of divisions in each coordinate direction n x , n y , and n z have to be determined case by case as the optimal values depend on the target volume size, the total number of particle entries in the dataset, and the complexity of flow solution in the vicinity of x M .
Each target subvolume now yields an associated subset s i,j,k ⊂ S containing a record of the particles that hit the corresponding subvolume V i,j,k centered at x V i,j,k = x M + dx i,j,k , where dx i,j,k is the displacement from the exact measurement location x M to the center of the subvolume V i,j,k .The obtained subsets can be independently postprocessed to generate sectional flux footprints f i,j,k , utilizing an estimator similar to Kurbanmuradov et al. (1999) (see also Rannik et al., 2000), but modified to approximate the footprint by computing the probability with which a fluid parcel released from a continuous source at x f = (x, y, h + Z o ) will lie within V i,j,k at any given time.Discretizing the source area (i.e., footprint grid) by x f = ( x f , y f , 0 ) with x f = y f = 2 m, the estimator reads which has an implicit dependence on the vicinity of x V i,j,k through the spatial confinement of s i,j,k .In Eq. ( 6) N i,j,k denotes the number of particle entries within the subset s i,j,k collected over a sufficiently long time period, and is the vertical velocity deviation of particle l from the spatially averaged mean flow value evaluated over the subvolume V i,j,k .Equation ( 7) relates to the coordinate rotation of the EC sensor, which eliminates the effect of w from the vertical flux evaluation by aligning the sensor with the mean wind (Aubinet et al., 2012, p. 76).Here, the evaluation of W l T proves particularly problematic due to the approximations associated with the use of w i,j,k and, therefore, it is a subject of further discussion in Sect.2.3.1.Finally, the function I = I (X l o , x f , x f ), which is responsible for distributing the hits on to the footprint grid based on the particles' coordinate of origin, is given as follows: www.geosci-model-dev.net/10/4187/2017/Geosci.Model Dev., 10, 4187-4205, 2017  The evaluation procedure (6) closely resembles that of Rannik et al. (2003), with the exception that here it is assumed that each particle is represented only once in each subset s i,j,k .
The individual sectional footprints are typically evaluated from subsets that contain an insufficient number of particle data entries needed to obtain a converged footprint distribution.(Hellsten et al., 2015) showed that, in an urbanlike environment, ∼ 10 6 particle hits are required to attain an adequately converged footprint distribution while ∼ 10 5 particle entries is sufficient to reveal the characteristic shape of the near-field distribution.In the piecewise postprocessing approach, the sectional footprint contributions may be constructed from an arbitrarily small dataset, but since the methodology substantially benefits from the ability to inspect and compare individual f i,j,k distributions, it is desirable to work with subsets s i,j,k containing more than 10 5 entries.To facilitate the postprocessing procedure, each f i,j,k should be individually stored as a stand-alone two-dimensional scalar field (i.e., raster map) that can be projected onto the threedimensional topography model of the LES domain to permit descriptive visualizations in the urban setting.The value of the denominator D i,j,k = N i,j,k x f y f featured in Eq. ( 6) has to be stored together with the footprint distribution because the assembly of the final footprint is carried out by computing where K is the set of all i, j, k combinations which have been selected via spatial sensitivity analysis (covered in Sect.2.3.2).

Coordinate rotation via far-field correction
The piecewise processing of the footprint carries an inherent difficulty that arises in situations where the mean flow displays strong gradients within the target volume.This is evidently present in the considered case study featuring an EC sensor mounted close to the top of a building.The difficulty relates to the evaluation of w i,j,k which is used in the footprint evaluation to extract the fluctuating velocity components about the mean value within its corresponding subvolume (as shown above).The initially implemented piecewise processing approach naturally involved utilizing LES to obtain the (45 min time-averaged) mean velocity distribution w from within the target box volume and evaluating the spatial average w i,j,k for each subvolume V i,j,k .However, with the Hotel Torni case study it became evident that this approach gave rise to a systematic negative bias in the footprints, which becomes immediately apparent in the far-field distributions.This outcome persisted until the discretization of the target volume was refined according to n x = 8, n y = 20, and n z = 12 such that the subvolumes corresponded with the uniform 1 m resolution of the LES grid.This involved generating n x × n y × n z = 1920 independent f i,j,k contributions via Eq.( 6), where the w i,j,k values were now directly obtained from their corresponding LES grid cells.Figure 5 illustrates the effect of target volume discretization by comparing crosswind integrated footprints f y at different levels of refinement.The comparison reveals how the negative bias in the far field and the reduction in near-field magnitude immediately emerge with coarser discretizations.It should also be noted that a targeted refinement in the z direction, while using coarser horizontal resolution in effort to generate thin subvolumes that approximate planes, does not remedy the situation because the mean flow gradients around the sensor site are significant in all directions.Such finite planes or one-cell-high grid layers are conventionally used as targets when Lagrangian-particle-based methods are utilized to evaluate footprints under heterogeneous conditions with undisturbed sensor sites (e.g., Steinfeld et al., 2008;Hellsten et al., 2015;Glazunov et al., 2016).Unfortunately, at the required level of target volume discretization, the excessive number of individual f i,j,k contributions causes the postprocessing to become highly tedious.Since the proposed LES-LS methodology is founded on the premise that the size of the original target box around the sensor site can be chosen arbitrarily, the postprocessing effort must entail a procedure that enables the exclusion of those f i,j,k contributions that are deemed unfit for the final assembly.However, this selection operation becomes overly laborsome to manage when the number of subvolumes becomes large (viz.values exceeding 10 2 ) and particularly when the individual sectional footprints inadequately converge and thereby become uninformative when examined independently.For instance, in this case study, the required level of discretization gives rise to f i,j,k contributions that are generated from ca. 10 4 particle entries, which is a decidedly insufficient amount even for generating informative approximations for the near-field distributions.For these reasons, it is deemed unacceptable that the evaluation of urban footprints solely relies on the established piecewise postprocessing method.In response, this paper introduces an augmented coordinate rotation technique, labeled far-field correction, which incorporates well into the proposed piecewise postprocessing strategy and brings significant savings in the associated data manipulation efforts.This alternative technique allows much coarser target volume discretization to be employed in the assembly of the final footprint without unacceptably compromising the result.
The method has a prerequisite that the deficiently obtained footprint (for instance, obtained via insufficient target box discretization) must exhibit a properly leveled off far field, because the approach fundamentally relies on the following simple assertion: if the footprint distribution plateaus in the far field, this asymptote can be amended to become the zero reference level, which deviates from the "correct" asymptote by a negligibly small offset.Accepting this assertion and the associated approximation paves the way for a corrective coordinate rotation scheme which can be laid out by first classifying the data contributing to the far-field footprint via sub-sets r i,j,k ⊂ s i,j,k which are defined as the sets of particle entries whose X o fall into the outermost portion of the domain Here X min o = min l∈{1,...,N rp } (X l o ) is the farthest upstream coordinate where particles are seeded (thus, farthest away from x M ) and β specifies the remotest percentage of the footprint across which the mean value of f y no longer changes, that is, ∂ f y ∂x ≈ 0 when averaging over the length of the far field.(The β value is case-specific, but a typical range is expected to fall between 10 and 20.)With the help of the far-field datasets r i,j,k , the fluctuating vertical velocity component, used in Eq. ( 6) and previously defined by Eq. ( 7), can now be evaluated as where defines the far-field-corrected mean vertical velocity, which is obtained by scaling the initially obtained value by a coefficient c i,j,k to satisfy the criterion that the particle entries in each r i,j,k do not contribute to the corresponding f i,j,k .This becomes a simple one-dimensional optimization problem in which the objective is to minimize J = β f i,j,k dx , where β represents the far-field domain, by the means of controlling c i,j,k .Thus, this technique bears resemblance to a planar fit method (Wilczak et al., 2001).Because the control variable here is a single scalar, a rudimentary implementation of an iterative gradient decent search algorithm suffices (see, for instance, Nocedal and Wright, 2006).
Table 2 displays selected diagnostic data obtained from an application of this far-field correction technique to the Hotel Torni footprint case study.The data indicate that, when the mean vertical velocity values are initially obtained from the LES solution, the c i,j,k scaling coefficients concentrate near the mean value of 0.9.The range of individual values naturally depends on the magnitude of the starting value w i,j,k in Eq. ( 11) which, in turn, depends on the chosen discretization of the target volume.But it is important to emphasize that, although the far-field correction method is guaranteed to yield a physically justifiable asymptotic behavior for the footprint, the combined effect of the correction method and the target volume discretization on the final footprint result cannot be inferred from Table 2.
This realignment of the coordinate rotation plane within a larger subvolume, when examined in contrast with the reference technique where the coordinate rotation is performed at full LES resolution, alters how some of the individual particles contribute to the footprint.However, this discrepancy gives rise to an error that is distributed throughout the footprint domain.Therefore, the validity of the far-field correction approach hinges upon the magnitude of this distributed error and its sensitivity to the target box discretization.The sensitivity can be established by carrying out the selective assembly of the footprint result for different levels of target box discretizations.

Selective assembly of the final footprint
Since its conception it has been clear that the piecewise postprocessing approach must be endowed with the capacity to incorporate a sensitivity analysis phase into the final assembly of the footprint result.One of the driving motivators for developing the piecewise approach arose from the need to reduce the computational cost of collecting a large number of particle hits by an arbitrarily sized target volume around x M .However, the reduction can only be achieved by the piecewise postprocessing approach if the sectional footprint results are allowed to be inspected and combined in a partially converged state.This is an important stipulation without which the proposed postprocessing strategy fails to offer considerable computational savings.
Thus, the process of selectively assembling the final footprint result begins by first defining an inadequately converged initial footprint, which represents the desired preform at x M .This reference footprint, labeled f REF , should be constructed from at least 10 6 particle entries to facilitate a sufficiently informative evaluation of sensitivities.The selection process proceeds by iteratively introducing partial contributions f (l) that are independent from f REF and evaluating the sensitivity of the footprint distribution with respect to the selection of target box indices in K (see Eq. 9).The objective is to obtain a sufficiently converged footprint while minimizing the discrepancy between the constituent f i,j,k included in the final result.Thus, the selection process is quantitatively guided by the evaluation of "deltas" between f REF and f (l) , constructed from a partial set K (l) of target box indices (i.e., f (l)  = f REF −f (l) ) and utilizing a norm over a subdomain ⊂ encompassing only the near field (a fraction of the total LES footprint domain nearest to x M ) as a measure for the associated discrepancy.In this connection, it has been found most effective to define the extent of such that the integral over the near-field domain constitutes approximately half of the total integral of the footprint: utilizing identically normalized footprints for this evaluation.In this study the footprints are normalized to yield f dx = 1.The exclusion of the outer portion of the footprint domain allows the relevant deviations in the near field to be reflected in || f (l)  || 2, while avoiding the contamination due to poorly defined "deltas" in the weakly converged outer region.In this study, the nearest 30 % of the total length of the LES footprint domain is used to represent the near field as this yields for the normalized reference footprint f REF dx = 0.51.The search for the fitting contributions entails an iterative procedure, which is described || 2, values are depicted in Fig. 6.The process begins by setting at the 0th iteration K (0) = {I M , J M , K M } ⊂ K, where the indices correspond to the subvolume containing x M .The obtained footprint f (0)  = f I M ,J M ,K M is composed of ca. 4 × 10 5 particle entries, which does not meet the desired level of convergence to act as f REF .Thus, through a qualitative inspection, the original set is augmented 1) , which is chosen as the reference footprint.
The iterative process continues such that new candidate contributions f (l) are introduced incrementally in a radially outward progressing manner.This process is demonstrated in Fig. 6 where intermediate entries f (2) -f (6) introduce differently combined additions in y, x and z directions.For the sake of brevity, the example contributions combine a relatively large number of f i,j,k entries.The decision to include a candidate contribution in the final assembly is done according to a criteria || f (l)  || 2, ≤ || f || max , where the maximum allowable discrepancy || f || max must be determined according to the case-specific requirements.In this case study, the threshold was set to include f (3) || 2, .Naturally this threshold level can be varied to generate alternative footprint assemblies (with different levels of convergence), which allow, in the context of the considered footprint applications, the impact and uncertainty associated with these choices to be transparently monitored.
The obtained final result, which combines the earlier accepted additions, features 20/45 of all subvolume contributions.The obtained footprint also exhibits adequate convergence in the far field, having been constructed from ca. 8 × 10 6 particle entries.Subsequently, the lowest vertical (k = 1) plane and the farthest (i = 3) plane were completely excluded from K in the final assembly.This outcome indicates that the contributions with the largest deviations arise from V i,j,k that are either in contact with the tower structure or in its wake region.Therefore, this suggests that it is not advantageous to set up V T such that the building structure cuts into the volume.
As long as the individual subsets contain a sufficient number of particle data entries (> 10 5 ), as is required by the far-field correction approach, it is beneficial to discretize the target volume as finely as possible (by increasing n x , n y , and n z ) as it enables a more flexible and fine-tuned assembly and permits a more accurate coordinate rotation treatment.Depending on the total number of particles gathered during the simulation and the size of the target box, the maximum number of admissible subvolumes is expected to be ∼10 2 .At this scale, when the postprocessing techniques are implemented with appropriate automations, the labor cost is not significantly affected by the total number of subvolumes.However, when the standard coordinate rotation is applied and the target volume discretization is carried out in accordance with the LES grid resolution, the total number of subvolumes readily exceeds 10 3 (as in this example study n x × n y × n z = 1920) the selective assembly phase becomes prohibitively laborsome.But, given sufficient computational capacity, the far-field correction approach can be exploited to perform the selective assembly process to provide a description for an effective target volume V T,eff = i,j,k∈K V i,j,k , which can subsequently be reassembled from the finely resolved V i,j,k contributions.Such a result is depicted in Fig. 7 together with two footprints that are obtained through an identically guided selection process utilizing far-field correction (FFC) with different combinations of n x , n y , and n z .The comparison reveals that the differences between the three results are remarkably insignificant, indicating, first, that a significant part of the error contributions, introduced by the farfield correction method, have a compensating effect and, second, that the obtained footprint is not highly sensitive to the sensor placement despite the variable flow conditions around the sensor.This demonstrates the utility and robustness of the selective piecewise postprocessing approach.From here on the presented results correspond to the n x = 3, n y = 5, n z = 3 target volume discretization level.

Outline of the procedure
Taking into account the far-field correction procedure, the postprocessing procedure for evaluating a footprint from a LES-LS obtained dataset can be described in the following steps.|| 2, indicating discrepancy between f REF and f (l) are shown where applicable.Acceptable candidates are marked by and the rejected by ×.Note the use of short-hand notation, e.g., 1:3 = 1, 2, 3.
1. Split the original dataset S into n x × n y × n z number of subsets labeled s i,j,k according to a Cartesian division of the target volume V T .
2. Evaluate an approximate footprint in a piecewise manner by applying Eq. ( 6) for each subset s i,j,k and assemble the result according to Eq. ( 9) by selecting all i, j, k values.(Here it is possible to use inaccurate data for the evaluation of w i,j,k as the objective is only to identify the far field).
3. Inspect the approximate footprint result to identify the extent of the far field (by specifying β) where the footprint reaches an asymptotic level to a good approximation, and specify β for the purpose of constructing r i,j,k .
4. Evaluate the sectional footprints f i,j,k from corresponding s i,j,k subsets by applying Eq. ( 6) with w * i,j,k evaluated through far-field correction approach as follows: a. select initial guess for c o i,j,k and w o i,j,k , and utilizing the data from r i,j,k compute the initial sectional footprint f i,j,k = f i,j,k (c o i,j,k ) and the corresponding far-field integral J o = β f i,j,k (c o i,j,k )dx ; b. perturb the coefficient c i,j,k = c o i,j,k + dc (initially with a guessed perturbation dc) and, using w * i,j,k = c i,j,k w o i,j,k and the data from r i,j,k , compute f i,j,k = f i,j,k (c i,j,k ) and J = β f i,j,k (c i,j,k )dx ; c. exit the loop if J < ε, where ε specifies the tolerance; d. compute derivative dJ dc = (J −J o ) dc and specify a new perturbation from dc = −γ dJ dc , where γ > 0 is a scaling parameter which, in this context, has been a experimentally set to ensure that the minimization problems converge sufficiently; e. set J o = J , c o i,j,k = c i,j,k and return to step 4b.
5. Select the appropriate set K of i, j, k combinations employing sensitivity analysis procedure in Sect.2.3.2.
Geosci.Model Dev., 10, 4187-4205, 2017 www.geosci-model-dev.net/10/4187/2017/(1) the far-field correction (FFC) method and selective assembly or (2) the standard coordinate rotation while utilizing the uniform 1 m resolution of the LES grid in the target volume discretization.The subvolume contributions included in the n x , n y , n z = 8, 20, 12 result were selected to correspond with the effective target volume V T,eff = i,j,k∈K V i,j,k determined via the selective assembly for n x , n y , n z = 3, 5, 3.
It is noteworthy that in step 2 for the approximate footprint evaluation and in step 4a for the initialization of the optimization loop, the values for the mean vertical velocities w i,j,k do not have to be accurate.Therefore, the use of vertical velocity data from LES can be omitted altogether, which simplifies the case setup and data handling considerably.The approximate values can be obtained more simply, for instance, by evaluating the mean of incident vertical velocity value from particle data in each s i,j,k .

Result assessment
The proposed methodology, founded on high-resolution LES-LS analysis and a piecewise postprocessing approach, has been shown to be a reliable, robust, and accessible, although computationally expensive, approach to generating topography-sensitive footprints in real urban applications.Since the underlying motivation for this development effort sprung from the need to evaluate the potential error that may arise when analytical, closed-form footprint models are applied to urban flux measurements, this work also proposes a technique to approximate the magnitude of this error in the absence of field validation studies.This approach hinges on the assumption that, in a real urban application, a topography-sensitive footprint obtained through a highly resolved LES-LS analysis features a higher level of accuracy and a lower level of uncertainty than any available closedform footprint model.
The proposed assessment technique compares the obtained LES-LS footprint result to an analytical model, which belongs to the group of closed-form models that would otherwise be employed in similar studies, by applying the footprint distributions to the land cover classification (LC) dataset in Fig. 1 that is presented at the same resolution as the topography height.In the following demonstration the closed-form footprint model by Korman and Meixner (2001) (KM), which is widely utilized in the EC community (e.g., Christen et al., 2011;Kotthaus and Grimmond, 2012;Nordbo et al., 2013), is used as an example analytical model.This choice is subjective and implies no preference over other available footprint models (e.g., Kljun et al., 2015;Horst, 2001).The KM model parameters and their specific values are declared in Table 3.The mean wind speed and the standard deviation of the crosswind component are extracted from Hotel Torni's anemometer measurements gathered on 9 September 2012 during the same 30 min time frame that was used to specify the meteorological conditions for the LES simulation (see Sect. 2.2.2).
A preliminary comparison between the obtained LES-LS and KM footprint distributions, f LES and f KM respectively, in the considered Hotel Torni case study draws immediate attention to the apparent differences that become discernible from the juxtaposition displayed in Fig. 8.The shown distributions have been normalized to yield f dx = 1 to aid the comparison.The LES-LS-generated footprint exhibits a complex, unpredictable probability distribution and a more pronounced spatial confinement, lacking the gradual asymptotic behavior of analytical models.In particular, the crosswind diffusion of the system is clearly overpredicted by the KM model even when the measurement height is taken to be the height of the sensor above the ground level minus the displacement height of 14.9 m, according to Nordbo et al. (2013).This value does take into account the surrounding buildings, but the ground level at Hotel Torni is 15 m above the sea level, which is also represented in the source area.
www.geosci-model-dev.net/10/4187/2017/Geosci.Model Dev., 10, 4187-4205, 2017 This exhibits the difficulty in choosing one representative parameter value for an analytical model applied to a real urban setting.The most evident deviations occur in the near field, where the f LES exhibits strong local variations between building tops and street canyons.Moreover, examining the crosswind integrated footprints in Fig. 9 reveals how f y LES reacts abruptly to changes in the example urban landscape, leveling off to a shallow descending slope much earlier than the gradually declining curve of f y KM .Thus, the presented comparison in the context of this case study succeeds in laying bare the nontrivial nature of urban footprints and highlights the importance of utilizing a high-resolution LES-LS approach to examine complex urban EC measurement sites.

Virtual assessment technique
The comparative technique proposed for assessing the potential error, that may arise if urban measurements are interpreted with closed-form footprint models, exploits the land cover dataset under the assumption that the LC distribution conveys the inherent urban heterogeneity sufficiently.Under this premise, the LC distribution can be adopted as a model distribution of sources Q such that each eth land cover type is assigned a constant mean source strength Q e = const.Thus, under this simplification the description of a measurement η in Eq. ( 1) can be decomposed as follows: (13) where the constituents of η are given by Here, A e is the footprint-weighted surface area of the eth land cover type and Now it is convenient to define two measures that facilitate a meaningful comparison between different footprints: the fractional contribution to the measurement from each constituent which require that Q e are assigned for each land cover type, and the source-area fraction that provide an easy estimate of the footprint's coverage independent of source strength information (or assuming identical Q e for all e).For proper assessment, these two fractions should be inspected in tandem.
The comparison is carried out by extracting the area corresponding to the LES domain from the LC dataset, shown in Fig. 10, which has been modified to include the relevant streets in the vicinity of the footprint for the purpose of including the effect of traffic emissions into the demonstration.The obtained f LES and f KM footprints are then projected onto this raster map to compute the required integrals and fractions.
A pie-chart of source-area fractions a e from Eq. ( 18) for the Hotel Torni's flux footprint is demonstrated in Fig. 11, which provides an informative overview on the differences in source-area coverages.The far-field-corrected (n x , n y , n z = 3, 5, 3) and the highly resolved (n x , n y , n z = 8, 20, 12) piecewise assembled LES-LS footprints agree within 0.2 %.In this particular example, the analytical KM model gathers a significantly larger contribution from the far-field, which is reflected in the significantly higher coverage of water surface area.However, assigning each land cover type its corresponding -potentially fictional -mean source strength Q e and evaluating the fractional contributions r e from Eq. ( 17) provides means to carry out simplified virtual experiments concerning particular EC measurements.To demonstrate with an example, consider CO 2 flux measurements in a hypothetical situation where 95 % of the CO 2 emissions originate from traffic (i.e., from roads) and 5 % arise from other anthropogenic sources, which are emitted through ventilation outlets on the building roofs.For the sake of simplicity, the contribution from water area is considered negligible and  vegetation is considered to act as a uniformly distributed sink over the land, which does not influence the ratio of source contributions in the measurement.Utilizing an undefined reference source strength Q ref , the sources are expressed as Q e = λ e Q ref , where the weights satisfy e λ e = 1.Thus, in this example λ 0 = 0.05 for buildings and λ 6 = 0.95 for roads.For this contrived situation the fractional contributions obtained with the LES-LS footprint become r 0 = 17.8 % and r 6 = 82.2%, whereas applying the Korman-Meixner footprint yields r 0 = 16.6 % and r 6 = 83.4%.In this example, while the two footprints have distinctly different source-area fractions for buildings and roads, their ratios are close since (A 6 /A 0 ) LES = 1.1(A 6 /A 0 ) KM , as seen in Fig. 11, which is the reason for obtaining such comparable measurement decompositions.
Repeating the introduced assessment technique for multiple representative meteorological conditions paves the way for a numerical approach that allows the obtained urban flux measurements to be interpreted either differently or with improved confidence.Naturally, having access to real source strength distributions opens up the ability to utilize LES-LS footprints (or positively assessed analytical footprints) to carry out detailed emission inventories (e.g., Christen et al., 2011).

Summary and conclusions
The utility of the eddy-covariance method in measuring the exchanges of mass, heat, and momentum between the urban landscape and the overlying atmosphere largely depends on the ability to determine the effective source area, or footprint, of the measurement.In situations where the heterogeneity of the surface becomes relevant, like for urban landscapes, and the structures surrounding the measurement site can no longer be considered as a homogeneous layer of roughness elements, the use of analytical footprint models becomes highly suspect.In order to diminish the resulting uncertainties and to obtain the ability to assess the applicability of analytical models, the ability to evaluate complex footprints with high resolution becomes essential.This work presents a numerical methodology to generate topography-sensitive footprints for real urban EC flux mea-  surement sites.This methodology is based on high-resolution LES-LS analysis where the simulation domain features a detailed description of the urban topography and accounts for the entire vertical extent of the atmospheric boundary layer.The online-coupled LS model within the LES solver is employed to simulate a constant release of inert gas emissions from the potential upwind source area of the considered EC sensor.The necessary data for the footprint generation are obtained from the LES-LS analysis by setting up a finite target volume around the sensor location and, over a sufficiently long simulation period, gathering a record of particles that hit this target.To generate an estimate for the flux footprint, this dataset is subjected to a postprocessing procedure that involves a coordinate rotation step, which eliminates the effect of the mean flow on the flux evaluation.But, if the considered EC sensor is mounted on a building (instead of a conventional tower-like structure) in the vicinity of which strong mean flow gradients occur, standard postprocessing techniques fail to produce physically meaningful footprints unless the target volume size is reduced to correspond with the LES grid spacing.This inevitably leads to prohibitive computational costs.Therefore, this work introduces a robust piecewise postprocessing strategy, which facilitates the evaluation of footprints despite the added complexity.The piecewise approach involves splitting the original dataset into a series of subsets which are all independently postprocessed to yield incompletely converged intermediate footprint estimates.The splitting is done by applying Cartesian discretization to the target volume in order to generate a series of subvolumes that correspond to the subsets.However, to facilitate a sufficiently accurate coordinate rotation treatment in the presence of strong gradients, the size of these subvolumes must also be reduced to match the resolution of the LES grid.This causes their number, and hence the number of intermediate sectional footprints, to become excessive, motivating the development of a new approximate scheme labeled far-field correction, which enables the subvolume size to be increased and the postprocessing effort to be reduced significantly.In the piecewise postprocessing approach, the final, completely converged, footprint is eventually selectively assembled from the obtained set of intermediates.
The methodology is demonstrated in a real urban application where the objective is to compute a highly resolved topography-sensitive footprint for the Hotel Torni EC flux measurement sensor mounted on the roof of a tall building situated in the downtown area of Helsinki, Finland.The EC sensor's measurement height is 60 m above the ground level and 36 m above the surrounding mean building height (24 m).The meteorological conditions for the LES simulation were adopted from measurements on 9 September 2012 when southwesterly winds and a neutrally stratified boundary layer of 300 m height were recorded.A detailed topography map of Helsinki at 2 m resolution from Nordbo et al. (2015) was utilized to construct the topography model for the LES-LS domain.The resolution of the computational mesh was set at 1 m throughout the domain to ensure that the relevant turbulent structures, even at the level of street canyons, were captured.An arbitrarily sized target box for sampling the Lagrangian particle hits was set up around the sensor location, which collected ca.19×10 6 particle hits during 3 h of simulation time.The obtained dataset was subjected to the proposed piecewise postprocessing method, demonstrating the functionality of the approach under various user-selected specifications.The obtained footprint stood in stark contrast to gradual ellipse-shaped analytical footprints: the distribution exhibited strong adherence to the building block arrangement in the near field, where the weight distribution changed abruptly between roof tops and street canyons.In comparison to the Korman and Meixner (2001) model, the LES-LS footprint also exhibited stronger contribution from the near field, but more rapidly diminishing contribution from the far field.
This paper also introduces an accessible technique to employ the obtained high-resolution topography-sensitive urban footprint in estimating the potential error that may arise when an analytical footprint model is used to interpret urban EC measurements.The underlying stipulation for this method is that it does not require knowledge of real source strength distributions.Thus, it is proposed that a detailed land cover type classification (LC) dataset is utilized as a model source strength distribution map for the urban surroundings assuming that it reflects the heterogeneity of the urban conditions sufficiently.Projecting a footprint distribution result onto such a LC map enables the evaluation fractional contributions, which indicate how each land cover type is represented in the measurement.This procedure provides a comparative technique to assess the effective deviations between different footprints.The demonstrated comparison between the LES-LS and analytical KM footprints in the EC measurement setup in Helsinki revealed substantial differences in the fractional contributions when all land cover types are considered equally relevant.The technique can also be applied by considering only selected land cover types and assigning each of them a variable source strength.This approach is demonstrated through a simple example, which mimics a hypothetical CO 2 flux measurement, where the effective source area is limited to only roads and buildings.
The context of this paper is limited to laying out the new methodology for generating urban footprints and exploiting them in the assessment of analytical models.It is evident that changes in the meteorological and anthropogenic conditions will influence the results and a proper assessment of the applicability of analytical models at a given EC measurement site will require that these conditions are varied, necessitating numerous footprint evaluations.This paper lays the numerical groundwork for such future investigations.

Figure 1 .
Figure 1.Raster maps of topography height h (a) and land cover types LC (b) from Helsinki area.The rectangle in the bottom left corner is aligned with southwesterly wind and represents the area of interest for the footprint analysis.In the surface type classification each pixel (2 m) is categorized according to the following numbering: 0 = building, 1 = impervious (rock, paved, gravel), 2 = grass, 3 = low vegetation, 4 = high vegetation, and 5 = water.

Figure 2 .
Figure 2. Visualization of the topography height distribution underlying the LES domain.The particle release area is enveloped by a white dashed line.The size of the precursor domain is outlined in the top left corner.The location of the turbulence recycling plane is marked by a black dotted line at x rc .

Figure 3 .
Figure 3.A three-dimensional rendering of the urban topography near Hotel Torni (a) and a close-up featuring the target box (T) for particle capturing (b).

Figure 4 .
Figure 4. Example discretization of target volume V T into n x × n y × n z subvolumes.A coarse illustration with n x = 2, n y = 3, and n z = 2 is shown.

Figure 5 .
Figure5.Crosswind integrated distributions of piecewise postprocessed footprints obtained with different levels of target volume discretizations using w i,j,k from LES solution.Anomalous contributions from subvolumes immediately behind or in contact with the tower structure were omitted from the assembly (refer to Sect.2.3.2).

Figure 6 .
Figure 6.Illustration of the selective assembly of the final footprint for n x ×n y ×n z =3×5×3=45.Values of || f (l)|| 2, indicating discrepancy between f REF and f(l) are shown where applicable.Acceptable candidates are marked by and the rejected by ×.Note the use of short-hand notation, e.g., 1:3 = 1, 2, 3.

Figure 7 .
Figure 7.Comparison of identically normalized footprint distributions (b) and their crosswind integrations (a) obtained either by applying(1) the far-field correction (FFC) method and selective assembly or (2) the standard coordinate rotation while utilizing the uniform 1 m resolution of the LES grid in the target volume discretization.The subvolume contributions included in the n x , n y , n z = 8, 20, 12 result were selected to correspond with the effective target volume V T,eff = i,j,k∈K V i,j,k determined via the selective assembly for n x , n y , n z = 3, 5, 3.

Figure 8 .
Figure 8.Comparison of identically normalized LES-LS (b) and KM (a) footprint distributions merged with the urban topography model of Helsinki.The location of the EC sensor (Hotel Torni) building is indicated with a white circle.

Figure 9 .
Figure 9.Comparison of normalized, crosswind integrated LES-LS and KM footprints.A light blue dashed line indicates the start of urban topography and the gray dashed line marks the location of the EC sensor.
Figure 10.Raster map of land cover types, LC, within the LES domain.The original surface type classification data in Fig. 1 has been augmented by adding streets (LC = 6) to the relevant footprint area.

Figure 11 .
Figure 11.Comparison of source-area fractions a e resulting from applying f LES and f KM to the raster map of land cover types in Fig. 10.

Table 2 .
Diagnostic data from the application of far-field correction in the coordinate rotation.The farthest 15 % of the source area in the LES domain is considered (i.e., β = 15).

Table 3 .
Parameters used in the Korman and Meixner footprint model.