Assessing Forest Canopy Impacts on Smoke Concentrations Using a Coupled Numerical Model

Abstract

The impact of a forest canopy on smoke concentration is assessed by applying a numerical weather prediction model coupled with a Lagrangian particle dispersion model to two low-intensity wildland (prescribed) fires in the New Jersey Pine Barrens. A comparison with observations indicates that the coupled numerical model can reproduce some of the observed variations in surface smoke concentrations and plume heights. Model sensitivity analyses highlight the effect of the forest canopy on simulated meteorological conditions, smoke concentrations, and plume heights. The forest canopy decreases near-surface wind speed, increases buoyancy, and increases turbulent mixing. Sensitivities to the time of day, plant area density profiles, and fire heat fluxes are documented. Analyses of temporal variations in smoke concentrations indicate that the effect of the transition from a daytime to a nocturnal planetary boundary layer is weaker when sensible heat fluxes from the fires are stronger. The results illustrate the challenges in simulating meteorological conditions and smoke concentrations at scales where interactions between the fire, fuels, and atmosphere are critically important. The study demonstrates the potential for predictive tools to be developed and implemented that could help fire and air-quality managers assess local air-quality impacts during low-intensity wildland fires in forested environments.

1. Introduction

The air-quality impacts of smoke from low-intensity wildland fires undertaken for fuel management are a particular concern for fire and air-quality managers [1]. In Larkin et al. [2], it is demonstrated that smoke from wildland fires can have significant regional impacts, particularly when multiple fires occur in a single airshed during a narrow time window. For wildland fires in the wildland-urban interface (WUI), the potential for smoke to affect nearby communities is almost always a concern [3]. Since the health and safety of fire management personnel can be affected by the timing, magnitude, and duration of smoke concentrations in the immediate vicinity of a fire [4,5], the local effects of smoke must be considered for all wildland fires, even for small low-intensity fires. Predictive smoke dispersion models are important tools for assessing the potential for smoke concentrations to exceed thresholds for safe fire operations and for smoke to affect nearby communities and other smoke-sensitive locations.

Small, low-intensity wildland fires occur and are prescribed in a variety of environments, in fuels that include short grasses, mixed grasses and shrubs, coarse woody debris, and understory vegetation in a forest [6]. While predicting smoke concentrations in any environment can be difficult, predictions for low-intensity fires under a forest canopy can be particularly challenging due to variability within the vegetation structure of the physical and meteorological factors that affect smoke transport and diffusion [7,8,9,10]. In three studies by Kiefer et al. [9,11,12], the development and testing of a fine-scale numerical modeling system, known as ARPS-CANOPY is documented, and it is demonstrated that the modeling system can be used to predict meteorological conditions associated with low-intensity wildland fires. These studies establish the importance of accounting for forest canopy characteristics and their impacts on fire-induced meteorology, and highlight the potential for these factors to affect smoke concentrations. As a point of clarification, a canopy is defined within the context of this study as the entire vegetation layer, including the overstory.

While ARPS-CANOPY simulations of fire-induced meteorology have been investigated for different canopy types and fuel configurations, simulations of smoke concentrations from the model have yet to be evaluated using field measurements. To use ARPS-CANOPY for simulating smoke concentrations, the meteorological model must be coupled with a numerical smoke transport and diffusion model (see [13] for a review of smoke transport and diffusion models). For this study, ARPS-CANOPY is coupled with FLEXPART-WRF [14], a version of the FLEXPART Lagrangian particle dispersion model [15]. The decision to use a Lagrangian dispersion model, rather than an Eulerian model like the Community Multiscale Air Quality (CMAQ) modeling system [16], is based on the expectation that smoke concentration predictions in forested environments are more accurate when at least some dispersion processes are resolved, as they are in Lagrangian models, than when all dispersion processes are parameterized, as they are in Eulerian models. It has been established that FLEXPART-WRF is capable of simulating smoke dispersion at synoptic- to meso-alpha scales [17,18] as well as continent-scale transport [19]. In Solomos et al. [17], it is suggested that WRF-FLEXPART can be applied using meteorological simulations with a grid spacing as small as several meters.

The objectives of this study are to evaluate simulations of smoke concentrations from the coupled ARPS-CANOPY/FLEXPART-WRF (AC/F) model and to assess the sensitivity of AC/F simulations to characteristics of the forest canopy. AC/F simulations are compared to observations of smoke concentrations (PM_2.5) collected during two prescribed fires that occurred in the New Jersey Pine Barrens in 2011 and 2012. Section 2 describes the methods, including descriptions of the characteristics of the prescribed fires (e.g., fuel types, canopy structure), the observed smoke concentrations, and the coupled numerical model. Section 3 consists of an overview of the simulated and observed smoke concentrations followed by an assessment of the sensitivity of smoke concentrations, plume heights, and associated meteorological variables in the AC/F simulations to characteristics of the forest canopy. Section 4 presents conclusions and discusses the implications of this study for smoke management decision-making.

2. Methods

2.1. Prescribed Fire Cases and Field Measurements

