Representing Low-Intensity Fire Sensible Heat Output in a Mesoscale Atmospheric Model with a Canopy Submodel: A Case Study with ARPS-CANOPY (version 5.2.12)

Mesoscale models are a class of atmospheric numerical model designed to simulate atmospheric phenomena with horizontal scales of about 2-200 km, although they are also applied to microscale phenomena, with horizontal scales less than about 2 km. Mesoscale models are capable of simulating wildland fire impacts on atmospheric flows if combustion by-products (e.g., heat, smoke) are properly represented in the model. One of the primary challenges encountered in applying a mesoscale model to studies of fire-perturbed flows is the representation of the fire sensible heat source in the model. Two primary methods 5 have been implemented previously: turbulent sensible heat flux, either in the form of an exponentially-decaying vertical heat flux profile or surface heat flux; and soil temperature perturbation. In this study, the ARPS-CANOPY model, a version of the Advanced Regional Prediction System (ARPS) model with a canopy submodel, is utilized to simulate the turbulent atmosphere during a low-intensity operational prescribed fire in the New Jersey Pine Barrens. The study takes place in two phases: model assessment and model sensitivity. In the model assessment phase, analysis is limited to a single control simulation in which 10 the fire sensible heat source is represented as an exponentially-decaying vertical profile of turbulent sensible heat flux. In the model sensitivity phase, a series of simulations are conducted to explore the sensitivity of model-observation agreement to (i) the method used to represent the fire sensible heat source in the model and (ii) parameters controlling the magnitude and vertical distribution of the sensible heat source. In both phases, momentum and scalar fields are compared between the model simulations and data obtained from six flux towers located within and adjacent to the burn unit. The multi-dimensional 15 model assessment confirms that the model reproduces the background and fire-perturbed atmosphere as depicted by the tower observations, although the model underestimates the turbulent kinetic energy at the top of the canopy at several towers. The model sensitivity tests reveal that the best agreement with observations occurs when the fire sensible heat source is represented as a turbulent sensible heat flux profile, with surface heat flux magnitude corresponding to the peak 1-min mean observed heat flux averaged across the flux towers, and an e-folding extinction depth corresponding to the average canopy height in 20 the burn unit. The study findings provide useful guidance for improving the representation of the sensible heat released from low-intensity prescribed fires in mesoscale models.

