Exploration of the impact of nearby sources on urban atmospheric inversions using large eddy simulation

In contrast to the ‘bottom-up’ approach to quantifying CO₂ 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¹ or the high resolution inventory Hestia (Gurney et al. 2012), the ‘top-down’ approach involves measuring CO₂ 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 top-down 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 one-third 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 top-down 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 model-observation pairs used in the inversion. A number of issues may lead to systematic errors in the transport model used in the top-down 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 state-of-the-science 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 grid-cell-averaged (i.e., ‘mean’, or model-resolved) atmospheric fields such as momentum, temperature, and moisture, as well as grid-cell-averaged turbulent variances and fluxes. PBL turbulence is parameterized. WRF can also be run as a large eddy simulation (WRF-LES) (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 chemistry-adapted version of WRF, WRF-Chem (Grell et al. 2005), that can simulate the chemical reactions and transport of various species, including CO₂.

During the INdianapolis FLUX experiment (INFLUX), twelve (12) communication towers in and around Indianapolis, IN were instrumented with tower-based 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, building-level CO₂ emissions from the Hestia product, and tower-measured atmospheric CO₂ mole fractions (Lauvaux et al. 2016). The mesoscale inversion system produced adjustments to the Hestia CO₂ emissions of about 20% over eight months (mostly fall and winter), possibly due to an under-estimation of CO₂ 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 “near-field” 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 one-dimensional (vertical only) PBL parameterization predicts the grid-cell-averaged statistics of turbulent variances and fluxes through variance closure assumptions; tendencies of model-resolved 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₂ and other species with a gradient-diffusion 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 point-source 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 mid-PBL 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 WRF-Chem, with a baseline model grid spacing of 1 km. Our procedure is to drive the WRF-Chem baseline configuration for the case of 28 September 2013 and predict the transport of CO₂ 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 WRF-Chem and WRF-LES 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 point-source diffusion equation for an idealized daytime PBL, and which approximate the differences in near-field 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 over-predicts vertical diffusion in the near-source 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 1-km mesoscale simulation and the 111-m LES are due to horizontal resolution, or due to inherent limitations of the K_Z–based approach for near-field 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.

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₂ plume released in the INFLUX domain, focusing on the near-field effects. ‘Near-field’ 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 chemistry-adapted version of WRF, WRF-Chem (Grell et al. 2005), for which CO₂ was treated as a passive tracer with a surface emission source. Model-resolved 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 WRF-Chem, 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 WRF-Chem 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 near-field 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 re-classified based on the 2006 version of the 30-m National Land Cover Database (NLCD) (Homer et al. 2015; Sarmiento et al. 2017). Time-dependent 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 point-source plumes in typical daytime PBL conditions using the baseline mesoscale configuration of WRF-Chem, rather than to verify the observed plume structure for this particular day. For simplicity, WRF-Chem assumed that the power plant emission occurred at the surface (in reality, the lowest half-layer height, at approximately 7 m above ground level (AGL)).

A five-domain one-way 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 Mellor-Yamada-Nakanishi-Niino (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 1-km domain (Domain 3) considered the ‘baseline’ mesoscale simulation. For the reference simulation, the 333-m 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 111-m domain (Domain 5), WRF-LES is used (Moeng et al. 2007), with the 1.5-order TKE subgrid closure of Deardorff (1980).

Because all nesting is one-way, 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 1-km mesoscale configuration versus a 111-m LES configuration. However, a separate simulation was run using a 111-m 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.

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 S-1). 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.

WRF-LES is an Eulerian grid implementation of LES that can be used as a nested grid in both real-case and idealized simulations. One advantage of WRF-LES 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. WRF-LES 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 111-m horizontal grid spacing, which is consistent with many studies using WRF-LES (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 non-resolved 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 non-periodic, it takes a finite time for LES to spin-up realistic turbulent eddies, either from a non-turbulent initial state or from non-turbulent 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 spin-up problem.

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 steady-state for a few hours.

The simulation period was 1200 UTC 28 Sep 2013 through 0000 UTC 29 Sep 2013 (0700–1900 LST). No spin-up 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 32-km 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).

Current urban inversion systems attempt to reproduce time-averaged 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 time-averaged, tower-based 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 large-scale 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) 1-km mesoscale output.