Smoke concentrations associated with two low-intensity wildland fires are investigated in this study. The fires were components of a Joint Fire Science Program field project [20,21] that entailed the collection of meteorological and smoke-concentration observations from prescribed fires in the New Jersey Pine Barrens. The fires occurred on 20 March 2011 (the 2011 case) and 6 March 2012 (the 2012 case) and were managed by the New Jersey Forest Fire Service. The 2011 case, the meteorological characteristics of which were presented in [11] and [21], consisted of a 107 ha burn block centered on 39.8726° N latitude and 74.5013° W longitude. The 2012 case, which was previously documented in [22], consisted of a 97 ha burn block centered on 39.9141° N latitude and 74.6033° W longitude. The 2011 burn block has canopy vegetation consisting of pitch pine (Pinus rigida Mill.) and mixed oak (Quercus spp.) with maximum canopy heights between 15 and 18 m, and understory vegetation consisting of blueberry (Vaccinium spp.), huckleberry (Gaylussacia spp.), and scrub oak. The 2012 burn block has canopy vegetation consisting of mixed oak and scattered pitch and shortleaf pines (P. echinata Mill.) with maximum canopy heights between 20 and 23 m; blueberry and huckleberry dominate the understory vegetation. Average understory fuel loadings for the 2011 and 2012 burn blocks during the fires were 1485 g m⁻² and 1104 g m⁻², respectively. The 2011 case was initiated with a single line fire (Figure 1a) while the 2012 case was initiated with a series of roughly parallel line fires. The line fires in both cases were allowed to propagate against the prevailing wind (backing fires) and burning was confined to surface fuels except for brief episodes of crowning in the 2011 case [21]. Approximate line fire mean spread rates were 1.50 m min⁻¹ and 0.33 m min⁻¹ for the 2011 and 2012 cases, respectively. Estimated fire intensities for the 2011 and 2012 cases, based on the formulation of [23], were 325 kW m⁻¹ and 52 kW m⁻¹, respectively [22].

For both cases, a network of meteorological towers and surface monitoring sites was constructed within and in the vicinity of the burn blocks. Towers were placed in the interior and near the perimeter of each burn block to provide spatial coverage of near-surface atmospheric conditions in the vicinity of the fires. A control tower was also placed to the northwest and west of the 2011 and 2012 burn blocks, respectively, to collect unperturbed meteorological and air-quality data for characterizing ambient conditions within the vegetation layers. Instrumentation mounted at multiple levels on each tower provided high frequency (0.5 or 10 Hz) measurements of meteorological variables as well as carbon monoxide (CO) and carbon dioxide (CO₂) concentrations (from Dataram monitors). The monitoring networks for the experiments also included surface-based monitors that featured phased array Doppler Sonic Detection and Range (SODAR) measurements of wind speed and direction, measurements of near-surface PM_2.5 concentrations, and ceilometer measurements of plume heights and boundary-layer PM_2.5 concentrations (2012 case only). The PM_2.5 concentrations from the surface monitors and the ceilometer measurements of downwind plume heights are used in this study to investigate the performance of the coupled numerical model. A complete description of the monitoring network and instrument specifications for the 2011 and 2012 experiments appears in [20].

2.2. ARPS-CANOPY Numerical Weather Prediction Model

The Advanced Regional Prediction System (ARPS) [24,25] is used to simulate meteorological conditions from the regional scale to the mesoscale, and ARPS-CANOPY [9], a modified version of ARPS that includes a sub-grid canopy parameterization, is used to simulate conditions within and above the forest canopy. For this study, a series of five one-way nested simulations are executed using four outer domains (8100, 2700, 900, and 300 m horizontal grid spacing) and one inner domain (100 m horizontal grid spacing). For reference, the outermost domain covers the Northeastern United States from Virginia northward and from Eastern Ohio eastward, and the inner domain covers a 100 km² area of the New Jersey Pine Barrens surrounding the burn unit. With the exception of case-dependent inner domain grid-center coordinates, simulations are performed for the two cases with identical domain configurations. For the 2011 and 2012 cases, the outer domains are initialized at 0000 UTC 19 March 2011 and 0000 UTC 5 March 2012, respectively, and run for 60 h, and the inner domain is started at 1200 UTC 20 March 2011 and 1200 UTC 6 March 2012, respectively, and run for 12 h. Based on experimentation with different initial and boundary conditions for the outermost domain, two data sources are utilized: the North American Regional Reanalysis (NARR) [26] (2011 case) and the Global Forecast System (GFS) [27] 0.5-degree analysis (2012 case). Land use and terrain data are input from the U.S. Geological Survey (USGS) 1-km (outer domains) and 100-m (inner domain) datasets [28]. Grid stretching is applied along the vertical axis in all simulations, providing 2-m vertical grid spacing in the lowest 84 m of the atmosphere on the inner domain.

ARPS-CANOPY is utilized for the inner domain only, wherein the bulk effect of the forest canopy is represented by plant area density (one-sided leaf area per unit volume; m² m⁻³). Plant area density information is derived from canopy density measurements collected using LiDAR remote sensing [29]. For the 2011 case, the burn block (107 ha) is divided into 121 grid cells. Grid cells are grouped in clusters (hereafter referred to as burn zones) and a 15.5 kW m⁻² surface sensible heat flux, decaying exponentially up to a height of 45 m, is progressively switched on and off in each burn zone to imitate the southwest-northeast movement of the 2011 line fire through the burn unit. As described in [11], the magnitude of the surface sensible heat flux is based on peak 1-min mean vertical heat flux measurements at the 20 m tower before, during and after line fire passage. It is assumed that if a hypothetical 100 m × 100 m square (i.e., the dimensions an ARPS-CANOPY grid cell on the inner domain) was centered over the observed line fire, at most 10% of the square would be burning at any given time. The exponentially-decaying heat flux profile is employed following [30] to attain more realistic plume-average properties. For the 2012 case, the burn unit (97 ha) is divided into a total of 94 grid cells, and multiple line fires progress through the burn unit from east to west as an 861.8 W m⁻² surface sensible heat flux, decaying exponentially up to a height of 45 m. The surface sensible heat flux is progressively switched on and off in each burn zone to imitate the movement of the multiple north-south oriented line fires through the 2012 burn unit. A summary of ARPS-CANOPY fire evolution and canopy vertical structure in each burn unit is provided in Figure 1. Additional details regarding the ARPS-CANOPY model setup (e.g., parameterizations, model time steps, boundary conditions) can be found in [11].

