Associate Editor: Paul Palmer; The University of Edinburgh, UK
1. Introduction
In contrast to the ‘bottomup’ approach to quantifying CO_{2} emissions, which involves compiling and estimating the characteristics of known emission sources into a quantitative database such as the US Environmental Protection Agency (EPA) Greenhouse Gas (GHG) Inventory^{1} or the high resolution inventory Hestia (Gurney et al. 2012), the ‘topdown’ approach involves measuring CO_{2} concentrations in the atmosphere and then using an atmospheric transport model to infer emissions (e.g., Enting and Mansbridge 1989; Bréon et al. 2015). The topdown approach can be effective in cases for which a source has a known location (e.g., Davoine and Bocquet 2007), but in situ source magnitude characterization is unavailable or unreliable, in which case concentration measurements in the vicinity of the source can help provide the missing information. This scenario is often true in the urban setting, for which a small number of essentially point or line sources at fixed locations may produce a large fraction of the area’s total GHG emissions (e.g., in Indianapolis approximately onethird of the total emissions is from the Harding St. power plant). Except for the simplest cases when known analytical solutions exist (e.g., the Gaussian plume model (Seinfeld 1986)), an atmospheric transport model is required to relate emission rates to concentrations and vice versa.
The challenge with the topdown method is that the relation between emission magnitude and subsequent concentrations is not always straightforward, and prone to errors including representativeness error and model error (Gurney et al. 2002; Tarantola 2005; Gerbig et al. 2008). Representativeness error derives from a mismatch between modeled and observed fields due to unresolved scales (Kookhan and Bocquet 2012). Random errors can be reduced by increasing the sample size of modelobservation pairs used in the inversion. A number of issues may lead to systematic errors in the transport model used in the topdown method. Systematic errors emerge if the distance to a point source is on the order of or smaller than the resolution of a transport model grid (Wilson and Sawford 1996). Transport error can also be due to an erroneous prediction of the prevailing ‘mean’ component of the wind, incorrect simulation or parameterization of turbulent perturbations to that wind, and errors introduced by model physical parameterizations (e.g., Lin and Gerbig, 2005). Errors in the surface energy balance, temperature profiles, and planetary boundary layer (PBL) structure may all indirectly lead to transport error by causing errors in the mean or turbulent wind field (Lauvaux and Davis, 2014; Sarmiento et al. 2017; Deng et al. 2017).
We use the Weather Research and Forecasting (WRF) model (Skamarock et al. 2008) to simulate atmospheric transport. It is a stateofthescience atmospheric model capable of performing both idealized and realistic simulations of the atmosphere with a full suite of model physics, including atmospheric dynamics, cloud physics, radiation, and surface fluxes. For global to mesoscale applications (e.g., Seaman et al. 2012), WRF predicts gridcellaveraged (i.e., ‘mean’, or modelresolved) atmospheric fields such as momentum, temperature, and moisture, as well as gridcellaveraged turbulent variances and fluxes. PBL turbulence is parameterized. WRF can also be run as a large eddy simulation (WRFLES) (Moeng et al. 2007). The LES method was developed in the 1960s–1970s (Smagorinsky 1963; Lilly 1967; Deardorff 1972). In LES the largest turbulent eddies are explicitly simulated, while a subgrid closure is used to parameterize the turbulent diffusion at smaller scales, including the continuous removal of energy by the downscale turbulent cascade of energy. LES thus offers a more explicit and hopefully realistic simulation of turbulent processes in the PBL. There is also a chemistryadapted version of WRF, WRFChem (Grell et al. 2005), that can simulate the chemical reactions and transport of various species, including CO_{2}.
During the INdianapolis FLUX experiment (INFLUX), twelve (12) communication towers in and around Indianapolis, IN were instrumented with towerbased continuous GHG analyzers (Miles et al. 2017) to monitor the urban emissions from fossil fuel and biogenic sources (Davis et al. 2017). To infer the urban emissions of GHGs, a high resolution atmospheric inversion system was developed using the WRF mesoscale atmospheric model coupled to a Lagrangian model, buildinglevel CO_{2} emissions from the Hestia product, and towermeasured atmospheric CO_{2} mole fractions (Lauvaux et al. 2016). The mesoscale inversion system produced adjustments to the Hestia CO_{2} emissions of about 20% over eight months (mostly fall and winter), possibly due to an underestimation of CO_{2} sources in Hestia, but also possibly due to errors in the transport model. This urban inversion system relies on the use of a mesoscale model run at 1 km resolution to simulate the observed atmospheric mole fractions of GHGs. Large GHG sources are sometimes within a few kilometers of the tower observation points. The mesoscale model may have systematic errors in its description of turbulent dispersion from these “nearfield” sources.
The largest eddies in the daytime convective PBL tend to be about 1.5 times the depth of the PBL (Kaimal et al. 1976), which is typically on the order of 1 km. Resolving individual turbulent eddies for a region as large as an entire city like Indianapolis for months to years is not computationally tractable. In the vertical dimension, turbulent fluxes and other atmospheric and tracer fields vary significantly over even smaller scales within the PBL, on the order of hundreds of meters or smaller. The typical configuration for a mesoscale model like WRF for the convective PBL employs vertical grid spacing on the order of tens of meters or less, but horizontal grid spacing on the order of a kilometer or more (e.g., Seaman et al. 2012). A onedimensional (vertical only) PBL parameterization predicts the gridcellaveraged statistics of turbulent variances and fluxes through variance closure assumptions; tendencies of modelresolved fields due to vertical turbulent transport can then be determined from the vertical gradients of the resolved fields and the PBL scheme output. Most WRF PBL schemes predict turbulent kinetic energy (TKE) and combine this with a length scale on the order of the size of the largest turbulent eddies to generate a vertical diffusion coefficient (K_{Z}), which can then be used to predict the vertical turbulent transport of CO_{2} and other species with a gradientdiffusion formulation. Horizontal turbulent transport is not explicitly predicted by the PBL scheme because the horizontal grid spacing is considerably larger than the scale of the largest PBL eddies, although for numerical stability reasons some horizontal diffusion operating on the scale of the horizontal grid spacing is typically present.
These parameterizations of PBL turbulence are not designed to simulate pointsource dispersion at scales less than or similar to the horizontal grid resolution of the mesoscale model. In addition, using K_{Z} with a length scale on the order of the size of the largest turbulent eddies is not appropriate when the model horizontal grid spacing becomes less than the size of the largest eddies. With this relation of grid spacing to large eddy size we are entering the “terra incognita” regime for which the modeling of turbulence is challenging (Wyngaard 2004). Furthermore, as discussed later, any formulation of K_{Z} solely dependent on integrated properties of the PBL turbulence does not adequately describe vertical scalar diffusion in the near field, as outlined in the seminal paper by Taylor (1921). LES with reduced grid spacing is however a suitable tool for numerical simulation of turbulent dispersion in these situations. The largest PBL eddies are resolved, and in the midPBL are the ones that dominate the turbulent diffusion process (e.g., Moeng and Wyngaard, 1984). Thus, both vertical and horizontal turbulent transport of scalars can be simulated more explicitly.
Our objective is to analyze the potential errors and biases of plume concentrations near a point source (i.e., within a few km) simulated using a mesoscale configuration of WRFChem, with a baseline model grid spacing of 1 km. Our procedure is to drive the WRFChem baseline configuration for the case of 28 September 2013 and predict the transport of CO_{2} from a known point source in the INFLUX study domain. We will assess the errors and biases of the baseline simulation with a nested domain that features a combination of WRFChem and WRFLES run at 111 m horizontal grid spacing, as well as other combinations of horizontal grid resolution and turbulent physics for sensitivity testing purposes. We will also use two versions of a modified Gaussian plume equation that are analytical solutions to the steady state pointsource diffusion equation for an idealized daytime PBL, and which approximate the differences in nearfield dispersion between the mesoscale model parameterized diffusion and the resolved turbulent transport from the LES. In combination these tools will provide evidence that the WRF mesoscale configuration inherently overpredicts vertical diffusion in the nearsource region. Section 2 will describe in more detail the methods used, the WRF model configurations, case study, analysis plan, and theoretical analytical solutions. Section 3 will present the results from the baseline comparison between the mesoscale and LES configurations, and show to what extent differences between the 1km mesoscale simulation and the 111m LES are due to horizontal resolution, or due to inherent limitations of the K_{Z}–based approach for nearfield plume dispersion. Section 4 will discuss some of the implications for the generation of concentration footprint functions for inversions by mesoscale models. Discussion is in section 5, and section 6 will present the summary and conclusions.
2. Methods
a) WRF model configuration
As stated in the introduction, we use the Weather Research and Forecasting (WRF) model (Skamarock et al. 2008) to investigate the transport and dispersion of a CO_{2} plume released in the INFLUX domain, focusing on the nearfield effects. ‘Nearfield’ will be defined more precisely later, but in general we mean distances on the order of a few km. In all cases we used the chemistryadapted version of WRF, WRFChem (Grell et al. 2005), for which CO_{2} was treated as a passive tracer with a surface emission source. Modelresolved transport for tracers is predicted by the dynamic core of WRF in the same manner as other conserved variables in default WRF. For the vertical turbulent transport of tracers in the mesoscale configuration of WRFChem, the model PBL scheme is used to derive the vertical eddy diffusion coefficient for scalars (K_{Z}), which is assumed to be the same as the coefficient used for heat transport. Then K_{Z} is passed to the chemistry module for the vertical turbulent tracer transport prediction. An advantage of using WRFChem to predict tracer transport is that the chemistry is fully inline with the atmospheric flow field prediction, so both can mutually evolve each model timestep.
To study nearfield dispersion, we use WRF to model the evolution of the plume emitted from the location of the Harding St. power plant (39.71°N, 86.20°W), located in the southwestern sector of Indianapolis. The Noah land surface model was used for the simulation, with some tiles reclassified based on the 2006 version of the 30m National Land Cover Database (NLCD) (Homer et al. 2015; Sarmiento et al. 2017). Timedependent emission data from the plant in 2013 were supplied by data obtained from the EPA Clean Air Market Division (CAMD) Emission Tracking System/Continuous Emissions Monitoring system (ETS/CEMs) for electrical generation, also used in the Hestia database over Indianapolis (Gurney et al. 2012). We note here that the primary objective of this study is to gain insight into the biases of pointsource plumes in typical daytime PBL conditions using the baseline mesoscale configuration of WRFChem, rather than to verify the observed plume structure for this particular day. For simplicity, WRFChem assumed that the power plant emission occurred at the surface (in reality, the lowest halflayer height, at approximately 7 m above ground level (AGL)).
A fivedomain oneway nested grid configuration was used for the simulations, with horizontal grid spacing of 9 km, 3 km, 1 km, 333 m, and 111 m, respectively, as seen in Figure 1. For the coarsest three model domains, the MellorYamadaNakanishiNiino (MYNN) Level 2.5 PBL scheme was used (Nakanishi and Niino 2006). This scheme was used for previous INFLUX studies (Lauvaux et al. 2017; Deng et al. 2017) and has been shown to simulate realistic and fairly accurate PBL depth and thermodynamic profiles for this city (Deng et al. 2017; Sarmiento et al. 2017). Among the domains is the 1km domain (Domain 3) considered the ‘baseline’ mesoscale simulation. For the reference simulation, the 333m domain (Domain 4) is also run using the MYNN scheme, recognizing that with this horizontal grid spacing the ‘terra incognita’ regime is being entered, in which turbulent eddies are beginning to be resolved, and neither mesoscale nor LES may be completely adequate (Wyngaard 2004). For the 111m domain (Domain 5), WRFLES is used (Moeng et al. 2007), with the 1.5order TKE subgrid closure of Deardorff (1980).
Because all nesting is oneway, all coarser domains can be considered to be independent of the evolution of the finer domains, while the finer domains only receive information from the coarser domains through their lateral boundaries during the simulation. Thus, a comparison of the plumes on Domain 3 and Domain 5 is a fair comparison between the effects of a 1km mesoscale configuration versus a 111m LES configuration. However, a separate simulation was run using a 111m domain with the MYNN parameterization, for the purposes of performing a sensitivity test on the impact of turbulent physics alone. Model parameters for all domains are shown in Table 1.
Domain 1  Domain 2  Domain 3  Domain 4  Domain 5  