In the dry daytime convective boundary layer, the vertical scale of the largest turbulent eddies is on the order of the mixed-layer height z_i, and their velocity scale is on the order of the Deardorff convective velocity scale, $w∗≡(gHziθ0)1/3$ (Deardorff 1970; Young 1988). Here g is gravitational acceleration, θ₀ 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 well-mixed in the convective CBL, is z_i/w_*, and is usually on the order of 10³ 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 $x0≡Uziw∗$ 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 111-m 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 near-field plume from that of the turbulent physics scheme.

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 steady-state 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 well-known 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 $Ky=U2ddxσy2$ and $Kz=U2ddxσz2$⁠. 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 (tracer-following) 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 energy-containing (largest) eddies. In (3) and (4), V_Y and V_Z are the velocity components in the horizontal cross-wind 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 x-dependent 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 $σV2$ and $σW2$⁠, while the squared displacements on the LHS of (3) become $σY2$ and $σZ2$⁠. These $σY2$ and $σZ2$ expressions can be used in the plume equation (2) provided crosswind-dependencies of the velocity variances are neglected. The result is the long-time (or equivalently far-source) formula for the CPS steady-state plume under gradient diffusion:

with $KY=TLYσV2$ and $KZ=TLZσW2$ (e.g., Wyngaard and Weil 1991). In the daytime convective PBL, T_LZ scales as z_i/w_*, while $σW2$ scales as w_*²; 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/far-source 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/near-source plume behavior. For t << T_L the mean squared displacement of particles is given by:

Conceptually, in the long-time 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 short-time 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 near-source relations (valid for x << UT_L):

where formally $KY=xσV2/U$ and $KZ=xσW2/U$⁠. Thus the near-source plume equation has a reduced effective eddy diffusivity compared to the far-source plume, but one that increases with downwind distance towards its far-source value at x ≈ UT_L.