To explore the sensitivity of smoke concentrations to the presence of the forest canopy, two ARPS-CANOPY simulations are run for each case. Simulations are performed with-canopy (“control”; hereafter 2011CL and 2012CL) and without-canopy (“no-canopy”; 2011NC and 2012NC).

2.3. FLEXPART-WRF Lagrangian Particle Dispersion Model

FLEXPART-WRF [14] is a version of the FLEXPART Lagrangian particle dispersion model [15], originally developed for use with global models, that has been adapted for use with the limited area Weather Research and Forecasting (WRF) model [31]. As a Lagrangian dispersion model, FLEXPART-WRF tracks individual particles and at each time step calculates the future position of each particle based on a combination of the mean (i.e., transport) wind and a random wind component. The mean wind is obtained from horizontal (u, v) and vertical (w) wind components from the meteorological model. The user may choose to have the random wind component computed internally using the Hanna turbulence parameterization [32], or computed as a function of turbulent kinetic energy (TKE) read in from the meteorological model. For either option, the user must also choose whether to have the requisite variables [friction velocity, surface heat flux, and planetary boundary layer (PBL) height] computed internally in FLEXPART-WRF or read in from the meteorological model. While particles can be removed in FLEXPART-WRF by dry and wet deposition and by decay, only dry deposition is used in this study. Additional details regarding FLEXPART-WRF can be found in [14].

A routine was developed that transforms output variables from ARPS or ARPS-CANOPY to the variables that FLEXPART-WRF normally ingests from WRF output. The files generated by this routine contain ARPS or ARPS-CANOPY model output, but are otherwise indistinguishable from WRF output files. Note that although the particle dispersion model is referred to in this study by its proper name “FLEXPART-WRF”, the simulations presented in this study are driven by ARPS-CANOPY, not by WRF.

A total of four FLEXPART-WRF simulations are conducted, and are referred to using the same domains and naming convention as the ARPS-CANOPY simulations described in Section 2.1: 2011CL (2011 “with-canopy”), 2011NC (2011 “no-canopy”), 2012CL (2012 “with-canopy”), and 2012NC (2012 “no-canopy”). With the exception of the particle release information, which is derived from emission estimates specific to each case (2011 or 2012; detailed below), the FLEXPART-WRF namelist settings are identical for all four simulations. Note that the same particle release information (i.e., number of particles, total mass, release timing) is used for the “with-canopy” and “no-canopy” simulation pairs. Only one species (PM_2.5) is released in the simulations, and it is assumed that PM_2.5 concentrations are indicative of smoke concentrations in the context of this study [2]. Parameters controlling the FLEXPART-WRF grid are set such that the grid is identical to that of ARPS-CANOPY. The FLEXPART-WRF synchronization interval (i.e., particle time step) is 10 s, the particle position output interval is 10 s, and the averaging interval applied to particle concentration calculations is 300 s. Finally, following the recommendation of [14], namelist options are set such that the random component of particle motion is computed internally using the Hanna turbulence formulation [32]. The PBL parameters required for the internal turbulence calculation (friction velocity, surface sensible heat flux, and PBL height) are ingested from ARPS-CANOPY.

Particle release information required by FLEXPART-WRF is derived using emissions estimates from the Fire Emission Production Simulator (FEPS) [33]. FEPS requires the following inputs to derive emissions estimates: fire event information, fuel loading profiles, fuel moisture profiles, hourly cumulative area burned, and hourly flame-height wind speed. In this study, observed fuel loading profiles representative of available fuels in the litter layer and understory from the 2011 (1479 ± 388 g m⁻²) and 2012 (1149 ± 266 g m⁻²) cases are input into FEPS, and FEPS-computed consumption values are replaced with observed average consumption values of 696.2 g m⁻² (2011) and 507.3 g m⁻² (2012) documented in Heilman et al. [22]. The 2011 and 2012 AC/F simulations employ 222,849 and 119,533 particle tracers, respectively. The larger particle tracer count in the 2011 case is attributable to both larger estimated total consumption and to the larger 2011 burn block (107 ha compared to 97 ha for the 2012 case). The hourly flame-height wind speed is calculated by assuming a natural log profile to extrapolate the 1 m above ground level (AGL) wind speed from each simulation to an approximate flame height of 0.5 m AGL. The 12-h mean reductions from the 1 m AGL wind speed to the 0.5 m AGL wind speed are 0.15 m s⁻¹ and 0.21 m s⁻¹ in the 2011 and 2012 cases, respectively. Total area burned each hour is specified based on the parameterized fire spread in ARPS-CANOPY, as guided by observed fire spread rates during the experiments. Emissions estimates from FEPS in the form of PM_2.5, distributed over time and space, are used as inputs to FLEXPART-WRF. For reference, the total mass of PM_2.5 released in the 2011 simulations (2.23 × 10³ kg) was nearly twice the total mass of PM_2.5 released in the 2012 simulations (1.20 × 10³ kg).