grid spacing that allows for some scales of turbulence to be explicitly resolved is a departure from traditional mesoscale model applications in which most or all turbulence is parameterized (Michioka and Chow, 2008;Kiefer et al., 2013;Chow et al., 2019).
Mesoscale models are capable of simulating wildland fire impacts on atmospheric flows tens of meters to hundreds of kilome-60 ters away from the fire, if combustion by-products (e.g., heat, smoke) are properly represented in the model. Examples of such mesoscale models are the Clark coupled atmosphere-fire model (e.g., Clark et al., 1996Clark et al., , 2004Coen, 2005), Active Tracer High Resolution Atmospheric Model (ATHAM) (e.g., Trentmann et al., 2006;Luderer et al., 2006), Weather Research and Forecasting (WRF)-SFIRE (e.g., Mandel et al., 2011;Kochanski et al., 2016), Meso-Non-Hydrostatic (Meso-NH)/ForeFire (e.g., Filippi et al., 2013), Advanced Regional Prediction System (ARPS) with fire heat source (e.g., Kiefer et al., 2008Kiefer et al., , 2009Kiefer et al., , 2010, 65 ARPS plus forest canopy submodel (ARPS-CANOPY) with fire heat source (e.g., Kiefer et al., 2014Kiefer et al., , 2015Kiefer et al., , 2016Kiefer et al., , 2018, and High-Resolution Rapid Refresh (HRRR)-Smoke (e.g., Ahmadov et al., 2017Ahmadov et al., , 2021. They are positioned in between computational fluid dynamics models with combustion submodels, such as the High-Resolution Model for Strong Gradient Applications (HIGRAD)/FIRETEC (Linn et al., 2002), and global-regional models that do not include combustion heat or other byproducts, such as the Global Forecast System (GFS) (Wang et al., 2019). Mesoscale models are able to simulate a broader range of 70 atmospheric scales than computational fluid dynamics models due to the use of fewer limiting assumptions in the derivation of model prognostic equations and the inclusion of a broader suite of physical parameterizations. However, this comes at the expense of explicit simulation of fine-scale flows near the fireline, due to the use of grid spacing too coarse to resolve smallerscale turbulence, the use of numerical smoothing at scales near that of the model grid spacing, and the lack of an explicit combustion model (Coen, 2018). Mesoscale models must be provided synthesized information about the combustion sensible 75 heat source and forest overstory vegetation structure to accurately simulate mean and turbulent flows associated with wildland fires in forested environments.
The primary way that atmospheric flows are impacted by wildland fires is through heat and momentum exchanges between the fire and atmosphere (Jenkins et al., 2001;Kremens et al., 2012). Sensible heat flux has received the bulk of attention in previous studies of fire-induced atmospheric perturbations (e.g., Hiers et al., 2009;Heilman et al., 2015;Clark et al., 2020), 80 a focus that is supported by the conclusion of Luderer et al. (2009) that latent heat flux plays a much smaller role in the development of buoyant plumes above wildland fires than sensible heat flux. One of the primary challenges encountered in applying a mesoscale model to studies of fire-perturbed flows is the representation of the fire sensible heat source in the model.
Because combustion occurs at scales too small to be resolved on the mesoscale model grid, the fire sensible heat source in the model is entirely subgrid-scale. Critically, there is no standard method for representing the sensible heat released from 85 wildland fire combustion in a mesoscale model (Clark et al., 1996(Clark et al., , 2004Sun et al., 2006;Kochanski et al., 2013;Kiefer et al., 2014;Kartsios et al., 2017;Kochanski et al., 2018). This is true of two-way coupled atmosphere-fire models, such as WRF-SFIRE, in which fire-induced changes to the simulated wind field feed back on the fire (and vice-versa), and simpler one-way coupled models, such as ARPS-CANOPY with fire heat source, in which the fire perturbs the simulated wind field, but atmospheric perturbations do not feed back on the fire. Two primary methods have been implemented previously: turbulent 90 sensible heat flux (method 1), either in the form of a heat flux profile (method 1a, HFP; e.g., Sun et al., 2006;Mandel et al., 2011;Kiefer et al., 2014) or surface heat flux (method 1b, SHF; e.g., Trentmann et al., 2006;Sun et al., 2006;Kiefer et al., 2009Kiefer et al., , 2010Filippi et al., 2013;Kiefer et al., 2015Kiefer et al., , 2016Kiefer et al., , 2018; and soil temperature perturbation or "hotplate" (method 2, HP; e.g., Heilman and Fast, 1992;Kiefer et al., 2008;Filippi et al., 2013). Beyond the methods themselves, the poorly defined nature of the model parameters used in the fire sensible heat source methods serves as motivation for exploring the sensitivity 95 of model-observation agreement to these parameters.
In this study, ARPS and ARPS-CANOPY are utilized to simulate background and fire-perturbed atmospheric conditions during a low-intensity operational prescribed fire in the New Jersey Pine Barrens (NJPB; Heilman et al., 2021).  . This study is one part of a broader effort within the SERDP project to explore how atmospheric dynamics, including ambient, fire-induced, and forest canopy-induced turbulence regimes within and near the fire environment affect fire propagation, energy exchange, and fuel consumption. The objective of this study is to provide useful guidance for improving the representation of sensible heat released from low-intensity prescribed fires in mesoscale models, through a comparison of 105 model simulations with the combustion sensible heat source methods outlined above and a suite of observations collected during the fire. Although this study focuses exclusively on low-intensity fires, the study findings may have relevance to mesoscale model simulations of higher-intensity fires, including wildfires. The availability of an extensive observational dataset from a low-intensity prescribed fire, and the previous application of ARPS-CANOPY to low-intensity fires, motivates this focus.
The remainder of this manuscript is organized as follows: the prescribed fire field experiment is described briefly in Sect. 2; 110 an overview of mesoscale model fire sensible heat source methods is provided in Sect. 3; the numerical experiment methodology is discussed in Sect. 4, including descriptions of the model (4.1), model configuration and experiment design (4.2), and analysis methodology (4.3); results and discussion are presented in Sect. 5, including details of the model assessment (5.1) and model sensitivity study (5.2) phases; and the paper is concluded in Sect. 6.

115
The prescribed fire was conducted on 13 March 2019 in an 11.2-ha burn unit located in the NJPB at the Silas Little Experimental Forest (39.9156 o N, 74.5956 o W) in south-central New Jersey, USA ( Fig. 1) (Heilman et al., 2021). The burn unit was situated in a mixed oak -pine stand, with chestnut (Q. prinus L.), black (Q. velutina Lam.), white (Q. alba L.), and scarlet (Q. coccinea Münchh.) oaks, and shortleaf (P. echinata Mill.) and pitch (P. rigida Mill.) pines forming the forest overstory.
Instrumentation consisted of a network of "fire-tracker" and multi-spectral sensors to measure flame arrival times, intensity, 130 and radiative heat flux; and six instrumented flux towers (Fig. 1). Fuel loading, moisture contents, and consumption were estimated from a combination of LiDAR and harvest measurements. To spatially characterize the spread of the prescribed fire through the burn unit, an array of 68 fire-tracker sensors was installed at ground level throughout the plot in a grid with approximately 40-m spacing between sensors (Fig. 1). The sensors consisted of 1.5-mm diameter K-type thermocouples attached to Arduino Feather data loggers (2796, Adafruit, New York, NY, USA) with thermocouple amplifiers (269, Adafruit) and GPS 135 antennas (746, Adafruit) to time stamp and provide location for each sensor. The fire-trackers were buried such that the tip of the thermocouple protruded through the surface fuels. All fire-tracker data were logged at a frequency of 2 Hz.
The flux towers consisted of a control tower (TWR C ), located approximately 180 m north-northwest of the northern burn unit perimeter, to characterize the ambient wind velocities and temperatures during the experiment, and five in-situ towers (east: TWR E , flux: TWR F , south: TWR S , north: TWR N , and west: TWR W ), to measure turbulent fluctuations and sensible 140 (i.e., convective) heat flux at multiple heights above the flame fronts. On the flux towers, sonic anemometers (Model 81000V, R. M. Young, Inc., Traverse City, MI, USA) were mounted at 3, 10, and 19 m above ground level (AGL) (hereafter, referred to as the lower, middle, and upper tower levels), and aligned in the true-north direction to measure variations in the three-dimensional wind velocity components (west-east, south-north, and vertical) and air temperature; TWR S was instrumented with sonic anemometers at the lower and upper tower levels only. Data were recorded at 10 Hz using Campbell Scientific 145 CR3000 data loggers (Campbell Scientific, Inc., Logan, UT, USA). Additional temperature measurements were made using thermocouples (Omega SSRTC-GG-K-36, Omega Engineering, Inc., Stamford, CT, USA) mounted at 0.25, 0.5, 1.0, 2.5, 5.0, 10.0, and 15.0 m AGL (all towers except TWR S , where instruments were mounted at 0.5 m AGL and every meter from 1-10 m AGL). Also, a SOnic Detection And Ranging (SODAR) wind profiler (Remtech PA0, The Villages, FL, USA) was deployed about 80 m south of TWR C to characterize planetary boundary layer wind speed and direction up to about 400 m AGL.

150
At 14:57 Eastern Daylight Time (EDT) on 13 March 2019, personnel from the New Jersey Forest Fire Service ignited a line fire using drip torches along the westernmost boundary of the burn unit, proceeding from south to north, with the first detection by the fire-tracker array at 14:59 EDT. The line fire was allowed to spread with the wind (i.e., a head fire) from the southwestern boundary through the burn unit and beneath the in-situ towers in a generally northeastward direction until reaching the eastern boundary at approximately 16:30 EDT. Burning was confined to the surface and understory fuels; very little 155 overstory vegetation was burned. Active burning was completed at about 16:45 EDT, with smoldering continuing thereafter, and the final detection by the fire-tracker array was at 17:38 EDT.

Fire sensible heat source methods
As stated in Sect. 1, there are two primary ways in which the fire sensible heat source has been represented in previous mesoscale modeling studies: turbulent sensible heat flux (method 1a, HFP; method 1b, SHF), and hotplate (method 2, HP).

160
The choice of method determines to what degree fire-specific model parameterizations are used to represent vertical sensible heat transport from the fire by unresolved fire-induced turbulent eddies, and to what degree the mesoscale model's native (i.e., non-fire condition) land-surface and subgrid-scale turbulence parameterizations are relied on to perform this function. From HFP to SHF to HP, progressively more of the model's native parameterizations are engaged to represent vertical sensible heat transport from the fire.

165
With the profile form of the turbulent sensible heat flux approach (method 1a, HFP), an exponential decay of turbulent sensible heat flux (H) with height (z) is prescribed as a function of a surface value (H S ) and an extinction coefficient (K E ), Equation (1) first appeared in Sun et al. (2006), and is based on Beer's law (Pfeiffer and Liebhafsky, 1959). At the height above the ground where the fire-induced subgrid-scale turbulent sensible heat flux becomes negligible, the turbulent sensible 170 heat flux is computed by the model via small-eddy theory (heat flux is proportional to the local vertical gradient of potential temperature; Stull, 1988), as is standard in subgrid-scale turbulence parameterizations in mesoscale models (e.g., ARPS, WRF; Xue et al., 2000;Powers et al., 2017), where w and θ are vertical velocity and potential temperature, respectively, K is eddy diffusivity, and () and () indicate 175 perturbations and grid-volume averages, respectively. In ARPS-CANOPY, K is parameterized as a function of subgrid-scale turbulent kinetic energy (TKE) and a turbulent length scale, the latter a function of grid spacing and atmospheric stability, and separate K values are computed for horizontal and vertical turbulent mixing (Xue et al., 2000;Kiefer et al., 2013). With the surface flux form of the turbulent sensible heat flux approach (method 1b; SHF), the combustion sensible heat flux is implemented only at the surface level (H S ), and small-eddy theory is used at all grid points above the lowest atmospheric 180 level.
Both forms of the turbulent sensible heat flux method were examined by Sun et al. (2006) via a comparison of model simulations to data collected during the Meteotron experiment (Benech, 1976). The authors found that depositing all of the sensible heat released from the fire within the lowest model grid layer above the ground (i.e., SHF method) led to an unrealistic spike in near-surface buoyancy flux within the buoyant plume above the fire. They also found that prescribing an exponentially-185 decaying sensible heat flux profile with an arbitrary e-folding extinction depth of 50 m resulted in an underestimation of near-surface buoyancy flux and vertical velocity within the buoyant plume (for reference, e-folding refers to a reduction of heat flux to e −1 times the surface value). They suggested that an e-folding extinction depth tied to the density of soot above the fire might provide more realistic plume quantities, but acknowledged that the e-folding extinction depth is likely to vary depending on a number of parameters, including fire intensity, flame height, and the environment surrounding the fire (e.g.,
Because of the previous implementation of the HP approach in mesoscale model studies of wildland fires (e.g., Heilman and Fast, 1992;Kiefer et al., 2008), this method is also explored in this study. With this method, sensible heat released from 195 the fire is incorporated via a soil temperature perturbation from ambient conditions. Sensible heat is communicated to the atmosphere via a bulk aerodynamic formula, standard in mesoscale model land-surface parameterizations, in which the surface sensible heat flux (H S ) is a product of the wind speed at the lowest atmospheric level (M 0 ) and heat exchange represented by an exchange coefficient (C h ) and the potential temperature difference between the lowest atmospheric level and the ground (θ 0 − θ G ) (Stull, 1988), As with the SHF approach, turbulent sensible heat flux at all vertical grid levels above the lowest atmospheric level is computed via small-eddy theory [Eq.
(2)]. It is critical to point out that outside of the burn unit, and outside of the fire sensible heat source application period for within-burn unit grid cells, the model's native land-surface and subgrid-scale turbulence parameterizations are employed.

Model Description
For numerical simulations of the NJPB prescribed fire, ARPS and ARPS-CANOPY are utilized. ARPS is a three-dimensional, compressible, nonhydrostatic atmospheric model with a terrain-following coordinate system (Xue et al., 2000(Xue et al., , 2001. ARPS has been applied across a range of spatial scales, from studies of turbulent flows, using grid spacing as fine as O(1) m (Dupont 210 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). In ARPS the bulk effect of a vegetation canopy on the atmosphere is represented within a single layer, beneath the lowest model grid point (as is standard in mesoscale models, e.g., WRF; Powers et al., 2017). Hereafter, this model is referred to as original ARPS to distinguish it from ARPS-CANOPY.
ARPS-CANOPY is a version of the ARPS model with a canopy submodel in which the effects of vegetation elements (e.g.,

215
branches and leaves) on drag, turbulence production/dissipation, radiation transfer, and the surface energy budget are accounted for through modifications to the ARPS model equations (Kiefer et al., 2013). Such changes allow for explicit simulation of airflow through a multi-level forest canopy. The vertical profile of plant area density (Ap), defined as the one-sided area of plant material per unit volume, is utilized in ARPS-CANOPY to represent the bulk effects of the canopy on various atmospheric processes (e.g., drag and radiative heating/cooling). For a full description of ARPS-CANOPY, see Kiefer et al. (2013).

220
A 1.5-order local turbulence parameterization with a prognostic equation for subgrid-scale TKE is utilized in original ARPS and ARPS-CANOPY, with the addition of canopy source and sink terms in the momentum and subgrid-scale TKE equations in ARPS-CANOPY. In all simulations, the original ARPS anisotropic turbulence option is used in which both horizontal and vertical components of eddy viscosity and diffusivity are computed; this option is recommended when vertical grid spacing is considerably smaller than horizontal grid spacing (Xue et al., 2000). This subgrid-scale turbulence parameterization has been 225 tested extensively and has been found to produce the correct vertical structure of mean variables and turbulent statistics with and without a forest canopy (e.g., Dupont and Brunet, 2008;Kiefer et al., 2013Kiefer et al., , 2014. Furthermore, a land-surface model based on Noilhan and Planton (1989) and Pleim and Xiu (1995), and radiation physics following Chou (1990Chou ( , 1992 and Chou and Suarez (1994) are utilized [with the addition of a canopy heat source term and with canopy shading effects on the ground energy budget accounted for (ARPS-CANOPY only)]. Lastly, for simulations with horizontal grid spacing larger than 100 m, 230 a non-local turbulence parameterization based on Sun and Chang (1986) is also utilized, in addition to the local turbulence parameterization. Table 1. Nested domain summary. In the model column, "Ao" and "Ac" refer to original ARPS and ARPS-CANOPY, respectively. The order of dimensions in the grid and domain size columns is x, y, z.  The vegetation canopy is represented in ARPS-CANOPY as a height-varying Ap profile specified at each grid point. Whereas 245 a meaningful estimation of Ap is challenging at even the scale of single tree stands, it has been demonstrated that LiDARderived canopy height profiles can be utilized to characterize the three-dimensional canopy structure (e.g., canopy bulk density, kgm −3 ) at high horizontal and vertical resolutions on spatial scales O (10 km) (Skowronski et al., 2011;Kiefer et al., 2014).
Here, we have derived Ap on grids with 90-and 30-m (horizontal) and 2-m (vertical) grid spacing from previously acquired aerial LiDAR data (Skowronski et al., 2020;Warner et al., 2020) [D5 dataset shown in Fig  of fire front position every 30 s (Eric Mueller, personal communication). Using GIS software, the period of fire presence within each D5 grid cell was determined. The heat source was applied in a particular burn-unit grid cell if the fire front was located anywhere within the grid cell. The heat source was applied steadily during the grid-cell fire window, with the heat source ramped up and down at the beginning and end of the fire window (total ramp-up/ramp-down period set to 10% of the fire window). A summary of the fire sensible heat source is provided in Fig. 4, beginning with a depiction of the heat source 260 methods (and their vertical distributions) in Fig. 4a. Fireline positions at four times, corresponding to the times of peak 1-min mean turbulent sensible heat flux at towers TWR S , TWR E , TWR F , and TWR N are depicted in Figs. 4b-e, along with the 1-min mean heat fluxes measured at all six flux towers. The time evolution of the model fire heat source area is depicted in Fig. 4f.
Similar to the LiDAR-based Ap dataset, the fire-related data provides ARPS-CANOPY with important information about the spatiotemporal variability of the fire sensible heat source.

265
The study focuses on the D5 simulations and takes place in two phases: model assessment and model sensitivity. In the model assessment phase, analysis is limited to a control simulation (HFP 2.4L ; Fig. 4a, Table 2) in which the fire sensible heat source is represented as an exponentially-decaying turbulent sensible heat flux profile (method 1a, HFP). For reference, case names in this study consist of the heat source method acronym (e.g., HFP) and an alphanumeric subscript with the number denoting the heat source magnitude (e.g., 2.4 for 2.4 kWm −2 ), and the letter denoting the e-folding extinction depth ["L" for low (4.4 For the concentric circles in panels (b)-(e): instrument height increases with distance from center, white shading is used to denote values less than 1 kWm −2 , and dark gray shading is used to denote the absence of middle tower level measurements at TWRS. to the 68-sensor-median peak fire-tracker temperature deviation from pre-fire conditions.

Analysis methodology 295
In this study, ARPS-CANOPY simulations are assessed using sonic anemometer measurements of the background and fireperturbed atmosphere. Although the tower thermocouple and SODAR data have been analyzed (not shown), the sonic anemometer data are the exclusive focus of this study since unlike the former data sources, sonic anemometer measurements provide us the ability to evaluate model-simulated variances (i.e., TKE) as well as mean quantities. For observational data and model output, 1-min means are computed from instantaneous time series (10 Hz sampling frequency for observations; 1 Hz output 300 frequency for model). For observations, quality control and assurance is performed in a four-step procedure. First, tower data are subjected to a despiking and filtering routine to remove erroneous data and values exceeding six standard deviations from running 1-h means . Second, a double-rotation tilt-correction routine is applied to the vertical windvelocity component data to correct the raw data for systematic errors originating from errors in the physical leveling of sonic anemometers in the field (Wilczak et al., 2001). Third, means computed from fewer than 90% of possible instantaneous obser-305 vations (i.e., fewer than 540/600 values) are masked. Finally, manual quality control is performed, resulting in the masking of lower tower level temperature at TWR C and TWR F , and upper tower level wind direction at TWR F .
In the model assessment phase, box-whisker plot vertical profiles, time series, vertical cross-sections, and three-dimensional surface plots of 1-min mean temperature, wind speed, wind direction, vertical velocity, and TKE, are compared between the control simulation (HFP 2.4L ) and observations at all six towers. The multi-dimensional approach taken in this study phase 310 allows for a more comprehensive assessment of simulated and observed variables than point comparisons alone. The approach provides helpful context for the assessment of differences between the model simulation and observations at the individual tower locations; for example, contoured plots can help identify spatial displacement errors (i.e., correct magnitude but incorrect location; Brown et al., 2011). As in Kiefer et al. (2014), the assessment of the control simulation in this phase is qualitative.
However, the multi-dimensional nature of the analysis and greater number of flux towers in this study distinguishes this ARPS-

315
CANOPY assessment from that of Kiefer et al. (2014).
In the model sensitivity phase, box-whisker plot vertical profiles and time series of 1-min mean temperature, wind speed, vertical velocity, and TKE, from all D5 simulations (  Tables 3-5 to Tables S1-S3 in supplemental material). The focus on TKE prediction is justified by the important role that turbulence plays in affecting fire behavior (Banerjee et al., 2020) and in controlling fire effects such as tree mortality (e.g., Kavanagh et al., 2010) and smoke dispersion (e.g., Charney et al., 2019). For the model sensitivity analysis, wind direction is omitted to enable larger figure panels that can accommodate the simultaneous plotting of all seven sensitivity tests ( Table 2). The evaluation of the D5 simulations in this phase is part qualitative and part quantitative. First, a 325 qualitative evaluation of box-whisker plot vertical profiles and time series is conducted. Second, a quantitative evaluation of the model simulations is performed using three statistics: mean difference, expressed as a percentage of the range of observed values, root-mean-square difference, also expressed as a percentage of the range of observed values, and index of agreement (Willmott, 1981), wind speed and vertical velocity appear to be simulated with highest fidelity, and temperature and TKE simulated with lowest fidelity. The model captures the vertical wind speed profile shape, including the positive wind shear between the middle and upper tower level sonic anemometers, and captures the increase in vertical velocity magnitude with height. However, the model 370 exhibits a negative temperature bias, even outside the burn unit at TWR C , and the model considerably underestimates TKE at the upper tower level, in particular at towers TWR E , TWR F , and TWR W .
Proceeding to the time series of simulated and observed variables (Fig. 6), the focus of the analysis shifts to an assessment of temporal variability and overall agreement in the timing of fire-induced variations (relative to the time of peak 1-min mean observed turbulent sensible heat flux). In examining Fig. 6, it is critical to keep in mind that the timing and duration of enhanced 375 values of sensible heat flux ("fire window") in each ARPS-CANOPY grid cell is informed by the fire-tracker data, not the flux tower data. The fire window within each model grid cell is generally longer than the fire window observed at the towers, and the model and tower-observed fire windows may begin and end at different times. Examining Fig. 6, the model is found to generally capture the temporal variability and range of wind speed, wind direction, and to a lesser degree vertical velocity, but not temperature and TKE. For the latter two variables, the magnitudes are considerably underestimated, and the timing of the To provide context for the assessment of ARPS-CANOPY at the flux towers, vertical cross-sections are presented in Fig. 7, along the same northwest-southeast oriented axis used earlier for the Ap vertical cross-section (Fig. 3b-c), and 3D surface plots are presented in Fig. 8. These multi-dimensional plots are intended to help identify displacement errors that can result when 390 the model correctly simulates the magnitude, but not the location, of fire-induced perturbations. Errors in both magnitude and location complicate the model assessment process (see Fig. 6.3 in Brown et al., 2011). Analysis begins with the vertical crosssections in Fig. 7: from northwest to southeast, the approximately 600-m-long cross-section axis intersects TWR C , TWR F , TWR W , and TWR S . Examining Fig. 7 as a whole, evidence of displacement error is found in most panels, but in particular the temperature, wind speed, and TKE panels. An example of spatial displacement error is seen in the left panel of the second row Assessment of of the control simulation concludes with an examination of horizontal slices of 1-min mean variables at the time of peak heat flux at TWR W (15:38 EDT), rendered in 3D space (Fig. 8) (Fig. 8b), and to a lesser degree, TKE (Fig. 8e), are found downwind (northeast) of the fireline at the 3 and 9 m model grid levels. For the other variables, and at the 19 m model grid level for all variables, it is difficult to discern any relationship between the simulated variables and the fireline. Linear structures in the wind speed (Fig. 8b), vertical velocity (Fig. 8d), and TKE ( Fig. 8e) panels, most noticeable at 19 m AGL, are suggestive of planetary boundary layer structures aligned with the mean wind (oriented southwest-northeast). Regarding TKE specifically, Fig. 8e reveals 1-min mean simulated TKE values of 410 5-7 m 2 s −2 at 19 m AGL within and adjacent to the burn unit at a moment when observed upper tower level TKE values at the five in-situ towers are between 1.8 and 9.2 m 2 s −2 . Although this comparison is encouraging, it is important to note that this does not necessarily imply that the model is correctly capturing the processes yielding the observed TKE values at the towers.
What Figs. 5-8 do suggest is that the model is reproducing the overall background and fire-perturbed atmosphere depicted by the tower observations. The shortcoming in accurately simulating upper tower level TKE values at TWR E , TWR F , and TWR W 415 motivates the model sensitivity tests that follow.

Model Sensitivity
In this study phase, box-whisker plot vertical profiles and time series of simulated temperature, wind speed, vertical velocity, and TKE from all seven D5 simulations are examined at TWR W (Figs. 9-10), along with model evaluation statistics (Tables   3-5). Recall from Sect. 4.3 that this tower is chosen for the model sensitivity analysis as it is one of three towers (TWR E , 420 TWR F , and TWR W ) where the control simulation (HFP 2.4L ) was found to considerably underestimate TKE at the upper tower level (Figs 5-8). The reader is reminded that this evaluation of simulations with different fire sensible heat source methods is intended to ultimately provide guidance for improving the representation of the sensible heat released from low-intensity prescribed fires in mesoscale models. The relevance of temperature, wind speed, and vertical velocity to predictions of, for example, smoke dispersion and fire behavior, motivates the examination of these variables in addition to TKE.

425
Beginning the model sensitivity evaluation with box-whisker plot vertical profiles of temperature (Fig. 9a), all simulations correctly depict decreasing median and range of temperature, with height. However, a broad range of lower tower level temperatures (and consequently, buoyancy) during the 40-m analysis period is found among the simulations. The NF 0 , HFP 2.4L , HFP 4.8L , and HP 400 simulations yield temperatures in relatively narrow ranges of 8. 4-10.8, 9.2-15.8, 8.9-20.2, and 9.2-19.4 o C, respectively, whereas the HFP 24L and SHF 24 simulations yield temperatures in the much broader ranges of 9.3-49.4 o C 430 and 9.8-63.8 o C, respectively. The simulation that most closely agrees with the observed 11.3 -28.1 o C range at the lower tower level is HFP 24H , with a range of 9.1-24.2 o C. Proceeding to wind speed (Fig. 9b)  Spacing between x-and y-axis tick marks is 90 m. Manual quality control of tower observations results in the exclusion of three data points: lower tower level T at TWRC and TWRF, and upper tower level WD at TWRF (denoted by black circles). Proceeding to vertical velocity (Fig. 9c), it appears that all simulations capture the observed increase in vertical velocity median and range with height, including the simulations with no or weak sensible heat sources (NF 0 , HFP 2.4L , HFP 4.8L , and the seven simulations are found to cluster in two groups: the three simulations with 100% of the peak 1-min mean observed turbulent sensible heat flux implemented (HFP 24L , HFP 24H , and SHF 24 ) form a "strong" wind speed group (wind speed∼1.5-3 ms −1 ), with the other simulations clustered in a "weak" wind speed group (wind speed∼0.5-2.5 ms −1 ). It is worth noting that despite differences in the vertical distribution of the sensible heat source in the "strong" wind speed group, wind speed differs little between the three simulations.

470
Examining vertical velocity in Fig. 10, only modest sensitivity to the sensible heat source method and parameters is noted.
The most noticeable aspect of the vertical velocity panels is the general lack of negative values at the lower tower level.    Table 5. As in Table 3, but for index of agreement (Willmott, 1981

Summary and Conclusions
In this study, we have examined different methods used to represent the sensible heat release from low-intensity prescribed fires in mesoscale models. Such fires are prevalent during prescribed fire operations in the eastern US and can impact the health and safety of both fire personnel and the general public. A series of simulations conducted with the ARPS and ARPS-CANOPY 500 models have been evaluated using observations collected during a low-intensity prescribed fire in the NJPB. Although this study focused exclusively on low-intensity fires, the study findings may have relevance to mesoscale model simulations of higher-intensity fires, including wildfires.
A two-phase model assessment and sensitivity test approach was utilized in which a control simulation (HFP 2.4L ) was examined first, followed by a comparison of simulations with different sensible heat source methods and parameters. The 505 multi-dimensional model assessment confirmed that the model reproduced the background and fire-perturbed atmosphere as depicted by measurements made at the six towers. However, the model was found to underestimate the upper tower level TKE values at some of the towers. The model sensitivity tests revealed that the best agreement with observations occurred when the fire sensible heat release was represented as a turbulent sensible heat flux profile (method 1a) with 100% of the peak 1-min mean observed heat flux averaged across the in-situ flux towers and an e-folding extinction depth of 18 m (100% of the average 510 canopy height in the burn unit) (HFP 24H ). The poorest agreement was found when (i) 100% of the peak 1-min mean observed heat flux was concentrated at the surface or in the lowest few grid levels (SHF 24 and HFP 24L ), or (ii) the peak 1-min mean observed heat flux value was heavily diluted (10-20% of peak 1-min mean observed value) before implementation in the model (HFP 2.4L and HFP 4.8L ).
The findings of this study have provided useful insight into the representation of the sensible heat source from a low-515 intensity fire in a mesoscale model. The findings suggest that methods that rely more heavily on the mesoscale model's native land-surface and subgrid-scale turbulence parameterizations to move heat vertically through the near-surface atmosphere (SHF and HFP with small e-folding extinction depth), can lead to a concentration of heat near the surface, excessively large nearsurface temperatures, and consequently an overestimation (underestimation) of near-surface (forest overstory) fire-perturbed wind speeds and TKE. On the other hand, implementing a turbulent sensible heat flux in the model corresponding to 10-20% 520 of the peak 1-min mean observed value (i.e., heavily diluting the point heat flux measurement), as in Kiefer et al. (2014), can lead to an overall underestimation of TKE near or just above the top of the forest overstory.
As with any study, limitations must be kept in mind when considering study findings and any potential applications. First, the study findings derive from simulations performed with a single mesoscale model, ARPS-CANOPY, and it's parent model, ARPS. It is expected that the findings will be applicable to other mesoscale models, but this will need be confirmed by the 525 broader model community, particularly for two-way coupled atmosphere-fire models like WRF-SFIRE. Although efforts to couple a fire behavior model with ARPS-CANOPY are underway, the model at this time is still one-way coupled. Second, the model simulations were evaluated with observations from a single low-intensity fire experiment. In future work, model simulations will need to be evaluated using observations from fires spanning different intensities, in different forest environments, and with different synoptic-mesoscale conditions. Third, this study examined only two extinction depths (4.4 and 18 m), with 530 a more complete examination of the extinction depth parameter space left to future work. Finally, the choice of e-folding extinction depth used in a particular mesoscale model simulation may ultimately depend on the intended model application. An e-folding extinction depth that yields satisfactory agreement with flux tower measurements of mean TKE may prove too large or too small in the context of smoke dispersion or fire behavior predictions (e.g., Kartsios et al., 2017;Kochanski et al., 2018).