The most fundamental difference in near-source (close to the source with respect to x_o) eddy diffusivity in comparison to far-source (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 near-source region, and in fact will systematically overestimate near-source vertical diffusion, assuming that it correctly predicts the far-source value of $KZ=TLZσW2$⁠. WRF-LES, however, does have the capability to capture the correct near-source 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 y-integrated form of these equations for simplicity. The details are found in the supplementary material (Supplementary Material Text S-1), but it can be shown that for the far-source version of the Gaussian plume the concentration surfaces near the edge of the plume ascend nearly vertically near the source. For the near-source plume, however, the concentration surfaces possess a finite slope near the source (in fact, they approach a linear slope).

Though the Gaussian plume equation (2) is usually a good description of far-source 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 far-source behavior in the vertical is not characterized by continuous spreading, but by the PBL becoming increasingly well-mixed. 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, $σW2$⁠, 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 near-source 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 near-source and far-source behavior.

For the 2-D 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 near-source effects can be included in f(x). This separability assumption leads to a Sturm-Liouville problem that may be solved given a particular form of g(z). However, the resultant procedures do not lead to readily analyzed closed-form expressions except in specialized cases.

Here, we consider the full 3-D case of (1). We assume that K_y has the form $Ky(x)=Kyfarf(x/UTL)≡Kyfar[1−exp(−xUTL)]$⁠, which smoothly merges the Taylor near-source and far-source diffusivity limits. We assume that K_z has this same x-dependence but also has a z-dependence in the separable form $Kz(x,z)=4Kzfar maxz^(1−z^)f(x^)$⁠, where $z^≡z/zi$ and $x^≡x/UTL$⁠. 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 mid-PBL peak and smaller values both at the surface and at z_i.

Given these functions, a unique regular analytic solution exists to the steady-state CPS equation, as shown by Nieuwstadt (1980) and Otte and Wyngaard (1996). More details appear in the supplementary material (Supplementary Material Text S-2), 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₀; thus $x^=x/cx0$⁠. Finally $q(x^)$ is the integral of the Taylor diffusivity modification function $f(x^)$⁠:

The y-integrated concentration equation is:

One convenience of this form is that we can test the impact of neglecting the near-source Taylor behavior in the Legendre analytical solutions by assigning $f(x^)=1$ and $q(x^)=x^$ instead of using (9). Henceforth we will call the solution that neglects the near-source Taylor behavior as the ‘far-source’ 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 far-source and blended Legendre solutions, respectively.

We begin by comparing the general PBL structures as simulated by the 1-km mesoscale and 111-m LES configurations around the 1800 UTC = 1300 LST reference time. The MYNN PBL parameterization of the 1-km configuration predicts a mean boundary layer height of 1400 m in the domain, which is consistent with the well-defined top of the model mixed layer as seen in a potential temperature cross-section 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 cross-section shows z_i to be quite similar to that in the 1-km 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 time-averaged surface sensible heat flux, H = 250 W m^–2. These yield a convective velocity scale of $w∗≡(gHziρcpθ0)1/3≈2.2 m 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≡−uρ∗3cpθ0(gkH)$ is the Monin-Obukhov 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 surface-generated 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 shear-driven phenomena like convective rolls appear (Lemone et al. 1973).

We proceed on the hypothesis that mixed-layer 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., Monin-Obukhov 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 $x0≡Uziw∗≈3.8 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 1-km mesoscale model simulation to the 111-m 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 well-mixed 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₀ of the source, and the maximum is narrower and more elongated in the along-wind 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 333-m mesoscale plume from Grid 4. In contrast to the 1-km mesoscale plume, the width of the 333-m mesoscale plume is now comparable to that of the 111-m LES plume, and the maximum concentrations are much closer in magnitude. Figure 6 also shows that temporal averaging of the 333-m 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₀ to decrease in magnitude.

Additional differences between the two plumes cannot be explained by spatial resolution. Neither version of the 333-m mesoscale plume in Figure 6 captures the highly elongated shape of the LES plume maximum in the along-wind 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 1-km plume maximum concentration has decreased by about a factor of two relative to the 95-m value, but otherwise the plume structure looks almost identical at the two levels. This feature is shared with the 333-m mesoscale plume (Figure 8). The LES plume maximum concentrations, by contrast, are reduced by about a factor of five relative to their 95-m AGL values; furthermore, the LES concentrations at this level are nearly zero within about a kilometer downwind of the source.

We now consider cross-sections 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 tank-model 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 cross-plume averaging distance of 14 km is used. Figure 9 compares integrated concentrations between the 1-km mesoscale domain and the 333-m mesoscale domain at 1800 UTC, normalized by the expected well-mixed value of Q/(Uz_i) (here, approximately 18 ppmv-m). At least within 2x₀ of the source, the cross-sections look quite similar, showing that the horizontal resolution difference has little impact on the near-field vertical plume distribution of cross-plume-integrated concentration (whereas, in contrast, horizontal resolution did have considerable impact on the non-integrated concentration field).

A comparison between the 333-m mesoscale cross-section with that from the 111-m LES is shown in Figure 10. The cross-sections were temporally averaged and performed over a 4-km cross-plume 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₀. Then from about 1.5x₀ to 2x₀ 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 large-scale plume behavior in both models is realistic.

The most notable difference between the 333-m mesoscale and 111-m LES plumes in Figure 10 appears within about 0.2z_i of the surface and 1.5x₀ 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₀. In the mesoscale plume, the concentration contours in the near-source region all resemble parabolas with no near-surface confinement, and the decrease of near-surface integrated concentration with down-plume distance is much more rapid.

In Figure 11, we again compare the 333-m mesoscale and the 111-m LES cross sections, but zoomed in more on the near-source region x < x₀. 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 near-source 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 near-source 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₀ 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 111-m 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 far-source and blended solutions, shown in Figure 13. Deriving the near-field 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 S-3. It is shown that for the far-source plume the isopleths of integrated concentration ascend as z = Cx, while for the blended plume the isopleths ascend as z = Cx². This accounts for the concave-up isopleths in the blended solution as well as the slopes becoming horizontal near the source, instead of remaining constant as in the far-source solution. It is also shown in Supplementary Material Text S-3 why in the far-source 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 far-source and blended solutions, respectively, suggest that the lack of inclusion of the Taylor near-source effect is the likely source of the integrated concentration differences we see between the mesoscale and LES plumes, at least in the near field.

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 mesoscale-derived 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 steady-state 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, $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 far-source/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 far-source (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 far-source plume, which we saw in the top-view 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 far-source plumes; this again corresponds to the LES vs. mesoscale plume behavior. We see that the analytical far-source and blended solutions have similar maximum concentrations, which shows that the Taylor near-source effect is not primarily responsible for the differences in maximum plume concentrations between the 111-m LES and 1-km mesoscale plumes in WRF. This is consistent with the previous discussion of the 333-m 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 1-km mesoscale and 111-m LES plumes.

If we use the ratio of the 100-m blended to far-source footprints in Figure 14 to estimate mesoscale footprint biases, the broader and more diffuse far-source 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₀), however, the virtual absence of the blended footprint would imply extremely large negative source biases from the use of the far-source 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 far-source 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 near-source and far-source regimes, it is helpful to return to (8) and consider separately the cases of $x^>>2$ and $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 far-source and blended solutions only differ in the form of the scaled horizontal coordinate, $q(x^)$⁠, so the blended, ‘true’ footprint can be derived from the far-source footprint through a horizontal transformation. For $x^>>2$⁠, we have $q(x^)=x^$ for the far source solution, and $q(x^)≈x^−1$ for the blended solution. At these distances the blended plume solution is simply that of the far-source solution at a non-dimensional distance of unity, or dimensional distance of cx₀, closer to the source. Using this ‘displaced difference’ to produce a more accurate depiction of far-field plume behavior is essentially equivalent to the ‘handicapped time’ usage in the study of Raupach (1989). Note that while for $x^>>2$ the blended diffusivity function asymptotes to the far-source diffusivity function, the blended $q(x^)$ continues to lag the far-source $q(x^)$ by one unit – essentially this is because the blended plume ‘falls behind’ the far-source plume in the near-field, and can never fully catch up.

Additionally for $x^>>2$⁠, it can be seen that all the Legendre modes except the zeroth-order uniform mode have become small – so both the blended and far-source plumes are well-mixed in the vertical, and have equivalent crosswind-integrated concentrations. However, the plume is wider in the cross-wind direction at $x^$ than at $x^−1$⁠, showing that the far-source/mesoscale plume is wider than the corresponding blended/LES plume. It can be shown that for large $x^,fLES/fmeso$ is greater than unity, leading to positive inferred source biases from the use of the far-source plume, whenever $(y −ys)2<σyLES2$⁠, where σ_yLES is the standard deviation of the blended/LES plume in the y direction at $x^$⁠. For $(y −ys)2<σyLES2$⁠, 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 $x^<<2$⁠, we can again use a horizontal transformation to convert the far-source/mesoscale plume into the equivalent blended/LES plume, but this time the transformation is $q(x^)=x^−1+e−x^$⁠, which has the approximate form of $x^2/2$ for small $x^$⁠. Additionally, the plumes are no longer vertically well-mixed. The blended/LES plume is both narrower and more closely confined to the surface than the far-source/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.

The modeling and theoretical analyses suggest that the 1-km 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 well-mixed, 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 near-source region, while underestimates of emission strength would occur far away from the plume centerline, and at elevated receptors in the near-source region.

As noted previously, running WRF-LES 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 near-field 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 near-field behavior and micrometeorology. LES could be conducted for a range of stability conditions, and these corrections then be applied to the near-field 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 monthly-integrated 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 far-source/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 counter-intuitive 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 near-zero concentration values in the LES simulation, but near-maximum concentrations in the mesoscale simulation. The ‘concentration gap’ zone in the blended/LES plume is maximized at high wind speeds, because x₀ 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 far-source/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 near-zero 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₂ peaks are systematically under-estimated 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 near-zero 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 top-down inversion results.

We showed that a 1-km 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 turbulence-resolving 111-m 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 ensemble-averaged 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 (crosswise-integrated 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 near-source 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 near-source 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 non-analytical 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. WRF-Chem-LES 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.

All model results can be made available upon request from the corresponding author.

The Hestia inventory product for Indianapolis and the hourly CAMD emission data were provided to us by Kevin Gurney and Jian-Hua Huang at Arizona State University. Daniel Sarmiento adapted and provided a high-resolution urban fraction land cover product for the Indianapolis region that was used in these model simulations.

This research is supported by the National Institute of Standards and Technology (Project # 70NANB10H245), by the National Aeronautics and Space Administration (Project #NNX15A137G), and by the National Oceanic and Atmospheric Administration (Project # NA13OAR4310076).

The authors have no competing interests to declare.

• Contributed to conception and design: BG, TL, KD

• Contributed to modeling results: BG, AD, TL

• Contributed to analysis and interpretation of modeling results: BG, TL, KD

• Drafted and/or revised article: BG, TL, KD, AD