3. Results and Discussion

3.1. Smoke Plume Overview and AC/F Assessment

AC/F simulated smoke plumes in the 2011 and 2012 control simulations (2011CL and 2012CL) are shown in Figure 2; Figure 3, respectively. The full-domain three-dimensional particle plots in the 2011 case (Figure 2, left panels) indicate an initially shallow plume at 10:32:30 EDT that develops into a plume with a depth of 1–2 km by 15:07:30 EDT. The 2011CL plume is tilted from vertical due to the mean flow, but it is more upright than the plume in the 2012CL simulation (cf. Figure 2 and Figure 3).

Comparing the three-dimensional view and the top-down view of the AC/F domain (Figure 2, left and center panels) further illustrates the evolution of the 2011CL simulated smoke plume. During the morning and early afternoon, a northeast surface wind transports smoke particles toward the southwest (10:32:30 and 11:42:30 EDT), and as the particles ascend through the deepening PBL, winds from the north transport particles toward the south (13:07:30 EDT). As the afternoon progresses (15:07:30 and 16:42:30 EDT), surface winds veer to the southeast, transporting ascending particles toward the northwest until they reach the upper portion of the PBL and lower free atmosphere, where winds from the northwest predominate. Particles are then transported toward the southeast across a broad swath of the domain. A closer view of particles in the vicinity of the burn unit, color-coded by elevation (Figure 2, right panels), shows that throughout the simulation, particles rise above the lowest 18 m of the atmosphere (the average canopy height across burn unit; Figure 1c) before or shortly after leaving the burn unit, and most particles outside the burn unit are above 18 m AGL.

As noted previously, the three-dimensional view of the simulated smoke plume in the 2012CL simulation (Figure 3, left panels) suggests a plume that is less upright than the 2011CL smoke plume. Examination of the three-dimensional and top down views (Figure 3, left and center panels) reveals mean winds from the northwest during the morning into the mid-afternoon (10:37:30 to 14:42:30 EDT) evolving into winds from the southwest by early evening (17:32:30 EDT). The top down view shows a “cone” of smoke downwind of the burn unit increasing in width with distance from the burn unit, indicative of both the lack of vertical directional wind shear and the tilt of the plume in 2012CL (compare to 2011CL; cf. Figure 2 and Figure 3).

Analyses that are zoomed into the vicinity of the 2012 burn unit (Figure 3, right panels) indicate that the temporal distribution of particles above and below 18 m can be separated into two phases. At 10:37:30 and 12:02:30 EDT, smoke particles released from the burn unit travel more than 1 km downstream before rising above 18 m. From 13:47:30 EDT onward, most smoke particles released from the fire rise above 18 m AGL before leaving the burn unit, which is similar to the particle evolution in the 2011CL simulation (Figure 2). Additionally, the zoomed-in particle plots highlight differences associated with the single line fire progressing from southwest to northeast in the 2011CL simulation and the multiple line fires that burn simultaneously in the 2012CL simulation.

Evaluations of AC/F-simulated smoke concentrations (Figure 2; Figure 3) are conducted in two ways; first via a comparison of AC/F PM_2.5 concentrations to concentrations measured by the Dataram monitors, and second via a comparison of AC/F estimated plume heights to ceilometer estimated plume heights. Time series of AC/F-simulated PM_2.5 concentrations in the lowest grid layer [0–2 m AGL] and measured PM_2.5 concentrations appear in Figure 4. This assessment considers both the timing and magnitude of PM_2.5 concentrations impacting the Datarams. Simulated concentrations are presented as 9-point mean, minimum, and maximum values centered on the model grid point that most closely corresponds to the location of each Dataram. The 9-point analysis methodology partially accounts for the inability of the model to resolve smoke gradients at the same scale as the grid spacing and also provides information about the spatial heterogeneity of smoke concentrations in the simulation. During the 2011 case (Figure 4a,b), only two of the Datarams (DR2 and DR3 in Figure 2) were impacted by smoke from the fire. At DR2, the period of elevated concentrations associated with the smoke plume starts at approximately the same time in the AC/F simulation and in the observed time series. Regarding concentration magnitudes, the 9-point maximum concentrations are comparable to observed values between 11:00 and 13:00 EDT, but greatly exceed observed values later in the day, whereas 9-point mean values are similar to the observed values from about 13:00 EDT onward. At DR3, elevated concentrations occur approximately one hour earlier in the simulation than they do in the observations. Regarding concentration magnitudes, 9-point mean values are generally most similar to the observed values, with 9-point maximum values approximately three times as large as the observed values. The inconsistent nature of the comparison of the 9-point mean and maximum values and the observations illustrates the challenge of assessing smoke simulations at one or two points.

Comparisons between the time series in Figure 4a,b and the top-down particle plots in Figure 2 (right panels) help explain the differences in timing at DR3. The timing differences between the AC/F simulation and the observations are attributed to uncertainties in the representation of the burn unit in AC/F. DR3 was in reality located at the southern tip of the burn unit [21], immediately east of the initial ignition line (ignition: 09:55 EDT), but due to the somewhat coarse representation of the burn unit in AC/F, the closest grid point to the instrument coordinates is located inside the simulated burn unit (Figure 2). Thus, while it took until 10:44 EDT for PM_2.5 concentrations at DR3 to exceed ambient values, simulated PM_2.5 concentrations in the 9-point zone exceed ambient by 09:57 EDT, only two minutes after ignition. The comparison between the AC/F and Dataram concentration time series during the 2011 case illustrates the difficulty inherent in making smoke predictions at such fine spatial and temporal scales: small-scale details of the fire, fuels, and meteorological conditions become critically important for accurate smoke simulations, but are not straightforward to represent with fidelity in the model.

