In this study, the Flexible Particle (FLEXPART)-WRF, a Lagrangian particle dispersion model, is employed to simulate pollutant dispersion in and near the Lehigh Gap, a gap in a prominent ridgeline in eastern Pennsylvania. FLEXPART-WRF is used to evaluate the diagnostic value of the ventilation index (VI), an index that describes the potential for smoke or other pollutants to ventilate away from a source, for indicating dispersion potential in complex terrain. Little is known about the effectiveness of the ventilation index in diagnosing dispersion potential in complex terrain. The modeling approach used in this study is to release a dense cloud of particles across a portion of the model domain and evaluate particle behavior and VI diagnostic value in areas of the domain with differing terrain characteristics. Although both horizontal and vertical dispersion are examined, the study focuses primarily on horizontal dispersion, assessed quantitatively by calculating horizontal residence time (HRT) within a 1-km-radius circle surrounding the particle release point. Analysis of HRT across the domain reveals horizontal dispersion patterns that are influenced by the ridgeline and the Lehigh Gap. Comparison of VI and HRT in different areas of the domain reveals a robust relationship windward of the ridgeline and a weak relationship leeward of the ridgeline and in the vicinity of the Lehigh Gap. The results of this study suggest that VI users should consider whether they are windward or leeward of topographic features, and highlight the need for an alternative metric that better takes into account the influence of the terrain on dispersion.
The ventilation index (VI), also known as the clearing index, is an index that describes the potential for smoke or other pollutants to ventilate away from a source (Ferguson 2001). VI is a standard part of fire weather forecasts issued by National Weather Service (NWS) forecast offices across the United States and is used by fire managers as a diagnostic tool to assess the potential for dispersion of smoke from wildland fires. In this study, the term “dispersion” is used to describe the movement of a pollutant away from its release point, through some combination of mean transport and turbulent diffusion, and should not be interpreted as implying turbulent processes alone. Although there is no standard VI formulation, it is typically computed as the product of a mixing height and a wind speed averaged within the mixed layer (i.e., transport wind speed).
VI is closely related to the box model, a type of dispersion model that is used as a screening tool for estimating pollutant concentration within a valley or basin (Goodrick et al. 2013). In a box model, pollutants are emitted uniformly across the bottom of a box with horizontal dimensions similar to the dimensions of a valley or basin, and pollutant concentration varies as smoke is emitted from the source and transported out of the valley or basin. Both VI and box models are built on the assumption that smoke is instantaneously mixed vertically through the depth of the mixed layer, and is transported laterally away from the source by a uniformly mixed wind profile above a surface-layer logarithmic wind profile. Although this assumption is most valid when applied on time scales on the order of hours to a day, VI is often applied on hourly or shorter time scales by fire weather forecasters and fire managers. Furthermore, the assumption that the pollutant is distributed vertically through the mixed layer may be violated near a point or area source, and may not be met until some distance downwind of the source (Goodrick et al. 2013).
In practice, a mixing height and mixed-layer-mean wind speed, measured at a single point in time and space, is used to compute VI (or drive a box model), and the conditions at the point are assumed to be representative of conditions over a broad area (on the order of 100 km) and over a length of time during which the instantaneous mixing assumption is likely to be valid (on the order of hours to a day) (Goodrick et al. 2013). An assumption of horizontal homogeneity and stationarity of mixed-layer depth and wind speed is not unreasonable over flat terrain during portions of the day (Stull 1988), but is more questionable in complex terrain where topography influences surface heating and induces wind regimes at a variety of spatial and temporal scales (Gorski and Farnsworth 2000). In complex terrain, one or more of the following may influence winds: terrain-forced flows, such as flow through gaps or around obstacles; and diurnal mountain winds, including slope winds, along-valley winds, cross-valley winds, and mountain-plain winds (Whiteman 2000). It is unclear whether VI, subject to the aforementioned limiting assumptions (instantaneous vertical mixing, time scales on the order of hours to a day, horizontal homogeneity and stationarity of wind and mixed-layer height), can be applied in complex terrain and over time scales of 1 h or less.
To date, the only effort to evaluate VI we are aware of is a study that was part of the Southeastern Smoke Project (Achtemeier et al. 2010). During the 5-yr project, fire activity data, concentrations of particulate matter with diameters ≤2.5 μm (PM2.5), and weather data were collected during 56 prescribed fires. The resulting comparison of smoke concentration and VI showed weak to no correspondence, although Goodrick et al. (2013) point out that the Achtemeier et al. (2010) comparison was made within a short distance of the burn site, and the relationship between VI and pollutant concentration may improve as one moves farther away from the vicinity of the source and the aforementioned assumptions about vertical mixing are progressively better satisfied. To our knowledge, validation of VI in complex terrain, either via a field campaign or a model study, has not been performed previously. As VI is a standard element of NWS fire weather forecasts across the United States, and is used by fire managers across the United States including areas of complex terrain (Ferguson 2001), a study of the relationship, if any, between VI and pollutant (e.g., smoke) concentration in and around complex terrain is potentially of great value.
In this study, the Flexible Particle Weather Research and Forecasting (FLEXPART-WRF) Model (Brioude et al. 2013), a Lagrangian particle dispersion model, is employed to simulate pollutant dispersion within a model domain centered on the Lehigh Gap, a water gap in eastern Pennsylvania. The Lehigh Gap is located where the Lehigh River flows through an opening in Blue Mountain, a prominent ridgeline along the eastern edge of the Appalachian Mountains. The modest nature of the terrain features in eastern Pennsylvaniaserves as a preliminary test of the diagnostic value of VI in complex terrain. The rationale for this initial study is that if VI performs poorly in an area like the Lehigh Gap, with relatively small terrain obstructions [O(400) m elevation change from ridge to valley], and a relatively simple terrain configuration, it is unlikely that it would perform better in areas such as the Intermountain West, with higher mountains, deeper valleys, and more complex configurations of basins, valleys, and gaps. The modeling approach used in this study is to release a dense cloud of particles across the model domain, compute residence time (RT) within a given horizontal or vertical distance from each particle’s release point, and compare VI and RT in areas of the domain with different terrain characteristics. RT is used in this study as a metric to quantify particle dispersion: longer (shorter) RT is equated with weaker (stronger) dispersion. Calculation of RT thus facilitates an examination of the statistical relationship (if any) between the diagnostic dispersion index, VI, and particle dispersion.
The setting (Lehigh Gap) and date (14 April 2013) chosen for this study correspond to that of a prescribed fire event at the Lehigh Gap Nature Center, a 756-acre nature preserve located immediately west of the gap. Approximately 400 acres of the nature preserve is a restored U.S. Environmental Protection Agency (EPA) Superfund site at the location of a former zinc smelter (in operation until 1980). The Superfund program was established by an act of the U.S. Congress to clean up sites contaminated by industrial activity (a description of the Superfund program may be found at https://www.epa.gov/superfund/what-superfund, and a description of the restoration of the Lehigh Gap site may be found at https://cumulis.epa.gov/supercpad/SiteProfiles/index.cfm?fuseaction=second.redevelop&id=0300624). It is important to keep in mind that the release of particles in this study is not meant to emulate the emission of smoke from a fire. Instead, the particle release strategy utilized in this study is designed to help clarify the impact of complex terrain on particle dispersion in the absence of heat or moisture release from a fire, and to facilitate an evaluation of the VI as a diagnostic tool for indicating dispersion potential in complex terrain. The choice to use the setting and date of a prescribed fire for this study is rooted in part from the desire to examine smoke dispersion in complex terrain during actual prescribed fire weather conditions, and to take advantage of the available meteorological measurements.
The remainder of this manuscript is organized as follows. A description of the Lehigh Gap study area and a brief summary of the 14 April 2013 prescribed fire event are provided in section 2, a summary of the numerical models used in this study and the experiment design are presented in section 3, results and discussion of the experiments are presented in section 4, and the paper is concluded in section 5.
2. Location description and fire event summary
As stated previously, the setting for the FLEXPART-WRF study is the Lehigh Gap in eastern Pennsylvania. The study area is displayed in Fig. 1, beginning with a broad view of surface elevation across the northeastern United States (Fig. 1a), indicating the location of the Lehigh Gap in eastern Pennsylvania, along the eastern edge of the Appalachian Mountains. Zooming in to an approximately 50-km2 area of eastern Pennsylvania (Fig. 1b), a segment of the 410-km-long southwest–northeast-oriented ridgeline known as Blue Mountain becomes visible. Also visible is a portion of the Lehigh River, a 175-km-long tributary of the Delaware River, intersecting Blue Mountain at the center of the image. Zooming in finally to the FLEXPART-WRF domain (Fig. 1c), it becomes clear that the Lehigh Gap exists where the Lehigh River flows approximately north–south through a gap in Blue Mountain. The ridgeline in the vicinity of the Lehigh Gap is approximately 420–450 m above mean sea level with the surface elevation at river level approximately 115 m above mean sea level where it flows through the gap.
On 14 April 2013, 9.6 acres of grassland along the lower slope of Blue Mountain, about 1.5 km northwest of the Lehigh Gap, was burned over the 1 h 48 min period beginning at 1312 eastern daylight time (EDT = UTC − 4) (burn unit indicated by red triangle in Fig. 1c). Vegetation at the burn site consisted primarily of prairie grasses, planted as part of the Superfund restoration, with scattered defoliated trees and areas of exposed shale also present. Outside of the Superfund site itself, vegetation in the area primarily consisted of oak and chestnut trees. The 14 April 2013 prescribed fire was part of a monitoring study at the Lehigh Gap Nature Center designed to inform future fire management at the site. Uncertainty regarding the release of heavy metal content and the transport of smoke emissions during prescribed burning of the prairie grass was the primary motivation for the study. Three instrumented flux towers collected data prior to, during, and following the burn: two in situ towers located within the burn block (not shown) and a control tower located approximately 0.5 km west of the burn block (red circle in Fig. 1c). The 10-m-tall control tower measured three-dimensional wind, air temperature, relative humidity, net radiation, pressure, and soil temperature at 10-Hz frequency. Data collected at this tower are used in section 3a to assess the simulated meteorological variables used to drive FLEXPART-WRF.
3. Model description and numerical experiment design
a. ARPS model configuration, parameterization, and evaluation
In this study, the meteorological fields required by FLEXPART-WRF originate from simulations conducted using the Advanced Regional Prediction System (ARPS), version 5.2.7 (Xue et al. 2000, 2001). ARPS is a three-dimensional, compressible, nonhydrostatic atmospheric model with a terrain-following coordinate system. ARPS has been applied across a range of spatial scales, from studies of in-canopy turbulence, using grid spacing as fine as O(1) m (Dupont and Brunet 2008), to studies of mesoscale and synoptic-scale phenomena, utilizing grid spacing of O(1–10) km (e.g., Xue et al. 2001; Parker and Johnson 2004; Michioka and Chow 2008). Furthermore, ARPS has been applied in complex terrain in a number of previous studies, including Smith and Skyllingstad (2005), Michioka and Chow (2008), and Kiefer and Zhong (2011). The choice to use ARPS in this study is based primarily on the established utility of ARPS by the authors for simulating flows in complex terrain as well as the implementation of terrain-shadowing effects in ARPS by Colette et al. (2003).
In this study, ARPS is run in one-way nested mode with four domains, beginning with an outermost domain over the northeastern United States, with 2.7-km horizontal grid spacing, and nesting down through intermediate 900- and 300-m horizontal grid spacing domains, to the innermost domain, with 100-m horizontal grid spacing. Initial and boundary conditions for the outermost domain are provided by the 12-km North American Mesoscale Forecast System (NAM) model (Rogers et al. 2009). For the outer three domains, U.S. Geological Survey (USGS) terrain and land-cover data (30 arc sec; http://www.caps.ou.edu/ARPS/arpsdown.html) interpolated to the model grid are utilized, whereas for the innermost domain, Shuttle Radar Topography Mission (SRTM) terrain data (3 arc sec; https://dds.cr.usgs.gov/srtm/version2_1/SRTM3/) interpolated to the model grid are implemented. SRTM terrain data are used in place of the USGS data as the latter data source was found to misrepresent the ridgeline elevation immediately west of the gap (not shown). In all domains, a 1.5-order subgrid-scale turbulence closure scheme with a prognostic equation for the turbulent kinetic energy is utilized, as well as a land surface and vegetation model based on Noilhan and Planton (1989) and Pleim and Xiu (1995) and radiation physics following Chou (1990, 1992) and Chou and Suarez (1994); in all but the innermost domain, a nonlocal turbulence parameterization based on Sun and Chang (1986) is also applied. Effects of topographic shading on radiative fluxes are accounted for following Colette et al. (2003). All four domains are initialized at 0800 EDT and run for 12 h; simulated variables are output from the innermost grid at 1-min frequency and are subsequently used to drive FLEXPART-WRF. However, the first 30 min of the innermost grid simulation are considered spinup time and are not utilized for any other purpose. The decision to use a 30-min spinup time is based on a consideration of the spatial scale of simulated topographic flows in and around the Lehigh Gap; examples of previous ARPS studies with comparable spinup time periods include Xue et al. (2003) (30-min spinup time) and Snook et al. (2015) (1-h spinup time).
Before proceeding, an evaluation of the innermost domain simulation is conducted in order to lend confidence to the use of meteorological output from ARPS to drive FLEXPART-WRF and evaluate VI in complex terrain. The purpose of this comparison is not to evaluate the accuracy of the ARPS simulation per se, but rather to consider if the ARPS-simulated surface winds and temperatures near the burn site are representative of conditions for which a prescribed fire could be conducted at the Lehigh Gap. Although a number of variables were measured at the control tower, this analysis is limited to wind speed, wind direction, and air temperature interpolated from the 3 and 10 m above ground level (AGL) instrument heights to 5 m AGL, the height of the lowest atmospheric grid level in ARPS. In Fig. 2, time series of 10-min-mean variables are compared between the tower observations and the ARPS-simulated values extracted at the nearest grid point to the tower. Comparing the observed and simulated time series of wind speed and direction first (Figs. 2a,b), note that the model underestimates wind speed, exhibits a westerly directional bias (outside of the period from 1200 to 1230 EDT, when a transient, likely terrain-induced eddy perturbs the ARPS-simulated wind direction), and exhibits weaker subhourly fluctuations than the observations. Comparing the observed and simulated temperature time series (Fig. 2c), bias is found to vary within ±0.6°C. Averaged over the 4-h time series, 5 m AGL wind speed, direction, and temperature observed (simulated) at the tower are 5.2 m s−1 (3.7 m s−1), 302.8° (275.6°), and 10.4°C (10.4°C), respectively. Efforts to improve the comparison of simulated and observed time series include the aforementioned replacement of USGS terrain data with SRTM data and modification of roughness lengths in the ARPS land surface model to more closely match vegetation types observed in and around the burn unit. The persistence of model errors following implementation of these modifications suggests that higher-resolution terrain and land-use data may be necessary for an improved comparison (left to future work).
Despite mean differences of −1.6 m s−1 for wind speed and −27.3° for wind direction, the weather conditions simulated by ARPS meet the weather criteria specified in the burn plan for the Lehigh Gap Nature Center on 14 April 2013 [S. Henry (burn boss) 2012, personal communication]: effective wind speed (combining the effect of midflame wind speed and terrain slope on fire spread) between 2 and 9 mph (0.9 and 4.0 m s−1), wind direction between west and north (270° and 360°), with northwest preferred, and air temperature between 30° and 90°F (−1.1° and 32.2°C). For reference, effective wind speed is a fire spread parameter used in fire behavior models (including those applied operationally) that is determined using either empirical or physics-based equations; for this brief model evaluation exercise, an empirical lookup table is used to estimate this parameter (see National Wildfire Coordinating Group 2017, p. 103). The assertion that ARPS-simulated wind speed falls within the range of effective wind speed specified in the burn plan is based on a log-wind adjustment of the 5 m AGL wind speed down to an approximate midflame height (1 m AGL), yielding a midflame wind speed of 2.5 m s−1. Estimating an average slope angle of 30% between the ridgeline and river, and assuming long-grass fuels, the effective wind speed for fire spreading upslope is about 3 m s−1; this value is well within the range specified in the burn plan. Thus, it is concluded that the weather conditions simulated by ARPS fall within the range of conditions for which burns are conducted at the Lehigh Gap Nature Center, allowing us to proceed with confidence to the FLEXPART-WRF experiments. This assessment is encouraging: given the complex nature of the topography in and around the Lehigh Gap, errors in the simulated wind speed and direction might plausibly have yielded simulated weather conditions outside the acceptable range specified in the burn plan (e.g., wind direction less than 270°, or effective wind speed greater than 4 m s−1).
b. FLEXPART-WRF Model configuration, parameterization, and experiment design
FLEXPART-WRF is a version of the Lagrangian particle dispersion model FLEXPART (Stohl et al. 2005), originally developed for use with global weather models, that has been adapted for use with the WRF Model (Powers et al. 2017); this study uses FLEXPART-WRF, version 3.3. Although we refer to the model by its proper name FLEXPART-WRF, the model in this study is actually driven by output from the ARPS model. The reader is reminded that the choice to use ARPS in this study is based primarily on the established utility of ARPS for simulating flows in complex terrain.
FLEXPART-WRF uses a Lagrangian reference frame, tracking individual particles and updating particle position at every time step based on a mean wind component and a turbulent wind component. The mean wind component is computed from the horizontal and vertical wind components ingested into FLEXPART-WRF from the meteorological model, with users given the option to use instantaneous winds, time-averaged winds, or instantaneous winds with a diagnostic calculation of vertical velocity; since ARPS does not output time-averaged winds (as in WRF versions 3.3 or higher), the instantaneous wind option is utilized in this study. For the turbulent wind component, users have two options: FLEXPART-WRF can parameterize this wind component internally using the Hanna turbulence scheme (Hanna 1982), as in the original FLEXPART model, or, FLEXPART-WRF can parameterize it using the turbulent kinetic energy (TKE) variable ingested from the meteorological model. As Brioude et al. (2013) concluded that the ingested TKE option violates the well-mixed criterion (a homogeneously mixed tracer in the PBL should stay homogeneously mixed), we follow their recommendation and use the internal Hanna turbulence formulation instead, which is a function of three planetary boundary layer (PBL) parameters: mixed-layer height, friction velocity, and surface heat flux. As with the turbulent wind component, users can choose to have FLEXPART-WRF compute the PBL parameters internally, or have FLEXPART-WRF use the PBL parameters ingested from the meteorological model; in this study, we choose the latter option to ensure consistency in mixing height calculations between ARPS and FLEXPART-WRF.
The FLEXPART-WRF simulation is initialized at 0830 EDT (30 min after the ARPS initialization), and is run for a total of 7 h, ending at 1530 EDT. The FLEXPART-WRF grid is consistent with the ARPS grid, with the same horizontal and vertical grid spacing, domain dimensions, and lower-left- and upper-right-corner coordinates. ARPS variables are read in once every minute from netCDF files that are generated from ARPS2FLEX, a routine that transforms the output from ARPS into the suite of variables that FLEXPART-WRF is designed to read in from WRF output. Three-dimensional particle position is updated every 10 s, and written out to ASCII files at the same interval; particle concentration is also computed but is not presented here. A total of 250 000 particles are released during the 7-h FLEXPART-WRF simulation, across an approximately 7.5 km × 7.5 km zone centered northwest of the Lehigh Gap (Fig. 1c). Particles are released uniformly across the zone, at a constant particle release rate, with initial particle location specified randomly. Particle diameter is restricted to the range of 2.5 μm or smaller (i.e., PM2.5), with a mean particle diameter of 0.5 μm, and a particle density of 1.36 g cm−3 (representative of oak; Levin et al. 2010). The use of nonzero particle density implies deposition processes are engaged in FLEXPART-WRF; however, in this study only dry deposition is applied. It is important to keep in mind that the particle release strategy employed in this study is idealized in nature and is not intended to represent actual emissions during the 14 April 2013 prescribed fire.
c. Analysis methodology
Analysis of the relationship between VI and particle dispersion in the Lehigh Gap vicinity proceeds in three stages. First, an overview of particle behavior is presented by analyzing spatial patterns of particle dispersion across the area encompassing the Lehigh Gap and Blue Mountain. Second, a statistical comparison of VI and RT is performed in different areas of the domain to assess the strength of the relationship between VI and particle dispersion, and the sensitivity of this relationship to the underlying terrain. Third, an analysis of meteorological conditions across the domain is presented, including an evaluation of the components of VI, mixing height and transport wind speed, and alternative dispersion metrics.
As a preliminary step in the analysis of particle dispersion, RT is computed for each of the 250 000 particles released during the FLEXPART-WRF simulation. For each particle, horizontal RT (HRT) is computed within a 1-km-radius circle centered on the particle release point, and vertical RT (VRT) is computed within a 100-m-deep layer above the particle release level; additional radii and depths are computed but are not shown here. In this study, longer (shorter) RT is equated with weaker (stronger) particle dispersion, and the terms are used interchangeably. Note that although VI does not indicate the preferred direction of dispersion (horizontal versus vertical), we feel it is a worthwhile endeavor to consider the relationship of VI to horizontal and vertical dispersion separately, in order to examine if VI is better correlated with one or the other. VI is computed as the product of a mixing height (MH), defined in this study as the height above the surface where virtual potential temperature first exceeds the lowest grid-level value, plus a 1-K perturbation (added to ensure that the bases of weak inversions are not misidentified as the mixing height), and a transport wind speed, defined here as the numerical average wind speed within the mixed layer (WSML). It is worth pointing out that the methods chosen for calculating MH and WSML are consistent with current NWS forecast practices (National Weather Service 2018). However, the 1 K increment used in this study, determined to be the smallest value that consistently bypassed weak inversions, yielding spatially consistent MH, is larger than the increment that is currently implemented at the NWS (0.5 K) (National Weather Service 2018).
4. Results and discussion
a. Particle behavior overview
Figure 3 shows the initial particle position of all particles released during the simulation, with the HRT and VRT of each particle (Figs. 3a and 3b, respectively) indicated by marker color. Note that the particle plots in Fig. 3 depict spatially inhomogeneous patterns of HRT and VRT. A comparison of RT and contoured surface elevation (dark gray contours in Fig. 3) suggests a relationship between RT and release location, relative to Blue Mountain and the Lehigh Gap, that is different for horizontal and vertical motion. Regarding horizontal motion (Fig. 3a), HRT is generally longer for particles released on the leeward (i.e., southeast facing) slope of Blue Mountain, compared to particles released on the windward (i.e., northwest facing) slope, and is also longer in an area of the hilly terrain 2–5 km northwest of the Lehigh Gap (recall the westerly mean wind direction in Fig. 2b). The longer HRT in the latter area appears to be related to stagnant near-surface flow simulated by ARPS during the morning (not shown). Overall, HRT is longest in the area 0–2 km southwest of the Lehigh Gap, on the leeward slope of Blue Mountain, and is shortest in the area 0–1 km northeast of the Lehigh Gap, on the windward slope of Blue Mountain. Regarding vertical motion (Fig. 3b), VRT is generally shorter for particles released along the ridgeline than for particles released in the surrounding lower elevation areas, but also appears to be shorter for particles released 1–3 km south of the Lehigh Gap.
Taken together, the HRT and VRT analyses in Fig. 3 depict horizontal and vertical dispersion patterns that are influenced by the terrain features of Blue Mountain and the Lehigh Gap. Horizontal dispersion is weaker on the leeward slope of Blue Mountain, especially southwest of the gap, and stronger on the windward slope of Blue Mountain; vertical dispersion is strongest along the ridgeline of Blue Mountain. Examination of terrain-induced meteorological phenomena and the relationship of particle dispersion to the three-dimensional wind field is analyzed in section 4c. Last, for the purpose of assessing VI diagnostic value in areas of the domain with relatively good ventilation and areas with relatively poor ventilation, and both in situ and upwind of Blue Mountain and the Lehigh Gap, five 1 km × 1 km analysis zones are defined [Fig. 3: upwind (ZUW), windward ridge (ZWR), leeward ridge (ZLR), windward gap (ZWG), and leeward gap (ZLG)].
With the analysis zones defined, a comparison of PM2.5 “plumes” from each zone is presented in Fig. 4. The term “plume” is used with care here, as the “plume” for each analysis zone is visualized by isolating all particles released from that zone; recall that particles are released throughout the 7.5 km × 7.5 km release zone depicted in Fig. 1c and Fig. 3. In Fig. 4, a general pattern of particle transport by the prevailing west-to-northwest wind is evident. Furthermore, as the particles are mixed to greater heights with time (Fig. 4, left panels), the axes of all five “plumes” reorient from generally west-northwest–east-southeast to predominately northwest–southeast (Fig. 4, center panels). In addition to a depiction of the general behavior of the particles in the left and center panels, the zoomed-in view in the right panels provides evidence of areas of stronger and weaker dispersion in the vicinity of Blue Mountain and the Lehigh Gap. Comparing ZWR and ZLR (green and dark blue colors, respectively), considerably larger particle counts are noted in the latter zone, consistent with weaker particle dispersion in the lee of Blue Mountain. Furthermore, relatively high particle counts are seen in the southeastern two-thirds of ZLG; note the dearth of particles in the portion of the zone that extends northwest of the ridgeline. Consistent with Fig. 3, a pattern of stronger (weaker) dispersion on the windward (leeward) slope of Blue Mountain is depicted in Fig. 4, with weakest dispersion noted southwest of the Lehigh Gap.
b. Ensemble-mean RT analysis
To assess the strength of the relationship between VI and dispersion, mean VI is now compared to ensemble-mean HRT and VRT. For each analysis zone and for each 10-min time block, mean VI is computed by spatially averaging VI across the zone, and temporally averaging the resulting quantity over the 10 values within the time block (1-min ARPS output), and ensemble-mean RT is computed by identifying the particles released during the time block, and averaging RT for that subset of particles. The outcome of this procedure is a series of scatterplots for mean VI and ensemble-mean RT; the process yields a total of 42 data points per panel (420-min simulation, 10-min time blocks), with each data point representing approximately 100 particles. As a method of assessing the strength of the VI–RT relationship, linear regression lines are fit to the data points, and corresponding statistics are computed: normalized slope S, computed as the slope divided by the plot aspect ratio , and coefficient of determination R2.
Beginning with HRT in Fig. 5, R2 values from 0.79 to –0.80 and S values from −0.63 to −0.67 are found in the two zones northwest of the Lehigh Gap and windward (i.e., northwest) of the ridgeline (ZUW and ZWR). The S values indicate a 63%–67% reduction in ensemble-mean HRT with a doubling of mean VI. Ensemble-mean HRT in both zones varies between about 2.5 and 8.5 min, with standard deviation values within each bin generally on the order of 2 min or less, outside of a few points corresponding to the very beginning of the simulation. Analysis of HRT and VI in an additional zone west-southwest of ZUW (not shown) reveals steeper slopes (i.e., greater sensitivity of HRT to changes in VI) in the aforementioned hilly area northwest of the Lehigh Gap where HRT is longer during the morning (Fig. 3a) because of stagnant near-surface flow. Nevertheless, the mean VI–ensemble-mean HRT relationship within ZUW and ZWR is arguably robust (R2 ~ 0.8 and S ~ −0.65). It is worth noting that this relationship was examined using 10-min time blocks, suggesting that despite the aforementioned limiting assumptions (e.g., horizontal homogeneity, stationarity) likely being violated, VI has diagnostic value even on time scales of an hour or less.
In ZWG, despite some similarity in appearance to the scatterplots in the other zones windward of the ridgeline (ZUW and ZWR), in terms of degree of scatter and standard deviation values, the R2 and S values of 0.34 and −0.15, respectively, indicate a weaker relationship between mean VI and HRT. It is also worth noting that ensemble-mean HRT varies across a narrower range than in ZUW and ZWR. In the two zones leeward (i.e., southeast) of the ridgeline (ZLR and ZLG), the R2 values of 0.13 and 0.26 are smaller than though comparable to those of ZWG, whereas the S values of −0.32 and −0.48 actually indicate steeper regression slopes than in ZWG. However, a visual comparison of the scatterplots in ZLR and ZLG shows greater spread among points and larger standard deviation values within each bin, compared to the scatterplots in the zones windward of the ridgeline (ZUW, ZWR, and ZWG). To summarize, the scatterplots indicate a higher degree of temporal and spatial variability of HRT within the zones leeward of the ridgeline, relative to zones windward of the ridgeline, and the smaller R2 values suggest a weaker relationship between mean VI and ensemble-mean HRT in the zones leeward of the ridgeline as well.
In an effort to compare the diagnostic value of VI for indicating horizontal dispersion potential versus vertical dispersion potential, the scatterplots of mean VI versus ensemble-mean HRT are now compared to scatterplots of mean VI versus ensemble-mean VRT (cf. Figs. 5 and 6). In the three zones windward of the ridgeline (ZUW, ZWR, and ZWG), the spread between points and the standard deviation values are consistently larger for VRT, and in ZUW and ZWR, R2 and S values are also smaller for VRT, indicating a weaker overall relationship between VI and VRT than between VI and HRT. In the two zones leeward of the ridgeline (ZLR and ZLG), the spread between points and standard deviation values are similar or smaller for VRT, whereas R2 values are comparable for HRT and VRT; however, S values are smaller for VRT than HRT. Taken as a whole, these results suggest that in ZUW and ZWR, where the mean VI–ensemble-mean RT relationship is stronger for horizontal motion than vertical motion, VI may have greater diagnostic value for horizontal dispersion than vertical dispersion. However, in ZLR and ZLG, where the mean VI–ensemble-mean RT relationship is somewhat weaker for horizontal motion than vertical motion, VI may have marginally weaker diagnostic value for horizontal dispersion compared to vertical dispersion. The diagnostic value of VI for horizontal versus vertical dispersion potential is a topic that merits further study; however, in the interest of brevity we choose to focus on horizontal dispersion in the remainder of the manuscript, and leave additional analysis of vertical dispersion to future work.
c. Meteorological analysis
An analysis of meteorological conditions in and surrounding the Lehigh Gap can provide insight into the relationship between VI and particle dispersion revealed in section 4b. This analysis proceeds in three stages: first, HRT is compared to the components of VI, MH, and WSML, in order to provide some context for the relationship of VI and horizontal dispersion evidenced in Fig. 5; second, vertical profiles of 1-h-mean virtual potential temperature and wind speed are examined to provide insight into the characteristics and temporal evolution of the PBL affecting both VI calculation and particle dispersion; third, plan-view maps of VI and its components are examined in order to assess the spatial variability of the PBL and the impact of the terrain on PBL flows, and thus particle dispersion.
The analysis begins with an examination of HRT and its relationship to the two components of VI, MH, and WSML, following the same ensemble-mean analysis procedure outlined at the beginning of section 4b. Motivated by the limited diagnostic value of VI indicated by Fig. 5 in ZWG, ZLR, and ZLG, and hypothesizing that particle dispersion in these zones is influenced more strongly by local winds than the average winds in the mixed layer, we also consider wind speed averaged across alternative, shallower layers. Scatterplots, along with linear regression lines and corresponding statistics (R2 and S), are presented in Fig. 7. Examining mean MH scatterplots first (Fig. 7, second column), it is clear that the relationship of ensemble-mean HRT and mean MH is generally weak and inconsistent across the zones. In ZWG and ZLG, S is negative, indicating decreasing ensemble-mean HRT (i.e., stronger dispersion) with increasing mean MH. This relationship between dispersion and MH is intuitive, as deeper mixing results in particles being mixed to greater heights where they are subject to stronger horizontal wind speeds. However, in ZUW, ZWR, and ZLR, S is positive, indicating increasing ensemble-mean HRT (i.e., weaker dispersion) with increasing mean MH. This relationship is counterintuitive and requires further analysis to explain the underlying processes. However, R2 values are at or below 0.35 in all five zones, even in ZWG and ZLG where the relationship is intuitive, indicating a weak to nonexistent relationship of horizontal dispersion and MH everywhere.
It must be kept in mind that the selection of a case with a consistently well-developed mixed layer throughout the day (see MH x-axis limits in Fig. 7), as opposed to a case with a more typical diurnal PBL evolution, in which a nocturnal boundary layer evolves into a mature mixed layer, reduces confidence in the broader applicability of our results. The synoptic setup during the previous night was not conducive to formation of a stable boundary layer, with weak cold-air advection at 850 hPa promoting the maintenance of a deep residual layer, despite an approaching area of surface high pressure promoting the development of a shallow surface inversion (not shown). The well-mixed nature of the atmosphere may help explain the apparent discrepancy between the finding of this study that VI has some diagnostic value for indicating smoke dispersion potential (see ZUW and ZWR in Fig. 5), and the suggestion by Achtemeier et al. (2010) that it does not.
Proceeding to the other component of VI, WSML (Fig. 7, column three), the S values from −0.49 to −1.20 for mean WSML and HRT are found to be larger than those for mean VI and HRT (S = −0.17 to −0.74; first column in Fig. 7). However, the difference in R2 values between the mean WSML and mean VI scatterplots is inconsistent between zones, with some zones exhibiting larger R2 values for mean WSML (ZWR and ZLR) and some showing larger R2 values for mean VI (ZUW, ZWG, and ZLG). Comparing scatterplots of mean WSML to scatterplots of wind speed averaged in the 0–100 and 0–10 m AGL layers (WS100 and WS10, respectively; compare last three columns in Fig. 7), we find that R2 values in all zones are larger for WS100 and WS10 than WSML. However, the largest percent increase in R2 values over WSML are found in ZWG, ZLR, and ZLG. It is in these three zones that the S values are also consistently larger for WS100 and WS10 than WSML. Thus, the greatest advantage in diagnostic value of the shallow-layer wind speed averages over the mixed-layer wind speed average is found in the three zones where VI performs least well (ZWG, ZLR, and ZLG; Fig. 7, first column).
Next, vertical profiles of ARPS 1-h-mean virtual potential temperature and wind speed, averaged across each analysis zone, are presented in Fig. 8; note that the midpoint of the 1-h period is referred to in the discussion that follows (e.g., 0900 EDT for 0830–0930 EDT). Beginning with virtual potential temperature (Fig. 8, top row), an archetypal mixed layer is found at all hours and across all zones. A mixed layer is evidenced by the 1000–1500-m-deep constant virtual potential temperature layer, sandwiched between a superadiabatic layer at the surface and the stably stratified free atmosphere above. Across all zones except ZWG, MH (again, defined as the height where virtual potential temperature first exceeds the surface value, plus a 1-K perturbation) peaks at about 1500–1700 m AGL between 1000 and 1100 EDT, and decreases somewhat thereafter. The lower MH during the afternoon appears to result from brief episodes of warming in the free atmosphere associated with terrain-induced gravity waves (not shown). In zone ZWG, MH follows a more typical diurnal pattern, gradually increasing during the simulation and peaking in the midafternoon.
Examining mean wind profiles next (Fig. 8; bottom panels), increasing mean wind speed with time is evident in all zones, along with at least some semblance of a constant wind speed layer. In ZUW and ZWR, archetypal mixed-layer-mean wind profiles are present by the early to midafternoon (1200–1500 EDT), with a logarithmic wind profile in approximately the lowest 10% of the mixed layer, and nearly constant wind speed through the depth of the mixed layer. In the remaining three zones (ZLR, ZWG, and ZLG), the mean wind profiles are more complex. The bottom of the constant wind speed layer in these zones is displaced to a height of 250–300 m AGL, and in ZLR, two shear layers are present, a surface-based one and an elevated one located between about 200–250 m AGL.
The values of mean WSML and WS100 printed in each panel offer additional information about the impact of the terrain on the wind field. Note that mean WSML is similar across the five zones, with respect to both magnitude and evolution, varying modestly between 0900 and 1100 EDT (~8 ±1 m s−1), increasing steadily between 1100 and 1300 EDT, then oscillating weakly from 1300 to 1500 EDT (~10.5 ±1 m s−1). In contrast to mean WSML, mean WS100 is not uniform across the zones, in terms of either magnitude or evolution. First, mean WS100 in ZWG is consistently 1–2 m s−1 stronger than in the other zones, reflecting the exposed nature of the zone immediately northeast of the Lehigh Gap. Second, although the evolution of WS100 in the remaining four zones is similar from 0900 to 1100 EDT, the evolution during the remainder of the simulation is noticeably different in the two leeward zones (ZLR and ZLG) compared to the two zones on the opposite side of Blue Mountain (ZUW and ZWR). In ZUW and ZWR, mean WS100 increases from 4.5–5 m s−1 at 1100 EDT to 7.4–8 m s−1 at 1300 EDT and then oscillates weakly thereafter. However, in the two leeward zones, mean WS100 increases from 4.4–4.8 m s−1 at 1100 EDT to 6–7.4 m s−1 at 1300 EDT, but then decreases to 4.3–4.9 m s−1 at 1500 EDT. Taken as a whole, it appears that the mean wind profiles in the vicinity of the Lehigh Gap deviate from the classic boundary layer wind profile (Stull 1988), that consists of a deep constant wind speed layer above a near-surface logarithmic layer, violating an assumption underpinning both VI and the related box dispersion model (Goodrick et al. 2013).
Analysis of the meteorological conditions impacting dispersion now proceeds to an examination of plan-view maps of contoured 7-h-mean VI, MH, WSML, and WS100, along with 7-h-mean horizontal wind vectors (Fig. 9). Note that the four variables presented in Fig. 9 correspond to the x-axis variables in the first four columns of Fig. 7 (WS10 is omitted from Fig. 9 because Fig. 7 showed similar statistics for WS10 and WS100). Examining mean VI first (Fig. 9a), one finds a spatially inhomogeneous field, with mean VI varying by ±15% around the release zone mean (1.4 × 104 m m s−1). Three areas (northwest of ZWR, southwest of ZLR, and southeast of ZWG) exhibit mean VI approximately 15% higher than the release zone mean, and two areas (one larger area southeast of ZLG and one smaller area mostly inside ZWG) exhibit mean VI about 15% lower than the release zone mean.
Examining the components of VI (MH and WSML; Figs. 9b,c), the aforementioned three areas of higher mean VI appear to be spatially correlated with three areas of deeper mean MH, and the area of lower mean VI south of the Lehigh Gap appears to be spatially correlated with an area of weaker WSML (all comparisons made to the release zone mean). Furthermore, MH (Fig. 9b) is shallower over and on the windward slope of Blue Mountain (MH ~ 1200–1400 m AGL), where surface virtual potential temperature is lower, and MH is deeper over approximately the southeast third of the release zone and in the lower elevations north of Blue Mountain (MH ~ 1400–1700 m AGL), where surface virtual potential temperatures are higher; differences in surface virtual potential temperature are primarily due to differences in vegetation cover across the domain (temperature and land surface fields not shown). Although at different scales of topography, the shallower (deeper) MH over the mountain peaks (valleys) noted here is consistent with findings from the Mesoscale Alpine Programme (MAP) reported in Rotach and Zardi (2007); however, it should be noted that the nature of MH over complex terrain is an open research question (De Wekker and Kossman 2015). WSML is weaker south of the Lehigh Gap (WSML ~ 7–8 m s−1) and is stronger along Blue Mountain and in an area northeast of ZUW (WSML ~ 10–11 m s−1), but is generally more spatially homogeneous than MH (cf. Figs. 9b,c); for reference, the release zone mean WSML is 9.3 m s−1. It is worth noting that the areas of shallower MH and stronger WSML along Blue Mountain overlap and appear to largely cancel each other out; similarly, the area of weaker WSML south of the Lehigh Gap appears to partially overlap with a small area of deeper MH south-southeast of ZLG (cf. Figs. 9a–c).
Next, we examine WS100 (Fig. 9d), as Fig. 7 suggests that ensemble-mean HRT is more strongly related to WS100 than WSML, especially in the zones where VI performed most poorly (e.g., ZLG). Examining Fig. 9d, note the correspondence between WS100 and the underlying terrain. Stronger (weaker) WS100 is simulated on the windward (leeward) side of Blue Mountain, and the domainwide strongest (weakest) WS100 is found northeast (southwest) of the Lehigh Gap. Across approximately the northwest half of the release zone, patches of stronger and weaker WS100 are found, with ZUW straddling two adjacent streaks. Last, the inset panel in Fig. 9d shows that near the center of ZLG, the wind field is not only weak in magnitude, but spatially variable in direction as well. The wind field described here is consistent with flows observed in the wakes of terrain and other flow obstacles (Whiteman 2000).
Finally, we briefly consider correlation coefficients between the 2D array of 7-h-mean HRT, constructed by averaging HRT over all particles released in each model grid cell over the 7-h release period, and the 2D arrays plotted in Fig. 9. Correlation coefficients are as follows: VI, −0.13; MH, 0.49; WSML, −0.50; WS100, −0.71. The magnitude and sign of the correlation coefficients are broadly consistent with the R2 and S values shown in the first four columns of Fig. 7, respectively, suggesting that for the well-mixed conditions on 14 April 2013, the variable with the greatest diagnostic value for indicating horizontal dispersion potential in the Lehigh Gap area is WS100 (correlation coefficient for HRT and WS10 is −0.39; not shown).
5. Summary and conclusions
In this study, a Lagrangian particle dispersion model (FLEXPART-WRF) was employed to simulate pollutant dispersion in and near the Lehigh Gap, a water gap in Blue Mountain in eastern Pennsylvania. The modeling approach used in this study involved the release of a dense cloud of particles over the course of a 7-h FLEXPART-WRF simulation, and the evaluation of particle behavior and VI diagnostic value in areas of the domain with differing terrain characteristics. The setting and date chosen for the study (14 April 2013) corresponds to that of a prescribed fire; however, the fire was not represented in ARPS or FLEXPART-WRF, either explicitly or implicitly. The choice to use the setting and date of a prescribed fire for this study is rooted in part from the desire to examine smoke dispersion during actual prescribed fire weather conditions. In other words, the purpose of this study was to examine the impact of complex terrain on particle dispersion and the utility of VI as a diagnostic tool in complex terrain, in the absence of heat or moisture release from a fire. Before proceeding to the analysis of particle dispersion and VI diagnostic value, ARPS-simulated 10-m wind speed and direction, and 2-m temperature, were compared to control tower observations and the weather conditions specified in the burn plan. This exercise verified that the ARPS-simulated weather conditions in the Lehigh Gap area on 14 April 2013 are conditions in which fire managers in the area actually burn.
The primary metrics used to quantify particle dispersion were RT within a 1-km-radius circle surrounding the release point (HRT) and RT within a 100-m layer above the release level (VRT), with shorter (longer) RT implying stronger (weaker) particle dispersion. As a preliminary step, RT-shaded markers were plotted at each particle’s release location and a contoured map of surface elevation was overlaid. The resulting figure depicted horizontal and vertical dispersion patterns that were influenced by the terrain features of Blue Mountain and the Lehigh Gap. For the purpose of assessing VI diagnostic value in areas of the domain with both relatively good and poor ventilation, and both in situ and upwind of Blue Mountain and the Lehigh Gap, five 1 km × 1 km analysis zones were defined. The primary means through which VI was evaluated were scatterplots of ensemble-mean RT and time- and zone-averaged VI, with linear regression lines and associated R2 and slope statistics used to assess the strength of the relationships. The scatterplot analysis showed a robust relationship between VI and RT windward of Blue Mountain and away from the Lehigh Gap, and a weak relationship between VI and RT leeward of Blue Mountain and in the vicinity of the Lehigh Gap.
The remainder of the manuscript was primarily devoted to an assessment of the underlying meteorological conditions driving the spatial patterns of dispersion seen in Fig. 3 and the VI–HRT relationships exhibited across the different analysis zones in Fig. 5. First, ensemble-mean HRT was compared in a series of scatterplots to mean VI, MH, and WSML, as well as alternative WS calculations, WS100 and WS10 In the zones where VI performed least well (ZWG, ZLR, and ZLG), the relationship between HRT and WS100 (or WS10) was consistently stronger than the relationship between HRT and WSML, and much stronger than the relationship between HRT and VI. Examining mean profiles of virtual potential temperature and wind speed next, the mean wind profiles in the vicinity of the Lehigh Gap were shown to deviate from the classic boundary layer wind profile, that consists of a deep constant wind speed layer above a near-surface logarithmic layer, violating an assumption underpinning both VI and the related box dispersion model. Plan view maps of contoured 7-h-mean VI, MH, WSML, and WS100 were subsequently examined and compared to the HRT-shaded particle position plot in Fig. 3a. A spatially heterogeneous pattern of mean VI was noted, and mainly attributed to spatial variations in mean MH; WSML was the most spatially homogeneous of the four quantities examined, with WS100 the most sensitive to the underlying terrain. Analysis of correlation coefficients computed between 2D arrays of gridcell-average HRT and corresponding variables (e.g., mean VI, WSML) revealed that the quantity with the strongest relationship to HRT was WS100 (correlation coefficient of −0.71).
Taken as a whole, the results of this study suggest that VI may have diagnostic value in some areas of complex terrain, especially upwind of terrain obstructions like Blue Mountain, even on time scales of an hour or less. However, VI may fail as an indicator of dispersion potential in the lee of terrain obstructions, and in the vicinity of gaps or passes. In these areas, VI fails to capture important terrain-induced flow modifications because of its reliance on MH and WSML, quantities that exhibit only modest sensitivity to the underlying terrain. It is important to keep in mind that the terrain in the vicinity of Blue Mountain and the Lehigh Gap is comparatively shallow and simple compared to terrain in other regions of the world, such as the western United States or the European Alps. Nevertheless, the results of this study suggest that users of VI should consider whether they are windward or leeward of topographic features, and highlight the need for an alternative metric that better takes into account the influence of the terrain on PBL flows, and thus the potential for pollutant dispersal. Although WS100 was found in this study to correlate best with HRT, highlighting the influence of surface-layer wind speed on pollutant transport, a potentially more useful metric would be one that, like VI, combines a transport wind speed component and a turbulent diffusion component. Development of such a metric, along with analysis of VI and pollutant dispersion in areas with more complex terrain, like the front range of Colorado, is left to future work. Also left to future work is the examination of cases with a greater range of VI, especially cases where the PBL undergoes a more typical transition during the day (nocturnal PBL transitioning to a well-mixed PBL). A follow-up study of VI and smoke dispersion in southwest Colorado, combining field campaigns and numerical modeling, is underway at the time of this writing; the results of the new study are expected to further elucidate the value of VI in complex terrain.
Nevertheless, the current study constitutes a step forward in our understanding of dispersion and the diagnostic value of VI in complex terrain. Furthermore, it is anticipated that the analysis framework developed in this study (e.g., ensemble-mean RT analysis, meteorological analysis) will prove useful in future studies of VI, regardless of terrain complexity.
Support for this research was provided by the USDA Forest Service via Research Joint Venture Agreements 13-JV-11242306-055 and 11-JV-11242306-058. This work was supported partially by the USDA National Institute of Food and Agriculture, Hatch project 1010691. Data from this study are freely available for download at http://eams2.usfs.msu.edu/Pub_Data/outgoing/kiefer/lehighgap/. The color map used in Figs. 3, 8, and 9 was developed by Cynthia Brewer at the Pennsylvania State University (http://colorbrewer2.org/). Finally, comments and suggestions from three anonymous reviewers were helpful in revising the manuscript and are greatly appreciated.
Northern Research Station, U.S. Forest Service, Morgantown, West Virginia.