Grid Spacing (km)  9 km  3 km  1 km  333 m  111 m 
Horizontal Dimensions (grid points)  101 × 101  100 × 100  88 × 88  151 × 151  151 × 151 
Vertical Levels  60  60  60  60  60 
Turbulent Physics  Mesoscale MYNN*  Mesoscale MYNN  Mesoscale MYNN  Mesoscale MYNN  LES Deardorff (Mesoscale MYNN) 
Land Surface Model  Noah  Noah  Noah  Noah  Noah 
The vertical grid spacing is the same for all nested domains, with 60 vertical levels (corresponding to 59 vertical layers). The vertical grid spacing is variable, increasing 10–25% per model layer per level. The midpoint of the first model layer is approximately 7 m above the surface, and there are eight model layers within about 200 m of the surface (refer to Supplementary Table S1). Ideally an isotropic grid would be used for LES, but we use the selected vertical grid configuration because it has been extensively used in previous INFLUX WRF simulations (Lauvaux et al. 2016; Deng et al. 2017) and allows for reasonably fast and stable simulations for our required temporal and spatial domains.
WRFLES is an Eulerian grid implementation of LES that can be used as a nested grid in both realcase and idealized simulations. One advantage of WRFLES is that it can be run as a fully heterogeneous nested grid receiving realistic lateral boundary conditions from parent domains and using the same physics suites as the parent domains, except for the representation of turbulence. WRFLES thus can be applied to realistic case study applications.
Our particular configuration was chosen based on specified computing resources to balance the competing concerns of sufficient horizontal resolution and sufficient horizontal areal coverage, both of which can be important for LES applications. We use 111m horizontal grid spacing, which is consistent with many studies using WRFLES (Moeng et al. 2007; Wang and Feingold 2009), but is coarse for representing the details of the PBL turbulence spectrum, since the effective horizontal resolution in WRF is around 7Δx (Skamarock 2004). It is also well known that the LES method can break down near the lower boundary, where the size of eddy structures are reduced and inherently nonresolved by the grid (e.g., Khanna and Brasseur 1998). So in this study we focus on the impacts of the largest PBL eddies on turbulent transport above the lowest few model levels, and over horizontal distances greater than several hundred meters.
Since by necessity the realistic lateral boundary conditions are nonperiodic, it takes a finite time for LES to spinup realistic turbulent eddies, either from a nonturbulent initial state or from nonturbulent lateral boundary conditions. Without periodic lateral boundary conditions, regions too near the inflow boundaries might never develop realistic eddy structures (Gaudet et al. 2012). The horizontal extent of the LES domain, about 17 km, is marginal for this purpose. In this study we focus only on plume behavior in local afternoon, by which time a mature daytime PBL has had a chance to develop, and we have the plume release near the center of Grid 5, far from the lateral boundaries, in order to mitigate the turbulent spinup problem.
b) Case study
The case of 28 September 2013 was selected for simulation. While observational verification of the meteorological details of a particular case is not a focus of this study, we did want to simulate the evolution of a simple daytime PBL. This case was chosen because it is free of precipitation, major synoptic weather changes, and low clouds (none within 6 km of the ground). By the afternoon fairly steady southwest flow was present, allowing the plume from the Harding St. power plant to become approximately steadystate for a few hours.
The simulation period was 1200 UTC 28 Sep 2013 through 0000 UTC 29 Sep 2013 (0700–1900 LST). No spinup of turbulent eddies in the LES domain was performed; as mentioned above, the gradual development of eddy structures from random thermal perturbations was sufficient because analysis was only performed on the mature convective PBL during local afternoon. No data assimilation was used to constrain the meteorological model fields within the model domains, but the 32km North American Regional Reanalysis (NARR) meteorological product was used to provide initial and lateral boundary conditions for the domains, as described in Deng et al. (2017).
c) Analysis plan
Current urban inversion systems attempt to reproduce timeaveraged GHG concentrations (e.g., Lauvaux et al. 2016; Staufer et al. 2016). Specifically, the INFLUX urban inversion system (Lauvaux et al. 2016) attempts to match timeaveraged, towerbased GHG concentrations with concentrations simulated using prior flux estimates and a Lagrangian footprint model that utilizes mean wind and turbulence fields from a mesoscale model with 1 km horizontal grid spacing. While the mesoscale simulations are intended to represent the mean concentration, an average over multiple eddies, the output from the LES contains explicit realizations of turbulent eddies and should be averaged over multiple model output times for comparison. The comparison between LES and mesoscale simulations is easiest to interpret when largescale forcing of the flow is not changing, and our LES resolves the largest fraction of the TKE spectrum when the PBL depth is the largest. For both reasons, we focus on the period from 1700–1900 UTC (1200–1400 LST) when the wind direction was relatively steady with respect to time and the boundary layer was deep. The LES output (every 10 minutes) was averaged over this period, and compared to the 1800 UTC (1300 LST) 1km mesoscale output.
In the dry daytime convective boundary layer, the vertical scale of the largest turbulent eddies is on the order of the mixedlayer height z_{i}, and their velocity scale is on the order of the Deardorff convective velocity scale, ${w}_{\ast}\text{\hspace{0.17em}}\equiv \text{\hspace{0.17em}\hspace{0.17em}}{\left({\scriptscriptstyle \frac{\mathit{\text{gH}}{z}_{i}}{{\theta}_{0}}}\right)}^{1/3}$ (Deardorff 1970; Young 1988). Here g is gravitational acceleration, θ_{0} is a reference potential temperature, and H is the magnitude of the surface sensible heat flux. If H and z_{i} are consistent between the mesoscale and LES configurations, we can achieve a controlled comparison of the simulated PBLs. We impose the same atmospheric initial conditions and land surface model in the mesoscale and LES configurations to ensure this consistency and comparability.
The timescale of the largest turbulent eddies in the convective PBL, which governs how quickly tracers become vertically wellmixed in the convective CBL, is z_{i}/w_{*}, and is usually on the order of 10^{3} seconds. So analysis of the vertical and downwind properties of the model plumes will be expressed in terms of the vertical distance z_{i}, and the characteristic downwind distance ${x}_{0}\equiv \text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}{\scriptscriptstyle \frac{U{z}_{i}}{{w}_{\ast}}}$ that tracers move in an eddy turnover time.
In order to better distinguish plume structure sensitivity due to differences in horizontal resolution from differences in model physics, we ran a sensitivity experiment, the same as the baseline experiment but this time with the 111m domain run using the mesoscale PBL closure. Normally this would be considered too fine a resolution for use of a PBL closure in the convective boundary layer (there would be the risk of double counting vertical turbulent transport if turbulent eddy structures are resolved in conjunction with the transport predicted by the PBL scheme) – but here we specifically want to isolate the effects of resolution on the nearfield plume from that of the turbulent physics scheme.
d) Gaussian plume nearsource analytical solutions
As an independent test of the plausibility of the WRF plume structures, we can use known analytical solutions of plume structure in certain idealized cases. In particular, suppose there is a tracer continuous point source (CPS) at the surface of strength Q (mass per unit time), in horizontally homogeneous meteorology with a mean constant horizontal wind U. The boundary conditions are that the tracer concentration c goes to zero as the crosswind coordinate, y, goes to ±∞ and the vertical coordinate z goes to +∞; at z = 0 is an impermeable barrier. Let us further assume that turbulent transport occurs in both y and z according to gradient diffusion, with eddy diffusivities K_{y} and K_{z}; in the downwind dimension x advection by the mean wind dominates turbulent transport (‘slender plume approximation’). Under these conditions the governing equation for the steadystate concentration is:
Here δ is the delta function. This equation may be solved by making the transformation t = x/U, which can be thought of as converting downwind distance to the time since release for a tracer element; by the slender plume approximation the two are related with proportionality factor U for all tracer elements. When K_{y} and K_{z} are arbitrary functions of x only, the solution is the wellknown Gaussian plume solution (Seinfeld 1986; Arya 1995):
for which the standard deviations of the concentration in the y and z dimensions, σ_{y}(x) and σ_{z}(x), are related to the eddy diffusivities by ${K}_{y}=\text{\hspace{0.17em}\hspace{0.17em}}{\scriptscriptstyle \frac{U}{2}}{\scriptscriptstyle \frac{d}{\mathit{\text{dx}}}}{\sigma}_{y}^{2}$ and ${K}_{z}=\text{\hspace{0.17em}\hspace{0.17em}}{\scriptscriptstyle \frac{U}{2}}{\scriptscriptstyle \frac{d}{\mathit{\text{dx}}}}{\sigma}_{z}^{2}$ . Despite all of the assumptions, this is a passable approximation to the plume behavior in the mature afternoon PBL of our case study, provided expressions for the diffusivities and/or σ can be found.
Taylor (1921) first derived expressions for the mean squared displacement of tracer particles released into the turbulent PBL by using a Lagrangian (tracerfollowing) framework, showing that the mean squared displacement of particles in dimensions Y and Z, at a time t since release, can be given by:
when t is large compared to T_{L}, the Lagrangian time scale of the turbulence, defined as (e.g., Dosio et al. 2005):
The Lagrangian time scale should be on the order of the turnover time of the energycontaining (largest) eddies. In (3) and (4), V_{Y} and V_{Z} are the velocity components in the horizontal crosswind and vertical dimensions, respectively; the primes indicate that these are perturbations from the mean wind. The brackets denote statistical averages performed over all particles with time t since their release. We can use X = Ut to transform (3) to an xdependent form. If then the bracketed averages are taken over sufficiently large samples (as would be the case for long time averages at stationary sensors such as a tower), the velocity variances on the RHS of (3) approach the total velocity variance in each dimension, denoted by ${\sigma}_{V}^{2}$ and ${\sigma}_{W}^{2}$ , while the squared displacements on the LHS of (3) become ${\sigma}_{Y}^{2}$ and ${\sigma}_{Z}^{2}$ . These ${\sigma}_{Y}^{2}$ and ${\sigma}_{Z}^{2}$ expressions can be used in the plume equation (2) provided crosswinddependencies of the velocity variances are neglected. The result is the longtime (or equivalently farsource) formula for the CPS steadystate plume under gradient diffusion:
with ${K}_{Y}=\text{\hspace{0.17em}\hspace{0.17em}}{T}_{\mathit{\text{LY}}}{\sigma}_{V}^{2}$ and ${K}_{Z}\text{\hspace{0.17em}\hspace{0.17em}}=\text{\hspace{0.17em}\hspace{0.17em}}{T}_{\mathit{\text{LZ}}}{\sigma}_{W}^{2}$ (e.g., Wyngaard and Weil 1991). In the daytime convective PBL, T_{LZ} scales as z_{i}/w_{*}, while ${\sigma}_{W}^{2}$ scales as w_{*}^{2}; thus K_{Z} is proportional to w_{*}z_{i} and to first order is solely a function of the surface sensible heat flux and PBL height. In this long time/farsource framework K_{Z} is solely a function of PBL turbulence statistics.
However, as noted in multiple studies (e.g., Seinfeld 1986; Arya 1995; Degrazia et al. 2001; Moreira et al. 2005; Kumar and Sharan 2010), Taylor (1921) found an alternate formula appropriate to short time/nearsource plume behavior. For t << T_{L} the mean squared displacement of particles is given by:
Conceptually, in the longtime regime, tracers experience a series of largely uncorrelated velocity increments, and as in the ‘random walk’ scenario mean squared displacement grows linearly with time. In the shorttime regime, however, tracer elements largely maintain their initial velocity perturbation, and the plume lateral extent is mainly a function of the velocity variance at the release point; mean squared displacement thus grows as the square of time. If we again use a large sample in the bracketed statistical average, and substitute into (2), we get the nearsource relations (valid for x << UT_{L}):
where formally ${K}_{Y}\text{\hspace{0.17em}\hspace{0.17em}}=\text{\hspace{0.17em}\hspace{0.17em}}x{\sigma}_{V}^{2}\text{\hspace{0.17em}}/\text{\hspace{0.17em}}U$ and ${K}_{Z}=x{\sigma}_{W}^{2}\text{\hspace{0.17em}}/\text{\hspace{0.17em}}U$ . Thus the nearsource plume equation has a reduced effective eddy diffusivity compared to the farsource plume, but one that increases with downwind distance towards its farsource value at x ≈ UT_{L}.
The most fundamental difference in nearsource (close to the source with respect to x_{o}) eddy diffusivity in comparison to farsource (far from the source with respect to x_{o}) eddy diffusivity is that it is no longer solely a function of PBL properties, but is also a function of x, which cannot be defined independently of source location. Arguably the whole concept of eddy diffusivity breaks down (Wyngaard 2010). In particular, this means that a mesoscale PBL parameterization that predicts K_{Z} from PBL statistics will not be correct in the nearsource region, and in fact will systematically overestimate nearsource vertical diffusion, assuming that it correctly predicts the farsource value of ${K}_{Z}\text{\hspace{0.17em}\hspace{0.17em}}=\text{\hspace{0.17em}\hspace{0.17em}}{T}_{\mathit{\text{LZ}}}{\sigma}_{W}^{2}$ . WRFLES, however, does have the capability to capture the correct nearsource plume behavior, because turbulent transport of scalars by the resolved PBL eddies is through advection by the resolved wind field rather than through an eddy diffusivity closure.
One consequence of the different behavior of (5) and (7) can be deduced by considering the yintegrated form of these equations for simplicity. The details are found in the supplementary material (Supplementary Material Text S1), but it can be shown that for the farsource version of the Gaussian plume the concentration surfaces near the edge of the plume ascend nearly vertically near the source. For the nearsource plume, however, the concentration surfaces possess a finite slope near the source (in fact, they approach a linear slope).
e) Alternate convective PBL analytical plume solution
Though the Gaussian plume equation (2) is usually a good description of farsource horizontal diffusion (e.g., Seinfeld 1986), is has inadequacies for vertical diffusion, even in idealized cases. Most importantly, the turbulent eddies in the PBL cannot transport tracers above the PBL top (apart from processes like convective venting), so the farsource behavior in the vertical is not characterized by continuous spreading, but by the PBL becoming increasingly wellmixed. This effect can be incorporated in the Gaussian plume model by imposing a reflective boundary condition at z_{i}; however, the resultant CPS solution becomes an infinite series of cosine functions. Second, the Gaussian plume model assumes no direct dependence of velocity variances on z, which is not typically a good assumption. In the daytime PBL the horizontal velocity variances are relatively independent of height, except near the lower boundary and PBL top (e.g., Wyngaard and Weil 1991); however, the vertical velocity variance, ${\sigma}_{W}^{2}$ , in the daytime convective PBL has a maximum value around 0.4z_{i} (Stull 1988; Moeng et al. 2007), and varies continuously as a function of altitude. In practice, applications of the Gaussian plume tend to use heuristic functions of downwind distance to take into account the different stability properties of different types of boundary layers (e.g., Venkatram 1996).
The WRF mesoscale PBL scheme does have the ability to predict the vertical dependence of K_{Z} and the presence of little or no turbulent transport above z_{i}, so in this respect mesoscale WRF should be superior to the Gaussian plume model. However, it is still subject to the limitation of not accounting for nearsource behavior. So for an independent comparison to the WRF mesoscale plumes, we would like an analytical solution that also includes the observed vertically dependent properties of the daytime PBL, but unlike the WRF mesoscale configuration includes both the correct nearsource and farsource behavior.
For the 2D case of c = c(x, z) and zero flux conditions at z = 0 and z_{i}, it has been shown (Wortmann et al. 2005; Kumar and Sharan 2010) that exact solutions to the CPS equation can be found when K_{z}(x, z) = f(x) g(z), for which nearsource effects can be included in f(x). This separability assumption leads to a SturmLiouville problem that may be solved given a particular form of g(z). However, the resultant procedures do not lead to readily analyzed closedform expressions except in specialized cases.
Here, we consider the full 3D case of (1). We assume that K_{y} has the form ${K}_{y}(x)\text{\hspace{0.17em}\hspace{0.17em}}=\text{\hspace{0.17em}\hspace{0.17em}}{K}_{\mathit{\text{yfar}}}f(x\text{\hspace{0.17em}}/\text{\hspace{0.17em}}U{T}_{L})\text{\hspace{0.17em}\hspace{0.17em}}\equiv \text{\hspace{0.17em}\hspace{0.17em}}{K}_{\mathit{\text{yfar}}}[1\mathrm{exp}({\scriptscriptstyle \frac{x}{U{T}_{{}_{L}}}})]$ , which smoothly merges the Taylor nearsource and farsource diffusivity limits. We assume that K_{z} has this same xdependence but also has a zdependence in the separable form ${K}_{z}(x,\text{\hspace{0.17em}\hspace{0.17em}}z)\text{\hspace{0.17em}\hspace{0.17em}}=\text{\hspace{0.17em}\hspace{0.17em}}4{K}_{\mathit{\text{zfar}}\hspace{0.17em}\text{max}}\widehat{z}(1\widehat{z})f(\widehat{x})$ , where $\widehat{z}\text{\hspace{0.17em}\hspace{0.17em}}\equiv \text{\hspace{0.17em}\hspace{0.17em}}z\text{\hspace{0.17em}}/\text{\hspace{0.17em}}{z}_{i}$ and $\widehat{x}\text{\hspace{0.17em}\hspace{0.17em}}\equiv \text{\hspace{0.17em}\hspace{0.17em}}x\text{\hspace{0.17em}}/\text{\hspace{0.17em}}U{T}_{L}$ . While this is not the most precise fit to the vertical structure of K_{z} in the convective PBL (e.g., Degrazia et al. 2001), it does at least capture its midPBL peak and smaller values both at the surface and at z_{i}.
Given these functions, a unique regular analytic solution exists to the steadystate CPS equation, as shown by Nieuwstadt (1980) and Otte and Wyngaard (1996). More details appear in the supplementary material (Supplementary Material Text S2), but the solution is:
Here P_{n}( ) are the Legendre polynomials, and a = 0.4, b = 0.2, and c = 0.4, where c ≡ UT_{L}/x_{0}; thus $\widehat{x}=x\text{\hspace{0.17em}}/\text{\hspace{0.17em}}c{x}_{0}$ . Finally $q(\widehat{x})$ is the integral of the Taylor diffusivity modification function $f(\widehat{x})$ :
The yintegrated concentration equation is:
One convenience of this form is that we can test the impact of neglecting the nearsource Taylor behavior in the Legendre analytical solutions by assigning $f(\widehat{x})\text{\hspace{0.17em}\hspace{0.17em}}=\text{\hspace{0.17em}\hspace{0.17em}}1$ and $q(\widehat{x})\text{\hspace{0.17em}\hspace{0.17em}}=\text{\hspace{0.17em}\hspace{0.17em}}\widehat{x}$ instead of using (9). Henceforth we will call the solution that neglects the nearsource Taylor behavior as the ‘farsource’ solution, whereas the one that uses (9) will be called the ‘blended’ solution. We will use this equation in direct comparison with the WRF mesoscale and LES plumes to determine how closely they resemble the farsource and blended Legendre solutions, respectively.
3. Results
We begin by comparing the general PBL structures as simulated by the 1km mesoscale and 111m LES configurations around the 1800 UTC = 1300 LST reference time. The MYNN PBL parameterization of the 1km configuration predicts a mean boundary layer height of 1400 m in the domain, which is consistent with the welldefined top of the model mixed layer as seen in a potential temperature crosssection through the power plant location (Figure 2, left panel). The mean depth grows slightly with time in the next few hours, but not by more than about a hundred meters. Boundary layer height is not explicitly predicted by the LES, but a crosssection shows z_{i} to be quite similar to that in the 1km simulation (Figure 2b). Surface heat flux, H, peaks right around 1400 LST, and then steadily declines; its value is quite sensitive to the local land use type, but is very similar between the two configurations because of their common land surface model and surface characterization (Figure 3). We compute a convective velocity scale using z_{i} = 1400 m and an estimate of the area and timeaveraged surface sensible heat flux, H = 250 W m^{–2}. These yield a convective velocity scale of ${w}_{\ast}\equiv \text{\hspace{0.17em}\hspace{0.17em}}{({\scriptscriptstyle \frac{\mathit{\text{gH}}{z}_{i}}{\rho {c}_{p}{\theta}_{0}}})}^{1/3}\approx 2.2\text{m}{\text{s}}^{1}$ .
The mesoscale and LES configurations do show some differences in PBL wind profiles near the emission source. The LES exhibits lower surface wind speeds, though both are around 6 m s^{–1} above about 50 m AGL (Figure 4). Vertical wind shear does have the potential to impact convective PBL structure, because it is a source of TKE. A standard measure of the validity of the convective CBL assumption is –z_{i}/L, where $L\text{\hspace{0.17em}\hspace{0.17em}}\equiv \text{\hspace{0.17em}\hspace{0.17em}}{\scriptscriptstyle \raisebox{1ex}{$u{}_{\ast}{}^{3}\rho {c}_{p}{\theta}_{0}$}\!\left/ \!\raisebox{1ex}{$(\mathit{\text{gkH}})$}\right.}$ is the MoninObukhov length, with u_{*} being the surface friction velocity and k the von Karman constant (about 0.4). For positive heat flux –L can be considered the maximum height at which TKE production by surfacegenerated wind shear exceeds buoyant production. For the mesoscale configuration u_{*} averages about 0.6 m s^{–1} in the vicinity of the power plant, while it is closer to 0.5 m s^{–1} in the LES. Thus we get –L≈42 m and z_{i}/L≈33 in the LES, while they are about 73 m and 19 in the mesoscale configuration. These low magnitudes of –L are indeed considered to represent a convective as opposed to neutral PBL (e.g., Gryning et al. 2007), though they are near the threshold at which sheardriven phenomena like convective rolls appear (Lemone et al. 1973).
We proceed on the hypothesis that mixedlayer convective scaling with z_{i} = 1400 m and w_{*} = 2.2 m s^{–1} will adequately describe plume turbulent transport in the bulk of the PBL for both simulations, with the caveat that other scaling would likely need to be used to describe the details of the surface layer, where the wind shear is significant (e.g., MoninObukhov similarity theory). Thus, the findings here should be scalable to other cases of convective PBL plume transport characterized by these parameters. Combined with a value of U = 6 m s^{–1}, we obtain the characteristic downwind eddy turnover distance of ${x}_{0}\equiv \text{\hspace{0.17em}\hspace{0.17em}}{\scriptscriptstyle \frac{U{z}_{i}}{{w}_{\ast}}}\approx 3.8\text{km}$ . We will normalize our results by this distance to enable our findings to be generalized to other cases.
Figure 5 compares the power plant plume from the 1km mesoscale model simulation to the 111m LES plume at a model height of 95 m AGL. The concentration normalization corresponds to the expected concentration enhancement if the source is advected through a wellmixed box of height and width z_{i}, and is about 13 ppmv. The resolved wind fields and plume direction for the mesoscale and LES plumes are similar, and the region with the lowest contoured tracer levels (normalized values about 0.1) cover roughly the same horizontal extent. Clearly though the LES plume has about five times greater maximum concentrations within a distance of about x_{0} of the source, and the maximum is narrower and more elongated in the alongwind direction.
Many of the differences between the plume structures are due to the difference in horizontal spatial resolution. Figure 6 is at the same height and time as Figure 5, but shows the 333m mesoscale plume from Grid 4. In contrast to the 1km mesoscale plume, the width of the 333m mesoscale plume is now comparable to that of the 111m LES plume, and the maximum concentrations are much closer in magnitude. Figure 6 also shows that temporal averaging of the 333m mesoscale plume every ten minutes from 1700–1900 UTC causes little change in the maximum concentration in comparison to the 1800 UTC model output, although the averaging procedure does cause the concentrations beyond about x_{0} to decrease in magnitude.
Additional differences between the two plumes cannot be explained by spatial resolution. Neither version of the 333m mesoscale plume in Figure 6 captures the highly elongated shape of the LES plume maximum in the alongwind direction. Figure 7 shows the same comparison as in Figure 5, but at a height of 206 m AGL (i.e., eighth model vertical level). The 1km plume maximum concentration has decreased by about a factor of two relative to the 95m value, but otherwise the plume structure looks almost identical at the two levels. This feature is shared with the 333m mesoscale plume (Figure 8). The LES plume maximum concentrations, by contrast, are reduced by about a factor of five relative to their 95m AGL values; furthermore, the LES concentrations at this level are nearly zero within about a kilometer downwind of the source.
We now consider crosssections of the integrated plume concentration in the crosswind direction as a function of downwind distance and height. In addition to allowing us to compare our results to the tankmodel LES plumes of Willis and Deardorff (1976), this metric tends to filter out the impact of horizontal model resolution, showing instead the vertical behavior of tracer material as it moves downwind. The averaging is performed perpendicular to and symmetric about a line drawn along the highest plume concentrations, assuming that this corresponds to the mean transport direction; a total crossplume averaging distance of 14 km is used. Figure 9 compares integrated concentrations between the 1km mesoscale domain and the 333m mesoscale domain at 1800 UTC, normalized by the expected wellmixed value of Q/(Uz_{i}) (here, approximately 18 ppmvm). At least within 2x_{0} of the source, the crosssections look quite similar, showing that the horizontal resolution difference has little impact on the nearfield vertical plume distribution of crossplumeintegrated concentration (whereas, in contrast, horizontal resolution did have considerable impact on the nonintegrated concentration field).
A comparison between the 333m mesoscale crosssection with that from the 111m LES is shown in Figure 10. The crosssections were temporally averaged and performed over a 4km crossplume transect because of the small size of the LES domain. These two integrated plumes show a number of similarities. The low concentration contours for both ascend to z_{i} at approximately 0.8x_{0}. Then from about 1.5x_{0} to 2x_{0} a slight concentration maximum appears in the upper part of the PBL in both plumes. Both of these features are similar to those reported in the laboratory tank convection experiments of Willis and Deardorff (1976) (their Figure 7). These resemblances give us confidence that at least largescale plume behavior in both models is realistic.
The most notable difference between the 333m mesoscale and 111m LES plumes in Figure 10 appears within about 0.2z_{i} of the surface and 1.5x_{0} of the source. In the LES plume, normalized integrated concentrations above unity are wholly confined to this zone, and the isopleths are nearly horizontal before the elevated plume maximum begins to develop around x_{0}. In the mesoscale plume, the concentration contours in the nearsource region all resemble parabolas with no nearsurface confinement, and the decrease of nearsurface integrated concentration with downplume distance is much more rapid.
In Figure 11, we again compare the 333m mesoscale and the 111m LES cross sections, but zoomed in more on the nearsource region x < x_{0}. The mesoscale plume concentration isopleths initially ascend at a greater rate than those of the LES plume, but then level off and turn concave down. The LES isopleths ascend much more slowly with downwind distance in the nearsource region, and some isopleths even seem to be concave up. This is at least qualitatively in agreement with the predictions of the edge of the unbounded Gaussian plume without and with the Taylor nearsource correction, as is the fact that at higher levels the mesoscale configurations show substantial concentrations near the source, while in the LES the concentration is displaced downwind of the source at higher levels. In the lowest couple hundred meters, the LES integrated plume concentrations are at least twice those of the mesoscale configuration. Above this level, however, the mesoscale model has higher integrated concentrations than the LES model nearly everywhere within x_{0} of the source.
To confirm that the main differences of integrated plume cross sections are not due to resolution difference, we performed a sensitivity simulation in which Domain 5 was run using the MYNN PBL parameterization, like the other four domains, instead of in LES mode. This would normally not be a reasonable choice of physics for modeling turbulent transport at this scale (the turbulent fluxes being parameterized as subgrid by the PBL scheme would largely be carried by turbulent eddies large enough to be explicitly resolved), but we try it here to isolate the impact of resolution from turbulent physics. When we look at the 111m mesoscale integrated concentration field for this time period (Figure 12), we clearly see that the additional enhanced horizontal resolution has little impact on integrated plume concentrations or the parabolic shape of the isopleths.
We now compare the WRF model output to the Legendre farsource and blended solutions, shown in Figure 13. Deriving the nearfield ascent rate for both forms of the Legendre solution is more involved than for the Gaussian plume case, but an analysis valid for z << z_{i} is found in Supplementary Material Text S3. It is shown that for the farsource plume the isopleths of integrated concentration ascend as z = Cx, while for the blended plume the isopleths ascend as z = Cx^{2}. This accounts for the concaveup isopleths in the blended solution as well as the slopes becoming horizontal near the source, instead of remaining constant as in the farsource solution. It is also shown in Supplementary Material Text S3 why in the farsource solution the maximum height of integrated concentration contours all fall along the same line radiating from the origin. The close resemblances between the mesoscale and LES plumes with the farsource and blended solutions, respectively, suggest that the lack of inclusion of the Taylor nearsource effect is the likely source of the integrated concentration differences we see between the mesoscale and LES plumes, at least in the near field.
4. Implications for concentration footprints
Above we used LES to assess the potential concentration biases of a mesoscale plume for a particular case. Using LES to simulate plume dispersion for an entire city over time domains of months to years, however, is computationally prohibitive. Given some idealized meteorological assumptions, we can use the 3D concentration formula (8) to help determine the expected systematic biases in the use of mesoscalederived footprints for the flux inversion problem.
Stated generally, the relation between a 3D concentration field at space/time coordinates (x, y, z, t) and a surface point source at (x_{s}, y_{s}, 0, t_{s}), t > t_{s}, is:
where Q has units of kg s^{–1} and f ( ), the forward transport function, has units of s m^{–3}. We now assume that: 1) the mean and turbulent atmospheric flow and thermodynamic fields are horizontally homogeneous, so the transport function only depends on relative position between receptor and source; 2) the atmospheric fields and Q are considered to be stationary, at least between times t and t_{s}; and 3) we make the slender plume approximation, so t – t_{s} is linearly related to the downwind relative position, x – x_{s} via the constant mean wind speed U. Under these conditions, the concentration field is also stationary, and is linearly related to the source strength by:
Here the transport function is just proportional to the steadystate concentration field for a surface point source given a representation of turbulent diffusion, which can correspond to either the theoretical solutions found above or to numerical simulations of surface point source emissions.
By linear superposition, the concentration resulting from surface sources from all values of x_{s} and y_{s} (to which only upwind positions contribute according to our assumptions) is given by the integral:
Here, $\overline{Q}$ has units of kg m^{–2} s^{–1}, and is the source strength per unit area over the element bounded by dx_{s} dy_{s}. So, if we assume that all of the tracer originates from a surface source, the concentration field at a given height z is a convolution of the surface source fields with a surface footprint function – i.e., the Green’s function of the surface source inversion operator for a receptor at height z – when cast as a function of x_{s} and y_{s}. Source inversion involves inferring Q given c and the footprint function.
Thus, either our theoretical or simulated plume solutions directly give us the relation between concentration and source strength for our idealized meteorology. Any systematic biases in these solutions will translate to systematic biases in inferred source strength for given concentration measurements. In particular, if we assume that the footprint derived from the blended solution of (8) corresponds to the LES footprint (henceforth f_{LES}), and that both of these are the ‘true’ footprint, then the used of the farsource/mesoscale footprint, f_{meso}, for a single point source with a given receptor concentration will lead to an inferred source strength bias of f_{LES}/f_{meso}.
Figure 14 shows the upwind footprint functions based on (8) for the farsource (left) and blended formulas (bottom) at a height of 100 m given the meteorology of the WRF simulations (about z/z_{i} = 0.07). We focus on the concentration at 100 m since that is a common sampling altitude for the INFLUX tower network (Miles et al. 2017). We see the more conical shape of the blended plume/footprint versus the rounder shape for the farsource plume, which we saw in the topview plots of the LES plumes vs. the mesoscale model plumes (Figure 5). We also see a horizontal gap between the source and concentration contours in the blended plume, increasing with height (Figure 15), which does not appear in the farsource plumes; this again corresponds to the LES vs. mesoscale plume behavior. We see that the analytical farsource and blended solutions have similar maximum concentrations, which shows that the Taylor nearsource effect is not primarily responsible for the differences in maximum plume concentrations between the 111m LES and 1km mesoscale plumes in WRF. This is consistent with the previous discussion of the 333m mesoscale plumes shown in Figure 6, which suggested that horizontal resolution is the primary cause of the difference in overall maximum plume concentrations between the 1km mesoscale and 111m LES plumes.
If we use the ratio of the 100m blended to farsource footprints in Figure 14 to estimate mesoscale footprint biases, the broader and more diffuse farsource footprint suggests that for most upwind distances the mesoscale inferred source strength close to the plume centerline will be positively biased, while away from the centerline it will be negatively biased (Figure 16). Sufficiently near the source (within a few tenths of x_{0}), however, the virtual absence of the blended footprint would imply extremely large negative source biases from the use of the farsource footprint. More plausibly, since the blended solution suggests that surface sources are nearly undetectable in their immediate vicinity by elevated receptors, the attempt to use an elevated receptor with a farsource footprint to assign a strength to a nearby surface source should be viewed with caution.
To aid in interpreting these footprint functions, and to distinguish the different behavior in the nearsource and farsource regimes, it is helpful to return to (8) and consider separately the cases of $\widehat{x}>>2$ and $\widehat{x}<<2$ , whose separation corresponds to a dimensional distance of about 3 km in this study. Recall that because of the separable nature of the solution, the farsource and blended solutions only differ in the form of the scaled horizontal coordinate, $q(\widehat{x})$ , so the blended, ‘true’ footprint can be derived from the farsource footprint through a horizontal transformation. For $\widehat{x}>>2$ , we have $q(\widehat{x})=\widehat{x}$ for the far source solution, and $q(\widehat{x})\text{\hspace{0.17em}\hspace{0.17em}}\approx \widehat{x}1$ for the blended solution. At these distances the blended plume solution is simply that of the farsource solution at a nondimensional distance of unity, or dimensional distance of cx_{0}, closer to the source. Using this ‘displaced difference’ to produce a more accurate depiction of farfield plume behavior is essentially equivalent to the ‘handicapped time’ usage in the study of Raupach (1989). Note that while for $\widehat{x}>>2$ the blended diffusivity function asymptotes to the farsource diffusivity function, the blended $q(\widehat{x})$ continues to lag the farsource $q(\widehat{x})$ by one unit – essentially this is because the blended plume ‘falls behind’ the farsource plume in the nearfield, and can never fully catch up.
Additionally for $\widehat{x}>>2$ , it can be seen that all the Legendre modes except the zerothorder uniform mode have become small – so both the blended and farsource plumes are wellmixed in the vertical, and have equivalent crosswindintegrated concentrations. However, the plume is wider in the crosswind direction at $\widehat{x}$ than at $\widehat{x}1$ , showing that the farsource/mesoscale plume is wider than the corresponding blended/LES plume. It can be shown that for large $\widehat{x},{f}_{\mathit{\text{LES}}}/{f}_{\mathit{\text{meso}}}$ is greater than unity, leading to positive inferred source biases from the use of the farsource plume, whenever ${(\mathit{\text{y}}{y}_{s})}^{2}<{\sigma}_{\mathit{\text{yLES}}}^{2}$ , where σ_{yLES} is the standard deviation of the blended/LES plume in the y direction at $\widehat{x}$ . For ${(\mathit{\text{y}}{y}_{s})}^{2}<{\sigma}_{\mathit{\text{yLES}}}^{2}$ , we have f_{LES}/f_{meso} < 1 and thus negative mesoscale inferred source biases. Both signs of bias approach unity as downwind distance approaches infinity.
For the case of $\widehat{x}<<2$ , we can again use a horizontal transformation to convert the farsource/mesoscale plume into the equivalent blended/LES plume, but this time the transformation is $q(\widehat{x})\text{\hspace{0.17em}\hspace{0.17em}}=\text{\hspace{0.17em}\hspace{0.17em}}\widehat{x}1\text{\hspace{0.17em}\hspace{0.17em}}+\text{\hspace{0.17em}\hspace{0.17em}}{e}^{\widehat{x}}$ , which has the approximate form of ${\widehat{x}}^{2}\text{\hspace{0.17em}}/\text{\hspace{0.17em}}2$ for small $\widehat{x}$ . Additionally, the plumes are no longer vertically wellmixed. The blended/LES plume is both narrower and more closely confined to the surface than the farsource/mesoscale plume. Near the plume axis f_{LES}/f_{meso} is much greater than unity near the surface. Above the surface, the dependencies get more complicated. As can be seen in Figure 15, for elevated regions there is a zone in immediate proximity to the source horizontal location for which the ratio f_{LES}/f_{meso} approaches zero.
5. Discussion
The modeling and theoretical analyses suggest that the 1km mesoscale modeling configuration does a reasonable job at reproducing plume structures on the scale of the Indianapolis urban area. For scales of a few kilometers or less, plume structures are subject to model errors, some due to model resolution, and some inherent to the representation of turbulent physics. The first can be addressed by ensuring the model grid spacing is substantially less than the source/receptor distance and the scale of the features to resolve; some aspects of the second could be addressed through an appropriate scale transformation.
When used directly, the LES and analytical results suggest that mesoscale plumes, even for sufficient horizontal resolution, diffuse too quickly with downwind distance both in the horizontal and, until the plume is wellmixed, in the vertical. This causes mesoscale concentrations to be too low along the plume centerline near the surface at large downwind distances, but too large away from the plume centerline and aloft, especially at small downwind distances. Under conditions of stationary and horizontally homogeneous meteorology, we could reasonably suppose that predictions of plume dispersion made with a mesoscale model over an extended period of record would also be too diffuse. This would lead to overestimates of emission strength downwind of sources in inversion systems, especially at low levels in the nearsource region, while underestimates of emission strength would occur far away from the plume centerline, and at elevated receptors in the nearsource region.
As noted previously, running WRFLES for an entire city over time domains of months is currently computationally prohibitive, and a wide variety of meteorological and turbulent conditions can be expected in such a period. However, the differences between the mesoscale and LES plume concentrations are predominantly in the near field. Since it only takes tens of minutes for emissions to advect from a source through the nearfield domain, the assumption of stationary micrometeorological conditions during this time is often reasonable. In this case, mesoscale simulations could simply be corrected in the near field from stationary solutions that contain the correct nearfield behavior and micrometeorology. LES could be conducted for a range of stability conditions, and these corrections then be applied to the nearfield of the mesoscale plume dispersion. Or, the stationary and horizontally homogeneous plume solutions in (13) can be applied, as long as for each time of record the appropriate micrometeorological conditions (e.g., wind speed and direction) are used.
A demonstration of what might happen in this exercise using actual meteorology is shown in Figures 17 and 18. These figures show the prediction of monthlyintegrated plumes at a height of 100 m from Hestia industrial point sources in the Indianapolis region, using the actual frequencies of wind speeds and direction of the climatological September wind rose of Indianapolis International Airport (KIND). (For the sake of this exercise, boundary layer height and the convective velocity scale were assumed constant.) Figure 17 uses the farsource/mesoscale analytical solution to model the plume transport, while Figure 18 uses the blended/LES analytical solution. We can see for both plots that the highest concentrations are in nearly circular patches in the vicinity of the point sources, but higher values are found in the mesoscale simulation (Figure 17). This counterintuitive result is explained by the aggregation over various wind directions of the plumes shown in Figure 14 for somewhat elevated receptors near the source, with nearzero concentration values in the LES simulation, but nearmaximum concentrations in the mesoscale simulation. The ‘concentration gap’ zone in the blended/LES plume is maximized at high wind speeds, because x_{0} is proportional to wind speed. Conversely, the plume aggregation generates wider areas with medium concentration values (yellow contours in Figure 18) in the LES. The medium concentration values in Figures 17 and 18 correspond to locations too far from a point source to be affected by the elevated concentration gap in the LES. At these locations the blended/LES field has systematically higher concentrations than the farsource/mesoscale field, because the mesoscale plume is more diffuse and has lower concentrations at its peak.
Overall, the mesoscale plume entails high concentration values near the source location and broader plume structures with a rapid attenuation of concentrations from the source location to the edges of the plume (Figure 14 – left panel). By contrast, LES plumes show nearzero concentrations above 100 m until about 1 km downwind of the source, with higher concentrations in the center of narrower plumes extending further downwind (Figure 14 – right panel). Previous studies (Lauvaux et al. 2016) noted that observed CO_{2} peaks are systematically underestimated in modeled mesoscale plumes. This result suggests that differences between mesoscale and LES plumes as presented here (i.e., broader plumes with lower maximum) may correspond to the mismatch between modeled and observed concentrations. Lauvaux et al. 2016 also noted that nearzero concentrations were too frequent in the model results, which would be consistent with our results if many mesoscale model plumes were erroneously too diffuse to appear above the lowest concentration bin. Future work will address this question over longer time scales to quantify systematic model errors and their impact on topdown inversion results.
6. Summary and conclusions
We showed that a 1km nested mesoscale model baseline simulation did a reasonable job at predicting downwind concentrations from the Harding St. Power Plant for a dry daytime boundary layer based on meteorology from 28 September 2013, using a corresponding turbulenceresolving 111m LES as a benchmark, and suitable spatial and temporal averaging of turbulent eddies. Some of the existing differences, such as substantially lower maximum concentrations, could be attributed to insufficient horizontal resolution in the baseline simulation. However, sensitivity tests though showed that some of the differences, especially within a few km of the source, were inherent to the use of ensembleaveraged turbulent physics, instead of explicit turbulent eddies, in the vertical turbulent transport prediction of the baseline configuration. The use of analytical solutions supported the conclusion that our LES configuration could capture the essential features of the short temporal/spatial scales in the Taylor (1921) diffusion framework, while the PBL parameterization in the mesoscale configuration could not. Among the consequences are that the mesoscale plume rises too rapidly with downwind distance (crosswiseintegrated concentrations ascend roughly linearly instead of quadratically near the origin) and spreads too rapidly in the lateral direction for a given downwind distance. Some of the consequences for using mesoscale transport models in source inversion systems have been discussed, in particular consequences for using elevated receptors in the nearsource region; however, simple transformations to the output of a mesoscale transport model might prove effective in removing most of the systematic biases.
For the particular case of power plant emissions, we should note that, unlike the case in this simulation, the emission source is elevated, which would have an impact on the nearsource concentrations. The Harding St. power plant emissions occur at multiple stack heights, including 80 m and 172 m (Gurney, personal communication). While this would impact the specific heights of largest plume concentrations in the near field, the reduced vertical dispersion in the near field away from the emission height would still be present, and our theoretical equations can be extended to this situation. Positive buoyancy of the power plant plumes would be another factor affecting plume dispersion; here we only consider neutrally buoyant plume behavior.
In the future we would like to extend our modeling work to a variety of meteorological conditions and associated PBL characteristics (e.g., neutral, weakly stable, nocturnal stable) which can be compared with the extensive tower and other observations taken during the INFLUX campaign. As long as the meteorology is approximately horizontally homogeneous, the symmetry between forward model transport and backward footprint Green’s functions is preserved. For alternate types of boundary layers, the strict applicability of these theoretical solutions and the LES model configuration may be reduced (because the largest eddies may not be well resolved and/or may not extend across the whole depth of the boundary layer). However, it is still expected that the methods used here will prove useful in modified form in other conditions. The LES and analytical models can complement each other; an appropriate analytical model can be applied to an extended data record of cases that would be too computationally extensive for use of LES in its entirety; however, LES is capable of predicting large eddy transport in meteorological conditions for which analytical solutions only approximately apply. Performing plume LES over a suitable set of scalable meteorological conditions could be used to derive nonanalytical footprint functions that could then be applied to an extended set of meteorological cases on a ‘best match’ basis, without the need of explicitly performing LES over the whole period of record. WRFChemLES can also be used to model plume dispersion for individual cases where horizontally heterogeneous meteorology is important in model/observational comparisons, but this application is beyond the scope of the current manuscript.
Data Accessibility Statement
All model results can be made available upon request from the corresponding author.
Supplemental Files
The supplemental files for this article can be found as follows:

Text S1. Nearsource Gaussian plume slopes.
File Type:
File Size:
Download
DOI: https://doi.org/10.1525/elementa.247.s1 
Text S2. Legendre solution to convective PBL plume.
File Type:
File Size:
Download
DOI: https://doi.org/10.1525/elementa.247.s1 
Text S3. Slopes of crossplume integrated Legendre solutions near surface source.
File Type:
File Size:
Download
DOI: https://doi.org/10.1525/elementa.247.s1 
Table S1. Approximate height of lowest model layers above ground level.
File Type:
File Size:
Download
DOI: https://doi.org/10.1525/elementa.247.s1