During the 2012 case, the Datarams were positioned such that all four instruments were impacted by smoke from the fire (Figure 2 and Figure 4c–f). Taken in totality, the agreement between AC/F and the Datarams is better for the 2012 case than for the 2011 case. The overall timing of elevated concentrations at DR1 and DR2 is consistent between AC/F and the observations, with the 9-point mean values of simulated concentrations generally matching the observations more closely than the maximum values, with the exception of two brief periods at DR2, when the 9-point maximum values more closely match observations (Figure 4c,d). In the 2012CL simulation, as well as during the field experiment, the location of DR3 closest to the burn unit resulted in the highest absolute concentration values of the four sites. At DR3, the peak in observed concentration (Figure 4e) occurred approximately 30 min after the start of the fire (11:00 EDT), with reduced but still elevated concentrations occurring throughout the remainder of the time series. The 9-point mean time series exhibits weak variability throughout the day, whereas the 9-point maximum time series exhibits a peak at approximately the same time as the observed time series, before diminishing somewhat and varying weakly thereafter. Examining the 10:37:30 EDT panels in Figure 3, a single line fire is evident along the eastern perimeter of the burn unit, with surface winds from the northwest transporting smoke toward the southeast that directly impacts DR3, brushes DR1 and DR2, and completely bypasses DR4. The timing of elevated concentrations at DR4 in AC/F and the Dataram measurements is consistent (Figure 4f), with concentrations remaining low until mid-to-late-afternoon (16:00–18:00 EDT) after which higher values occur. The smoke plume alignment during the 2012 case (Figure 3, center and right panels) indicates that the arrival of elevated PM_2.5 concentrations in the 2012CL simulation at DR4 is the result of winds from the northwest backing to the southwest during the afternoon.

Ceilometer measurements during the 2012 case allow further evaluation and assessment of the 2012CL simulations. Before proceeding, three ceilometer details need to be addressed. First, the ceilometer was not deployed during the 2011 case, so plume height assessments are possible only for the 2012 case. Second, the ceilometer is a LiDAR instrument that measures backscatter, with concentrations derived using Equations (1) and (2) in [34]. Third, the ceilometer is capable of detecting as many as three smoke layers with heights up to 7.5 km with a resolution of 20 m. For the purpose of comparison to AC/F five-minute mean smoke concentrations, observed ceilometer concentrations are averaged over five-minute time blocks. To determine plume height based on ceilometer concentrations, a methodology was chosen in which the smallest concentration threshold that masked out background (i.e., ambient) concentrations in a time-height cross-section of five-minute mean ceilometer concentrations was considered to be the top of the smoke plume. This threshold, as indicated in Figure 5a, was determined to be 15 µg m⁻³; thus, plume height for the sake of this assessment is defined as the height of the 15 µg m⁻³ concentration contour in the ceilometer observations and in the AC/F simulations. The height of the ceilometer 15 µg m⁻³ contour is compared to the minimum, mean, and maximum height of the 15 µg m⁻³ value in the AC/F simulations for the nine-grid-point box centered on the ceilometer location, as well as the domain-maximum value (Figure 5b). Figure 5c,d shows the orientation of the simulated smoke plume to the ceilometer location at 13:00 and 15:25 EDT. Overall, the domain-maximum plume height most closely matches the observed temporal evolution of plume height indicated by the observations both in magnitude and in reproducing the observed increases and decreases. There are three periods (centered around 12:25, 14:15, and 15:25 EDT) when the domain maximum and the 9-point minimum, mean, and maximum values differ the least. These three periods correspond to periods when the plume is more upright (e.g., 15:25 EDT; Figure 5d), in contrast to the remainder of the simulation when the plume is less upright (e.g., 13:00 EDT; Figure 5c). It is noteworthy that all of the simulated plume height values are larger than the observed values between 15:00 and 16:00 EDT. This analysis highlights the challenges in simulating plume heights with sufficient accuracy to confidently assert that a model is able to reproduce observed plume behavior in a low-intensity fire under a forest canopy.

3.2. Smoke Plume Sensitivity to the Forest Canopy

The sensitivity of aspects of the simulated smoke concentrations (overall character of the smoke plume, PM_2.5 concentration, smoke particle residence time, and plume height) to the presence of the forest canopy (Figure 1c,d) is now investigated. The influence of the canopy is discussed separately for each fire case.

3.2.1. 2011 Case

Comparing with-canopy (2011CL) and no-canopy (2011NC) particle plots for the whole domain reveals that overall differences between the two simulations are minor (not shown). However, zoomed-in top-down panels (e.g., right panels in Figure 2) reveal differences that occur within 1–2 km of the fire. Comparing the near-fire smoke concentrations in 2011CL and 2011NC (Figure 6), it is apparent that the direction of smoke transport near the fire is largely insensitive to the presence of the canopy. However, fewer particles remain below 18 m after exiting the burn unit in 2011CL than in 2011NC (c.f. black and green dots, Figure 6, left and right panels). This result suggests that the near-surface plume in the 2011CL simulation is more upright than in the 2011NC simulation. In 2011CL, buoyant particles move more vertically than horizontally within a column positioned more or less above the fire, while in 2011NC particles move more horizontally than vertically within about 1–2 km of the fire. At distances greater than 2 km from the fire, buoyancy becomes dominant and the particles accelerate upward in both simulations.

The residence time of smoke particles within a range of horizontal radii and vertical distances extending from the release point of each particle is now investigated. An ensemble mean residence time is computed using individual particle residence times grouped into 15-min–wide bins based on the release time of each particle. Figure 7 presents scatter plots of ensemble mean residence times for 2011CL and 2011NC. The horizontal mean residence times (Figure 7a,b) for 500-m and 1500-m radii are on average 53% and 28% larger, respectively, in 2011CL than in 2011NC. Vertical mean residence times for 20 m and 100 m depths (Figure 7c,d) show a weaker and less consistent relationship than for horizontal mean residence time, with on average 6% shorter and 6% longer residence times for particles in the 20 m and 100 m deep layer, respectively, in 2011CL compared to 2011NC. Examination of intermediate depths (40, 60, and 80 m) reveals that the finding of shorter mean residence time in 2011CL only applies to the 20 m depth; for depths of 40 m or greater, mean residence time is longer in 2011CL than 2011NC. Since canopy heights within the ARPS-CANOPY domain are on average about 18 m (Figure 1c), this result suggests that processes within the forest canopy may locally enhance vertical dispersion through the canopy, while the canopy imparts a damping effect on vertical turbulent mixing within the surface layer as a whole. Nevertheless, based on the percent difference magnitudes (28–53% vs. 6%), it appears that the impact of the forest canopy on vertical dispersion is more limited than for horizontal dispersion. The limited effect of the forest canopy on vertical dispersion suggests that vertical mixing in this case is largely determined by PBL processes that are minimally influenced by the forest canopy. It should also be noted that there is more turbulence energy in the horizontal wind components than in the vertical wind component in the simulations (not shown), consistent with turbulence measurements during the fire [21], so the effect of turbulence on dispersion is anisotropic. It is not clear from these simulations what impact, if any, turbulence anisotropy has on the sensitivity of vertical dispersion to the presence of the canopy, but it is a complicating factor that should be investigated in future studies.

Time series of the ratio of the mean residence times in each simulation (Figure 7, inset panels) show that there are no clear trends or persistent changes in the ratios for either of the horizontal analyses or for the 100 m vertical analysis that can be attributed to the inclusion of the canopy. However, there is a change in the ratio after 17:55 EDT in the 20 m vertical residence time analysis (Figure 7c) that persists through the end of the time series. This time corresponds approximately to the onset of the nocturnal PBL in the simulation (local sunset: 19:10 EDT). It is therefore likely that a change from ratios slightly below 1 (slightly shorter 20 m vertical residence times in 2011CL) to ratios approaching 2 (longer 20 m vertical residence times in 2011CL) can also be attributed to changes in the near-surface PBL processes that alter the particle dispersion such that longer residence times occur in 2011CL than occur in 2011NC.

To help explain aspects of the smoke plume and particle behavior in the particle position and mean residence time analyses, vertical profiles of the meteorological variables that affect particle position and residence time are examined. Figure 8 presents vertical profiles of time- and space-averaged mean wind speed, potential temperature, and TKE in 2011CL and 2011NC; note that the TKE in Figure 8 is total TKE, (resolved + subgrid-scale), using the method described in [11] to calculate the resolved TKE. The meteorological variables are spatially averaged across all grid cells for which the fire heat source is applied during each five-minute period, and then temporally averaged over a one-hour period.

The influence of the forest canopy is evident in the reduction of near-surface mean wind speed (Figure 8, top row), where wind speeds in the lowest 18 m (mean canopy height) are about 5–6 m s⁻¹ in 2011NC, but only 1–2 m s⁻¹ in 2011CL. Weaker winds in the near-surface layer are consistent with the longer horizontal residence times noted for 2011CL in Figure 7. Recall from the average plant area density profile in Figure 1c that the densest vegetation in the 2011 burn unit was located in the understory, coincident with the fire, with a secondary maximum in vegetation density located at about 2/3 canopy height. Near-surface values of mean potential temperature are 6–7 K higher in 2011CL than in 2011NC. This result is consistent with the findings of Kiefer et al. [12], wherein reduced mixing due to canopy drag was shown to result in higher air temperatures inside the canopy. However, the colocation of the heat source and densest vegetation in the 2011 case (Figure 1c) suggests that the effect of vegetation on potential temperature may be greater here than in simulations presented in Kiefer et al. [12], wherein the densest vegetation was displaced about 12 m above the heat source (see their Figure 1). The resulting potential temperature profile yields greater surface buoyancy in 2011CL, consistent with the more upright smoke plume compared to 2011NC (not shown).

TKE profiles in Figure 8 show smaller TKE values below 20–25 m AGL in 2011CL than in 2011NC (compare green and brown lines). Note that although particles are more positively buoyant in 2011CL, their mean vertical residence time in the lowest 20–100 m of the atmosphere differs by at most 6% from their counterparts in 2011NC (Figure 7c,d). A possible explanation for this apparent discrepancy is that two effects of the canopy on vertical motion are cancelling each other out: canopy drag directly opposes vertical acceleration, as it does horizontal acceleration, but it also indirectly increases the particle buoyancy through decreased turbulent mixing and longer residence times. Two additional possibilities, as discussed in the context of Figure 7, are that PBL processes in these simulations are overwhelming canopy processes, yielding similar vertical dispersion patterns regardless of the presence or absence of the forest canopy, and that turbulence is dominated by horizontal wind components (i.e., is anisotropic).

Finally, Figure 9 presents the impact of the forest canopy on plume height, diagnosed from the height of the 15 µg m⁻³ domain-maximum five-minute mean concentration contour in the AC/F simulations (as in Section 3.1). The complex pattern of points in the scatter plot makes an assessment of the impact of the forest canopy on plume height challenging. The time series of the ratio of plume heights in the two simulations (Figure 9, inset panel) shows that the plume height ratio is more variable during the first half of the simulation (09:55–14:55 EDT) than during the second half (14:55–19:55 EDT). The mean difference between the 2011CL and 2011NC plume heights during the first half of the simulation is 60 m (~9% of the mean). A deeper plume in 2011CL is consistent with greater buoyancy resulting from higher surface potential temperatures (Figure 8), and longer mean residence times in the vicinity of the fire (Figure 7). The scatter plot shows that plume heights of about 1.6 km occur frequently in 2011CL (between 14:55 and 19:55 EDT, when mixing heights near the fire vary between about 1.5 and 2 km; not shown), whereas during the same time period, plume heights in 2011NC vary between about 0.8 and 1.8 km. Despite this variability, the plume height ratio is smaller during the second half of the simulation because the overall plume heights are greater. The mean difference between the 2011CL and 2011NC plume heights during the second half of the simulation is 101 m (~7% of the mean), again indicating a deeper plume with the forest canopy present.

3.2.2. 2012 Case

Top-down particle position plots zoomed in to the burn unit for the 2012 case are shown in Figure 10. A comparison between Figure 10 and the three-dimensional plots in Figure 3 (left panels) reveals two phases of smoke plume behavior. From the morning into the very early afternoon (Figure 10; 10:37:30 and 12:02:30 EDT), the plume is tilted in both 2012CL and 2012NC, with ‘streamers’ of smoke aligned with the northwest mean wind. Differences between the simulations are mostly limited to the precise orientation of these streamers to the mean wind direction. However, starting in the early afternoon (Figure 10; 13:47:30–17:32:30 EDT), greater differences in near-fire smoke concentrations in the two simulations are evident. A larger number of particles downwind of the burn unit remain at or below 18 m AGL in 2012NC than in 2012CL. Particles exiting the burn unit almost immediately rise above 18 m in 2012CL, while particles in 2012NC sometimes travel more than 1 km before rising above 18 m. This result is consistent with the canopy assessment of the 2011 case presented in the previous subsection.

Scatter plots of mean residence times for the 2012 case (Figure 11) reveal differences between the 2012CL and 2012NC simulations in the horizontal (Figure 11a,b) and vertical (Figure 11c,d) directions. For horizontal radii of 500 m and 1500 m, mean residence times are on average 93% and 44% longer, respectively, for 2012CL; the differences are larger than, but of a similar magnitude to, what was found in the 2011 case. Larger differences in the 2012 case may be related to the weaker fire heat source, which allows the influence of the canopy on horizontal dispersion to be more pronounced. Vertical mean residence times for 20 m and 100 m depths are on average 1% shorter and 7% longer, respectively, in the 2012CL case. Similar to the 2011 case, the finding of shorter mean residence time in 2012CL only applies to the 20 m depth; for depths of 40 m or greater, mean residence time is longer in 2012CL than 2012NC. Also similar to the 2011 case, the impact of the forest canopy on horizontal dispersion is greater than for vertical dispersion.

Time series of the ratio of the mean residence times in 2012CL and 2012NC (Figure 11, inset panels), as with the 2011 case (cf. Figure 7, inset panels) show that there are no clear trends or persistent changes in the ratios before 17:55 EDT. However, while only the 20 m vertical plot (Figure 7c inset) showed a change in the ratio after 17:55 EDT for the 2011 case, all of the time series indicate a change after 17:55 EDT for the 2012 case. As with the 2011 case, 17:55 EDT corresponds approximately to the onset of the nocturnal PBL in the simulation (local sunset: 18:56 EDT). The differences between the two cases are likely attributable to the weaker heat flux in the 2012 case. This result suggests that across a range of scales (10–1500 m), the effect on smoke dispersion of the transition from a daytime convective to a nocturnal stable PBL is weaker when sensible heat fluxes from the fires are stronger.

Time- and space-averaged mean wind speed, potential temperature, and TKE in 2012CL and 2012NC (Figure 12) reveal that although the impact of the canopy on mean wind speed is comparable for the 2011 and 2012 cases, the impact on mean potential temperature and mean TKE is smaller in 2012 (compare green and brown lines in Figure 8; Figure 12). Specifically, the reduction in mean wind speed due to the canopy is about 2–4 m s⁻¹ in the 2011 and 2012 cases, while the reduction in mean TKE due to the canopy is smaller in 2012 (0.5–2 m² s² reduction compared to a 3.5–5 m² s² reduction in 2011), and the increase in mean potential temperature due to the canopy is smaller in 2012 (1–1.5 K increase compared to a 2.5–4 K difference in 2011). Despite the differences in plant area density profiles in the two cases (Figure 1c,d), the impact of canopy drag on mean wind speed is comparable. However, the impact of the canopy on potential temperature, and consequently buoyancy, is weaker in the 2012 case than the 2011 case. This is partly the result of the weaker fire heat source in 2012, and partly due to the relative absence of vegetation near the ground where the fire heat source is strongest. The mitigating effect of the canopy vegetation on turbulent mixing near the ground is weaker in 2012 than 2011 (compare CL and NC lowest-level mean TKE in Figure 8; Figure 12).

Finally, an assessment of plume height in 2012CL and 2012NC, diagnosed from the height of the 15 µg m⁻³ domain-maximum five-minute mean concentration contour in the AC/F simulations, is presented in Figure 13. The plume height scatter plot indicates that the overall influence of the canopy on plume height in the 2012 case is minimal, outside of the period from 15:00 to 16:00 EDT when plume heights are greater in the 2012CL simulation due to a more upright plume (cf. Figure 13 to Figure 5b,d). The time series of plume height ratios (Figure 13 inset) suggest that the result is otherwise insensitive to time of day.

4. Conclusions

This study presents an evaluation of a coupled numerical model that is used to simulate smoke concentrations from low-intensity wildland fires in the presence of a forest canopy. A numerical weather prediction model that includes the effect of a forest canopy (ARPS-CANOPY) is coupled with a Lagrangian particle dispersion model (FLEXPART-WRF). The coupled numerical model is used to simulate smoke concentrations and plume heights from two low-intensity wildland (prescribed) fires that occurred in March 2011 and 2012 in the New Jersey Pine Barrens. The simulated smoke concentrations and plume heights are compared with smoke and meteorological observations collected during the two fires, and sensitivity analyses are performed to assess the impact of the forest canopy on smoke concentrations, plume heights, and the meteorological conditions that affect smoke transport and dispersion.

Results indicate that, in some circumstances, the coupled numerical model can reproduce the timing and magnitude of increases and decreases in surface smoke concentrations that appear in the observations. Differences between the simulated and observed smoke concentrations are attributed to the representation of the burn block and the fire in the coupled numerical model, and to the inability of the model to resolve small-scale smoke gradients. Overall, simulated smoke concentrations for the 2012 case more closely match the observations than do simulated smoke concentrations for the 2011 case. For the 2012 case, the coupled numerical model reproduces the trend in observed plume height and some of the temporal and spatial variations, and it was shown that the results are highly sensitive to details in the model formulation. Sensitivity analyses highlight the effect of the forest canopy on simulated meteorological conditions as well as the smoke concentrations and plume heights that depend on those conditions. Results indicate that the overall effect of the forest canopy is to decrease near-surface wind speed, increase buoyancy, and decrease turbulent mixing in the canopy. However, there are sensitivities to time of day, the extent to which the smoke column is oriented vertically, and differences in plant area density profiles and surface heat fluxes in the two cases. The smoke plumes in both cases are more upright in simulations that include a forest canopy than in simulations without a forest canopy, although the effect is more evident in the 2011 case than in the 2012 case due to differences in the characteristics of the forest canopy and in the magnitudes of the surface heat fluxes. The 2011 case, which has a stronger heat flux and more vegetation near the ground where the heat flux is strongest, exhibits an upright smoke plume that is dominated by the effects of buoyancy, which allows smoke to remain in the vicinity of the fire for a longer period of time. Conversely, the weaker heat flux and less vegetation near the ground in the 2012 case produces a smoke plume that is less upright and on which the impact of buoyancy is weaker than in the 2011 case. Temporal variations in smoke concentrations indicate that the effect on smoke dispersion of the transition from a daytime to a nocturnal PBL is weaker when sensible heat fluxes from the fires are stronger.

The differences highlighted in the comparison between simulated and observed smoke concentrations for the two cases illustrate the need for additional field studies that collect more comprehensive smoke concentration observations. For example, observations at additional locations for a given prescribed fire would enable one to determine whether the differences can be primarily attributed to errors in the magnitude of the simulated emissions or to spatial/temporal displacement of the simulated smoke plume. Observations from a larger number of prescribed fires would support a meaningful quantitative analysis of differences between observed and simulated smoke concentrations, and would motivate additional sensitivity studies to determine whether other model formulations reduce the differences. A rigorous quantitative validation of model performance using multiple prescribed fires and additional model formulations should be undertaken to address this uncertainty and to help inform the implementation of models such as AC/F for operational smoke management decisions. Future model validation studies should also address the impact of model spatial and temporal resolution on the results, to determine whether smoke managers should employ simulations on spatial scales similar to those presented here or if comparatively coarse resolution models would be sufficient. This study documents the challenges in predicting meteorological conditions and smoke concentrations under a forest canopy at scales that are relevant for assessing the local effects of smoke from low-intensity wildland fires. Comparisons between canopy and no-canopy simulations show that assessments of smoke exposure in the immediate vicinity of the fire are sensitive to whether the meteorological effects of the canopy are included. Since variations in smoke concentration on these scales directly affects the health of fire managers and also impacts other sensitive receptors in the vicinity of the fire, including the effect of the canopy is important because it could impact both prescribed fire planning and go-no-go decisions for initiating a fire. Accurate predictions of local smoke concentrations could help fire managers refine prescribed fire plans such that they can better anticipate unexpected impacts from smoke, and can lead to revised plans that enable prescribed fires to occur in conditions that would otherwise be considered marginal or out of prescription. However, the comparisons presented in this study between predicted and observed smoke concentrations highlight the difficulty in assessing smoke concentrations at fine spatial and temporal scales. These results suggest that accounting for fine-scale details of the fire, forest canopy, and fire-induced meteorological conditions is critically important for producing accurate predictions of smoke concentrations and plume heights. Additional testing and demonstration of predictive tools, such as AC/F, are planned to help transition the models into products that fire and air-quality managers can use to assess local air-quality impacts during low-intensity wildland fires in forested environments.