**In the Line of Fire: Consequences of Human-Ignited Wildfires to Homes in the U.S. (1992–2015)**

**Nathan Mietkiewicz 1,2,3, \* , Jennifer K. Balch 1,2,4 , Tania Schoennagel 4 , Stefan Leyk 1,4 , Lise A. St. Denis 1,2 and Bethany A. Bradley 5**


Received: 21 August 2020; Accepted: 5 September 2020; Published: 7 September 2020

**Abstract:** With climate-driven increases in wildfires in the western U.S., it is imperative to understand how the risk to homes is also changing nationwide. Here, we quantify the number of homes threatened, suppression costs, and ignition sources for 1.6 million wildfires in the United States (U.S.; 1992–2015). Human-caused wildfires accounted for 97% of the residential homes threatened (within 1 km of a wildfire) and nearly a third of suppression costs. This study illustrates how the wildland-urban interface (WUI), which accounts for only a small portion of U.S. land area (10%), acts as a major source of fires, almost exclusively human-started. Cumulatively (1992–2015), just over one million homes were within human-caused wildfire perimeters in the WUI, where communities are built within flammable vegetation. An additional 58.8 million homes were within one kilometer across the 24-year record. On an annual basis in the WUI (1999–2014), an average of 2.5 million homes (2.2–2.8 million, 95% confidence interval) were threatened by human-started wildfires (within the perimeter and up to 1-km away). The number of residential homes in the WUI grew by 32 million from 1990–2015. The convergence of warmer, drier conditions and greater development into flammable landscapes is leaving many communities vulnerable to human-caused wildfires. These areas are a high priority for policy and management efforts that aim to reduce human ignitions and promote resilience to future fires, particularly as the number of residential homes in the WUI grew across this record and are expected to continue to grow in coming years.

**Keywords:** WUI; fire; defensible space; prescribed fire; community vulnerability; fire suppression costs; Zillow

#### **1. Introduction**

Wildfire poses a direct threat to communities and people living in fire-prone areas [1,2]. Communities that meet or intermingle with undeveloped wildland vegetation, creating zones known as the wildland-urban interface (WUI; <10% of the U.S. land area), are at the greatest risk of wildfire across the U.S. [3–5]. Humans' relationship with wildfire in the WUI is a complex interaction between an increasing number of people living in flammable landscapes (increased number of at-risk communities and ignition sources) [6–9], a warmer and drier climate that is more conducive to fire [10–12], and increased fuels accumulated in some places due to years of fire suppression [13,14]. Elements of these interactions have been the focus of contemporary scientific debate and public discourse [15], but we still lack fine-scale estimates of homes threatened from wildfire and how WUI communities contribute to wildfire at a national scale.

The WUI has expanded in the U.S. by a third between 1990 and 2010 [7] and is projected to double by 2030 [3]. Currently, there are an estimated 2.5 million residential structures in fire-prone areas of the WUI in the coterminous U.S., equivalent to \$1.4 trillion in property value at risk [5]. Importantly, only ~15% of the WUI in the west is developed while the remaining 85% is available for development and deemed potential WUI expansion areas [16], highlighting the pressing need to implement strategies now that will reduce future wildfire risk in an expanding WUI. Moreover, exposure to wildfire smoke can have profound negative health effects [17]. An estimated 17.5 million people living in these fire-prone WUI areas [5] in the western US are expected to experience longer and more intense smoke waves, in which unhealthy particulate matter from wildfire smoke persists for more than two days [18]. Displacement from wildfires due to home loss or poor air quality are expected to increase as wildfire occurrence continues to rise [15,19].

Despite the call to adapt to increasing wildfires [15], we lack yearly, fine-scale estimates at a national scale of the number of homes actually threatened by wildfire and the role of human-related ignitions in starting fires near these communities. Previous work estimated that 286,000 homes were threatened by large wildfires (small wildfires, <400 ha, were not accounted for) between 2000–2010, based on decade-interval data within census blocks [7]. Also, it has been documented that 84% of the nation's wildfires were caused by people between 1992 and 2013 [8], but this work does not quantify the role of human-caused wildfires where people actually live. Further, it is problematic that there are no estimates of the wildfire threat to homes past 2010, as four of the top ten largest wildfire years (>32,000 km<sup>2</sup> ) occurred in 2011, 2012, 2015, and 2017 [20].

Human influence on the landscape is an important predictor of wildfire ignition likelihood [4,21–23]. For example, proximity to roads at local scales is positively correlated with wildfire ignition density [4,22], while fire density has a hump-shaped relationship with population density, peaking around ~10 people/km2 [24]. At regional to national scales, it has been shown that proximity to WUI areas influences the occurrence of large wildfires (>40 ha) [25]. Recent work by Syphard et al. [9] has further shown that human presence is more important in wildfire ignition than climate across the conterminous US, but burned area is primarily driven by climate conditions [10,26]. Despite the prominent role of humans in igniting fires, it remains unknown how the WUI acts as a source of human ignitions as well as how much fire damage the WUI potentially sustains at fine scales and on an annual basis. Here, we investigate how human ignitions have altered fire frequency, burned area, and seasonality in the WUI compared to wildland areas. We further assess the number of residential structures threatened and the associated suppression costs of human-started wildfires within the WUI, providing a comprehensive assessment of the consequences and costs of human-started wildfire in our most vulnerable communities.

To address these questions, we utilized over 1.6 million wildfire events (1992–2015) designated as human or lightning caused [27] and a novel spatio-temporal housing dataset of over 200 million housing records derived from Zillow's ZTRAX dataset [28] to provide a direct and improved estimate of the number of residential units threatened (defined based on being located within 1 km of the wildfire perimeter) by wildfire by year. We quantified wildfire costs from ~150,000 wildfire situation reports (1999–2014) that provide daily estimates of fire suppression costs [29]. Combining these three datasets we assessed the role of human-caused wildfires in the WUI (769,087 km<sup>2</sup> ; Interface WUI = 162,368 km<sup>2</sup> ; Intermix WUI = 606,719 km<sup>2</sup> ), very low-density housing (VLDH; 2,283,410 km<sup>2</sup> ), and wildland areas (2,565,320 km<sup>2</sup> ) [7]. We also provide an estimate of the increase in the number of homes in the WUI across the study time period.

#### *Hypotheses*

We hypothesized that human-started wildfires would be the dominant type near the WUI, while the VLDH would contain a mix of human- and lightning-started wildfires, and lightning-started wildfires would dominate the wildlands. Similarly, we expected that wildfires would be more expensive in the wildlands due to the larger wildfire footprint and the difficulty in navigating the terrain during suppression efforts. Conversely, we expected the cost of wildfires near the WUI to be less per wildfire, but more expensive proportional to wildfire size due to the more proximate threat to communities and people. We further hypothesized that the distance to human settlement and the median home density were strongly related to the prevalence of human versus lightning ignitions. Fire frequency and fire season length were expected to vary with the type of ignition source (i.e., lightning versus human) and the distance from the urban core boundary. Lastly, we hypothesized that there would be important regional differences and therefore, we refined some of our analysis based on reorganized level 1 ecoregions to capture the major eastern and western U.S. patterns.

#### **2. Datasets**

#### *2.1. The Wildland-Urban Interface*

A WUI spatial database [7] was used to assess the distribution of wildfire activity within the WUI, very low-density housing (VLDH), and wildland areas for 1990, 2000, and 2010. Each polygon within this dataset represents a U.S. census housing block group, which is defined as an interface, intermixed, or non-WUI category. Intermix WUI and Interface WUI both have census blocks that exceed 6.17 housing units per km<sup>2</sup> and have >50% wildland vegetation or <50% wildland vegetation, respectively. Interface WUI areas must also be located within 2.4 km of an area larger than 5 km<sup>2</sup> containing >75% wildland vegetation. We combined the interface and intermixed WUI zones to represent the total area of the WUI within the US. Very low-density housing areas were delineated based on census block groups that had <6.17 housing units per km<sup>2</sup> and wildland vegetation that is >50%. Wildlands were defined as census block groups housing density of 0 units per km<sup>2</sup> and wildland vegetation >50%.

Previous research has encompassed the WUI within a 2.4 km buffer based on the assumption that fire events outside of the WUI can produce firebrands that are likely to be transported ahead of a wildfire front and ignite a new wildfire in a WUI census block [30,31]. In this current work, we opted to not use this 2.4 km buffer and instead directly evaluate the wildfire activity in each of the main classes independently (e.g., Intermix WUI, Interface WUI, very-low density housing, and wildlands) [32]. In effect, the VLDH category represents a more spatially explicit and representative area between the WUI and "true" wildlands than the 2.4 km buffer zone. Additionally, we acknowledge that the wildland agricultural interface (WAI) and wildland industrial interface (WII) are important areas of wildfire risk near human settlement, the data product used in this study did not differentiate those locations from the WUI or VLDH [7].

#### *2.2. U.S. Forest Service Fire Program Analysis-Fire-Occurrence Database*

The U.S. Forest Service Fire Program Analysis-Fire-Occurrence Database (FPA-FOD) [27] documents the location, cause, discovery date across ~1.8 million wildfires from 1992–2015 from various sources, including U.S. federal, state, and local records of wildfires on public and private lands. We collapsed the US Forest Service-designated cause categories of equipment use, smoking, campfire, railroad, arson, debris burning, children, fireworks, power line, structure, and miscellaneous fires into a human-started wildfire class. All fires <0.001 acres or fires that had a missing/undefined cause were removed from the FPA-FOD dataset [8]. Wildfire size was evaluated by small (<4 km<sup>2</sup> ), large (4–500 km<sup>2</sup> ), and very large (>500 km<sup>2</sup> ) classes. There is acknowledged reporting bias over time (e.g., completeness issues) at the state level [33] making trend analysis difficult [8].

#### *2.3. Monitoring Trends in Burn Severity*

Annual burned area and large wildfire perimeters (1992–2015) were obtained through a combination of buffering the FPA-FOD points based on the size of the fire and from the Monitoring Trends in Burn Severity (MTBS) group [34] (see http://mtbs.gov/ for information on pre-processing, specific information on index calculation, and derived dataset available). MTBS leverages the extensive Landsat TM/ETM+/OLI record to produce spatially explicit fire perimeters for all large fires (>400 ha in the western US and >200 ha in the eastern US). These data are produced using the differenced Normalized Burn Ratio (dNBR), which uses pre- and post-fire spectral response to define full extents of the disturbance events. In addition to this automated process, and to ensure consistency and precision, manual digitization is performed at on-screen display scales between 1:24,000 and 1:50,000.

#### *2.4. Zillow Transaction and Assessment Dataset*

To measure residential settlement at very fine spatial scales and annually we used the Zillow Transaction and Assessment Dataset (ZTRAX), which contains unique data on housing transactions, home values, rental estimates, spatial location, home- and property-related information as well as built-year information, for existing homes and certain other properties across the United States. The ZTRAX data are obtained from a major large third-party provider as well as through an internal initiative, called County Direct. County Direct prioritizes counties based on different characteristics and supplements the third-party coverage by collecting data directly from county Assessor and Recorder's offices and represents a growing share of the ZTRAX dataset. The ZTRAX database over the study period includes public records and assessor data for approximately 200 million parcels in over 3100 counties in the U.S. (https://www.zillow.com/ZTRAX/). While this database is of unique nature covering most of the U.S. there are some issues related to data quality as described above, including spatial, temporal and thematic uncertainties that we assessed and measured in order to create accompanying uncertainty layers. The database is under continuous revision and will be updated regularly. The raw data and the created SQLite databases cannot be shared publicly per the established data share agreement, but there are recently published aggregations to 5-year intervals between 1810 and 2015 [28].

In this study, we measure residential land use at a given point in time using the construction year information within the ZTRAX database. Furthermore, we use the geographic information provided for each record, which represents approximate address point locations, to characterize residential land use at fine spatial and temporal resolution [28]. Using the raw ZTRAX point database, we selected only housing points that were defined as residential properties and had a defined built year to reduce temporal uncertainty. In recent research, data limitations including low spatial precision in census-based housing data [7], low classification accuracy in land cover data [35], and lack of small fires (e.g., MTBS [34]) resulted in limited accuracy of national assessments. Improved estimates of WUI-wildfire interactions are critical to assess the vulnerability of current and future housing development to wildfire risks and costs. By utilizing the novel ZTRAX data in conjunction with extensive federal fire records we were able to measure threatened homes with high spatial precision and on an annual basis, with future potential to continue to identify regions with elevated risk. These new fine-resolution data products provide unprecedented opportunities for such risk analysis at fine analytical scales, fundamentally changing the way patterns of human-fire interactions can be described.

#### *2.5. United States National Incident Command System Historical ICS-209 Reports*

We further leverage an expanded dataset mined from the United States National Incident Management System (NIMS) Incident Command System (ICS) Historical ICS-209 Reports acquired from the publicly available ICS-209 data on the Fire and Aviation Management website (https://fam. nwcg.gov/fam-web/), created by some of the co-authors. We performed systematic and reproducible

data cleaning protocols to access these data in a usable form [29]. The cleaned ICS-209 database, henceforth known as ICS-209-PLUS-WF, resulted in 20,353 unique incidents compiled from over 150,000 daily records between 1999–2014, containing information on total estimated fire suppression cost, total residential and commercial structures destroyed and threatened, total fatalities, and wildfire cause. Here, we defined a community threatened by wildfire when threatened or destroyed structures for a given event were >0. These data are a spatially explicit summary of 67% of the total estimated fire suppression costs from all wildfires that required an incident management team but represents <2% of the total wildfires that actually occurred from 1999–2014. The fire suppression field is systematically underestimated [29]. Therefore, all of our cost estimates are lower than true costs, but they are proportionally similar to the total accrued costs reported by the National Interagency Fire Center (NIFC) [20].

#### **3. Methods**

The details of the WUI classes were extracted based on the point location of wildfire ignition in the FPA-FOD database. The classes within the WUI layer [7] are based on statistics and classified at the census block group level, which vary widely in size from 0.001 km<sup>2</sup> to 3403 km<sup>2</sup> , with a mean of 0.843 km<sup>2</sup> across the conterminous U.S. The size of the census block groups from east to central to west vary from 0.0007 km<sup>2</sup> to 1231 km<sup>2</sup> (mean = 0.469 km<sup>2</sup> ), 0.001 km<sup>2</sup> to 1020 km<sup>2</sup> (mean = 0.937 km<sup>2</sup> ), 0.001 km<sup>2</sup> to 3403 km<sup>2</sup> (mean = 1.55 km<sup>2</sup> ), respectively. While the mean area of census block groups for the WUI, VLDH, and wildlands are 0.651 km<sup>2</sup> , 5.72 km<sup>2</sup> , and 1.57 km<sup>2</sup> , respectively, the classes in the WUI layer are interconnected along an urban-rural continuum characterized by large homogenous clusters of class-similar census block groups that tend to be on the order of 10 s to 1000 s km<sup>2</sup> in size. Thus, the census block groups allow for an appropriate minimum mapping unit to satisfy the assessment from Short [33]; "The FPA FOD should provide point locations of wildfires at least as precise as a PLSS section (2.6 km<sup>2</sup> grid)", while retaining the finer detail in class variation across space without having to aggregate further to an arbitrary pixel unit. Given these vector mapping units of the WUI layer at fine and coarse scales, we believe that our analysis is more than appropriate for quantifying human- and lightning-caused wildfires at the WUI, VLDH, and the wildlands.

Information contained in the WUI data was spatially joined to the FPA-FOD and ICS-209-PLUS-WF point database and temporally aligned such that the discovery wildfire year was matched with the appropriate WUI delineation between 1992–2015. This resulted in an FPA-FOD + WUI and ICS-209-PLUS-WF + WUI databases totaling 1,386,595 and 23,607 wildfire ignition points, respectively. Additional ancillary data layers were subsequently intersected (e.g., level 3 ecoregions, state boundaries), allowing us to further refine boundaries into general regions (i.e., west, central, southeast, etc.). A systematic 50-km grid was imposed to quantify the relative contributions of human- vs. lightning-started wildfire ignition within Intermix WUI, Interface WUI, VLDH, and wildlands. Using these point databases, we were then able to create summary statistics across all grouping combinations (i.e., WUI classes, ignition type, seasonality, fire size, regions, 50-km grid, etc.). We calculated the frequency of wildfire ignitions, mean and 95th confidence interval of fire size, interquartile range of discovery day-of-year. Seasonality was defined as the season that fire ignition was most prevalent across all classes.

For the FPA-FOD and ICS-209-PLUS-WF database, we estimated burned area footprint by developing a circular buffer around each wildfire ignition point, equating to the total burned area associated with the fire event. For each wildfire ignition point that was not associated with an MTBS (large fires) wildfire perimeter, we used this buffering logic to estimate the burned area footprint; for all fires that did fall within the wildfire perimeter and have an association with the MTBS fire, we used the MTBS wildfire polygon. In addition, and because of known inaccuracies in the FPA-FOD point locations [33], we buffered the MTBS polygons equivalent to 2.6 km<sup>2</sup> as suggested by Short et al. (2015), which identified only 16 MTBS wildfires outside those fire polygons, which we subsequently buffered to the fire size (see Accuracy Assessment for more details).

To create estimates of wildfire threat, buffered rings were created around the estimated wildfire polygons at distances of 250 m, 500 m, and 1 km from the wildfire perimeter edge. Three buffer distances were selected to better understand and quantify the estimated number of homes threatened than with a single metric because there are no national standard defining distances from a wildfire's edge that determine evacuations and highest home threats. Those definitions are made at the moment and can vary across space as the wildfire event unfolds. These boundary distances were chosen to best exemplify a gradient of moderate (1 km) to high (250 m) home threat by wildfires. Homes within 1 km of a wildfire are not only at higher threat of wildfire relative to homes further than 1 km from wildfire edge, but those homes will likely fall within an evacuation zone and directly impact homeowners during the event.

We extracted all residential homes from the cleaned ZTRAX database for all years leading up to and including the year of the wildfire event that fell within the wildfire perimeter and each of the buffered FPA-FOD and ICS-209-PLUS-WF rings. We further spatially intersected the WUI census block groups with these buffered polygons to estimate the total WUI class area that was burned/threatened by wildfire. ANOVAs with a pairwise comparison of the means using the Tukey HSD function were used to test for differences between fire size, ignition type, homes threatened, and wildfire suppression costs between the WUI, VLDH, and wildlands.

To test whether wildfire ignition type, timing, and frequency is a function of the distance from human settlement and home density, we calculated the distance of all wildfires from the urban core boundary, proxied by the high-density urban class in the WUI database, and subsequently binned the data by 10-km grid cells across the contiguous US. This allowed us to more easily generalize relationships while retaining the spatial integrity of the data. We also used the 10-km grid cell to estimate the median home density. We employed generalized additive models (GAMs), a common statistical tool in wildfire science [36,37], to explore the relationship between the frequency and IQR of human versus lightning-caused wildfire against median home density and distance the WUI edge.

#### *Accuracy Assessment*

We conducted a comparative analysis to estimate a baseline accuracy metric using two datasets, (1) the MTBS polygons that were found within the FPA-FOD database and (2) the buffered burned area estimate using the FPA-FOD points that had an associated MTBS polygon. We did this assessment based on homes found in the ZTRAX database within wildfire perimeters. We found 288,846 homes threatened within buffer-estimated wildfire polygons. The raw MTBS polygons reported 156,915 homes threatened within the burned perimeter. This equates to the buffered burned area overestimating the number of homes threatened within the wildfire polygon by 84%. The major source of these differences between the buffered points and the raw MTBS polygons are due to the irregular shape of the raw wildfire polygons that tend to be more elongated and follow the topography, vegetation distribution, and wind direction while the buffered points do not take into account the underlying land structure.

Of the 17,381 wildfires found in the MTBS database from 1992–2015, 8576 were found within the FPA-FOD, which equates to 0.6% of the total wildfires (1,700,188) found in the FPA-FOD, or 60% of all fires >200 ha (14,327) in the FPA-FOD. We found that 49% of the wildfires in the MTBS product between 1992 and 2015 had an associated FPA-FOD ignition point that was located within the wildfire burned perimeter. When we buffered the perimeters equivalent to 2.6 km<sup>2</sup> , we only found 16 additional FPA-FOD points associated with the MTBS product out of the 8864 polygons that were not found within the perimeters, suggesting that the remaining polygons that did not have an FPA-FOD ignition within the perimeter and beyond were not solely due to spatial inaccuracies but instead these fires were completely absent from the FPA-FOD product. It is known that there is a lack of overlap between the FPA-FOD when compared to satellite-based detections, which may relate to intentional burns, a lack of suppression efforts or reporting on small fires, or just the absence of actual events [38]. The additional 16 polygons that had an associated FPA-FOD ID within the 2.6 km<sup>2</sup> buffer had a median distance from the perimeter of 76 m, with over half of those points located with 135 m from the perimeter edge further suggesting that the 2.6 km<sup>2</sup> metric may be overestimating the inaccuracies of the FPA-FOD database. More refined assessments of this validation are needed for the wider fire community.

#### **4. Results**

#### *4.1. Percent of All Wildfires Originated in the WUI*

The WUI represents ~10% of the total conterminous U.S. land area in 2010 (ca. 7,780,091 km<sup>2</sup> ), yet 32% of all wildfires (N ~ 437,233 between 1992 and 2015) originated in the WUI. The remaining 29% and 39% of wildfire events started in very low-density housing (VLDH) and wildlands, respectively. Wildfires starting in the WUI burned a total of 16,412 km<sup>2</sup> , representing 4% of the total area burned in this record. Mean fire size in the WUI was significantly smaller (4 ha (3.3–4.3)) than in the VLDH (23.6 ha (21.8–25.9)) and wildland areas (58.3 ha (54.4–62.5); values in brackets indicate a 10,000 repeated bootstrapped 95th confidence interval around the mean). Further, across the record, the total number of residential homes in the WUI (both interface and intermix) increased by 145% with ~32 million new residential homes in the combined WUI between 1990 and 2015 (Figure S1).

#### *4.2. In the WUI, Humans Caused Nearly All Wildfires, Doubled Fire Season Length, and Increased the Number of Wildfires More than 20-Fold, Compared to Lightning-Started Wildfires*

Humans caused 97% of all wildfires (n = 424,700) in the wildland-urban interface, 85% of all wildfires (n = 459,054) in the very-low-density housing, and 59% of all wildfires (n = 242,047) in the wildlands between 1992 and 2015 (Figure 1A, Table 1). All fires < 0.001 acres or fires that had a missing or undefined cause were removed from the dataset prior to analysis [8]. In total, human-caused wildfires, irrespective of where they were ignited, burned 12,411 km<sup>2</sup> of the WUI (Figure 1B), 78,484 km<sup>2</sup> of the VLDH (Figure 2B,C), and 82,934 km<sup>2</sup> of the wildlands (Figure 2B,C). In the WUI and VLDH ~50% and ~30% of that area were burned by small fire events (< 4 km<sup>2</sup> ), while in the wildlands ~62% of that area was burned by large fire events (4–500 km<sup>2</sup> ). Human-caused large wildfires (4–500 km<sup>2</sup> ) burned 3311 km<sup>2</sup> of the WUI. Notably, there were no very-large wildfires (>500 km<sup>2</sup> ) that started in the WUI. Large, human-caused wildfires burned >4 times more of the WUI in the western U.S. (n = 123; 1549 km<sup>2</sup> area burned) than the eastern U.S. (n = 113; 1228 km<sup>2</sup> area burned), while small human-caused fires were frequent and burned >7 times more land area of the WUI in the eastern U.S. (n = 327,108; 7027 km<sup>2</sup> area burned) as compared to the west (n = 69,449; 1008 km<sup>2</sup> area burned) (Table S1). Wildfires in the WUI were predominantly caused by people across all regions but wildfires in the WUI, VLDH, and wildlands were most abundant in the southeast and western coastal areas of the U.S. (Figures 1A and 2A). Outside of the WUI, lightning-caused wildfires were most frequent and burned the most land area throughout the intermountain west, while human-caused wildfire dominated wildfire burned area and frequency along the Pacific coast and throughout the eastern U.S. (Figure 2). The eastern regions experienced the greatest growth in wildfire occurrence within the WUI, VLDH, and wildlands due to humans, by a factor of 5, 5, and 6, respectively, compared to the increase in wildfire occurrence due to lightning (Table S1).

Human-caused fires exceeded twice the length of the fire season in the WUI, VLDH, and wildlands, expanding it to 148, 164, and 141 days, respectively, compared to lightning wildfire season (WUI = 72 days, VLDH = 45 days, and wildlands = 42 days; Table 1, Figure 3). Human-caused wildfire ignitions were most prevalent during the spring months (March, April, May) in the east and the summer in the west (June, July, August) (Figure 4, Table S2). Human-caused wildfires in the WUI during the shoulder seasons caused the greatest area burned (spring = 4665 km<sup>2</sup> ; fall = 4641 km<sup>2</sup> ) and the highest number of fires (spring = 179,736; fall = 76,041) (Table S2a). While in the VLDH and wildlands, lightning-started wildfires in the wildlands during the summer months caused the greatest area burned (VLDH = 59,563 km<sup>2</sup> , wildlands = 135,614 km<sup>2</sup> ) and the highest number of fires in the spring for VLDH area and summer for wildland areas (VLDH = 178,206, wildlands = 134,309) (Table S2). The median discovery date for human-caused fires in the WUI was over 3-mo earlier (20

May) than lightning-caused fires (26 July), in the VLDH was 2-mo earlier (9 May) than lightning-started fires (26 July), and wildlands over 3-mo earlier (22 June) than lightning-started fires (27 July) (Table S1). Humans-caused fires in the WUI resulted in a longer wildfire season, as the median discovery date for autumn fires was 33 days later for human-caused fires than for lightning-caused fires (16 October and 13 September, respectively).

**Figure 1.** (**A**) The total number of wildfires (dot size) that originated in the wildland-urban interface stratified by the proportion of wildfires caused by humans, (**B**) the total fire size (dot size) as a percentage of WUI that was burned by human-caused wildfires, and (**C**) the total number of homes threatened in the WUI stratified by the proportion of wildfires caused by humans within each 50-km pixel between 1992 and 2015. Black lines represent state boundaries. For (**A**,**C**), reds indicate a greater proportion of human-caused wildfires, while blue indicates a greater proportion of lightning-caused wildfires. Note, only the buffer estimated burned area was used (**B**), not the 250/500/1000 m buffers.

#### *Fire***2020**, *3*, 50


**Table 1.** Key wildfire metrics within the wildland-urban interface (WUI), very low-density housing (VLDH), and wildlands across the conterminous U.S. (CONUS) and stratified by ignition between 1992 and 2015.

**Figure 2.** Spatial representation of fire effects in the very-low-density housing and wildlands. (**A**) The total number of wildfires (dot size) stratified by the proportion of wildfires started by humans, (**B**) the total fire size (dot size) as a percentage of each class that was estimated burned by human-started wildfires, (**C**) the total fire size (dot size) as a percentage of each class that was estimated burned by lightning-started wildfires, and (**D**) the total number of homes threatened stratified by the proportion of wildfires started by humans within each 25-km grid cell from 1992–2015. Black lines represent state boundaries. For (**A**) and (**D**), reds indicate a greater proportion of human-started wildfires, while blue indicates a greater proportion of lightning started wildfires.

Within the WUI, humans increased the number of wildfires, relative to the number of lightning-caused wildfires, 36-fold in the eastern (n = 318,163) and 23-fold in the western U.S. (n = 66,232). Human-caused ignitions at and near the WUI (WUI plus areas of very-low-density housing) resulted in 883,680 (64%) wildfire ignitions, 74,603 km<sup>2</sup> of area burned (20%). The effect of human ignitions on wildfire frequency (Figure S2) and season length (Figure S3) is clearly evident within at least 80 km from the high-density urban areas, compared to lightning ignitions. See Tables S3 and S4 for state and level 3 ecoregion analysis.

**Figure 3.** Frequency distributions of human- and lightning-caused wildfires by Julian discovery day of year stratified by wildfires that started either in the wildland-urban interface (WUI), Very low-density housing (VLDH), or wildlands within the eastern or western U.S. between 1992 and 2015.

**Figure 4.** Spatial representation of season that fire ignition was most prevalent across all classes stratified by human and lightning started wildfire in (**A**) WUI, (**B**) very-low-density housing, and (**C**) wildlands. Winter is defined as the months of December, January, February; Spring as March, April, May; Summer as June, July, August; Fall as September, October, December. Dot sizes indicate the total number of wildfire ignitions from 1992–2015 in a 50-km grid cell.

Overall, human-caused wildfire frequency peaks at median home densities associated with the WUI in the west (40 homes/10 km<sup>2</sup> ), central (34 homes/10 km<sup>2</sup> ), and southeastern (14 homes/10 km<sup>2</sup> ) U.S. (Figure 5), and in higher median home densities in the northeastern U.S. characteristic of urban zones (3100 homes/10 km<sup>2</sup> ). There are no clear trends for lightning-caused wildfire as a function of

median home density, but generally decreasing relationship with increasing home density. Across all regions, fire season length of human-caused wildfires peak at the transition between VLDH and WUI areas (Figure 5).

**Figure 5.** Regional relationships between human- and lightning-caused wildfire ignition frequency and the log median home density within 10-km pixel between two decades (1994–2004 and 2005–2015). Solid lines are based on the best fit Generalized Additive Model regressions with 95th confidence envelope. Dotted vertical lines indicate the division between the WUI and VLDH categories, assuming constant vegetation cover, where urban/WUI boundary equals log(741.3162) home density and the WUI/VLDH boundary equals log(6.17) home density.

#### *4.3. Suppression Costs were Significantly Higher When Protecting Homes*

Overall suppression costs/km<sup>2</sup> were more than ten times greater when any homes were at risk, regardless of ignition type, compared with when no homes were at risk (\$57,670/km<sup>2</sup> vs. \$4076/km<sup>2</sup> ). Additionally, in communities that were directly threatened by wildfires, the associated suppression costs proportionally increased with wildfire size, irrespective of ignition type (Figure 6).

The average cost of suppression in the WUI doubled (\$32.5 million to \$65.1 million per wildfire) and tripled in the wildlands (\$559,000 to \$1.7 million per wildfire) from 1999–2014. Suppression of human-caused wildfires in the WUI cost \$78,105/km<sup>2</sup> per year (\$58,202–\$98,682). Similarly, the per year total suppression costs of controlling human-caused wildfires in the WUI averaged \$26 million (\$18,875,066 – \$34,843,955) (Table S5). Human-caused wildfire costs also varied geographically, with suppression of WUI fires in the West (\$220,080/km<sup>2</sup> ) costing twelve times more than in the eastern U.S. (\$18,300/km<sup>2</sup> ). Similarly, VLDH and wildlands fires in the West (VLDH = (\$75,944/km<sup>2</sup> (\$52,615–\$102,362), wildlands = \$98,105/km<sup>2</sup> (\$72,142–\$127,172)) cost six times more than the eastern U.S. (VLDH = (\$9,756/km<sup>2</sup> (\$5,852–\$15,222), wildlands = \$21,225/km<sup>2</sup> (\$9,964–\$38,665)). Regardless of where the fire originated, the yearly average costs to suppress wildfires nearly doubled from \$809 million to \$1.55 billion between 1999 and 2014. Further, over half of the costs to fight wildfires were spent on wildfires that threatened residential homes (53% or \$8.8 billion; Table S6).

**Figure 6.** The log mean number of (**A**) fire suppression costs and (**B**) homes threatened per human- and lightning-started wildfire event initiated within the wildland-urban interface (WUI), Very low-density housing (VLDH), and wildlands, stratified by fire size in hectares between 1992 and 2015. Error bars indicate the 95th confidence interval around the mean of each group. Tukey's HSD pairwise comparison of the means is represented by differing letters and letter combinations indicating significant differences among groups (*p* < 0.0001).

#### *4.4. Sixty Million Residential Homes Cumulatively were Within One km of Human-Ignited Wildfires in the WUI*

During the study period (1992–2015), over one million homes cumulatively were within perimeters of human-caused wildfire in the WUI, with another 58.8 million homes within 1 km of a wildfire in the WUI (see methods for accuracy assessment). Regardless of ignition location, human-caused wildfires threatened 96% (74,243,133) of all homes threatened by wildfire between 1992 and 2015, while lightning-caused wildfire threatened 2,702,625 homes. Human-caused wildfires accounted for 97% (1,037,018), 93% (127,570), 91% (133,333) of total cumulative residential homes threatened within fire perimeters that ignited in the WUI (Figure 1C), VLDH, and wildlands (Figure 2D), respectively (Table 1).

Human-caused wildfires in theWUI threatened an annual average of 2,493,730 (2,213,141–2,787,175) homes between 1999–2014 (Table S5). Human-caused wildfires in the WUI threatened 30 times more residential homes per year (51,677 (42,692–63,871)) than lightning-caused wildfire (1640 (1021–2404)). Although human-caused fires account for less than half of the overall burned area, they were the main cause of homes threatened, with over 74 million residential structures cumulatively within 1 km of human-caused wildfires. Regardless of development class, small human-caused wildfires were more prevalent and consequently the greatest aggregate threat to homes, with 92% (73,738,356) of homes within wildfire perimeters or threatened by small fires (Table S1). The average number of homes threatened was proportional to the relative size of the wildfire area (Figure 6B). For instance, the smallest wildfires (<1 km<sup>2</sup> ) threatened on average 1.4 homes (1.3–1.5), while the largest wildfires (>500 km<sup>2</sup> ) threatened on average 551 (71.1–1350.0) homes.

#### **5. Discussion**

This study illustrates how theWUI, which accounts for only a small portion of U.S. land area (<10%), acts as a major source of wildfires, almost exclusively human-caused (Figure 1A). WUI-originated fires also account for a disproportionate number of homes threatened and suppression costs, relative to its area, and expand to burn large areas of surrounding very low-density housing areas (Table 1). Nearly a third of all wildfires between 1992 and 2015 started in the WUI, making these areas a high priority for policy and management aimed at reducing human ignitions as well as promoting resilience to future fires. Our findings directly connect two primary aims of the U.S. Forest Service's National Cohesive Wildland Fire Management Strategy [39], namely, protecting homes, communities and other values at risk; and managing human-caused ignitions. Furthermore, understanding that the WUI is a major source of human-caused ignitions, which pose significant threats to WUI communities, will be useful in designing more effective WUI-centric outreach programs to reduce the number of unwanted human-caused fires [40].

Human-caused ignitions increased the average fire season length two-fold and fire frequency more than two-fold relative to lightning ignitions within the WUI. Previous studies have found that early onset warming and drier conditions [41,42] from anthropogenic climate change [43] and human-caused wildfire [8] have increased the frequency and length of wildfire season across the U.S. This study found a more nuanced relationship: the vast majority of spring wildfires across the U.S. were almost exclusively anthropogenic and occurred within the WUI. By delineating where wildfires are occurring (i.e., the WUI) and at what time of the year the wildfire season is expanding (i.e., the shoulder seasons [8,44]), policy makers and government officials can more accurately assess budget allocations and risks associated with seasonal behaviors to reduce wildfire occurrence. We further found that wildfire in the WUI, where communities are at most risk, occurred throughout the entire year, solely by human-caused ignitions.

We documented that the presence of people enables human-related ignitions well beyond the WUI, but there are important regional differences (Figure 5 and Figures S2–S4)) in human-dominated fire frequency and fire season length. Human-caused wildfires are dominant from the WUI to at least 80 km from the urban core in the western U.S., and greater than 100 km in the northeast and central U.S. This result implies that the human imprint on wildfire extends far beyond the currently defined WUI boundary, that human ignitions may replace the role of lightning ignitions in some places, and that the footprint of WUI development may be underestimated, and more varied, than currently represented [7]. The very low-density housing class, which represent 29% of US land area, will be a critical transition zone for wildfires to burn into WUI areas due to their relatively close proximity to the WUI, yet are not accounted for in previous WUI analyses [3,5,7,31].

Changing climatic conditions have increased temperatures and promoted longer and more extreme precipitation events in some places [45,46] leading to an higher fuel loads, both live and dead, in many fire some systems [15]. Over the past several decades, greater fire activity and total burned area has occurred [10,11,19,43,47–51], which is predicted to increase in the U.S. in coming years [52–54]. This study found that the majority of human-caused wildfires were relatively small (<4 km<sup>2</sup> ) but were responsible for the vast majority of homes threatened (92%). Of critical importance, wildfires did not have to be large to threaten homes. In the future, with more extreme climatic [55] and weather conditions [49] these small wildfires, however, have a greater potential to grow into larger wildfires [47].

Our results highlight a critical overlap between human ignitions and infrastructure assets at risk from fire. Though the ZTRAX database underestimates the total number of homes across the U.S. by roughly 30% (86.7 million properties) compared to the U.S. Census (125.8 million homes) (more detailed evaluation of the ZTRAX data and the derived HISDAC-US surfaces is described in Uhl et al. [56]), the increased spatial and temporal precision allowed us to analyze the potential impact of all wildfire sizes (Figure 7), resulting in an estimated five times more homes threatened within wildfire perimeters than previously documented [7]. Moreover, from 1990 to 2015, there were an estimated ~32 million more homes built within the WUI (145% increase), demonstrating substantial growth over this period.

**Figure 7.** Examples of three wildfires that were human-caused illustrating the enhanced spatial resolution of homes threatened by wildfires using the ZTRAX database compared to the information contained within the census block groups via SILVIS. (**A**) The Topanga wildfire was an urban, human-caused wildfire near Santa Monica, California that started on 28 September 2005. The WUI dataset estimated 7595 homes, while the ZTRAX estimated 936 homes threatened within the fire perimeter. (**B**) The Black Forest wildfire was an arson-caused wildfire in the suburbs of Colorado Springs, CO that started on 11 June 2013. The WUI dataset estimated 1480 homes, while the ZTRAX estimated 859 homes threatened within the fire perimeter. (**C**) A rural, human-caused wildfire from debris burning in northeastern Tennessee that started on 2 March 2006. The WUI dataset estimated 554 homes, while the ZTRAX estimated 54 homes threatened within the fire perimeter.

We found that suppression costs are tightly linked to structures threatened, rather than simply being in the WUI, suggesting that implementing defensible spaces around homes remains a critical step in mitigating fire risk and suppression costs. While it has been suggested that the tripling in fire suppression expenses since the early 1990s has been partly due to increased wildfire risk within the WUI [16,54], we found that suppression costs from wildfires originating within the WUI were small (Table 1), likely because WUI fires tended to be small. WUI-originated fires may be likely detected, accessed, and suppressed earlier and more easily—likely resulting in smaller fires on average (mean fire size WUI = 5 ha vs. wildlands = 75 ha) without the need for an incident command response. Nonetheless, of the fires that needed an incident command response, the relatively low contribution of the WUI to fire suppression costs is counterbalanced by the significant threat of wildfires to structures in the WUI. Overall, while suppression costs within the WUI were low relative to total suppression costs across the U.S. when homes were threatened by wildfire, suppression costs were always high (Table S3).

We show that people are directly responsible for increasing wildfire risk where there are homes, highlighting the importance of policy efforts that aim to reduce ignitions in the WUI. There has been strong messaging around reducing the human-related ignitions that start forest fires, e.g., the Smokey Bear advertising campaign around campfires [57]). But it is critical to note that the primary sources of human ignitions in the WUI are debris burning, arson and equipment use; and outside the WUI the most prominent are railroads, roadways, power lines, and campsites [8,21]. There is a critical need to target messaging and efforts directed at reducing ignitions from these different and complex human activities in the WUI; we effectively need a suburban Smokey Bear. Further, given the year-round fire season—predominantly driven by human-related ignitions—it is also important to increase awareness about fire risk outside of the typical summer wildfire season, and support restrictions on debris burning in spring and fall months when fire danger is high, especially in the eastern U.S. In addition, new initiatives could focus on areas with much finer precision as an important result of employing novel settlement data products (Figure 7; [28,56]) in this study, which showed that risk across the WUI is much more varied than assumed and can be closely interlinked to the location of residential structures. Generally, more effective ways to communicate fire risk to the public and decrease ignitions near vulnerable communities is needed [58,59].

Approaches for reducing wildfire risk to communities have proposed ignition-resistant building materials and fuel-reduction treatments in and around homes in the WUI fire-prone areas of the western United States [30,31,60–62]. Community protection is a key priority in national firefighting strategies, but we should consider where fuel reduction strategies would be most effective—at the individual home level or focused within the WUI [15,63]. In conjunction with defensible space management, wildfire suppression that is more narrowly focused on lands surrounding communities that are most at-risk should be prioritized, allowing wildfire to burn in the wildlands where fire will have an ecological benefit and support future wildfire hazard mitigation and utilizing prescribed burns as a mechanism to mitigate extreme future fires [15,54,64,65].

#### *Limitations*

There are some limitations and observed biases in the utilized datasets that are important to acknowledge. First, the FPA-FOD database captures many small events; wildfires less than or equal to 100 ha account for 98% (1,675,166 out of 1,700,188 events) of the database. There is an implicit bias towards reporting of small fires near communities, as they are more likely to be reported than small events in the wildlands. These smaller events are also more likely to be human-caused, and because they are near human settlement, they are more likely to be shorter in duration and smaller because they are identified quickly and easier to access. Due to the smaller size, these fires are almost always overlooked when exploring wildfire presence on the U.S. landscape, particularly because one of the most commonly used fire perimeter datasets has a cut-off threshold for inclusion of 1000 acres in the west and 500 acres in the east (e.g., MTBS). Yet, there has been a longer selection bias in the wildfire literature excluding these smaller fires (e.g., [19,41,47,49,66,67]), which are responsible for the vast majority of homes threatened. Furthermore, smaller wildfires are much higher in frequency than large events, which may be a greater concern to WUI communities under a changing climate [43,44].

Second, while the FPA-FOD database is unique in that it contains critical information on the type and location of the wildfire ignition, it lacks final wildfire perimeter vector polygons. When perimeters were not available through the MTBS dataset, we used a circular buffer equal to the final size of the wildfire burn area. For FPA-FOD wildfire ignition points that had an associated MTBS wildfire ID we substituted that particular wildfire perimeter for the generalized buffer. These MTBS wildfire perimeters account for 83% (7693 out of 9260) of all large wildfires in the combined database (>400 ha). This addition of the 1 km (or 100 ha) buffer is large enough to capture the variation in wildfire perimeter due to topographic/weather/and or suppression efforts. When the wildfire was less than 100 ha, the 1-km buffer completely consumed the known perimeter, making the more preferable and detailed wildfire perimeter unnecessary for those calculations. Because wildfires less than or equal to 100 ha

account for 98% of the wildfire database, we are confident that our estimates for homes threatened within 1 km of an event are reliable. Further, it has been estimated that embers can travel up to 1 mile (1.6 km) or further from a fireline, a mechanism generally explored only for larger wildfires [68]. Since the majority of the wildfires in our database are smaller fires (<100 ha; n = 1,675,166) rather than larger fires (>400; n = 7693), we chose a more conservative estimate of distance from the fireline of 1 km rather than 1 mile.

Third, based on comparisons with the US census data for 1990, 2000, and 2010, the ZTRAX data underestimate housing counts by 30%. Therefore, the estimates provided here are likely underestimating the true threat to homes. ZTRAX housing property counts and county stats on housing units are strong, >95% in urban counties and >75% in rural counties during this time period [56].

#### **6. Conclusions**

We demonstrate the remarkable effect that humans have had on introducing wildfire to all landscapes of the conterminous US, most importantly within the WUI, through changes in wildfire frequency, seasonality, burned area, and associated costs. We found that human-caused wildfire was responsible for the vast majority of residential homes threatened between 1992 and 2015. People are starting almost all of the wildfires that threaten our homes. Further, the costs to protect homes from wildfires account for over half of the suppression budget. These results have two implications; first, there is a great need to determine how best to reduce human-related ignitions that result in damaging wildfires; and second, this provides greater justification for implementing fuel treatments and prescribed burns, where safe, to mitigate the risk and threat of future wildfires, particularly near and within the WUI. A better understanding of our own relationship with fire, how we promote it, and how we are vulnerable to it, ultimately will help us live more sustainably with fire.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2571-6255/3/3/50/s1.

**Author Contributions:** Conceptualization, N.M., J.K.B., S.L. and B.A.B.; Data curation, N.M. and L.A.S.D.; Formal analysis, N.M.; Funding acquisition, J.K.B.; Investigation, N.M. and T.S.; Methodology, N.M., J.K.B., T.S. and S.L.; Project administration, J.K.B., T.S. and B.A.B.; Resources, J.K.B., S.L. and L.A.S.D.; Supervision, J.K.B. and T.S.; Validation, N.M.; Visualization, N.M.; Writing—original draft, N.M., J.K.B., T.S., S.L., L.A.S.D. and B.A.B.; Writing—review & editing, N.M., J.K.B., T.S. and S.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by the Earth Lab through the University of Colorado, Boulder's Grand Challenge Initiative. Partial support was provided through the Humans, Disasters, and the Built Environment program of the National Science Foundation, Award Number 1924670 to the University of Colorado Boulder.

**Acknowledgments:** We thank A. Mahood and R. Chelsea Nagy for valuable discussions and J. Uhl for the original processing of the ZTRAX database. We also thank A. Park Williams, Todd Hawbaker, and Melanie Vanderhoof for their detailed manuscript review. Data and materials availability: The authors were provided access to the Zillow Transaction and Assessment Dataset (ZTRAX) through a data use agreement between the University of Colorado Boulder and Zillow Inc. More information on accessing the data can be found at http://www.zillow.com/ZTRAX. The results and opinions are those of the author(s) and do not reflect the position of Zillow Group. Support by Zillow Inc. is gratefully acknowledged. The code used to analyze the data can be found at https://github.com/ NateMietk/human-ignitions-wui. All data used, with the exception of raw ZTRAX, is publicly available; a gridded version of the ZTRAX data can be found here: https://dataverse.harvard.edu/dataverse/hisdacus.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Parcel-Level Risk Affects Wildfire Outcomes: Insights from Pre-Fire Rapid Assessment Data for Homes Destroyed in 2020 East Troublesome Fire**

**James R. Meldrum 1, \* , Christopher M. Barth 2 , Julia B. Goolsby 1 , Schelly K. Olson 3 , Adam C. Gosey 3 , James (Brad) White 3 , Hannah Brenkert-Smith 4 , Patricia A. Champ <sup>5</sup> and Jamie Gomez 6**


**Abstract:** Parcel-level risk (PLR) describes how wildfire risk varies from home to home based on characteristics that relate to likely fire behavior, the susceptibility of homes to fire, and the ability of firefighters to safely access properties. Here, we describe the WiRe Rapid Assessment (RA), a ¯ parcel-level rapid wildfire risk assessment tool designed to evaluate PLR with a small set of measures for all homes in a community. We investigate the relationship between 2019 WiRe RA data collected ¯ in the Columbine Lake community in Grand County, Colorado, and whether assessed homes were destroyed in the 2020 East Troublesome Fire. We find that the overall parcel-level risk scores, as well as many individual attributes, relate to the chance that a home was destroyed. We also find strong evidence of risk spillovers across neighboring properties. The results demonstrate that even coarsely measured RA data capture meaningful differences in wildfire risk across a community. The findings also demonstrate the importance of accounting for multiple aspects of PLR, including both hazards and susceptibility, when assessing the risk of wildfire to homes and communities. Finally, the results underscore that relatively small actions by residents before a fire can influence wildfire outcomes.

**Keywords:** parcel-level risk; risk assessment; post-fire analysis; risk mitigation; rapid assessment; natural hazards

#### **1. Introduction**

In October 2020, Grand County, Colorado, was severely affected by the 193,812-acre East Troublesome Fire. Fueled by drought, beetle-killed trees, and red flag weather conditions (i.e., high winds, dry fuels, and low relative humidity), the fire grew rapidly, including a spread of 87,093 acres in the 24 h starting on the afternoon of October 21. The fire resulted in two deaths and destroyed 366 homes. This paper investigates whether parcel-level rapid wildfire risk assessment data collected before the fire can help to explain why these homes were destroyed while others in the fire's path were not. In so doing, it more broadly investigates the relevance of parcel-level characteristics in determining wildfire risk to homes and people.

The East Troublesome Fire was one of many devastating wildfires in 2020. The general risk that wildfire presents to society is well-recognized, e.g., [1–3], as are the roles of development patterns [4–6] and climate change [7,8] in contributing to increasing risk.

**Citation:** Meldrum, J.R.; Barth, C.M.; Goolsby, J.B.; Olson, S.K.; Gosey, A.C.; White, J.; Brenkert-Smith, H.; Champ, P.A.; Gomez, J. Parcel-Level Risk Affects Wildfire Outcomes: Insights from Pre-Fire Rapid Assessment Data for Homes Destroyed in 2020 East Troublesome Fire. *Fire* **2022**, *5*, 24. https://doi.org/10.3390/fire5010024

Academic Editor: Travis Paveglio

Received: 22 November 2021 Accepted: 10 February 2022 Published: 12 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Reflecting broad concern, many organizations with interests in reducing wildfire risk conduct some version of wildfire risk assessment. The U.S. federal government supports two databases with national coverage that seek to describe wildfire risk at either the census tract (Federal Emergency Management Agency's National Risk Index, https: //hazards.fema.gov/nri/, accessed on 11 February 2022) or municipality scale (United States Department of Agriculture's Wildfire Risk to Communities, https://wildfirerisk.org/, accessed on 11 February 2022). State-wide assessments include the Oregon Wildfire Risk Explorer (https://oregonexplorer.info/topics/wildfire-risk?ptopic=62, accessed on 11 February 2022), and the Colorado Forest Atlas (https://coloradoforestatlas.org/, accessed on 11 February 2022). These assessments tend to follow the risk framework laid out by Finney [9], Thompson and Calkin [10], and Scott et al. [11] and packaged for a public audience at https://wildfirerisk.org. This framework, which is generally consistent with other established risk reduction frameworks [12,13], recognizes wildfire *risk* as the intersection of wildfire *hazard* and *vulnerability* to fire. Wildfire *hazard* encompasses the *likelihood* of a fire at a given location and the *intensity* of that fire if it occurs. *Vulnerability* encompasses both the *exposure* and *susceptibility* of people or assets to expected fire behavior. For example, people are exposed to hazards if they live in homes that are at risk, while their susceptibility depends in part on whether they would be able to safely evacuate during a wildfire. As another example, two homes might be similarly exposed to wildfire hazards but differ widely in susceptibility, and thus overall risk, if one is a wooden cabin and the other is a concrete bunker.

Here, we refine this general wildfire risk framework to focus on parcel-level wildfire risk to homes and residents. Specifically, we define *parcel-level risk* (PLR) as the combination of the local wildfire *hazard* posed to a residential parcel and the *vulnerabilities* of people and property to that hazard (gold boxes in Figure 1). PLR emphasizes intra-community heterogeneity in risk and is embedded within broader-scale contexts that also might influence risk to homes and residents, such as general social vulnerability or determinants of landscape-level hazards such as proximity to wildland vegetation (light green with dashed lines in Figure 1).

**Figure 1.** Conceptual model of parcel-level risk (PLR) that shows wildfire risk to homes and residents parsed into hazard and vulnerability (gold); actions that can reduce risk by addressing its components (light blue); and categories useful for a parcel-level wildfire risk assessment such as the WiRe RA ¯ (green). PLR is embedded within broader-scale contexts (light green with dashed lines) that influence risk to homes and residents but are typically not measured at scales sufficiently refined for capturing relevant parcel-level variation.

As described in the subsections below, four categories of action (blue in Figure 1) can address these aspects of risk at the parcel level. To summarize, fuel reduction near the structure and support for safe and effective fire suppression can reduce both the likelihood and intensity of wildfire on and near the property; planning and preparation before a wildfire can reduce the exposure of people and property as well as the susceptibility of exposed people and property; and structural hardening can reduce the vulnerability of homes and the people living in them. Thus, in the PLR conceptual model, risk varies at the level of individual parcels as a function not only of conditions associated with the possibility of wildfire and conditional fire behavior at that parcel, but also of conditions related to structural vulnerability and the potential for safe and effective response by fire suppression resources.

As such, PLR provides a framework for describing wildfire risk to homes and residents, built on a foundation of fire science that combines theoretical modeling, laboratory experiments, and post-fire investigations. Although embedded within the general framework of risk described above, PLR is motivated by the recognition that large scale wildfire risk assessments tend to have limited coverage of conditions that vary at the scale of individual properties and treat entire classes of assets, such as residential property, as responding uniformly to a hazard. In contrast, PLR offers an explicit focus on the heterogeneity of risk among the residential parcels that comprise a geographic community, consistent with fire science establishing differences in wildfire vulnerability [14] and recent studies emphasizing the influence of factors at multiple scales on the effectiveness of fire risk reduction strategies [15]. By design, PLR focuses on property conditions that can be influenced by residents or property owners.

The concept of PLR can be implemented by the *parcel-level wildfire risk assessment* of a relatively small set of attributes on residential properties. Specifically, four categories of observable attributes (green in Figure 1)—parcel-level hazard, defensible space, access, and structure—succinctly describe the status of the parcel as it relates to PLR. In practice, parcel-level wildfire risk assessments tend to take one of two forms, each filling a separate dimension of programmatic need: in-depth evaluations conducted on-site, or rapid assessments conducted from the road or sidewalk. An in-depth property evaluation is intended to provide detailed guidance (e.g., remove four tagged trees) for mitigating risk on a property and typically entails on-site engagement between a wildfire practitioner and a property owner. A rapid assessment, when conducted as a census of all properties within a community, can serve numerous purposes, including: identifying risk reduction priorities; generating material for outreach and education; tracking changes in vegetation and mitigation over time; informing suppression strategies; exploring gaps between perceived and assessed risk; and conducting post-fire analyses such as this paper. To these ends, some organizations with interests in wildfire risk mitigation, such as fire departments, state agencies, and non-governmental organizations, conduct parcel-level rapid wildfire risk assessments [16]. Although past research supports the relevance of detailed property characteristics assessed in an in-depth assessment [17], such as whether windows are double-paned and whether roofing includes a non-combustible fiberglass underlayment [18], such attributes are considered beyond the scope of parcel-level rapid wildfire risk assessment. Instead, rapid assessments typically measure relatively few attributes within each category, prioritizing those that can be evaluated rapidly and at a distance (i.e., from the road or sidewalk), and with only a small set of possible responses for each attribute.

Using a quick, systematic approach facilitates the assessment of PLR for all residential properties within a given community. However, one might question whether a small set of coarsely measured attributes, often measured imperfectly due to privacy and resource constraints, can be useful for representing a property's risk. This paper addresses that question. Specifically, we test whether parcel-level rapid wildfire risk assessment data, collected before the fire, explain the destruction of structures from a wildfire. To do so, we make the most of an unfortunate opportunity arising from the East Troublesome Fire. A

year before that fire, member organizations of Grand County Wildfire Council (GCWC) conducted a parcel-level rapid wildfire risk assessment for a select set of communities in Grand County, Colorado. They employed the WiRe RA, a parcel-level rapid wildfire risk ¯ assessment approach that has been developed into a standardized tool by the authors and partners. The WiRe RA was developed over more than 15 years of collaboration by a team ¯ of researchers and practitioners, referred to as the Wildfire Research (WiRe) Team, as part ¯ of a systematic data collection and integration approach (the WiRe Approach) intended ¯ to inform local wildfire risk education efforts and allow for the monitoring of community adaptation over time [19]. For all residential homes within the selected communities, the WiRe RA recorded data on parcel-level attributes pertaining to defensible space, structural ¯ hardening, property access and identification, and parcel-level hazard. One of the selected communities, Columbine Lake, lay in the path of the East Troublesome Fire's rapid spread. On the night of 21 October 2020, the fire destroyed 30 structures. (Note that thirty structures in the Columbine Lake community were reported destroyed in the East Troublesome Fire. Of these, 26 of these were included in the rapid assessment data from 2019. Of these 26 assessed but destroyed structures, three have at least one attribute marked as unobservable in the rapid assessment data, resulting in 23 completely assessed structures destroyed by the fire. None of these numbers include the many structures that sustained damage during the fire (e.g., burns, smoke damage) but were not considered complete losses) in that community).

Few studies investigate the relationship between pre-fire parcel-level attributes and the likelihood of a structure being destroyed in a fire [17]. However, such analysis can be valuable. While some parcel-level attributes are impossible to determine after the fact if structures are destroyed in a fire [20], post-fire analysis is uniquely positioned to capture the dynamics of an event, which include the complexity of fire behavior as well as decisions and constraints regarding fire suppression and response. Post-fire analyses using pre-fire parcel attributes are rare for at least three reasons. First, the probability of a wildfire spreading into any community in any given year is quite low, suggesting a small pool of potential study locations. Second, the data must be recent, because many fire-relevant attributes (e.g., vegetation around the home) change over time. Third, sample sizes (i.e., the number of homes exposed to the fire) tend to be low, a constraint we seek to overcome by estimating the effect of overall risk scores, or individual assessed attributes, on whether a structure was destroyed in a series of individual models, each of which accounts for spatial correlations in both dependent and independent variables.

Our results demonstrate that assessing PLR using a parcel-level rapid wildfire risk assessment helps to explain the outcomes for individual residential properties within the Columbine Lake community. These results underscore the importance of evaluating both wildfire hazard and vulnerability when considering the risk of wildfire to people, homes, and communities. They also demonstrate the utility of a parcel-level rapid wildfire risk assessment, such as the WiRe RA, in representing PLR and its heterogeneity across ¯ parcels. Furthermore, the results suggest risk interdependencies across parcels but also that relatively small actions taken by residents can affect home survivability during a wildfire event.

#### **2. Parcel-Level Wildfire Risk Assessment**

This section summarizes wildland fire science relating to the four categories of observable attributes and describes their implementation in our parcel-level rapid wildfire risk assessment, the WiRe RA. This section also notes which of these attributes can be directly ¯ influenced by the residents or owners of individual properties.

#### *2.1. Parcel-Level Hazard*

In the general model of wildfire risk, likelihood refers to the probability of wildfire at a location, and intensity refers to fire behavior conditional on wildfire. Three main factors generally influence wildfire behavior: weather, fuel type, and topography [9]. In the PLR conceptual model, likelihood refers to the probability of a wildfire in the vicinity being transmitted to the property, and intensity refers to fire behavior on and near the property, if that were to occur. So defined, both likelihood and intensity are influenced by variation in fuels and topography at the parcel level. For example, post-fire simulations have found different fuel types to predict home destruction [21]. Grasses typically support fast-moving, low-intensity fires, whereas denser fuels support higher-intensity flames; both can pose high hazards in specific conditions [22]. Post-fire simulation has found slope, relative to the direction of a spreading fire, to predict home destruction [21]. A steep topography can increase the rate of wildfire spread or cause erratic behavior [23–25], and rough terrain can impede firefighter access [26]. Thus, a steep topography places nearby homes at risk [27], especially during rapid, wind-driven fires when suppression is particularly challenging [20,28].

Furthermore, nearby structures can themselves become ignition sources. Post-fire studies have found that the spatial arrangement of buildings, including building density [26] and distance to nearby buildings [29,30], predicted housing damage and that close home proximity allowed wildfire to spread quickly between homes [31]. Fuels and structure proximity can also interact; one post-fire study found that low-density, rural developments were at higher risk from wildfire than high-density developments, with possible explanations including the associated increased exposure to flammable vegetation, limited road access, and complex terrain [27].

Thus, local geographic characteristics affect PLR by influencing both the likelihood and intensity of wildfire on a property. The WiRe RA operationalizes parcel-level hazard ¯ with field-identified attributes measuring the distance to hazardous topography (e.g., valleys, cliffs, chimneys), the slope of the ground near the structure, and the general type and density of fuels on and around the property (i.e., attributes 1 through 3 on Table 1). Although not assessed during the initial Grand County project, distance to the closest home has since been integrated into the WiRe RA, and it was calculated retroactively for Grand ¯ County properties using spatial data (i.e., attribute 4 on Table 1). Furthermore, although weather conditions can vary systematically at the parcel level, they are not separately assessed during the WiRe RA. ¯

**Table 1.** Description of attributes collected in the Wildfire Research (WiRe) Rapid Assessment and ¯ descriptive statistics for structures in the Columbine Lake community, Grand Lake, CO, with complete assessments, by whether structures were destroyed in the East Troublesome Fire or not. Shading depicts grouping of attributes into categories (i.e., 1–4: parcel-level hazard; 5–6: defensible space; 7–9: access; 10–12: structure).



#### **Table 1.** *Cont.*

#### *2.2. Defensible Space*

Wildfire likelihood and intensity are also affected by the presence and quality of defensible space around the structure. The phrase 'defensible space' originates from the reference to this zone as supporting safe and effective fire suppression activities around a structure. Indeed, defensible space has been found to reduce entrapment risk for firefighters, ostensibly increasing their effectiveness [32]. Given that fuels in the defensible space zone can pass flame to a structure through direct point of contact, radiative heat, or firebrands (airborne embers) [14], defensible space also provides passive protection in terms of reducing fuels available near the structure.

While firefighters' experience can confirm the value of a zone that supports safe fire suppression activities (i.e., the active mechanism), laboratory experiments support the value of reducing fuels that could transmit fire to the home (i.e., the passive mechanism). In experiments, structures burned upon direct contact with flames [33], but mock home structures survived radiant heat from an active crown fire that was at least 10 m away [34]. Empirical post-fire case studies also indicate the role of defensible space in reducing fire transmission from nearby vegetation. Analysis of vegetative conditions derived from pre-fire aerial imagery in southern Australia demonstrated a significant impact of 40 m defensible space on home survival, because it distanced the home from the higher radiative heat and ember density near burning vegetation; increased distance from upwind shrubs and trees up to 100 m away had a lesser but still significant effect on home survival [35]. Two other post-fire studies found that pre-fire vegetation conditions near the home, measured with aerial and satellite imagery, helped to determine structure damage during fires in

eastern Australia [36] and northern California [30]. Aerial imagery after fires in San Diego County, CA, suggested that up to 30 m of defensible space increased home survival, with little added benefit beyond that [22]. The same study suggested that only 30 to 40 percent of vegetation must be removed to achieve defensible space protection, provided the remaining vegetation is properly spaced and does not overhang the structure [22,27]. However, effectiveness can depend on specific fire conditions, as some studies found vegetation measures, including defensible space, to be relevant but less important predictors of home destruction than other factors such as topography in two wildfires in San Diego, CA and Boulder, CO [26].

Although these empirical studies tend to focus on the passive value of fuels reduction when interpretating their findings, the passive effectiveness of fuel reduction is confounded by active effectiveness in terms of how it affects decisions about suppression made dynamically during an event. Indeed, defensive actions have also strongly predicted structure damages [36]. Thus, results regarding defensible space in one incident may have limited transferability to other contexts, particularly if fire behavior and response differ.

To represent these dual mechanisms of risk reduction, the PLR conceptual model (Figure 1) links defensible space to fuels reduction near the structure and support for safe and effective fire suppression. The WiRe RA operationalizes defensible space in two at- ¯ tributes. The first, vegetative defensible space, measures the closest distance from the residence to overgrown, dense, or unmaintained vegetation (i.e., attribute 5 on Table 1). The second, other combustible materials, measures the closest distance from the residence to other combustible materials, such as wood piles, wicker furniture, propane tanks (i.e., attribute 6 on Table 1). The defensible space category is differentiated from parcel-level hazard not only because of the relationship with supporting safe and effective fire suppression near the structure but also because, in contrast to the parcel-level hazard category, the two attributes measuring defensible space reflect the outcomes of actions taken (or not) by the residents of a property. Thus, attributes in the defensible space category offer opportunities for residents to intervene in the wildfire risk to their homes or themselves.

#### *2.3. Access*

Access considerations make up the third category of the PLR conceptual model (Figure 1). Access considerations relate to two interventions in wildfire risk to a home and its residents: the ability of suppression resources to safely access the home during a wildfire, which can reduce the hazard, and the ability of residents to evacuate safely, which can reduce their exposure to the hazard.

Ingress/egress refers to available routes for connecting a property to a reasonably distant, safe location. Many fire-prone communities, including in Colorado, have just one path of ingress/egress [37], making firefighter access and resident evacuation difficult [38]. A study of various southern California wildfires found distance to major roads to be an important predictor of home destruction in low-density communities, perhaps due to firefighter response times [27]. Similarly, safe ingress/egress achieved via preparatory fuel clearing allowed firefighters to access burning properties in a southern Californian fire [32].

Driveway clearance, including width, length, and the presence of a turnaround, affects the ability for fire engines to enter a property—and rapidly exit if necessary. Ideally, the driveway or cleared space around it are wide enough for two vehicles to pass each other and for an emergency vehicle to turn around [39,40]. Driveway width requirements were found to increase firefighter safety in southern California, allowing personnel to return to burning houses behind the fire front [32].

Address visibility can also support safe and effective fire response. At night or in smoky conditions, visible addressing, GPS address data [41], and paper maps [32] can all help firefighters to orient themselves and reach the correct destination. During a fire, personnel may need to call in additional resources to specific locations; visible addressing allows firefighters to quickly locate homes and respond to an active fire. This can be particularly important during large fires when personnel come from other locations and

are not familiar with the area. Accordingly, community wildfire preparedness plans or codes often include standards pertaining to address visibility, such as reflectivity, font size, posting location, and combustibility. Despite this perceived importance, we know of no study that assesses the efficacy of address visibility in influencing the risk to homes during a wildfire.

The WiRe RA operationalizes access considerations into three attributes: ingress/egress, ¯ driveway clearance, and address visibility (i.e., attributes 7 through 9 in Table 1). In some cases, ingress/egress is determined at the community level; however, individual homes within communities with otherwise good ingress/egress options might have only one feasible route to the property due to, for example, cul-de-sacs or dead-end roads. Unless enforced by local codes and zoning, driveway clearance and address visibility will vary at the parcel level and are often only observable by visual observation of the property. These latter two attributes are also within the control of residents and/or property owners.

#### *2.4. Structure*

Structural characteristics make up the fourth category of the PLR conceptual model (Figure 1). Depending on materials, design, and condition, roofing, siding, and items attached to a home (e.g., decks, fences), all can be vulnerable to firebrands and can ignite the rest of a home [18]. However, some studies have found siding to be less important to wildfire mitigation than other home characteristics, because most siding materials can resist ignition, given adequate distance from vegetation or other combustibles [18,27]. Decks and fences can pass flames to the home or to other combustibles [18,42], even with "non-combustible" materials; debris accumulated in cracks or underneath a deck are more common sources of deck ignition than the flammability of the deck itself [43]. Post-fire analysis of two fast-moving southern California wildfires found attached wooden fencing to significantly predict home destruction, regardless of whether homes were in a high- or low-exposure environment [20]. Although analysis was constrained by a lack of pre-fire data, that study also suggested that risk is reduced by clearing debris from gutters, eaves, and roofs and by using fire-resistant or non-combustible materials wherever possible [20].

The WiRe RA operationalizes structural considerations with three attributes: roofing ¯ material, most vulnerable siding material, and the presence and combustibility of attachments (i.e., attribute 10 through 12 on Table 1). These attributes primarily relate to structural hardening and associated reductions in vulnerability to wildfire, but the interface between fuels and attachments also affects the localized hazard. As with access considerations, these attributes typically vary at the parcel level. In some cases, existing databases (e.g., county assessor records) include relevant information, but these attributes are often only discernable by visual observation of the property. Structural attributes can be changed by residents and/or property owners through renovation or replacing materials, although initially building with less combustible materials tends to be significantly less costly [44].

#### *2.5. Overall Risk*

Notwithstanding its emphasis on understanding and communicating the individual attributes, the WiRe RA enables calculating an overall parcel-level wildfire risk score. Each ¯ possible response for the thirteen attributes is assigned a point value based on subjective expert judgment following the fire science described above, with the weighted sum providing an overall wildfire risk score (see Tables 1 and 2). For example, out of 1000 possible points, a wood shake-shingle roof is worth 300, whereas insufficient driveway clearance is worth 10, reflecting the general understanding that a shake-shingle roof influences risk much more than driveway clearance does. This weighting system can be adjusted to specific programmatic needs and local contexts, and indeed, the authors have worked with numerous programs to adjust attributes and their weights accordingly. The "Assessor Reference Guide" (ARG) presented in Figure S1 in the Supplementary Materials provides rationale and additional considerations for attributes and response categories as collaboratively determined for the WiRe RA project in Grand County Colorado. ¯

**Table 2.** Description of summary scores calculated from the Wildfire Research (WiRe) Rapid Assess- ¯ ment and statistics (mean, Moran Test *p*-value) for structures in the Columbine Lake community, Grand Lake, CO with complete assessments, by whether structures were destroyed in the East Troublesome Fire or not. A Benjamini–Hochberg procedure [45] suggests that only *p*-values of *p* < 0.001 on this table (bolded) are significant after adjusting for multiple comparisons with an assumed 10% false discovery rate.


Similarly, weighted sums can be calculated for each of the four categories contributing to PLR shown on Figure 1 (i.e., parcel-level hazard, defensible space, access, and hazard) by summing the point values for included attributes. Attributes and point values have been developed iteratively through more than 15 years of collaborative partnerships between the WiRe Team and wildfire practitioners throughout the western United States. Although the ¯ WiRe RA aggregates individual attributes into overall risk using a simple weighted sum, ¯ the literature suggests that some attributes may interact. For example, flammable siding has been found to pose less risk with proper defensible space [18], and homes on steep slopes have been found to benefit from greater defensible space distances than homes on shallower slopes [22]. While the linear sum can be modified to represent such interactions, we believe that doing so in any rigorous way, including generalized across different fuel and fire behavior contexts, would require a substantially more nuanced understanding of the interactions than currently exists.

The WiRe RA is focused on identifying parcel-level heterogeneity within a community. ¯ Thus, it typically does not explicitly include information from broader-scale assessments, such as wildfire behavior modeling that relates to hazards at the landscape scale or social vulnerability assessments that relate to vulnerability at the community scale. However, we acknowledge that comprehensive assessment of the wildfire risk to homes and residents entails consideration of the broader-scale contexts in which the community, and thus the PLR, is embedded (see Figure 1, light green boxes with dashed borders). Depending on purposes, the assessment of risk might also benefit from additional considerations not addressed here, such as the potential for parcel-level variation in typical weather conditions or in weather conditions expected during extreme events.

Finally, we note that the concept of PLR expands upon the related concept of the home ignition zone (HIZ), which emphasizes that reducing the susceptibility of homes to wildfire, by appropriately managing the materials, design, and maintenance of the home in relation to the immediate surroundings, is key to reducing the risk of home loss [33,46]. While the HIZ concept focuses on susceptibility, PLR further encompasses localized wildfire hazard and exposure to provide a more comprehensive assessment of risk at the parcel level. This accommodates the inclusion of not only structural conditions and defensible space but also other localized influences on the wildfire's likelihood and intensity (i.e., parcel-level hazard) and conditions that interface with suppression and preparation (i.e., access). PLR also recognizes the role of defensible space in supporting safe and effective fire suppression, which relates more to the localized wildfire hazard than to the susceptibility of a home. Most importantly, these modifications recognize that the interventions possible for one actor, at one scale, can interact with the interventions possible for other actors and at other scales; for example, a resident's decision to widen their driveway can improve a fire department's ability to safely protect that resident's home.

#### **3. Materials and Methods**

#### *3.1. Study Context*

Columbine Lake is a small community approximately two miles away from downtown Grand Lake, CO, a gateway town on the western side of the Rocky Mountain National Park (see Figure 2). In October of 2020, this community was one of many in Grand County, CO, affected by the East Troublesome Fire. The fire was reported in the Arapaho National Forest on 14 October and spread to over 10,000 acres in the first three days. The fire then grew rapidly, covering 187,964 acres by the end of 23 October. The fire destroyed an estimated 366 residences and 214 outbuildings and commercial structures [47] and damaged many more. Despite the destruction of 26 out of 461 assessed structures in the community of Columbine Lake, the community received somewhat of a glancing blow from the fire, in comparison to other communities in the area that received the full brunt of the extreme fire conditions and consequently faced near total destruction of their structures. Note that the number of destroyed structures reported throughout this paper pertains to complete losses among structures with complete rapid assessment data from 2019; this number does not include numerous structures that sustained damage (e.g., burns, smoke damage) but were not considered complete losses, nor does it include structures destroyed by the fire but not covered by the 2019 rapid assessment. This total does include 3 structures with only partially complete assessments, as described in more detail below.

**Figure 2.** Map of East Troublesome Fire perimeter [48] and structures in the community of Columbine Lake, Grand Lake, CO that were assessed by Grand County Wildfire Council in 2019 and either destroyed (red) or not (blue) in the East Troublesome Fire. Inset shows the location within Colorado. Basemap image is the intellectual property of Esri and is used herein under license. Copyright © 2022 Esri and its licensors. All rights reserved.

ē

ē

ē

ē

ē ē

#### *3.2. Data*

Previously, GCWC had selected Columbine Lake as one of six communities to be the site of a detailed data collection effort in collaboration with the Wildfire Research (WiRe)¯ Center. (The Wildfire Research (WiRe) Center (https://wildfireresearchcenter.org, accessed ¯ on 11 February 2022) is a non-profit organization that works with wildfire practitioners to seek locally tailored pathways to create fire-adapted communities through the WiRe¯ Approach described above and by ref. [19].) As part of this effort, GCWC used the WiRe¯ RA tool to assess 461 homes in the community of Columbine Lake in the summer of 2019. Data collection was conducted by individuals from each of the five fire protection districts who had been trained to conduct the WiRe RA. The assessment was conducted from the ¯ roadside and supplemented using online imagery where necessary. As is common in roadside assessments, numerous attributes could not be clearly observed for all parcels and were marked as unobservable in the WiRe RA dataset. While no more than three ¯ parcels were marked as unobservable for 10 out of the 12 attributes, two attributes ("other combustibles" and "attachments") were frequently unobservable (see Table S1 for details). Omitting parcels with unobservable attributes truncates the dataset by 24 percent, down to a total of 352 fully assessed structures, 23 of which were destroyed in the fire. To avoid introducing unnecessary measurement error bias, we used this set of 352 structures with complete assessment data for our main analysis. However, because the large proportion of incomplete assessments also introduced gaps in the spatially lagged measures, Table S2 replicates our primary results for the full dataset with unobservable attributes coded as having the riskiest rating for that attribute. Because proportionally fewer destroyed structures had incomplete data versus structures not destroyed (12% vs. 24%), this coding was conservative with respect to identifying meaningful attributes for explaining structures destroyed in the East Troublesome Fire. The results shown in Table S2 are generally consistent with the main results presented below.

GCWC assessed twelve attributes in 2019: the eleven represented by attributes 1–3 and 5–12 in Table 1, plus an additional attribute describing the length of the driveway that was omitted from the analysis here due to a lack of variation. However, the WiRe¯ Team continues to revise the assessment tool as necessary based on new understandings. Based on revisions to the standard WiRe RA implemented in March 2020 (i.e., before the ¯ East Troublesome Fire), we supplemented the WiRe RA data collected by GCWC with ¯ an estimate of the distance between structures (attribute 4 in Table 1). We calculated the distance to the nearest home using point data for the locations of assessed structures. To adjust for point versus polygon data, 40 feet were subtracted from all calculated values (based on an assumed 20 feet from centroid to exterior wall for each of two homes) before converting to categorical ratings for inclusion in the overall risk score. We then calculated the overall risk score using the relative point values agreed upon by wildfire mitigation experts on the WiRe Team on March 2020, shown in Table 1. ¯

Table 1 provides descriptive statistics for data pertaining to structures with complete assessment data destroyed and not destroyed within the Columbine Lake community. After the East Troublesome Fire, GCWC identified which structures within the WiRe RA ¯ dataset for Columbine Lake had been destroyed, as shown on the map in Figure 2. For reference, we overlaid spatial data on the fire perimeter [48], along with a 50-foot buffer to accommodate imprecision in spatial data, upon a point layer representing the locations of all assessed structures. Although 21 out of 23 structures with complete assessments were identified as within the burn perimeter, we note that all structures in the community were within a quarter mile or less from this active burn perimeter and that the entire community was exposed to wind-driven ember showers and active suppression activity during the events of 21 October 2020. Further, we note that although the fire perimeter is influenced by fire behavior at broader, landscape-level scales, the location of the perimeter within the community is not exogenous, but rather is at least partially simultaneously determined by structure destruction, due to parcel-level conditions, suppression activity, and because structures themselves can act as a source of fuel for the fire. Accordingly, we considered

the entire community as being exposed to the fire and estimate our main models using the full dataset. To test the robustness of this decision, we replicated our main models with a subset of the data, constrained to the 116 structures within the burn perimeter; the results (shown in Table S3) are similar in direction to those from the preferred models but generally greater in magnitude while being considerably less precise.

#### *3.3. Empirical Analysis*

We focused our primary analysis on all structures within the community with complete WiRe RA data, modeling the influence of assessed variables and summary scores on ¯ whether a structure was lost in the East Troublesome Fire. Although, ideally, we would model the likelihood of a structure being destroyed as a function of all assessed attributes jointly, the small number of observations limited the degrees of freedom to do so. Instead, we estimated a series of separate models, each with one assessed attribute or summary score as the main independent variable.

Past research [49], practitioner experience, and the spatial patterns apparent in Figure 2 all suggested the possibility of spatial clustering of not only fire outcomes but also many of the assessed attributes. We tested and found strong evidence for spatial dependence in whether structures were destroyed (Moran test χ <sup>2</sup> = 187.96, *p* < 0.001) and in data for seven of the 12 assessed attributes (last column of Table 1). Thus, we estimated relationships using a series of spatial Durbin models [50], given by:

$$\begin{aligned} y &= \mathfrak{a} + \beta \mathfrak{X} + W \mathfrak{X} \mathfrak{y} + \rho Wy + \mathfrak{e} \\ \mathfrak{e} &\sim N\left(0, \sigma^2 I\_{\mathfrak{U}}\right) \end{aligned} \tag{1}$$

In all models, the dependent variable was defined as *y* = 1 if the assessed structure was destroyed and *y* = 0 otherwise (although the binary outcome variable suggests that a logit or probit specification would be more appropriate, the theory and implementation of a spatial Durbin logit/probit model is not well established. Regarding the main limitation of the linear probability model approach, Woodridge [51], (p. 455) notes, " . . . If the main purpose is to estimate the partial effect of [the independent variable] on the response probability, averaged across the distribution of [the independent variable], then the fact that some predicted values are outside the unit interval may not be very important." See Table S4 for a comparison of an ordinary versus logistic regression for all attributes in the dataset, without controlling for spatial effects). Each model included: a constant for the intercept (*α*); a slope (*β*) multiplied by a sole independent variable (*X*) (corresponding to the numerical score for one of the twelve assessed attributes shown in Table 1 or the five summary scores shown in Table 2); a spatial weights matrix (*W*), calculated as the inverse-distances among the centroids of all assessed structures; a spatial lag (*γ*) for the weighted independent variable *WX*; a spatial lag (*ρ*) for the weighted dependent variable (*Wy*); and an error term (*ǫ*) consisting of a normally distributed mean-zero disturbance with constant variance (*σ* 2 ) multiplied by the *n*-dimensional identity matrix (*In*), where *n* is equal to the number of observations.

We estimated models in Stata/SE 16.0 (any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government) using the *spregress* command, employing the generalized spatial two-stage least-squares estimator and treating errors as heteroskedastic. Because of the spatial relationships, parameter estimates for the spatial Durbin model do not have straightforward interpretations and only make intuitive sense in terms of combined effects. Accordingly, we reported and focused on the estimated impacts, as defined by LeSage and Pace [52], which summarize the average marginal effects of the independent variable on the reduced form mean and are calculated with the *estat impacts* command in Stata. The reported impacts included direct impacts, which report the average of the own-marginal effects of the independent variable; indirect impacts, which report the average of the spillover effects of the independent variable, and total impacts, which are the sum of the other two and thus of primary interest.

That is, direct impacts reflect the average of the effects of each property's attributes directly on whether that property's own structure was destroyed, ignoring any spillover effects across properties. Indirect impacts can be interpreted as reflecting the average of either the spillover from one property's attribute upon the outcomes for all its neighbors or the spillover from all of a property's neighbors' attributes upon the outcome for that property. As the sum of direct and indirect impacts, total impacts thus reflect the average overall effect of a property's attributes on whether that property's structure was destroyed, accounting for the interactions amongst neighbors.

#### **4. Results: Risk Assessment Data Help Explain Destroyed Structures**

In this section, we investigate the relationships between WiRe RA data and whether ¯ homes in the Columbine Lake community were destroyed by the East Troublesome Fire using the modeling approach described above. Tables 3 and 4 show the results for seventeen separate spatial Durbin models for structures in the Columbine Lake community with complete assessments.

Table 3 reports estimated parameters for each model. In all models, the coefficient on the spatial lag on the dependent variable (*Wy*) is large, positive, and strongly significant, consistent with the spatial clustering of structures that were lost and with the possible effect of structures acting as a source of fuel and increasing the hazard to their neighbors during the event. The strong significance of the spatial lag on the independent variable (*WX*) in nearly all models provides further support for estimating spatial effects with the spatial Durbin model.

**Table 3.** Coefficients estimated for 17 separate spatial Durbin models, one for each independent variable listed, of whether a structure was destroyed (y) as a function of the independent variable from the Wildfire Research (WiRe) risk assessment (RA) for structures with a complete assessment ¯ (*n* = 352; *y* = 0.07). Shading depicts grouping of attributes into categories (i.e., 1–4: parcel-level hazard; 5–6: defensible space; 7–9: access; 10–12: structure). A Benjamini–Hochberg procedure [45] suggests that all *p*-values of *p* = 0.077 or less on this table (bolded) are significant after adjusting for multiple comparisons with an assumed 10% false discovery rate.


**Table 4.** Impacts estimated from 17 separate spatial Durbin models, one for each independent variable listed, of whether a structure was destroyed (y) as a function of the independent variable from the Wildfire Research (WiRe) risk assessment (RA) for structures with a complete assessment ( ¯ *n* = 352; *y* = 0.07). Shading depicts grouping of attributes into categories (i.e., 1–4: parcel-level hazard; 5–6: defensible space; 7–9: access; 10–12: structure). A Benjamini–Hochberg procedure [45] suggests that all *p*-values of *p* = 0.021 or less on this table (bolded) are significant after adjusting for multiple comparisons with an assumed 10% false discovery rate.


That said, as described above, coefficients on the independent variable (*X*) and its spatial lag (*WX*) cannot be separately interpreted in an intuitive way. Thus, we focus on the average impacts shown in Table 4 as our main results. Each row reports the average total impact, the average direct impact, and the average indirect impact, with the average estimated impact, standard error, and *p*-value shown for each measure. Impacts pertain to marginal effects, meaning that they describe the average change in the probability of the structure being destroyed for a 1-point change in the independent variable, calculated at the reduced form mean. As shown in Table 1, attribute ratings are scored anywhere from 5 to 300 points, depending on the attribute, such that higher points always pertain to higher expected risk.

Table 4 shows that eight of the twelve individual attributes measured by the WiRe RA ¯ are estimated as having positive and significant total impacts on the likelihood of a structure being destroyed. Strikingly, none of these variables are estimated as having significant direct impacts upon a given structure, whereas all of the assessed attributes with significant total estimates are estimated as having significant indirect impacts. Indirect impacts can be interpreted in two ways: the perspective *from* an observation, which relates to how a change in a single parcel influences the risk to all other structures, or the perspective *to* an observation, which relates to how changes in all other parcels influence a single parcel [52]. From either perspective, the results imply strong spillovers in risk across properties, in which the hazardous conditions on individual parcels play a large role in determining the risk to their neighbors as well as to themselves through the interactions of fire risk between structures.

Given the small sample size, the linear probability model specification, and the role of expert-assigned weights in determining parameter magnitudes, we focus primarily on the direction and significance of the estimated average total impacts for interpretation. Among the four variables for the parcel-level hazard (1–4), adjacent fuels and the distance to nearest home are estimated as positive and significant. Both defensible space attributes, pertaining to vegetation (5) or other combustibles near the home (6), have positive total

effects on risk. Driveway clearance (8) and address visibility (9), two out of three attributes related to access, are positive and significant, as are both siding (11) and attachments (12) for the structure category. In addition, all five summary measures (13–17), including the overall risk score (variable 17), are also found to have a positive and significant total impact. The mean estimated impact of 0.0512 from the (scaled) overall risk score variable suggests an approximate 5% increase in the likelihood of a structure being lost for every 100-point increase in risk, with the result coming from the combined direct and indirect impacts. For further intuition, Figure 3 depicts the distributions of overall wildfire risk scores for structures destroyed versus not destroyed within the entire community (a) as well as within the fire perimeter (b). The figure suggests a meaningful relationship in both cases, with destroyed structures on average having higher overall risk scores.

−

**Figure 3.** Comparison of the distribution of overall risk score for structures with complete assessments that were destroyed (orange) versus not destroyed (blue) in the East Troublesome Fire, for the entire community of Columbine Lake, Grand Lake, CO (**a**) and within the fire perimeter (**b**).

Finally, Tables S5 and S6 replicate the main results of Table 4 for all available indicators or for binary indicators equal to 0 for attributes rated at the lowest risk level and equal to 1 otherwise, respectively. These results do not rely on the relative weights assigned for calculating the summary risk scores. Although nuanced differences exist, these replications demonstrate the general robustness of the implied linear response to the expert-judgment score values assigned to attribute levels. In particular, the results for the binary measures for the measured attributes, shown in Table S6, are quite strong and generally analogous to the main results of Table 4.

#### **5. Discussion**

The results suggest that the WiRe RA provides strong explanatory power for whether ¯ structures in the Columbine Lake community were destroyed in the East Troublesome Fire. The overall wildfire risk score from the WiRe RA strongly relates to whether a structure was ¯ destroyed. Because the summary risk score from the WiRe RA is explicitly constructed to ¯ represent overall PLR, this suggests that PLR meaningfully describes variation in wildfire

risk across homes within a community. However, it is also notable that the distributions for overall risk scores overlap substantially for structures destroyed versus not destroyed, demonstrating that although the WiRe RA data provide insights for whether structures ¯ were destroyed by the fire or not, the results are far from deterministic. Some structures with relatively very low risk scores were destroyed in the fire, while other structures with relatively high risk scores were spared, and this is true both within the fire perimeter and across the entire community.

Given the additional attribute- or category-specific results discussed above, and the fact that results pertaining to the overall risk score rely on the subjective expert judgment of relative attribute weights, we maintain that a single summary measure of overall risk does not tell the entire story, but rather, it can mask important—and actionable—information. All four category-specific summaries have significant total impacts on the likelihood that a structure was destroyed in the East Troublesome Fire, suggesting that each of these four categories represents a relevant component of PLR. This suggests numerous points of entry for reducing PLR, including common recommendations related to defensible space and structural hardening as well as access considerations. Most of the individual attributes measured in the WiRe RA relate to whether structures were destroyed or not, ¯ offering nuanced insights. For example, one of the most common recommendations for reducing wildfire risk to a home is the maintenance of a defensible space around the home; however, details of recommendations vary. Recent guidance, e.g., [53], often emphasizes the importance of the nearest zone, within approximately 5 feet of the structure. Past research from other contexts found no measurable differences in risk reduction beyond about 100 feet [22] or 130 feet [35] from the home. Here, although our main results generally find increasing risk as the distance to dense or unmaintained vegetation decreases, the results estimated with separate indicator variables suggest similar increases in risk for any measured amount of defensible space less than the maximum observation of 150 feet. In other words, our results suggest that effective defensible space possibly extends to at least 150 feet from the structure. Our results also suggest the importance of removing other combustible items, including but not limited to flammable furniture, woodpiles, and propane tanks, from the vicinity of the home. Indeed, our binary-indicator model suggests a similar increase in the likelihood of a home being destroyed due to having less than 30 feet of distance between any such items and the structure as compared to having less than 150 feet of distance from maintained vegetation.

We also find strong evidence of the interdependence of risk among structures, in which the decisions made by one's neighbors can affect the risk to one's own property [49,54]. Evidence includes the significant impact of the proximity to the nearest home attribute, which is consistent with other research that found the distance between structures to be an important component of wildfire risk to homes, e.g., [29,31]. Evidence also comes from the spatial correlations in whether structures were destroyed, demonstrated robustly across all estimated models, which suggests that burned structures tend to be clustered. Further evidence for risk interdependencies comes from the significantly estimated average indirect impacts for many attributes, which signify risk spillovers in terms of the attributes' influence on properties. This all underscores the importance of considering homes and other structures as not only *recipients* of fire, but also as *drivers* of wildfire behavior. These results also suggest important limitations for the applicability of landscape-level hazard models to a populated built environment, due in part to the practice of masking "urban" or "developed" pixels as "unburnable" terrain in the fuel models that underly such modeling [55,56].

In addition to missing heterogeneity in parcel-level hazard or the local fuel characteristics represented by the defensible space category, landscape-level models that treat all residential property the same neglect the heterogeneity in vulnerability that is driven by differences in access or structure characteristics. Our results are consistent with the expectations set by the literature reviewed above in demonstrating the value of structural hardening for reducing risk, such as having less vulnerable siding materials and avoiding

combustible attachments such as decks, fences, or porches. We also find that two of the access-related attributes, driveway clearance and address visibility, have significant impacts on the likelihood that a structure was destroyed.

Because it cannot be assessed with post-fire data, the role of address visibility is worth considering further. Our models suggest that whether an address was visibly posted or not is significantly related to the likelihood that a structure was destroyed in the fire. As discussed above, it is quite plausible that address visibility supports a safe and effective response to a given property, and the strength of this result suggests value from further investigation. However, operational memory suggests that address visibility was not particularly important in the conditions of this fire, and thus this attribute might serve as a proxy for variation in risk reduction effort not otherwise reflected in the coarse measurements of the rapid assessment. For example, residents who increase their address visibility might also be more likely to remove debris from porches and gutters. Either way, this result suggests that small actions taken by residents can make a real difference in their wildfire risk. Further, these results reflect the potential importance of the complex dynamics of a wildfire event, such as the interactions between parcel-level characteristics and suppression decisions, in determining the realized wildfire risk.

Finally, decisions about wildfire risk mitigation actions are influenced by many factors that often interact in complex ways, including but not limited to risk perceptions, costs, and perceived effectiveness, e.g., [44,57–60]. Because rapid assessments represent parcel-level heterogeneity in risk within a community, rapid assessments can inform the development and implementation of programs that encourage or support these decisions. For example, an organization might decide to use rapid assessment data to prioritize properties for more time-intensive detailed site assessments or for more costly cost-share incentive programs. Of course, comprehensive parcel-scale assessments may face challenges for any given community in terms of practicality, cost-effectiveness, and coverage. Practitioners who are interested in utilizing such approaches must consider the tradeoffs involved in implementing such assessments and whether the benefits provided, in terms of the potential uses, warrant the financial and time costs to their organizations.

#### **6. Conclusions**

The WiRe RA conducted by GCWC in the summer of 2019 captured elements of ¯ the risk of a home being destroyed during the East Troublesome Fire of October 2020. Despite small sample sizes and the coarse measurements underlying the WiRe RA, the ¯ overall risk score constructed from the weighted sum of assessed attributes meaningfully depicts heterogeneity in the risk that the East Troublesome Fire presented to homes in the Columbine Lake community. Furthermore, seven out of the twelve attributes assessed in the original WiRe RA, including attributes from each of the four categories (i.e., parcel- ¯ level hazard, defensible space, access, and structure), help to explain which homes were destroyed. The results also support incorporating proximity to the nearest home into parcellevel wildfire risk assessment, as that attribute is also strongly related to the likelihood of whether a home was destroyed.

Thus, the parcel-level conditions described by the PLR conceptual model are important for understanding and communicating about wildfire risk to homes. These results are consistent with other recent works, e.g., [14,33,46], in arguing for the importance of conditions within the home ignition zone for understanding wildfire risk to homes. That said, this paper reports on an empirical analysis of a unique event; we do not claim that the attributes found to be meaningful here are the most important for reducing risks, nor that these results should supersede those found from other contexts, whether in terms of geographic, social, or fire conditions. In terms of generalization, the details of which specific attributes mattered in this case study are less important than the basic fact that some property-level attributes were found important. Since parcel-level characteristics as measured by the WiRe RA tool matter here, there is reason to believe they might matter ¯ in the next fire. These findings suggest that having this type of data for all parcels in a

community could help community programs to identify where to start, or how best to prioritize resources, for reducing wildfire risk to homes and their residents.

Our results underscore that parcel-level variation is an important consideration for understanding risk to homes for homeowners, communities, fire managers, and policy makers. Indeed, our results suggest that parcel-level WiRe RA data offer insights not ¯ available from landscape-level hazard modeling outputs. In other words, landscapelevel hazard information cannot replace detailed parcel-level information in terms of understanding wildfire risk to a home. This suggests the importance of working toward the greater integration of parcel-level variation in risk-related characteristics (i.e., about PLR) into existing landscape-level wildfire risk assessments. Given the robustness of the results on the distance between structures here and in the past literature, a promising first step could be more consistently including such a measure—assessed at the level of individual parcels, rather than simply as an average housing density measure—into broader-scale assessments.

The results also demonstrate the usefulness of rapidly collected, coarse data in representing a property's wildfire risk. Measurement for the WiRe RA is imperfect, based on ¯ professional judgment with at most four levels for a given attribute, and conducted from the roadside to accommodate privacy concerns and resource constraints. The overall risk score is based on subjective expert judgment pertaining to the relative weights of different attributes. It is reasonable to ask whether such imperfect data are useful for planning or prioritization; our results suggest that they can be. As our results are robust to the use of binary versus expanded indicator variables for each attribute, future research might explore the potential for efficiency gains through further simplifying current approaches to parcel-level rapid wildfire risk assessment, such as employing binary rather than three- or four-leveled indicators.

That said, the lack of relationship found here of course does not imply that a variable might not be critically important in a different context. For example, numerous studies support the importance of noncombustible roofing materials for reducing risk despite the lack of observed importance in our case study—but only two assessed structures in our dataset had wooden roofs; one of them burned down and the other did not. Our estimation approach cannot overcome the limited support for many of the possible measurement levels in the rapid assessment. Most of our estimated average effect sizes are relatively low, reflecting the fact that numerous other considerations not represented by the data here also contribute to whether a home may be destroyed. As parcel-level rapid assessment data become more widely collected, there will likely be more opportunities to investigate the effectiveness of this type of data in representing wildfire risk to homes; we hope such unfortunate intersections of past data collection and hazardous events can be seized upon to reduce the social costs of wildfires.

Overall, this study finds that resident or homeowner actions matter for reducing wildfire risk to the home. While the resident of a property might have limited capacity to change conditions underlying a hazardous burn probability or conditional flame length rating, that same resident likely can influence many of the attributes measured by the WiRe RA. One of the many influential attributes in our model—address visibility—is also ¯ an easy, low-cost attribute to change. Furthermore, action at the parcel and community scales constitutes a critical piece of the story for assessing the full wildfire risk to homes, even though landscape-level wildfire hazard assessments might not be able to capture these efforts. As such, assessed risk levels should not be taken as given at the scale of an individual property; rather, residents can reduce their own risk, and that of their neighbors, by focusing on aspects of PLR over which they have some control. We hope such information empowers wildfire risk practitioners, and the residents they serve, to take meaningful actions to reduce the wildfire risk to homes and communities.

**Supplementary Materials:** The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/fire5010024/s1, Table S1: Number of structures for which attributes were marked as "unobservable" in the Wildfire Research Risk Assessment (WiRe RA); Table S2: ¯ Replication of Table 4 using full dataset with unobservable attributes for structures coded as having the riskiest rating for that attribute instead of being coded as missing and dropped from subsequent analysis; Table S3: Replication of Table 4 constrained to structures within the burn perimeter; Table S4: Comparison of ordinary versus logistic regression for all assessed attributes (jointly modeled), without considering spatial effects; Table S5: Replication of Table 4 with full indicators for the levels of each attribute of a structure, rather than using the numerical score as an implied linear measure; Table S6: Replication of Table 4 with binary indicators for each attribute of a structure, coded to 0 for the lowest-risk level and 1 for all others, rather than using the numerical score as an implied linear measure; Figure S1: Assessor Reference Guide (ARG) developed for Grand County WiRe RA. ¯

**Author Contributions:** Conceptualization, analysis, visualization: J.R.M.; methodology: J.R.M., C.M.B., S.K.O., A.C.G., H.B.-S., P.A.C., J.G.; writing—original draft preparation: J.R.M., J.B.G.; writing—review and editing: all. All authors have read and agreed to the published version of the manuscript. The findings and conclusions in this article are those of the authors; they do not necessarily represent the views of the USDA and should not be construed to represent USDA agency determination or policy.

**Funding:** This research was funded, in part, by USDA Forest Service, State and Private Forestry. Rapid assessment data and data on destroyed structures were funded by Grand County Wildfire Council with the support of the Fire Districts of Grand County.

**Data Availability Statement:** Rapid assessment data and data on destroyed structures are owned by Grand County Wildfire Council (https://bewildfireready.org/, accessed on 11 February 2022) with the support of the Fire Districts of Grand County and are available upon request; other data are available from NIFC [47] as described in the text.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results. The Wildfire Research (WiRe) Center (https: ¯ //wildfireresearchcenter.org, accessed on 11 February 2022) is a nonprofit organization that works with wildfire practitioners and collaborated with Grand County Wildfire Council for data collection; the author Brenkert-Smith is a member of the Wildfire Research (WiRe) Center Board, and the authors ¯ Meldrum, Champ, and Barth are on the Wildfire Research (WiRe) Advisory Committee. ¯

#### **References**


## *Article* **Mitigating Source Water Risks with Improved Wildfire Containment**

#### **Benjamin M. Gannon 1,2, \*, Yu Wei <sup>1</sup> and Matthew P. Thompson 3**


Received: 21 July 2020; Accepted: 20 August 2020; Published: 21 August 2020

**Abstract:** In many fire-prone watersheds, wildfire threatens surface drinking water sources with eroded contaminants. We evaluated the potential to mitigate the risk of degraded water quality by limiting fire sizes and contaminant loads with a containment network of manager-developed Potential fire Operational Delineations (PODs) using wildfire risk transmission methods to partition the effects of stochastically simulated wildfires to within and out of POD burning. We assessed water impacts with two metrics—total sediment load and frequency of exceeding turbidity limits for treatment—using a linked fire-erosion-sediment transport model. We found that improved fire containment could reduce wildfire risk to the water source by 13.0 to 55.3% depending on impact measure and post-fire rainfall. Containment based on PODs had greater potential in our study system to reduce total sediment load than it did to avoid degraded water quality. After containment, most turbidity exceedances originated from less than 20% of the PODs, suggesting strategic investments to further compartmentalize these areas could improve the effectiveness of the containment network. Similarly, risk transmission varied across the POD boundaries, indicating that efforts to increase containment probability with fuels reduction would have a disproportionate effect if prioritized along high transmission boundaries.

**Keywords:** water supply; erosion; wildfire containment; Potential fire Operational Delineations; Monte Carlo simulation; transmission risk

#### **1. Introduction**

Improved wildfire containment is an attractive strategy to mitigate the risk of degrading water quality beyond limits for treatment because of the potential to limit fire sizes and impacts to tolerable levels without the need to completely exclude fire from the landscape. Recent efforts to make containment planning more proactive, focus on zoning the landscape into fire management units called Potential fire Operational Delineations (PODs) using existing high probability control features such as roads, rivers, and fuel transitions [1,2]. Beyond the inherent value of engaging managers in the process to identify and critique potential control features, the resulting POD areas become relevant spatial units for pre-fire analysis of endogenous and transmitted wildfire risk to inform response strategies that are appropriate for the predicted direction and magnitude of fire effects to water supplies and other natural resources and human assets [2]. While there has been substantial progress engaging managers in the bottom up approach to develop and employ PODs and their associated response strategies [2–8], less attention has been paid to evaluating the risk mitigation effectiveness of containing wildfire within these units and what functional improvements should be made to the size and spatial arrangement

of the containers to maximize their protection benefit for water supplies and other values, such as wildlife, that depend on the scale of fire activity.

Wildfire is often harmful to water quality because reductions in surface cover and infiltration cause increases in surface runoff and erosion that can mobilize and transport contaminants into surface drinking water sources [9–12]. While the specific contaminants and concentrations of concern may vary by watershed and water system [12–14], water quality degradation generally becomes problematic when large quantities of sediment are mobilized by intense rainfall causing contaminant concentrations to exceed thresholds for effective water treatment (e.g., [15]). Post-fire sediment loads are influenced by fire size and burn severity, topography, soil properties, and rainfall intensity [10,16,17]. Previous efforts to account for fire effects on watersheds and water supplies account for some of these factors [18,19], but the use of relative fire effects measures makes it difficult to evaluate whether a given fire will degrade water quality. This shortcoming has been addressed in recent years with increasing use of spatially explicit erosion and sediment transport models to make quantitative predictions of sediment yield from modeled wildfires (e.g., [20–24]). Sediment yield models have been widely used to examine the risk mitigation effectiveness of area-wide fuel treatments meant to reduce burn severity [23–26] but they have not yet been used to evaluate the performance of fire containment strategies to reduce area burned.

Some water systems have discrete features, such as terminal reservoirs, that could be targeted for protection within a single POD, but many municipal watersheds in the western USA are hundreds to thousands of square kilometers in size and therefore require some level of internal compartmentalization to protect water supplies. In theory, the size and spatial arrangement of PODs could be designed to mitigate the risk of water quality degradation by both containing fires with potential for large growth and subsequent contaminant loads near their ignition sources and ensuring that within-POD burning does not result in adverse consequence. Managers consider both values at risk and presence of control features when delineating PODs, which often results in smaller PODs near developed areas and larger PODs in the backcountry [2,4]. However, it is not clear that the size and configuration of manager-delineated PODs will reduce risk of wildfire-related water quality degradation. Several attempts have been made to automate the processes of identifying suitable control features and aggregating them into PODs [27,28] using roads, streams, watershed boundaries, and spatial models of suppression difficulty and potential for control [29–31], but data-driven approaches have yet to inform the desired size and spatial configuration of PODs to mitigate a particular risk.

Recognizing the importance of fire size, location, and burn severity for watershed response, several previous studies have employed Monte Carlo wildfire simulation to characterize watershed exposure and water supply risk [19,32–34]. Their results suggest that most risk to water supplies is associated with a small subset of total fire activity. Moreover, the source locations of damaging wildfires tend to cluster in certain parts of the landscape, which implies containment benefits will depend strongly on location. Simulated fire ignition locations and fire extents can be intersected with relevant management units to partition fire impacts from burning within the unit of origin and transmission to the surrounding landscape [2,35,36]. Analyzing risk transmission across a network of PODs could help to identify locations with high source risk that would benefit from investment in activities to improve containment probability, such as roadside fuels reduction. Areas with fuels conducive to fast fire spread tend to transmit the most fire [36], which will result in high water supply risk when adjacent areas have high erosion potential and/or short transport paths to water supplies. Analysis of water supply risk from self-burning could also identify high risk PODs that would benefit from further compartmentalization.

The goal of this study is to provide a proof of concept model to evaluate the effectiveness of a containment network at mitigating risk of source water quality degradation. The general approach should also be relevant for assessing risk to other resources that depend on disturbance size. We utilized Monte Carlo wildfire simulation, erosion, and sediment transport modeling to quantify the potential water supply impacts from a set of simulated wildfires with and without containment. We analyzed risk and risk mitigation with two measures of water supply impact—total sediment load and frequency of exceeding turbidity limits for treatment—to highlight how considering the scale-dependent effects of wildfire changes the perceived mitigation value of fire containment. Risk transmission was analyzed to identify possible improvements to the containment network with measures of transmitted risk highlighting those PODs and POD boundaries that could benefit from activities to improve containment probability and measures of self-burning indicating areas in need of further compartmentalization.

#### **2. Materials and Methods**

#### *2.1. Evaluation Framework*

The evaluation framework was designed to contrast the water quality impacts of uncontained wildfires and wildfires contained within the POD of origin in terms of total sediment load and average post-storm suspended sediment (Figure 1). Total sediment load is similar to the commonly used net value change measures in risk assessment [37,38] insomuch as more is interpreted as bad and any marginal reduction decreases risk. However, using change in total sediment load as a measure of risk has the potential to falsely assign mitigation benefit to containment when either the load from the uncontained wildfire is already below a meaningful threshold of water quality degradation or containment reduces erosion but the resulting load is still above the treatment threshold. Average post-storm suspended sediment concentration is used here to estimate whether fires will degrade water quality beyond limits for water treatment and whether degradation outcomes change with containment. This measure of risk better approximates the threshold-dependent nature of water quality degradation owing to the size of the receiving waterbody and the water system sensitivity to contaminants.

**Figure 1.** The evaluation framework focuses on total sediment yield and average post-storm suspended sediment as measures of water quality degradation risk. Variable inputs are in light grey. Stochastically simulated wildfire perimeters were combined with estimates of burn severity to model post-fire erosion and sediment transport to the water supply both with and without containment. Sediment yield was converted to average post-storm suspended sediment concentration using the receiving waterbody volume and the annual frequency of sediment generating storms.

Our evaluation framework focuses on the key uncertainties in wildfire-water quality degradation risk related to the extent of the watershed burned and post-fire rainfall (Figure 1). As further described in following sections, many plausible wildfire perimeters were simulated with the Monte Carlo ignition and spread model RANDIG [35,39], which were then clipped to their POD of origin to approximate a strategy of improved containment. Post-fire erosion was then simulated for each perimeter using crown fire activity predicted with FlamMap 5.0 [40] as a proxy for burn severity to modify the cover and soil variables in the Revised Universal Soil Loss Equation (RUSLE) [41]. We accounted for uncertainty in post-fire rainfall by modeling erosion for three rainfall scenarios ranging from common to extreme. We estimated annual sediment loads to the water supply based on the predicted proportion of sediment transported off hillslopes and through channels using Sediment Delivery Ratio (SDR) models [42,43]. Post-storm suspended sediment concentrations were estimated by assuming average storm sediment loads are diluted in the mean daily flow volume of the river during the May to October thunderstorm season, which is associated with most post-fire erosion and water quality degradation in the study region [15,16]. All analyses were completed with R version 3.5.3 [44] except where noted otherwise.

#### *2.2. Study Area*

The study area encompasses 3021 km<sup>2</sup> of the Front Range Mountains in Colorado, USA (Figure 2). The Front Range has a history of large and severe fires that have caused extreme erosion, reservoir sedimentation, and water quality degradation [15,45–48]. The names of the focal municipal water supply and other geographic features within the study area are withheld for security reasons. The extent of the study area was defined to include the contributing area to a municipal pipeline diversion (1254 km<sup>2</sup> ) and a network of PODs developed by the local National Forest and their partnering state and local fire management agencies (an additional 1767 km<sup>2</sup> ). PODs that intersected a 5 km buffer around the watershed were included to analyze fires that spread into the watershed from nearby areas. Elevation ranges from 1559 to 4135 m above sea level across the study area. The climate is continental with warm dry summers and cold winters. Most erosion in this region results from intense rainstorms during the summer and early fall [16,49]. The study area is primarily forest (71.7%, most of which is dominated by conifers) and the remainder is a mix of shrubland (9.0%), sparsely vegetated alpine (8.9%), and grassland (8.7%) [50]. Land ownership is split between the USDA Forest Service (55.3%), private (18.8%), National Park Service (18.1%), local government (6.4%), and state (1.4%).

#### *2.3. Potential Fire Operational Delineations*

PODs were developed by fire and resource management specialists from the local National Forest and external fire management partners from other federal, state, and local agencies. The PODs range in size from 502 to 23,672 ha with a mean of 4316 ha and a median of 3516 ha. The PODs tend to be smallest near human settlements due to both the increased presence of control features and greater need for fire containment around communities. PODs larger than 10,000 ha are clustered in the higher elevation, western portion of the of the study area where much of the land is publicly owned and the transportation network is sparse. PODs also tend to be large along the major river canyon that runs west to east across the study area (Figure 2) due to limited presence of high probability control features other than the river and highway in the canyon bottom.

The rugged topography, rocky soils, and dense forests of the study area are major constraints on firefighter and equipment accessibility and operability. Accordingly, managers preferentially chose roads as the control features to bound PODs; of the 1386 km of POD edge, 985 km are roads (71.0%), 167 km are trails (12.0%), 150 km are ridges (10.8%), 46 km are streams (3.3%), and the remaining 40 km are fuel transitions, lakes/reservoirs, or lacking defined control features (2.9%). Many of the trails and ridges selected as control features are in barren or sparsely vegetated areas of the alpine, so roads make up an even larger proportion of POD edges in the fuel types where wildfire transmission is a concern. Numerous observational studies have documented that roads benefit fire control by serving as hard fire breaks that either stop fires passively or in combination with suppression firing or holding

activities [30,51–53]. The frequent use of roads in this POD network suggests containment probability should be high along most boundaries under low to moderate fire weather and many boundaries have potential for containment under more extreme conditions with well-coordinated suppression tactics.

**Figure 2.** Map of the study area featuring the focal watershed and PODs that intersect a five km buffer around the watershed. Landcover is from LANDFIRE [50]. Barren is sparsely vegetated alpine. The inset maps the location of the study area in the USA.

#### *2.4. Fire Occurrence*

We used the Monte Carlo fire simulation program RANDIG, which is a command-line version of the FlamMap minimum travel time module [39], to model a plausible set of 5000 large fire growth events across the study area. The inputs to RANDIG include raster surfaces of fuels, topography, and ignition density, and a set of fire scenarios describing the fuel moisture, wind speed, wind direction, spot probability, and burn duration for the simulations and their probabilities of occurrence. The intent of our model parameterization is to approximate the distribution of potential area burned during the initial growth period of large fires owing to variation in wind direction and wind speed. We focused on the early growth period of fires to align with the desire to contain most fires before they leave the POD of origin. Modeling fire growth over longer periods would increase fire size and thus the avoided area burned and water quality impacts but would also introduce greater uncertainty about final fire extent as more potential containment features are encountered and weather conditions are likely to moderate.

Raster fuels and topography data representing landscape conditions circa 2014 were acquired from LANDFIRE [50] including canopy cover, canopy bulk density, canopy base height, canopy height, surface fire behavior fuel model [54], elevation, slope, and aspect. Fuels were adjusted in lodgepole pine (*Pinus contorta* var. *latifolia*) forests by lowering the canopy base height by 20% and changing the fire behavior fuel model to high load conifer litter (TL5 from [54]) to better match recent observations of extreme fire behavior in these forests [55]. The other spatial input is a raster surface of ignition density, which influences the relative probability of fire ignition across the modeling domain. Spatial point

locations of historical fires from Short [56] were generalized into a raster surface of ignition density using a kernel density function with a search distance of 10 km in ArcGIS 10.3 [57].

Fuel moisture, wind speed, and wind direction for the fire scenarios (Table 1) were informed by data from a Remote Automated Weather Station [58] located in the northern half of the study area at 2500 m above sea level. Fuel moisture and wind speed percentiles were calculated with FireFamilyPlus 4.1 [59] and wind speed was converted from a 10-min to 1-min average based on Crosby and Chandler [60]. Most large fires in this region occur in early summer during drought years or in the fall when fuel moisture is extremely low. These conditions were approximated using the historical 3rd percentile fire season (1 April–31 October) fuel moistures, which are 2, 3, 6, 30, and 60 percent for the 1-h, 10-h, 100-h, herbaceous, and woody fuels, respectively. Fuel moisture was held constant across all wind scenarios because it exhibits little meaningful variation below the 10th percentile. Wind scenarios were designed to approximate the joint probability distribution of wind speeds and directions that are problematic for fire growth. The 50th, 90th, and 97th percentiles of 1-min average wind speeds are 19.3, 33.8, and 43.5 kph (at 6 m). We generalized these into three levels of wind speed (16.1, 32.2, and 48.3 kph) and their associated spotting probabilities (0.02, 0.05, and 0.10) that we assigned relative probabilities of occurrence of 0.90, 0.07, and 0.03. Previous large fires in this landscape are associated with strong westerly winds and our analysis of the historical record found that 74.1% of all winds greater than or equal to 16.1 kph were from the northwest, west, or southwest, which have relative probabilities of occurrence equal to 0.29, 0.48, and 0.23. We combined the three levels of wind speed and spotting probabilities with the three variations of wind direction into a total of nine fire scenarios (Table 1). Burn duration was set to four hours for all scenarios, which was determined by incrementally adjusting burn duration in 30 min time steps until the largest simulated fire was within ±5% of 20,000 ha, which we judge as a reasonable upper bound for fire size during a single burn period in this landscape based on other fires in the region [46].


**Table 1.** Fire scenarios used to simulate fires in RANDIG. Burn duration was set to 240 min and fuel moisture was held constant at the 3rd percentile of the historical record.

#### *2.5. Fire Behavior and Severity*

Crown fire activity [61] was modeled as a proxy for burn severity with FlamMap 5.0 [40] by mapping surface fire, passive crown fire, and active crown fire to low, moderate, and high severity, respectively. Crown fire activity is commonly used to estimate burn severity for watershed modeling [24,33,62] because it captures the trend of increasing fire intensity along the gradient of surface to active crown fire behavior. Fuel moisture was set to the same 3rd percentile fuel moisture described in the fire occurrence section. The same topography and modified fuels rasters were also used as the landscape inputs to FlamMap. To simplify the analysis, we modeled burn severity for the middle wind speed scenario (32.2 kph at 6 m) and used the wind blowing uphill option to represent a consistent worst-case scenario for all aspects.

#### *2.6. Post-Fire Watershed Response*

Post-fire erosion and sediment transport to the water diversion was predicted with a system of coupled hillslope erosion, hillslope sediment transport, and channel sediment transport models (Figure 1) that has been calibrated to make reasonable predictions of post-fire sediment yields within the study region [24]. The NHDPlus raster and watershed network products [63] were used to represent the topological connections between upland sediment sources and the water diversion point via sub-catchment drainage paths to the flowline network and the series of intervening flowlines between each catchment and the diversion. First, gross hillslope erosion was modeled for each fire with a raster Geographic Information System implementation [64] of RUSLE [41]. Sediment transport to streams was predicted using an empirical model of post-fire hillslope sediment delivery ratio from the western USA [42] to estimate the proportion of sediment generated in each pixel that makes it to the flowline network. Third, the total sediment from each catchment was routed down the flowline network to the diversion point using a simple model of channel sediment delivery ratio [43] adapted for the channel types in the study area.

#### 2.6.1. Hillslope Erosion

RUSLE predicts gross erosion (Mg ha−<sup>1</sup> year−<sup>1</sup> ) as the product of factors for rainfall erosivity (R), soil erodibility (K), length and slope (LS), cover (C), and support practices (P) [41]. Rainfall erosivity is calculated as the product of storm maximum rainfall intensity and kinetic energy per unit area [41]. First year post-fire erosion was modeled at three levels of May to October rainfall erosivity—403, 887, 5168 MJ mm ha−<sup>1</sup> h <sup>−</sup>1—representing the 2, 10, and 100-year recurrence interval rainfall erosivity (hereafter "rainfall erosivity") for the regional climate [65,66]. The May through October period was selected because most post-fire erosion in this region occurs in response to high intensity summer rainfall [16]. LS was calculated from a 30 m resolution digital elevation model [63] following the methods of Winchell et al. [67] with a maximum limit on flow accumulation of 0.9-ha imposed to approximate the original hillslope length guidance in Renard et al. [41]. Baseline K came from the Soil Survey Geographic Database where available and the State Soil Geographic Database to fill missing data [68]. Post-fire erosion was simulated by modifying the K and C factors based on wildfire extent and burn severity [24,69]. No support practices were considered to model the unmitigated erosion hazard. Baseline erosion is not a major concern for water quality, so we focused our assessment on the post-fire increase in erosion. First-year post-fire increase in erosion (A) was calculated with Equation (1) for each level of rainfall erosivity.

$$\mathbf{A} = \mathbf{R} \times \mathbf{L} \mathbf{S} \times [(\mathbf{K}\_{\mathbf{b}} \times \mathbf{C}\_{\mathbf{b}}) - (\mathbf{K} \times \mathbf{C})].\tag{1}$$

The subscript b indicates the burned condition for K and C factors. We limited hillslope erosion predictions to 100 Mg ha−<sup>1</sup> year−<sup>1</sup> based on the maximum observed values reported in the study region [49].

#### 2.6.2. Hillslope Sediment Transport

An empirical model of post-wildfire hillslope sediment delivery ratio (hSDR) from the western USA [42] was used to estimate the proportion of sediment generated in each pixel that makes it to the stream network. The NHDPlus flowlines were first extended to include all pixels with a contributing area greater than 10.8 ha [70] to better approximate the extent of the post-fire channel network. Post-fire hSDR was then estimated with the annual length ratio model from Wagenbrenner and Robichaud [42]. We applied this model to predict hSDR as a function of the flow path length from each pixel to the nearest stream channel as the "catchment length" and the flow path length across the pixel as the "plot length" (Equation (2)). Flow path length to the nearest channel was calculated from a 30 m digital elevation model [63] in ArcGIS 10.3 [57]. We doubled the predicted hSDR to account for under-sampling of suspended sediment in the model training data and to roughly calibrate our net sediment yield predictions to the small catchment yields from the Hayman Fire in Colorado [42]. This increased the maximum hSDR from 0.27 to 0.54 for areas near streams and it increased the minimum hSDR from 0.05 to 0.10 for locations furthest from streams. We later compare our modeled gross and net hillslope

sediment yields to relevant field observations in the discussion to demonstrate that this assumption is reasonable. Channel pixels were assigned hSDR of 1.

$$\log(\text{hSDR}) = -0.56 - 0.0094 \times \text{(flow path length to channel/flow path length across pixel)}, \quad \text{(2)}$$

The first-year mass of sediment (Mg) delivered from a catchment to the stream network (TS) was calculated as the sumproduct of the post-fire hillslope erosion (A), the pixel area, and hSDR for all burned pixels (N) in the catchment (Equation (3)).

$$\text{TS} = \text{SUM(A}\_{\text{i}} \times 0.09 \text{ ha/pixel} \times \text{hSDR}\_{\text{i}}) \text{i} = 1 \text{ to } \text{i} = \text{N}\_{\text{i}} \tag{3}$$

#### 2.6.3. Channel Sediment Transport

Sediment was routed through the NHDPlus flowline network to the diversion by adapting the channel sediment delivery ratio (cSDR) model of Frickel et al. [43] to the channel types in the study watershed [24]. In montane streams of this region, sediment retention is generally highest in low order channels because of high roughness and limited transport capacity and very low in the high order channels with high transport capacity [45]. Observations of post-fire sediment transport in a similar watershed in Wyoming suggest transport of fine sediments in suspension should be very efficient in high order channels even during base flow conditions [71]. These trends are approximated in our model by assigning cSDRs of 0.75, 0.80, 0.85 and 0.95 per 10 km of stream length to 1st, 2nd, 3rd, and 4th or higher-order streams, respectively. Sediment retention in lakes and reservoirs was accounted for by assigning as a cSDR of 0.05 to the terminal flowline in each waterbody. The annual mass of fire-related sediment (Mg) delivered to the water diversion (TD) was calculated as the sum of sediment delivered to streams for all upstream catchments multiplied by the product of cSDRs for the intervening flowlines (Equation (4)).

$$\text{TD} = \text{SUM(TS}\_{\text{\%}} \times [\text{PRODUCT}(\text{cSDR}\_{\text{k}}) | \text{k} = 1 \text{ to k} = \text{P}]) \text{j} = 1 \text{ to j} = \text{O},\tag{4}$$

The subscript j is the index for the O upstream catchments and the subscript k is the index for the P intervening flowlines between catchment j and the pipeline diversion.

#### *2.7. Water Supply Impacts*

The first metric of water supply impact is the total wildfire related sediment delivered to the diversion (Mg). The second metric is the per-fire average post-storm suspended sediment concentration (SSC). Wilson et al. [66] found that a threshold rainfall intensity of 7 mm h−<sup>1</sup> best predicts when post-fire hillslope erosion will occur in this region. This intensity is exceeded on average four times per year in the study watershed. We make the simplifying assumption that the first-year post-fire sediment load from the coupled erosion and sediment transport model is divided equally among four storms. We estimate that 35% of the hillslope erosion predicted by RUSLE is part of the fine-grained inorganic and organic components that contribute to suspended sediment based on observations of soil particle sizes generated from post-fire hillslope erosion and transported in suspension after summer thunderstorms in the region [71,72]. Post-fire water quality is usually degraded for short periods (hours to days) following rainstorms in this region [48,73], so we calculate post-storm suspended sediment concentrations using the average storm load of fine sediment and the daily flow volume past the diversion point, which averages 1.48 × 10<sup>9</sup> L per day for the May to October period (gage-adjusted estimates from [63]). Suspended sediment concentration is rarely monitored directly, so limits for treatment are more commonly expressed in turbidity. For this analysis, we use the high end of 100 Nephelometric Turbidity Units (NTU) reported in the literature [15,74] to be conservative in our judgement of exceeding limits for treatment. A conversion equation (Equation (5)) developed from

post-fire monitoring of the Fourmile Canyon Fire was used to predict turbidity (NTU) from SSC (mg L−<sup>1</sup> ) [15].

$$\text{NTU} = (\text{SSC} - \text{2.84}) / 1.166,\tag{5}$$

#### *2.8. Containment E*ff*ectiveness Evaluation and Prioritization*

To quantify the effectiveness of containment, we focused on the difference between the total water impact measures with and without containment including watershed area burned, sediment delivered to the diversion, and number of turbidity threshold exceedances. The difference between impact measures for the uncontained and contained scenarios is the avoided transmitted risk [36]. Total sediment load is a continuous value whereas turbidity exceedance is a binary outcome. Water quality degradation was only considered transmitted when the outcome changed from below 100 NTU for within POD burning to above 100 NTU for the entire fire footprint. To prioritize improvements along the potential control lines that bound PODs, we calculated risk transmission across the POD edges based on their proportional engagement with the fires that originate in their respective PODs; that is, the outcomes associated with fire spreading to the surrounding landscape were divided among the lines based on their intersected length. It is anticipated that the primary mitigation action would be fuels reduction along the control lines, so transmission risk was normalized by length to compare the relative benefit of hardening control lines.

#### **3. Results**

#### *3.1. Fire Occurrence*

Historical fire ignitions from the FOD [56] were concentrated in the lower and middle portions of the focal watershed and along the southern boundary of the study area (Figure 3a) reflecting both variation in fire season length and human use of the landscape. The 5000 wildfires simulated with RANDIG ranged in size from 0.09 to 20,868 ha with a mean of 1961 ha and a median of 1469 ha. We selected the 3040 fires that burned at least part of the focal watershed for further analysis. Their size distribution did not vary substantially from that of the full simulation set. The excluded fires either did not grow large enough to intercept the focal watershed, or the predominant wind direction caused them to spread away from it. The middle and lower portions of the watershed are predicted to burn most frequently due to both the greater ignition density and the presence of fuel types that promote faster spread (Figure 3b). The high elevations in the western half of the study area are predicted to burn infrequently due to low ignition density and sparse fuels. The southeast corner of the study area near the water diversion has low burn probability because the fuels have not yet recovered from a recent wildfire.

#### *3.2. Fire Behavior and Severity*

Crown fire activity is predicted to vary across the watershed due to differences in fuels and topography (Figure 4a). A notable portion of the alpine and some recently burned areas are mapped as non-burnable cover types (13.7%). Surface, passive crown, and active crown fire are predicted on 25.9%, 39.3%, and 21.1% of the watershed area respectively, which we use as proxies for low, moderate, and high burn severity. This translates to predictions of low severity effects in grass and shrub fuel types and moderate or high severity effects in most forests. High severity effects are most common in forests with high horizontal and vertical continuity on steep slopes. Our prediction that approximately 60% of the watershed should burn at moderate or high severity is in line with the observed severity of recent large wildfires in Colorado [75].

**Figure 3.** (**a**) Fire Occurrence Database (FOD) records of historical ignitions and interpolated surface of relative ignition density used in the RANDIG simulations. (**b**) Burn probability from the simulated fires that intercept the study watershed.

**Figure 4.** (**a**) Predicted burn severity using crown fire activity categories of surface, passive crown, and active crown fire as proxies for low, moderate, and high severity fire. (**b**) Predicted post-fire erosion with 2-year rainfall erosivity. (**c**) Combined Sediment Delivery Ratio (SDR) accounting for both hillslope and channel transport. (**d**) Predicted sediment delivery to the water supply diversion with 2-year rainfall erosivity.

#### *3.3. Watershed Response*

Like burn severity, the magnitudes of post-fire erosion and sediment transport vary widely across the watershed owing to variation in topography, soils, and proximity to the diversion. Figure 4 illustrates this for the 2-year rainfall erosivity. The greatest sediment hazard is associated with steep terrain near the major channels that is predicted to burn at moderate or high severity. Post-fire erosion and sediment transport potential is generally low in the flatter terrain in the northeast quadrant of the watershed, the high mountains above major waterbodies, and the recently burned areas. The spatial distribution of sediment hazard is similar for 10-year and 100-year rainfall erosivity, but the absolute magnitude increases considerably. Table 2 summarizes the distribution of predicted erosion, sediment delivery to streams, and sediment delivery to the diversion for the 3040 simulated wildfires that burned in the watershed. The predicted mean post-fire gross erosion for the simulated wildfires is 12.3, 20.4, and 46.4 Mg ha−<sup>1</sup> for the 2, 10, and 100-year rainfall erosivity, respectively. Much of this sediment should be retained in the watershed, especially where waterbodies interrupt sediment transport (Figure 4c), so delivery to the diversion averages only 4.2, 7.0, and 15.9 Mg ha−<sup>1</sup> for the 2, 10, and 100-year rainfall erosivity, respectively.


**Table 2.** Summary statistics of first-year post fire erosion, sediment delivery to streams, and sediment delivery to the water supply diversion (div.) in Mg ha−<sup>1</sup> by rainfall erosivity for the simulated wildfires that burned into the watershed. These are total sediment yields including the coarse and fine fractions.

#### *3.4. Avoided Watershed Area Burned*

For improved containment at POD boundaries to avoid water supply impacts, the target fires must leave the POD of origin under unmanaged conditions. Of the 3040 simulated wildfires that burned at least part of the focal watershed, 2351 of them (77.3%) burned at least some area outside the origin POD. Fires occasionally burned more than ten PODs, but of the fires that burned more than one POD, most burned between two and five PODs (77.9%). This suggests that most fire transmission during the initial burn period is between a POD and its adjacent neighbors, but some rare events may burn across multiple POD boundaries.

Containing all fires within their POD of origin would reduce the average watershed area burned from 1361 to 562 ha per fire, a 58.7% reduction (Table 3). The distributions of watershed area burned for the contained and uncontained scenarios are shown in Figure 5a. Containing large fires has the greatest potential to avoid watershed area burned; the 1396 fires that burned more than 1000 ha account for 93.8% of the avoided area burned. Containment in the POD of origin would eliminate fires that burn more than 10,000 ha in the watershed, which numbered 26 (0.9%) in the uncontained scenario. The percentage of fires burning greater than 5000 ha would be reduced from 4.0 to 0.2. Watershed area burned by fires that originate from PODs that are wholly or mostly outside the watershed should be reduced to negligible levels under the containment scenario, but these PODs account for only a small fraction of area burned when fires are allowed to grow freely (Figure 6a). Most fires start in the central and eastern portion of the watershed (Figure 3) and the predominant west winds means that PODs in the lower 2/3rds of the watershed are the source of fires that burn the greatest area (Figure 6a). All else equal, larger PODs are larger sources of fire because they have more ignitions. Containment reduced watershed area burned from fires that ignited in 61 of the 70 PODs, but some of the largest PODs still have substantial watershed area burned with containment (Figure 6a) because fires have room to grow large before encountering a potential control feature.


**Table 3.** Summary of water supply impacts across all fires by containment scenario and rainfall erosivity. A turbidity threshold of 100 NTU was used to compute the number of exceedances.

**Figure 5.** Summary of containment effects on distribution of fire-level indicators of water supply risk by rainfall erosivity including: (**a**) watershed area burned, (**b**) first-year post-fire sediment to the diversion, and (**c**) first-year post-fire average post-storm turbidity (vertical black line marks the 100 NTU threshold for treatment).

#### *3.5. Avoided Sediment*

Containment reduced the total sediment load to the pipeline diversion by 50.4–55.3% depending on rainfall erosivity from an average of 6.1–23.2 thousand Mg per fire to an average of 3.1–10.4 thousand Mg per fire (Table 3). The distributions of sediment delivered to the diversion for the contained and uncontained scenarios are shown in Figure 5b. Sediment loads vary across several orders of magnitude due to differences in fire size, erosion and sediment transport potential, and post-fire rainfall. The effect of containment on sediment load is roughly equivalent to reducing rainfall erosivity one level (Figure 5b). The spatial distribution of sediment source risk is similar to that of watershed area burned (Figure 6b). PODs that are partially or wholly outside the watershed are a minimal risk to water supplies after containment, but fire activity in the larger PODs situated in the middle of the watershed is still expected to produce large sediment loads.

**Figure 6.** Spatial summary of containment effects on distribution of POD-level indicators of water supply risk for the 2-year rainfall erosivity including: (**a**) watershed area burned, (**b**) first-year post-fire sediment to the diversion, and (**c**) frequency of turbidity exceedances for fires that originate within each POD.

#### *3.6. Avoided Water Quality Degradation*

Containment effects on water quality degradation were less substantial than for watershed area burned and total sediment to the diversion (Table 3; Figure 5c); turbidity exceedances were reduced by 33.5, 21.3, and 13.0 percent for the 2, 10, and 100-year rainfall erosivity, respectively. With containment, 36.5, 49.4, and 63.2 percent of fires are predicted to exceed the 100 NTU threshold for the 2, 10, and 100-year rainfall erosivity, respectively. Most fires that caused turbidity to exceed limits for treatment originated in the large PODs in the middle of the watershed (Figure 6c). The three PODs with the most turbidity exceedances are all larger than 10,000 ha. Containment only reduced the number of turbidity exceedances from these PODs from 640 to 568 (an 11.3% reduction) for the 2-year rainfall erosivity, and containment offered almost no mitigation benefit (1.0% fewer exceedances) for these PODs under the most extreme rainfall scenario. In contrast, containment reduced turbidity exceedances by more than 50% in 33 of the 70 PODs under median rainfall conditions. These PODs range in size from 502 to 14,153 ha with a mean of 3548 ha. Many of these PODs are mostly or wholly outside the watershed, but some are smaller PODs inside the watershed.

#### *3.7. Prioritizing POD Network Improvements*

The limited effect of containment on turbidity exceedances highlights the need to break up the three large PODs with high source risk in the middle of the watershed (Figure 6c). These three PODs are also the top priorities for further compartmentalization based on watershed area burned and total sediment load from self-burning. With containment, an additional eight PODs were the source of 20 or more turbidity exceedances under median rainfall conditions. Cumulatively, these top 11 PODs account for 91.4% of the fires that degraded water quality in the contained scenario, so efforts to further reduce fire sizes in these PODs should have high benefit.

Prioritizing improvements along the potential control lines that bound PODs can be informed with measures of risk transmission (Figure 7). Total sediment to the diversion was transmitted at the highest rates along POD edges in the middle portion of the watershed (Figure 7a) where there is high potential for fires to spread into erosion prone terrain near the diversion (Figures 3b and 4). In contrast, transmitted water quality degradation was more concentrated along the edges associated with the smaller PODs in the north central portion of the watershed (Figure 7b). Transmission risk was also high for several control lines in the eastern half of the watershed that are nearly perpendicular to the dominant wind direction. Mitigation priorities differed depending on which metric of transmission risk was used (Figures 7 and 8). The two metrics both identify a similar order of priorities (Spearman's ρ = 0.89) but they have moderate disagreement about the magnitudes of potential risk mitigation (Pearson's R = 0.71), especially for the highest-ranking edges (Figure 8). Most notably, few of the POD edges associated with the three large PODs that are the source of most turbidity exceedances (Figure 6c) are high priorities for mitigation because containment at these locations infrequently changes the water quality outcome despite the potential to avoid large quantities of sediment.

**Figure 7.** Total transmitted risk for all fires from (**a**) sediment to diversion and (**b**) turbidity exceedances normalized to edge length in kilometers for 2-year rainfall erosivity.

**Figure 8.** Edge transmission risk comparison for 2-year rainfall erosivity. Edges ranked in the top 20 using either metric are colored red.

#### **4. Discussion**

This proof of concept analysis demonstrates the potential for improved early containment of large fires to lower watershed area burned by 58.7% and to reduce risk to source water between 13.0% and 55.3% depending on impact metric considered. Proportional reductions in total sediment load to the diversion ranged between 50.4% and 55.3%, but the potential to avoid exceeding turbidity limits for treatment was notably lower—varying between 33.5% and 13.0% reduction for the 2- and 100-year rainfall erosivity, respectively (Table 3). The contrasting response of our water impact metrics to increasing rainfall erosivity (Table 3) reveals that avoiding large quantities of sediment may not translate to avoiding degraded water quality if the residual sediment load is still large. The sources of water supply risk and potential mitigation benefits of fire containment varied widely across the POD network (Figure 6) suggesting the potential to further improve mitigation effectiveness with targeted divisions to reduce the size of PODs with high risk from self-burning and fuels reduction to improve containment probability along high transmission boundaries (Figure 7).

Our analysis built on previous studies of wildfire-water supply risk and wildfire risk transmission to estimate the avoided water supply impacts from improved fire containment within pre-identified PODs. Omi [18] approached this issue from the related perspective of fuel break construction and maintenance in California using estimates of avoided area burned and a relative damage index to value fuel break benefits. Monte Carlo wildfire simulation and watershed effects analyses capture similar information on exposure and impacts with the added benefit of associating fire outcomes with their ignition locations and final extents [32,33]. A recent effort to zone the study landscape into PODs provided the operationally relevant fire containers used to estimate avoided water supply impacts using risk transmission methods [35,36] as suggested by Davis [76] to estimate the area saved from burning after encountering a control feature. The avoided area burned and sediment load measures we modeled are similar to the impact metrics used to value the benefit of containment in previous studies, but our evaluation of water quality degradation provided a unique opportunity to evaluate whether the size and spatial arrangement of the PODs are appropriate to mitigate a scale-dependent risk. Our results suggest POD-based containment could meaningfully reduce risk of exceeding turbidity limits for treatment (Table 3), but the large percentage of unmitigated risk implies that the containment network could be more effective with smaller PODs.

Our estimates of avoided impacts are premised on the simplifying assumption that all fires are contained within their POD of origin, which is likely realistic for many of our modeled fires but optimistic for the most extreme fires in the region [35,46]. We chose not to address the probability of containment in this study because existing models focus on characteristics of the control features, surrounding fuels and topography, and fire behavior [30,77,78] but do not explicitly consider the effects of suppression [79]. Managers in this landscape primarily identified roads as control features because they aid firefighter access [29] and suppression firing [51]. It is also anticipated that proactively identifying control features and response strategies will lead to timely and well-coordinated tactics that increase the probability of containment. For example, extensive pre-season planning has been credited with improving the strategic use of suppression firing and aerial retardant drops to contain fire in PODs during extreme weather [80]. We did not account for suppression firing in this study, which can sometimes substantially increase area burned [81] and thus would dampen the contrast between our containment scenarios. However, managers ideally use backing fire to minimize adverse effects [80]. Improved modeling of suppression actions and effects would help to refine our estimates of risk mitigation.

The post-fire erosion and sediment transport modeling used here has several limitations that are important to acknowledge. First, the linked fire and erosion model system (Figure 1) is subject to multiple data, model, and model linkage uncertainties that have potential for prediction error as discussed extensively in previous publications [24,25]. Recent work has shown that water quality at the basin scale is sometimes minimally impacted despite modeled increases in hillslope erosion [82], emphasizing the need to test and refine erosion and sediment transport models with empirical observations at multiple scales [83]. Most of our predicted first-year post-fire hillslope erosion yields for the 2-year and 10-year rainfall erosivity scenarios (Table 2) are close to the study-wide means of 9.5–22.2 Mg ha−<sup>1</sup> and the range of individual hillslope observations of 0.1–38.2 Mg ha−<sup>1</sup> from previous fires in the region exposed to moderate rainfall [11,17,47,84]. Many of these studies had hillslope sediment fences fill and overtop, so the reported yields are usually interpreted as a lower bound estimate of the true erosion rate. For the 100-year rainfall erosivity, only the top decile of modeled fires exceed the 72 Mg ha−<sup>1</sup> of rill and interrill erosion reported in the first year after the Buffalo Creek Fire in response to similarly extreme rainfall (converted from volume estimates of [45] using bulk density of 1.6 Mg m−<sup>3</sup> ). Despite doubling the efficiency of hillslope transport in this study, only the net sediment delivery to streams for the upper decile of fires with 10-year rainfall erosivity and the upper half of fires with 100-year rainfall erosivity (Table 2) approach the small catchment sediment yields of 22.0–38.6 Mg ha−<sup>1</sup> observed in the first two years after the Hayman Fire [85,86]. This seems reasonable given the larger size of most catchments in this study. After our rough calibration, our combined hillslope and channel SDR values (Figure 4c) are close to SDR values estimated with similar travel time methods [87,88]. None of the simulated fires at any rainfall level (Table 2) are predicted to deliver sediment to the diversion at a rate close to the whole watershed sediment yield of 52.5 Mg ha−<sup>1</sup> for the first year of the Buffalo Creek Fire [45], likely because we did not account for channel erosion.

Our water degradation analysis also layers on additional assumptions that the annual suspended sediment load is evenly divided among the annual average of four sediment-generating storms and the storm sediment load is evenly mixed in the average daily flow volume of the river during the thunderstorm season. Despite these approximations, the resulting turbidities—which averaged 309, 516, and 1181 NTU for the 2, 10, and 100-year rainfall erosivity, respectively—align well with common observations in the region of post-fire turbidities between 100 and 1000 NTU and occasional observations >1000 NTU [15,48,89]. The assumption that storm load is an equal division of annual load does not account for the substantial intra-annual variability in storm characteristics [15,83,90], seasonal trends in runoff and erosion [91], nor the interannual variability in the frequency of storms

with sufficient intensity to cause erosion [34,66]. Similarly, unaccounted for variability in daily flow volume should influence the vulnerability of the water source. Given these simplifications, we have more confidence in our contrasts of containment benefits across scenarios than we do in our absolute estimates of degradation risk. Our analysis also focused exclusively on the acute periods of severe water quality degradation after rainstorms in the first year after fire, so it is unclear if containing fires to smaller sizes will avoid elevated carbon, nitrogen, phosphorus, manganese, and suspended solids concentrations that may persist for years after fires in Colorado [15,89], increasing treatment complexity and cost and raising concerns about the formation of disinfection byproducts [74,92]. Similar water quality responses and treatment challenges have been observed after wildfires in Canada, Australia, and Europe [12,13].

Despite uncertainties in the precise magnitude of risk reduction, improved containment appears promising compared to other mitigation strategies. We found that limiting fires to their POD of origin should reduce the total sediment load from wildfire between 50.4 and 55.3% (Table 3). Previous assessments of landscape-scale fuel treatments in the western USA predict long-term sediment reduction of 19% [25] and up to 34% reduction in sediment costs [24]. Salis et al. [93] project that treating 15% of a landscape in Sardinia, Italy would only reduce average sediment yield 4–12%, but their treatment scenarios were not prioritized to avoid erosion. Based on the narrowest contrast in these figures (34% for fuel treatment and 50.4% for containment), POD-based containment should compare favorably to landscape scale fuels reduction as long as the containment failure rate is less than 32%. Furthermore, compartmentalizing fire in small units of the landscape has the potential to avoid disrupting multi-source water systems by limiting fire impacts to a single source. The benefit of containing individual wildfires should vary widely (Figure 5), as fire encounters with control features and associated impacts beyond the POD of origin depend strongly on where the fire ignites.

We also demonstrated how risk transmission metrics could inform improvements to the POD network, which should be relevant to fire, land, and water managers engaged in spatial fire planning. The small number of PODs with high risk from self-burning are high priorities for further compartmentalization, which could require improving firefighter access and/or reducing fuels. Fine scale analyses of risk factors and containment opportunities would benefit these efforts. If further divisions are not feasible or practical (e.g., because of wilderness or wildlife habitat concerns), these PODs could be candidates for fuels reduction with prescribed or managed fire. It is also valuable for water managers to identify areas that are not conducive to proactive risk mitigation, so they can plan how to best respond to the anticipated effects of future fires. As previously discussed, we did not estimate the probability of containing wildfire at POD boundaries and how containment probability would change with fuels reduction, but managers are interested in identifying potential control lines in need of improvement to support safe and effective fire response. Measures of transmission risk across the POD edges (Figure 7) highlight where these efforts should be targeted to maximize their benefit. However, priorities differed depending on the water supply effects measure used (Figure 8); most notably, there is greater potential to avoid degradation by improving containment probability around the smaller PODs. Further analyses are needed to evaluate if fuel conditions around these POD edges necessitate treatment for firefighting effectiveness and safety.

The style of Monte Carlo exposure and effects analyses we present should also be useful for evaluating fire protection strategies for other high value resources and assets that depend on the scale of disturbance. For example, most ecological concerns relate to the area and spatial pattern of high severity effects on vegetation and the resulting consequences for wildlife and reforestation by dispersal-limited species (e.g., [94–96]). If intolerable levels of fire exposure or effects can be defined for ecological values, similar methods could be used to assess the protection value of POD-based containment. Wildfire impacts to homes and other values in the wildland-urban interface (WUI) are almost always negative, but consequences often become disastrous when the area and assets affected by fire overwhelm firefighting resources [97]. Managers intuitively design smaller PODs in the WUI [2], but it has not been tested whether these PODs are appropriately sized to avert WUI disasters—i.e., whether asset exposure for most fires is below the fire protection capacity. Similarly, wildfire impacts to transportation networks may cross thresholds of concern for evacuation when traffic exceeds the capacity of the available routes. Explicitly defining performance objectives for these and other fire protection concerns could help to tailor POD size and spatial arrangement in future fire planning efforts.

#### **5. Conclusions**

Improved wildfire containment has potential to meaningfully reduce wildfire risk to water supplies, but these effects are scale dependent. In our test cases, approximately 75% of fires intersected potential control features and, if these fires were contained within their POD of origin, watershed area burned would be reduced by 58.7%, total sediment load to the diversion would be reduced between 50.4 and 55.3%, and water quality degradation beyond limits for treatment would be reduced between 13.0 and 33.5%. Risk mitigation was higher for total sediment load than water quality degradation because containment did not always change water quality outcomes. Moreover, priorities to improve the network design by modifying the size of the PODs or improving containment probability along their edges differ depending on the effects measure used. This highlights the importance of properly defining water supply impacts for wildfire risk assessment and mitigation effectiveness studies. Similar analyses could be applied to other scale-dependent resources at risk of wildfire to inform containment network design.

**Author Contributions:** Conceptualization, B.M.G., Y.W. and M.P.T.; methodology, B.M.G.; software, B.M.G.; formal analysis, B.M.G.; investigation, B.M.G.; data curation, B.M.G.; writing—original draft preparation, B.M.G.; writing—review and editing, B.M.G., Y.W. and M.P.T.; visualization, B.M.G.; supervision, Y.W.; project administration, Y.W.; funding acquisition, Y.W. and M.P.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by joint venture agreement 19-JV-11221636-170 between the USDA Forest Service Rocky Mountain Research Station and Colorado State University; cost share agreement number 17-CS-11021000-032 between the USDA Forest Service, Arapaho and Roosevelt National Forests and Pawnee National Grassland and the Colorado Forest Restoration Institute at Colorado State University; and agreement number 19-DG-11031600-062 between the USDA Forest Service, Southwestern Region and the Colorado Forest Restoration Institute at Colorado State University.

**Acknowledgments:** The authors thank Codie Wilson for sharing the rainfall data used in the analysis.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

## **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Assessing Potential Safety Zone Suitability Using a New Online Mapping Tool**

**Michael J. Campbell 1, \* , Philip E. Dennison 1 , Matthew P. Thompson <sup>2</sup> and Bret W. Butler 3**

	- matthew.p.thompson@usda.gov

**Abstract:** Safety zones (SZs) are critical tools that can be used by wildland firefighters to avoid injury or fatality when engaging a fire. Effective SZs provide safe separation distance (SSD) from surrounding flames, ensuring that a fire's heat cannot cause burn injury to firefighters within the SZ. Evaluating SSD on the ground can be challenging, and underestimating SSD can be fatal. We introduce a new online tool for mapping SSD based on vegetation height, terrain slope, wind speed, and burning condition: the Safe Separation Distance Evaluator (SSDE). It allows users to draw a potential SZ polygon and estimate SSD and the extent to which that SZ polygon may be suitable, given the local landscape, weather, and fire conditions. We begin by describing the algorithm that underlies SSDE. Given the importance of vegetation height for assessing SSD, we then describe an analysis that compares LANDFIRE Existing Vegetation Height and a recent Global Ecosystem Dynamics Investigation (GEDI) and Landsat 8 Operational Land Imager (OLI) satellite image-driven forest height dataset to vegetation heights derived from airborne lidar data in three areas of the Western US. This analysis revealed that both LANDFIRE and GEDI/Landsat tended to underestimate vegetation heights, which translates into an underestimation of SSD. To rectify this underestimation, we performed a bias-correction procedure that adjusted vegetation heights to more closely resemble those of the lidar data. SSDE is a tool that can provide valuable safety information to wildland fire personnel who are charged with the critical responsibility of protecting the public and landscapes from increasingly intense and frequent fires in a changing climate. However, as it is based on data that possess inherent uncertainty, it is essential that all SZ polygons evaluated using SSDE are validated on the ground prior to use.

**Keywords:** firefighter safety; safe separation distance; safety zones; LCES; Google Earth Engine; lidar; LANDFIRE; Landsat; GEDI

#### **1. Introduction**

Wildland firefighters are tasked with a wide variety of fire management duties, many of which place them in close proximity to flames. One of the primary tasks is the removal of fuels and construction of containment lines in order to limit the potential damage to lives, property, and other critical resources [1–3]. Particularly when engaged in a direct attack, whereby firefighters may be working within a few meters or less of the flaming front, the potential risk for safety incidents is elevated [4]. Sudden or unexpected changes in fire behavior can have devastating effects to vulnerable fire personnel on the ground [5]. Events such as the Yarnell Hill fire in 2013, which claimed the lives of 19 firefighters and the South Canyon fire, which resulted in 14 firefighter fatalities, demonstrate the tragedy that can occur in the wildland fire profession [6–8]. Beyond these well-known, high-fatality events, there is an additional and significant background level of mortality that occurs among on-duty wildland firefighters [9]. The causes of death are varied, and include heart

**Citation:** Campbell, M.J.; Dennison, P.E.; Thompson, M.P.; Butler, B.W. Assessing Potential Safety Zone Suitability Using a New Online Mapping Tool. *Fire* **2022**, *5*, 5. https://doi.org/10.3390/fire5010005

Academic Editor: James R. Meldrum

Received: 12 November 2021 Accepted: 5 January 2022 Published: 7 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

attacks, vehicular and aircraft accidents, falling trees, and smoke inhalation, to name a few. Between 1990 and 2016, there were 480 wildland firefighter fatalities, nearly one fifth (19%) of which were due to burnovers or entrapments. Burnovers are events in which flames overcome fire personnel, either on the ground or while in a vehicle, and entrapments are when firefighters become trapped by surrounding flames, unable to evacuate to safety [10]. As wildland fires increase in frequency, extent, and intensity, wildland firefighters may be put at heightened risk while working in the increasingly complex fire environment [11–15]. Passage of the Infrastructure Investment and Job Act in the US developing a requirement for federal agencies to "develop and adhere to recommendations for mitigation strategies for wildland firefighters to minimize exposure due to line-of-duty environmental hazards" further underscores the importance and continued policy relevance of wildland firefighter safety (see H.R. 3684 §40803(d)(5)(A)).

To avoid injury or fatality, firefighters employ a range of safety measures in anticipation of and throughout fire management operations. Safety is often considered the first priority of wildland firefighters, as evident in the 10 Standard Firefighting Orders and the 18 Watch Out Situations, colloquially referred to as the "10s and 18s" and known by all fire crews in the US [16,17]. One of the most important safety protocols is the interconnected and interdependent system of lookouts, communications, escape routes, and safety zones (LCES) [18]. First formalized by Paul Gleason in 1991, and later modified to include anchor points by Thorburn and Alexander in 2001 [19], proper implementation and continued reevaluation of LCES is essential to firefighter survival. Lookouts are members of a fire crew well-trained in the interpretation of fire behavior and weather conditions who are placed at a vantage point in the fire environment that ensures continued visual observation of the fire crew, the fire itself, and the surrounding landscape. Communications ensure that critical information is conveyed in a clear, timely, and orderly manner between various resources deployed on a fire, including ground crews, lookouts, vehicular and aerial assets, and incident command. Escape routes are pre-defined evacuation pathways that enable crews to retreat to a safety zone or other safe area. And finally, safety zones (SZs), which are the primary focus of this paper, are areas on the landscape that firefighters can retreat to in dangerous situations in order to avoid bodily harm. SZs are large, open areas with little or no flammable material contained within. They can take a variety of forms, including naturally sparsely or unvegetated areas and areas that have had fuel removed either mechanically, or through burning.

The effectiveness of SZs is defined by the extent to which they can provide safe separation distance (SSD) from surrounding or nearby flames [4,5,20–22]. SSD is defined as the distance one must maintain from fire to avoid burn injury, without deployment of a fire shelter [5,22]. The current operational guideline in the US for defining SSD emerged from the work of Butler and Cohen in 1998, who used radiative heat transfer modeling to determine that firefighters should maintain a distance of at least 4 times the height of the proximal flames [22]. Although this guideline has since been widely adopted [23], the research that underlies it is based solely on one heat transfer mechanism: radiation. Heat transfer by convection is also a major—sometimes dominant—force, particularly in the presence of steep slopes and high winds [24–26]. In the presence of such convective heat, particularly if a fire crew is upslope and/or downwind of flames, SSD will increase [5,21]. Thus, the four times flame height rule is likely insufficient in these conditions. Recent work by Butler et al. has sought to update this guideline with the inclusion of a "slope-wind factor", which adds a multiplicative term to the SSD equation to account for the effects of convective heat transfer [5,21,26–28]. In addition, given that SZs should be designated prior to, rather than during, the presence of flames, the four times flame height rule requires firefighters to predict how tall the flames might eventually be, which is a challenging endeavor. Accordingly, the newly proposed guidelines assume that, in a crown fire, flame height is approximately equal to twice the vegetation height [28]. As a result, the new SSD equation is defined as:

$$\text{SSD} = 8 \times \text{VH} \times \Delta,\tag{1}$$

where VH is vegetation height and ∆ is the slope-wind factor. Butler recently defined these slope-wind factors seen in Table 1, based not only on slope and wind speeds, but also on the burning conditions, as dictated by fuel conditions (e.g., moisture) and weather (e.g., relative humidity) [28].



Although guidelines for use on the ground are valuable, they still require the firefighters themselves to make the calculation of SSD on the ground while engaged in other fire management activities. This requires the ability to accurately estimate vegetation height and terrain slope and anticipate wind speed and fire intensity. Moreover, even if these difficult interpretations and predictions can be made, an even more challenging endeavor is to identify an area on the ground cleared of vegetation that provides the calculated SSD in all directions. In order to reduce the subjectivity of this process and the potential for interpretation error, Campbell et al. sought to develop a robust, geospatially driven approach for identifying and assessing the suitability of SZs in advance of a fire using high-resolution airborne lidar, which provides detailed terrain and vegetation height information [20]. Their method automatically identified open areas, determined the height of surrounding trees, and calculated a score for the relative suitability of potential safety zones within the open areas based on tree height and SSD. The potential for broad applicability of their approach, however, is limited by three main factors: (1) the lack of widely available airborne lidar data, (2) the temporal relevancy of and lack of scheduled updates to available lidar data, and (3) the fact that only existing clearings could be assessed. With respect to the first limitation, at present, high-quality, publicly accessible airborne lidar data is particularly sparse in the Western US, where fires are most frequent, large, and intense. Regarding the second limitation, vegetation structure is dynamic in fire-prone environments, where SZ assessment is most important. Thus, having outdated lidar data potentially limits accurate SSD calculation. Lastly, although existing open areas may be viable SZs, firefighters often use areas that have already burned ("black" areas) or create SZs through mechanical fuel removal.

To resolve the limitations of this previous work and to improve wildland firefighter safety, we are introducing a new, interactive, web-based, open-access mapping tool for estimating SSD and evaluating potential SZ effectiveness through geospatial analysis. Instead of relying on lidar, this tool uses LANDFIRE Existing Vegetation Height data, which is both nationally available in the contiguous US and is updated every few years. Additionally, instead of only assessing SSD-driven suitability on clearings that already exist, this tool allows users to draw their own SZ polygon to evaluate the potential suitability of a SZ in any environment. Since LANDFIRE vegetation heights may not be as accurate as airborne lidar, given that it is a modeled product driven by satellite imagery, it is important to quantify the effects of differing sources of vegetation height data on SSD evaluation. Accordingly, the primary objectives of this study are to: (1) introduce and describe the algorithm that underlies a new tool for calculating SSD and analyzing SZ suitability; (2) compare SSD and

potential SZ suitability using different sources of vegetation height data; (3) demonstrate an example case study of the tool in a simulated wildland firefighting situation.

#### **2. Materials and Methods**

#### *2.1. Algorithm Description*

The Safe Separation Distance Evaluator (SSDE) algorithm is built and applied in Google Earth Engine (GEE), a cloud-based platform for processing and analyzing GIS and remotely sensed data, using the JavaScript application programing interface [29]. GEE was selected for few reasons: (1) it enables the production of user-facing applications that can be widely accessed by anyone with an internet connection; (2) it hosts an immense catalog of geospatial data, including datasets necessary for the analysis of SZ suitability; (3) its cloud computing capabilities provide for rapid execution of complex geospatial functions, allowing users to quickly assess SZ suitability. Accordingly, all data processing described in this section is conducted using the GEE.

The SSDE evaluates SSD in two primary ways (Figure 1). The first is per-pixel SSD, which is a representation of how far one must be from that pixel (e.g., in meters) in order to avoid burn injury (Figure 1a). This is calculated at the individual pixel level across an entire area of interest based on the vegetation height and terrain slope within each pixel, and user-defined wind speed and burn condition classes. It provides a landscape-scale view of SSD and can be used to aid in the delineation of potential SZ polygons. However, it is perhaps more important to evaluate SSD at the level of the SZ polygon, as this can help fire personnel determine the suitability of a potential SZ. Accordingly, the second way that SSDE evaluates SSD is through the analysis of proportional SSD (pSSD) within potential SZ polygons (Figure 1b). pSSD quantifies the extent to which a potential SZ polygon provides SSD from surrounding vegetation/flames, considering the average per-pixel SSD contained within a series of segments (or clusters of contiguous pixels) around the SZ polygon. Measured in percent, a pSSD of 100% or greater for a given pixel would mean that, factoring in vegetation height surrounding the polygon, slope, wind speed, and burn condition, the pixel's location should provide sufficient SSD, should fire personnel opt to use this location as a SZ. Conversely, a pixel with a pSSD of less than 100% would indicate that firefighters located within that pixel may risk injury from burning vegetation outside the boundary of the polygon. A detailed description of the computation of both SSD and pSSD follows.

**Figure 1.** Conceptual depiction of the two different ways that SSDE calculates SSD, including a per-pixel SSD, representing the distance one must maintain from a pixel in order to avoid burn injury (**a**), and within-SZ polygon pSSD, representing the relative extent to which a potential SZ polygon provides SSD from surrounding flames, based on segment-level mean SSD (**b**).

Per-pixel SSD is calculated as follows. The algorithm relies on two primary input datasets, both of which are available on a nationwide basis in the US. The first is LANDFIRE

Existing Vegetation Height, version 2.0, which is a 30 m spatial resolution raster dataset where each pixel represents the average height of the dominant vegetation within that pixel (Figure 2b) [30,31]. The second is a Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM), also at 30 m spatial resolution [32]. Using the SRTM DEM, terrain slope in percent is calculated (Figure 2c). This slope raster is then converted into a slopewind factor dataset, based on user-defined wind speed and burning condition categories (Figure 2d; Table 1). Using Equation (1), SSD is calculated from the slope-wind factor and vegetation height data, producing a spatially exhaustive representation of the distance one must maintain from each pixel in the image dataset in order to avoid injury (Figure 2e). As discussed further in Sections 2.3 and 3.3, there is a user-selected bias correction option that is built into SSDE to account for underestimation of tree heights in the LANDFIRE data.

**Figure 2.** Illustration of the process for calculating per-pixel SSD in an example focus area in central New Mexico (**a**), including the vegetation height data from LANDFIRE (**b**), the SRTM-derived slope data (**c**), the slope-wind factor under moderate burning conditions and moderate wind speeds (**d**), and the resulting SSD raster (**e**). Legend labels contain parenthetical indication relevant sub-figure. Note that this figure does not illustrate the SZ selection process in a real or simulated wildland fire; it merely conveys the calculation of per-pixel SSD in a landscape in which fires could occur.

The per-pixel, landscape-wide layers representing SSD, vegetation height, and slope, as well as the high-resolution Google imagery base map, can be used to identify potential SZ on the landscape. In a realistic wildland fire scenario, this identification process would likely be undertaken by a GIS specialist working with an incident management team, who has specific knowledge of the current fire extent, projected changes in fire behavior, and crew assignments. However, given the open-access nature of SSDE, this process can likewise be undertaken by anyone with an interest in wildland firefighter safety. To calculate pSSD at the SZ polygon level, the user can define a potential SZ using the polygon drawing tools in SSDE, guided by the conditions both within the polygon and surrounding the polygon. The best SZs are those that contain no flammable material within, naturally or otherwise, so ideally this SZ polygon would be drawn in an area with low fuel loading, such as short or sparse grasses or litter. Alternatively, a SZ polygon could also be drawn

in an area that has recently burned, or an area that would be targeted for fuel removal to create a SZ.

pSSD within a SZ is dependent upon the slope and vegetation height of the surrounding landscape. As discussed in the Introduction, heat transfer from flames generally increases with increasing vegetation height and terrain slope. Accordingly, a SZ in the midst of steep terrain and tall vegetation will require a larger SSD than a SZ in the midst of flat terrain and short vegetation. If we assume that the SZ itself contains little or no flammable material, then the primary concern for SZ evaluation is the area surrounding the SZ. Accordingly, to calculate pSSD within the SZ polygon, slope and vegetation height need to be evaluated within a "buffer" surrounding a SZ [20,33]. In an environment with short vegetation, a relatively small buffer surrounding a SZ needs to be evaluated, because flame heights will tend to be low and more distant vegetation will not affect the suitability of a SZ given its limited capacity for transferring heat over long distances (Figure 3a). Conversely, in the presence of tall trees and steep slopes, larger buffer areas are needed to account for the effects of more distant fuel and flames (Figure 3b). For example, a stand of 20 m tall forest in a low-slope environment with moderate wind speed and burning conditions would require an SSD of 320 m—that is, to maintain safety, firefighters would need to find a SZ that provides at least 320 m of separation from this stand. However, if the SZ evaluation procedure only included an analysis of vegetation within a small (e.g., 20 m) buffer around the SZ, then the potentially dangerous effects of the heat emitting from a crown fire in this stand could be missed, putting the firefighters at risk.

**Figure 3.** Relationship between the safe separation distance (SSD) defined by vegetation surrounding the safety zone and the size of the buffer needed around the safety zone to ensure that all relevant surrounding vegetation is considered in the safety zone analysis. With shorter vegetation, a smaller buffer is needed (**a**), whereas with taller vegetation, a larger buffer is needed (**b**). Note that the relationship between vegetation height and SSD is not drawn to scale.

Accordingly, to ensure that an appropriately sized buffer area is evaluated around a potential SZ polygon, the buffer size must be defined by the highest per-pixel SSD in the surrounding area. To do this, the largest theoretically possible SSD in the US (SSDUSmax) is first calculated as follows:

$$\text{SSD}\_{\text{USmax}} = 8 \times \text{VH}\_{\text{USmax}} \times \Delta\_{\text{max}} \tag{2}$$

where VHUSmax is the height of the tallest tree in the US (116 m) and ∆max is the maximum slope-wind factor from Table 1 (10). As a result, SSDUSmax is equal to 9280 m. While this is purely a theoretical SSD, it provides a useful initial buffer distance. The potential SZ polygon is buffered by this SSDUSmax, and then the highest SSD pixel value within the initial buffer (SSDSZmax) is then identified. This value is then used to create the final buffer around a designated SZ polygon, which will serve as the basis of SSD evaluation.

The buffer area is then used to clip out a subset of the SSD raster (Figure 4a). Given that each pixel possesses its own SSD (based on its unique combination of vegetation height and slope-wind factor), it is important to consider the variation in SSD that is found throughout the area surrounding the SZ polygon. Vegetation may be taller or slopes may be steeper on one edge of a SZ polygon as compared to the other. Accounting for this spatial variability at the pixel level would greatly increase computation time, so the SSDE algorithm applies an image segmentation function on the clipped SSD raster to generate a series of segments (i.e., clusters of similar pixels) representing areas with relatively homogeneous SSD values (Figure 4a). From each segment, Euclidean distance is calculated on a continuous basis as a 30 m spatial resolution raster where each pixel represents the distance from that segment. Each segment also has a mean SSD value, aggregated from the SSD pixel values within the segment. Accordingly, a comparison between the Euclidean distance raster and the segment's mean SSD value should provide insight into whether pixels in the SZ polygon are within or beyond the SSD. To determine this on a relative basis, the Euclidean distance is then divided by the segment's mean SSD, resulting in a dataset containing pixels representing pSSD, where pixel values under 1 (or 100%) represent areas that would be unsafe if the fuel within that segment were burning under crown fire conditions (Figure 4b). Values over 1 (or 100%) represent areas that would be safe, according to the SSD guideline. This process is repeated for every segment surrounding the SZ polygon, producing a series of pSSD rasters. To reduce these rasters down to a single representation of pSSD within the SZ polygon, the minimum pixel value (representing the "worst-case" pSSD for all segments surrounding the polygon) is extracted among all of the rasters (Figure 4c). This pSSD raster enables the identification of the safest area within the potential SZ polygon, defined as the location with the highest pSSD (Figure 4d). Lastly, a threshold is applied to the pSSD to distinguish between safe areas, or areas that are equal to or greater than 100% of pSSD, and unsafe areas, or areas that are less than 100% of pSSD (Figure 4d). The size of the resulting safe areas is important for wildland firefighters, as it provides a sense for how many resources (e.g., firefighters, trucks, dozers, engines) can utilize the SZ at once. It has been estimated that approximately 5 m<sup>2</sup> (50 ft<sup>2</sup> ) is recommended for each firefighter, and 28 m<sup>2</sup> (300 ft<sup>2</sup> ) for each piece of heavy equipment [33,34]. To enable users to interact with the resulting data outside of GEE, there are options for downloading all of the resulting data layers, including the SSD raster, pSSD raster, safe/unsafe raster, the SZ polygon vector feature, and the safest point vector feature.

#### *2.2. Vegetation Height Analysis*

Given that LANDFIRE Existing Vegetation Height is a modeled product based on multispectral satellite imagery, it can be assumed that there is inherent inaccuracy in its vegetation height estimates. Conversely, airborne lidar data, which result from laser pulses emitted from an airborne platform reflecting off the ground and aboveground surfaces, provide direct measurements of vegetation height. Thus, airborne lidar is considered a reliable source of data for deriving accurate and precise canopy height models in vegetated environments [35]. In fact, the models that derive the LANDFIRE Existing Vegetation Height product are trained in part using airborne lidar data, in addition to ground-level plot measurements [36]. Inaccuracy in vegetation height results in inaccuracy in SSD, and inaccuracy in SSD could have significant effects on the safety of wildland firefighters. Accordingly, to quantify these effects, we compared LANDFIRE-driven SSD to lidar-driven SSD.

**Figure 4.** Illustration of the process for calculating within-safety zone (SZ) safe separation distance (SSD) in the same example focus area in central New Mexico shown in Figure 1, including the manual input of a SZ polygon, the buffering of that SZ polygon, and the generation of SSD segments (**a**), the calculation of pSSD from a single example segment (**b**), the combined minimum of pSSD from all segments (**c**), and derivative safety products, including the safest point and a binary classification of safe versus unsafe within the SZ polygon (**d**). Legend labels contain parenthetical indication relevant sub-figure(s).

Although LANDFIRE is one of the most widely used, nationwide vegetation height products in the US, particularly among wildland fire scientists, it is not the only one. A recent study by Potapov et al. introduced a new global forest height product driven by a combination of Global Ecosystem Dynamics Investigation (GEDI) data and Landsat 8 OLI imagery [37,38]. GEDI is a spaceborne lidar instrument that produces large footprint (~25 m diameter) full waveform surface elevation data. Distinct from airborne lidar, however, which produces a dense point cloud of measurements that can be used to interpolate spatially exhaustive, high-resolution height models, GEDI's sampling design is such that large gaps exist between successive and adjacent footprints, limiting the capacity to produce spatially contiguous height models. However, Potapov et al. used GEDI data to train a model to predict forest height using Landsat imagery on a global scale. Similarly to LANDFIRE, this GEDI/Landsat hybrid product is a modeled product and, as such, possesses inherent uncertainty. In the interest of basing our algorithm on the best and most reliable vegetation height data, we also compared the GEDI/Landsat-derived SSD to airborne lidar-derived SSD.

To compare LANDFIRE, GEDI/Landsat, and airborne lidar, we selected three areas in fire-prone regions of the Western US. These areas were selected to capture a range of terrain and vegetation height conditions. Importantly, each of the three areas was required to have recently collected, freely available airborne lidar data and no recent major disturbances. The 3 areas selected, each 20 × 20 km in size, are shown in Figure 5. The first area is in western New Mexico (Figure 5a), ranging in elevation from 1659 to 3032 m (mean = 2253 m), ranging in slope from 0 to 153% (mean = 26%), and is dominated by ponderosa pine forests and piñon-juniper woodlands. According to the airborne lidar data, vegetation heights range from 0 to 41 m (mean = 13 m). The second area is in western Wyoming (Figure 5b), ranging in elevation from 1741 to 3038 m (mean = 2284 m), ranging in slope from 0 to 179% (mean = 35%), and is dominated by Douglas fir and lodgepole pine forests and sagebrush shrublands. According to the airborne lidar data, vegetation heights range from 0 to 40 m (mean = 15 m). The third area is in northern California (Figure 5c), ranging in elevation from 137 to 1399 m (mean = 729 m), ranging in slope from 0 to 195% (mean = 30%), and is dominated by Douglas fir forests, black oak woodlands, and ruderal grasslands. According to the airborne lidar data, vegetation heights range from 0 to 72 m (mean = 22 m).

Airborne lidar datasets for each of these 3 areas were acquired from the USGS 3D Elevation Program (3DEP) [39]. The New Mexico data came from the USGS 3DEP dataset titled "NM\_SouthCentral\_2018\_D19", collected in 2018 with an average point density of 6.0 pts/m<sup>2</sup> . The Wyoming data came from the USGS 3DEP dataset titled "WY\_Southwest\_2020\_D20", collected in 2020 with an average point density of 7.6 pts/m<sup>2</sup> . The California data from the USGS 3DEP dataset titled "USGS\_LPC\_CA\_NoCAL\_Wildfires \_B4\_2018", collected in 2018 with an average point density of 10.1 pts/m<sup>2</sup> . It should be noted that there are some temporal discrepancies between the airborne lidar data (2018–2020), LANDFIRE data (2016), and GEDI/Landsat data (2019). To avoid issues associated with major vegetation height changes that may have occurred between these time frames, these areas were selected specifically to avoid containing any fires that had occurred between 2015 and 2021. To match the training and validation data used in the development of LANDFIRE and Potapov et al.'s height products, airborne lidar canopy height models were derived at a 30 m spatial resolution using the 90th percentile of aboveground point returns within each pixel. All of the airborne lidar data processing was conducted using LAStools [40].

To assess SSD uncertainty, the algorithm described in Section 2.1 was applied to each of the three vegetation height datasets (LANDFIRE, GEDI/Landsat, and airborne lidar). However, to enable the flexibility needed for this comparative analysis, the algorithm was implemented in Python with heavy reliance on the Esri ArcPy library, rather than GEE. Given that the GEDI/Landsat product does not include vegetation heights lower than 3 m (as it is considered a "forest" height product), LANDFIRE vegetation heights were used in areas lower than 3 m in height for the GEDI/Landsat analysis. For all three datasets, the same terrain slope was used, derived from a USGS 3DEP 30 m digital elevation model. In each of the three test areas, SZ polygons were automatically generated, using a random shape generator created for this study. Briefly, it starts with a randomly located point, and then creates a shape based on a randomly defined number of equally spaced vertices that are placed at random distances from the center point. The distances are selected from a normal distribution with a mean between 100 and 1000 m and a standard deviation of 0.25 times the mean. The resulting shape is then smoothed to produce a final SZ polygon, as can be seen in Figure 5a–c. In all, 25 SZ polygons were generated for each test area (75 total). Since each polygon was evaluated independently, randomly placed polygons could overlap. For each SZ polygon, every unique combination of burning condition and wind speed was evaluated. As a result, the algorithm was run 2025 times (3 vegetation height datasets × 75 SZs × 3 burning conditions × 3 wind speeds).

**Figure 5.** Three test areas used to assess SSD uncertainty, including a site in western New Mexico (**a**), western Wyoming (**b**), and northern California (**c**), as well as a locator map illustrating their broader geographic context (**d**). The LANDFIRE SSD shown represents low burning conditions and light wind. The SZs represent randomly generated polygons used for testing.

Four statistical comparisons were made between the results from the three different vegetation height datasets. The first was a per-pixel SSD comparison. For this, a set of 25 random points was generated within each combination of vegetation height dataset, burning condition, and wind speed (675 points total), and pixel values were extracted from SSD rasters. The second was the maximum pSSD found within each of the randomly generated SZ polygon for each combination of vegetation height dataset, burning condition, and wind speed. The third was the total safe area (the areas meeting or exceeding SSD) within the SZ polygon for each combination of vegetation height dataset, burning condition, and wind speed. All three of these comparisons were made using ordinary least squares regression between LANDFIRE and airborne lidar, GEDI/Landsat and airborne lidar, and LANDFIRE and GEDI/Landsat. The fourth and final statistical comparison was aimed at determining the geographic divergence between the locations identified as being the

safest point within each SZ. For this, the horizontal distance between safe points was calculated between LANDFIRE and airborne lidar, GEDI/Landsat and airborne lidar, and LANDFIRE and GEDI/Landsat, from which histograms and associated descriptive statistics were derived.

#### *2.3. LANDFIRE Bias Correction*

As mentioned in Section 2.1 and discussed in greater detail in the Section 3.2, the per-pixel SSD comparison revealed that, in comparison to the airborne lidar-based SSD values, the LANDFIRE and GEDI/Landsat SSD values were underestimated. To ensure that the algorithm is not overestimating the relative safety of drawn SZs, a bias correction procedure was developed. Although both LANDFIRE and GEDI/Landsat featured a similar trend in bias, the bias correction was only tested on the LANDFIRE data, given the fact that LANDFIRE was ultimately selected as the vegetation height dataset used in SSDE, the justification for which is further discussed in Section 4. To correct for underestimation in SSD, the linear regression model resulting from the per-pixel SSD comparison was used to adjust LANDFIRE SSD pixel values through a linear transformation. To test the extent to which this bias correction procedure increased the agreement between LANDFIRE- and lidar-derived SSD values and SZ suitability metrics, the same four statistical comparisons previously described were conducted: (1) per-pixel SSD linear regression; (2) within-SZ maximum pSSD linear regression; (3) within-SZ total safe area linear regression; (4) safest point distance histogram comparison. Lastly, to enable users to select between using the raw LANDFIRE data or the bias-corrected LANDFIRE data as a basis of SSD calculation, a selection option was added to SSDE.

#### *2.4. Use Case Demonstation*

To provide an example use case for SSDE, we simulated a situation in which wildland firefighters are tasked with selecting a SZ during an active wildfire. To do this, we created a fictional fire perimeter in a ponderosa pine forest in an area of northern New Mexico that has a variety of potential SZs in the form of open grassy meadows. Two of these potential SZs are delineated and analyzed using SSDE. To test the viability of these potential SZ under different weather and fire conditions, they are each evaluated for two sets of extremes: (1) light wind and low burn condition; (2) high wind and extreme burn condition.

#### **3. Results**

#### *3.1. GEE Application*

The SSDE was developed in GEE and a free, open-access, web-based application can be viewed at https://firesafetygis.users.earthengine.app/view/ssde (accessed on 3 January 2022). Screen captures of the interface can be seen in Figures 6–8. Figure 6 represents the view that a user would have when zoomed to an area of interest, including an instructional panel on the left that explains how to use the tool (Figure 6a), the main map interface with a Google imagery base map (Figure 6b), and the polygon drawing tool that can be used to digitize a potential SZ polygon (Figure 6c).

The user then has the option to select the burning condition, wind speed, and whether or not to correct for vegetation height bias (Figure 7a). Following these selections, the user generates a per-pixel SSD map (Figure 7b). The user also has the option to download the resulting raster image (Figure 7a).

The user can then manually draw in a potential SZ polygon and evaluate a variety of suitability metrics by selecting 'Calculate SZ SSD' (Figure 8a). As a result, three new spatial layers will be added to the map, each within the extent of the SZ polygon, including the pSSD (Figure 8b), as well as a point feature representing the safest location and a safe (pixels greater than or equal to SSD) and unsafe (pixels less than SSD) raster. These layers can be toggled on and off in the 'Layers' menu (Figure 8c). The layers can also be individually downloaded (Figure 8a). Lastly, a series of tabular metrics are also computed and displayed to aid in the evaluation of SZ suitability (Figure 8d). The user can evaluate

and compare multiple different wind speeds, burning conditions, and whether or not bias is corrected for.

**Figure 6.** A screen capture of the Wildland Fire Safe Separation Distance Evaluator Google Earth Engine web application, including an instructional panel on the left that explains how to use the tool (**a**), the main map interface with a Google imagery base map (**b**), and the polygon drawing tool that can be used to digitize a potential SZ polygon (**c**).

**Figure 7.** A screen capture of the Wildland Fire Safe Separation Distance Evaluator Google Earth Engine web application, including the user-selected options for generating and downloading a SSD map (**a**), and an example of the resulting map (**b**).

**Figure 8.** A screen capture of the Wildland Fire Safe Separation Distance Evaluator Google Earth Engine web application, including the widgets used for calculating SZ suitability and downloading the resulting layers (**a**), a layer representing pSSD (**b**), the 'Layers' menu for toggling on and off layer visibility (**c**), and a widget with several quantitative suitability metrics for two different wind speed-burning condition combinations (**d**).

#### *3.2. Vegetation Height Analysis*

The results of the per-pixel SSD uncertainty analysis can be seen in Figure 9. When compared to airborne lidar-derived SSD, which we assume to be the most accurate given the fact that airborne lidar vegetation heights are directly-measured rather than modeled, it becomes clear that both LANDFIRE- and GEDI/Landsat-derived SSD values, on average, are underestimated (Figure 9a,b). Since all of the other variables (slope, wind speed, burning condition) were held constant, this means LANDFIRE and GEDI/Landsat vegetation height products underestimated vegetation heights, as taller vegetation results in higher SSD. Per-pixel SSD values are still fairly strongly correlated (R<sup>2</sup> LANDFIRE-lidar = 0.63; R 2 GEDI/Landsat-lidar = 0.68), suggesting that the modeled products trend towards similar results as the airborne lidar-derived per-pixel SSD, but the regression slopes of less than 1 indicate that both LANDFIRE and GEDI/Landsat will underestimate per-pixel SSD, particularly on the high end. This further points towards the fact that LANDFIRE and GEDI/Landsat fail to capture the structure of very tall vegetation well, though this failure appears more pronounced in the LANDFIRE data, given the lower regression line slope. Of particular note is the presence of a number of sample points where the LANDFIRE and GEDI/Landsat per-pixel SSD values approach zero but the airborne lidar per-pixel SSD values are greater than zero. A visual assessment of the three vegetation height datasets in comparison to aerial imagery revealed that, at least in some cases, this can be attributed to the presence of standing dead vegetation such as in a post-fire environment. The airborne lidar captures this standing dead vegetation since the laser pulses interact with any aboveground objects. However, the LANDFIRE and GEDI/Landsat, both of which rely heavily on spectral data tend to suggest there is no vegetation in these areas, given the fact that the typical vegetation reflectance signal (e.g., high near infrared reflectance, low visible wavelength reflectance) is no longer present after a disturbance event. The comparison between LANDFIRE and GEDI/Landsat directly further suggests that LANDFIRE tends to

underestimate vegetation heights to a greater extent than GEDI/Landsat, though they are strongly correlated (R<sup>2</sup> LANDFIRE-GEDI/Landsat = 0.74) (Figure 9c).

**Figure 9.** Results of the per-pixel SSD comparison between LANDFIRE and airborne lidar (**a**), GEDI/Landsat and airborne lidar (**b**), and LANDFIRE and GEDI/Landsat (**c**). The red line represents the ordinary least-squares regression model between the *x* and *y* variables, whereas the black line represents the 1:1 line, where *x* is equal to *y*.

The results of the SZ polygon-level maximum pSSD uncertainty analysis can be seen in Figure 10. The results follow a similar trend from the per-pixel SSD comparison, which makes sense given that within-SZ pSSD is driven by the pixel-level SSD values surrounding the SZ polygon. Both LANDFIRE and GEDI/Landsat overestimate the maximum pSSD within SZ polygons as compared to airborne lidar (Figure 10a,b). This can once again be attributed to vegetation height underestimation from these modeled products as compared to airborne lidar data. In addition, there is a much greater correlation at the SZ polygon level than at the pixel level. For example, at the pixel level, LANDFIRE and lidar SSD have a coefficient of determination (R<sup>2</sup> ) of 0.63, whereas at the SZ polygon level, LANDFIRE and lidar maximum pSSD have an R<sup>2</sup> of 0.92. This is likely due in part to the segmentation procedure that aggregates adjacent pixels with similar SSD values and computes a mean, which drives the pSSD calculation, thus reducing the pixel-level noise and producing a stand-level estimate of vegetation height. This high correlation also suggests that SZ polygon size is likely the major driver of within-SZ pSSD, outweighing the effects of differing vegetation height estimates. When comparing between the 2 modeled products, LANDFIRE and GEDI/Landsat produce very similar results with a high correlation (R<sup>2</sup> LANDFIRE-GEDI/Landsat = 0.88) and a regression line slope of close to 1 (0.94) (Figure 10c). Again, even though LANDFIRE tended to produce slightly lower per-pixel SSD values due to lower vegetation height estimation (Figure 9c), when pixels are aggregated at the segment level and when those segment-level pSSD estimates are further aggregated at the SZ polygon level, subtle differences in vegetation height bear little effect on pSSD calculation.

The results of the SZ polygon-level total safe area uncertainty analysis can be seen in Figure 11. These results mirror those seen in Figure 9, which makes sense given that safe area is a direct product of the within-SZ pSSD calculation. Both LANDFIRE and GEDI/Landsat tend to overestimate total safe area within SZ polygons as compared to airborne lidar (Figure 11a,b). While such an overestimation is potentially problematic, the most acutely problematic disagreement between the modeled products and the airborne lidar is the large number of SZ polygons where lidar indicates there is 0 ha safe area and the modeled products indicate there is greater than 0 ha safe area—in some cases as much as 100 ha of safe area. However, it is important to recognize that since the safe versus unsafe classification is based on a defined pSSD threshold, a SZ polygon that reaches 99% of SSD would still be mapped as unsafe, so small differences in vegetation heights can have large effects on these results. Once again, the LANDFIRE and GEDI/Landsat safe

area estimates are both highly correlated and relatively unbiased when compared to one another (Figure 11c).

**Figure 10.** Results of the safety zone polygon-level maximum pSSD comparison between LANDFIRE and airborne lidar (**a**), GEDI/Landsat and airborne lidar (**b**), and LANDFIRE and GEDI/Landsat (**c**). The red line represents the ordinary least-squares regression model between the *x* and *y* variables, whereas the black line represents the 1:1 line, where *x* is equal to *y*.

**Figure 11.** Results of the safety zone-level total safe area comparison between LANDFIRE and airborne lidar (**a**), GEDI/Landsat and airborne lidar (**b**), and LANDFIRE and GEDI/Landsat (**c**). The red line represents the ordinary least-squares regression model between the *x* and *y* variables, whereas the black line represents the 1:1 line, where *x* is equal to *y*.

The results of the safest point geographic displacement analysis can be seen in Figure 12. LANDFIRE and GEDI/Landsat result in safest points being on average between about 100 and 150 m from the lidar-derived safest point location (Figure 12a,b). Given the 30 m spatial resolution of the analysis, that means, on average, the safest points are within 3–5 pixels of one another. This suggests that the geography of the safest point is perhaps not as sensitive to the vegetation height as it is to the geometry of the SZ polygon. However, there are exceptions, given the right-skewed tail of the distributions. In a select few instances, safest points were found to be 800 m apart or more. From a firefighter safety perspective, this is a significant deviation, as the safest point is likely the target gathering point within the SZ for the crew and equipment. The LANDFIRE and GEDI/Landsat safest points are, on average, closer to one another than they each are to the airborne lidar-derived safest points (Figure 12c).

**Figure 12.** Results of the safety zone polygon-level total safest point location comparison between LANDFIRE and airborne lidar (**a**), GEDI/Landsat and airborne lidar (**b**), and LANDFIRE and GEDI/Landsat (**c**).

#### *3.3. LANDFIRE Bias Correction*

The results of the LANDFIRE bias correction procedure can be seen in Figure 13. The bias correction resulted in SSD and SZ suitability metrics that align more closely with those derived from airborne lidar. Specifically, in all three regression analyses, including the per-pixel SSD (Figure 13a), the within-SZ polygon maximum pSSD (Figure 13b), and the within-SZ polygon total safe area (Figure 13c), the regression slope was much closer to one than the un-corrected data. The bias correction bore little effect on the safest point location comparison, which makes sense given the fact that vegetation heights were all increased linearly, meaning the placement of the safest point would not have changed significantly (Figure 13d).

**Figure 13.** Results of the comparison between bias-corrected LANDFIRE data and airborne lidar data, including per-pixel SSD (**a**), within-SZ polygon maximum pSSD (**b**), within-SZ polygon total safe area (**c**), and distance between modeled safest points (**d**).

#### *3.4. Use Case Demonstration*

The results of the case study demonstration can be seen in Figure 14. The two potential SZ polygons varied in shape, size, and proximity to the fire perimeter. In a direct attack situation, a fire crew would be working in close proximity to the perimeter of the fire. Thus, having a SZ that is in close proximity to the crew is the best way to ensure that they can access the SZ quickly via a pre-planned escape route. Based on proximity alone, SZ 1 is superior to SZ 2 (Figure 14a). Indeed, provided that wind speeds were light and burn conditions were low, SZ 1 would provide sufficient SSD for the crew, providing a maximum within-SZ pSSD of 115%, and a total safe area of 0.36 ha (Figure 14b). However, under high winds and extreme burn conditions, this SZ would no longer be suitable (Figure 14c). In comparison, SZ 2 is further away, but due to its larger size, provides greater SSD, particularly under light winds and low burn conditions (Figure 14d). Even though this polygon is much larger, it still does not meet the SSD guideline under high winds and extreme burn conditions, providing a maximum within-SZ pSSD of 74% (Figure 14e). Armed with this information, however, the crew could opt to expand SZ 2 to increase its suitability through mechanical fuel removal or controlled burning.

#### **4. Discussion**

We envision the SSDE being of broad interest to the wildland fire community, from fire scientists to incident management personnel to wildland firefighters. Its open-access nature allows anyone to explore, examine, and interact with the concepts of SZs and SSD. Even if not used in an operational context, there is great value in being able to quickly and easily examine the conditions that define potential SZ suitability on a broad spatial scale. Wildland firefighters designate SZs on a daily basis as a part of their fire management duties. Built into this designation process is an inherent degree of subjectivity that can result in differences in the interpretation of SZ suitability between and among crews. By using an objective tool for SZ suitability analysis that can be broadly applied in the US, fire crews

across the country can increase the consistency and reliability of the SZ evaluation process. However, given that this is a web-based platform requiring an internet connection, that also does not translate well to a mobile environment, we do not envision this as a real-time decision-making tool that firefighters could use on the ground. Instead, operational use of SSDE could be at the incident command level, for daily or more frequent evaluation of potential SZs for crews working on a fire. Given the dynamic and fast-paced nature of fire management, it is essential to be able to make rapid assessments of SZ suitability, particularly as fire conditions change. For example, cross-referencing near-real time data representing the current fire perimeter with SSDE can enable the evaluation of whether or not previously burned areas can provide SSD from nearby unburned fuel. Additionally, our SZ suitability driver analysis revealed the strong influence of wind on maximum within-SZ pSSD. This highlights the need to continually re-evaluate SZ suitability not only as the fire evolves, but as local weather conditions change as well.

There is relatively little research on the actual utilization of SZ in wildland firefighting. As discussed in the Introduction, SZ are largely designated ad hoc by firefighters, with little support from GIS and remote sensing data. As a result, there are essentially no publicly available datasets that represent SZ that were used by firefighters. If they were available, we could perform more empirically driven studies of SZ effectiveness, but as it stands, we are limited to using theoretical models to predict safety zone effectiveness based on fire physics [21,22,25–27]. While we do not know how many firefighters have suffered injuries or fatalities due specifically to the selection of inadequate SZ, we do know that burnovers and entrapments occur on an annual basis, and that these injuries and fatalities would be much less likely to occur if escape routes and SZ were properly employed. Every wildland fire safety incident is different from the last, so there is no systematic way to determine the extent to which a tool such as SSDE would have mitigated the incident. However, with the increased infusion of reliable geospatial data and modeling into wildland firefighter safety, we hope to minimize future incidents. In future research, it would be beneficial to compare the results of SSDE to actual SZ used by firefighters. This could be accomplished through the GPS collection of SZ in the field or perhaps through simulation of SZ selection in a virtual environment, to minimize risk to fire personnel while still gaining valuable, realistic insight into SZ selection [41–43].

The SSD guideline upon which SSDE is based (Equation (1), Table 1) can be quite restrictive, particularly in landscapes featuring rugged terrain and tall trees. For example, in a worst-case scenario (high winds, extreme burning conditions), trees that were 25 m tall and on slopes that were steep (>40%) would produce an SSD equal to 2 km, making SZ designation impractical, if not impossible. However, even if no SZ is found to meet guidelines, SSDE still can provide an analysis of which safety zones are closest to meeting the SSD requirement. SSDE enables evaluation of newly burned potential SZs or existing clearings that nearly meet the SSD requirement and could be enlarged through prescribed burning or mechanical fuel removal. For example, if a potential SZ polygon delineated within an existing clearing was determined through the use of SSDE to only reach a maximum pSSD of 90%, SSDE could be further used to evaluate the relative suitability of one or more potential SZ expansion options. These options could be informed by the aerial imagery, terrain, vegetation height, and SSD data built into SSDE, which could aid in the determination of the areas that might be best suited for fuel removal (i.e., removing sparse, short shrubs on flat slopes rather than dense, tall trees on steep slopes to improve the SZ).

Our comparative analysis between SZ suitability metrics derived from LANDFIRE, GEDI/Landsat, and airborne lidar revealed that the former two 30 m products tended to underestimate vegetation height and thus overestimate the suitability of potential SZs. However, given that airborne lidar data are not widely available in the Western US, and even where they are available may not reflect current vegetation structural conditions, the SSDE had to be built using either LANDFIRE or GEDI/Landsat to enable broad application. If one product or the other demonstrated a significantly closer alignment with the lidarderived suitability metrics, then we would have built SSDE using the superior product, yet

no such clear winner emerged, with LANDFIRE marginally outperforming GEDI/Landsat at the SZ level and GEDI/Landsat marginally outperforming LANDFIRE at the pixel level. Thus, we were left to rely on other factors in choosing which dataset to base the SSDE algorithm on. We chose LANDFIRE for the two primary reasons. First, LANDFIRE is a well-established interagency US government-funded program that has a lengthy history and is already well-known to the wildland fire community in the US. Secondly, LANDFIRE data are scheduled to be continually updated to reflect changes in vegetation structure over time. This is particularly important given the dynamic nature of vegetation in disturbanceprone regions of the US. SSDE could potentially be expanded to other regions of the world if adequate vegetation height data becomes available. The GEDI/Landsat forest height data used in this study are globally available but does not include vegetation below 3 m in height; thus, additional height information would be needed to fill in this lower range. Alternatively, regionally specific height models could be developed where airborne lidar is more common.

Although other variables were found to be more important predictors of maximum within-SZ pSSD than vegetation height (e.g., wind speed and slope; Appendix A), we only performed an uncertainty analysis on the vegetation height source data. The slope data in our study were derived from SRTM DEMs, which are widely regarded to be an accurate source of elevation (and, in turn, slope) data [44,45]. That being said, a previous comparison between SRTM and other DEM source datasets have revealed that SRTMderived slope estimates tend to be comparatively high [46]. However, given the relatively coarse categorization of slope classes used to calculate SSD in our model, we believe the results would be minimally affected by using different sources of elevation data. With respect to the other variables (wind speed, burn condition, SZ geometry), these are all user-defined in SSDE. That is not to say that user definition does not come with inherent uncertainty, but the uncertainty analysis can, and indeed should, be conducted by the user of SSDE. For example, by simulating different wind speed and burn conditions, the user can gain a valuable sense of the variability of SZ suitability. There is also inherent uncertainty in the fire physics research that has defined the SSD equation (Equation (1)); however, testing that uncertainty is beyond the scope of this study. The SSDE platform is sufficiently flexible to enable future updates to the fundamental SSD equation that underlies the algorithm.

This study represents an important contribution to a growing body of literature surrounding the application of GIS and remote sensing to wildland firefighter safety. This study builds upon previous work by Dennison et al. and Campbell et al. [20,33] using airborne lidar to identify and evaluate SZ suitability but scales the application to the contiguous US. Escape routes are a second component of LCES that can be modeled, including variation in travel rates with slope [47,48], impacts of lidar-derived vegetation density and surface roughness on travel rates [49], and a landscape scale Escape Route Index [50]. Travel rates used in geospatial modeling can simulate the evacuation of a fire crew in comparison to predicted fire behavior in order to assess potential trigger points on the landscape that a fire crew could use [51]. Beyond SZs and escape routes, there is an emergence of other geospatially driven approaches for fire management that have implications for firefighter safety, such as the mapping of snag hazards [52], the spatially explicit estimation of suppression difficulty [2,53], the prediction of potential control locations [3,54], and the mapping of potential wildland fire operational delineations [55,56]. These tools are now commonly used in the US to enhance situational awareness and improve the quality of strategic decision-making on large and complex wildfire incidents [57]. Taken together, these tools stand to greatly improve the efficiency and safety with which wildland firefighters can engage in their life-saving work.

#### **5. Conclusions**

SZs are one of the most important safety tools available to wildland firefighters. The difference between a good SZ (one that provides sufficient SSD) and a bad SZ (one that does not), can be the difference between life and death. Although guidelines exist for the assessment of SSD, errors in estimating vegetation height, distances, and SSD based on anticipated slope, wind, and fire conditions could result in firefighter injury or fatality. In this study we have introduced a broadly applicable methodology and associated openaccess web mapping tool for the assessment of SZ suitability based on vegetation height, terrain slope, wind speed, and burning conditions. The tool enables users to delineate their own potential SZ anywhere in the contiguous US and rapidly compute a variety of suitability metrics that are driven by free and publicly available datasets. SSDE allows for those with wide-ranging expertise to explore and examine the complexities of potential SZ suitability, and how the conditions that affect SZ suitability vary over space with vegetation height and slope.

As with all geospatially driven tools, however, it is important to consider the effects that data inaccuracy and uncertainty may have on the tool's results. This importance is greatly amplified when the tool is designed to ensure the safety of others. To that end, we compared two different modeled vegetation height datasets to reference data derived from airborne lidar data with respect to their relative agreement on the array of SZ suitability metrics produced by the tool. Although correlations between the two datasets and airborne lidar were relatively high across all metrics, they tended to underestimate vegetation heights, resulting in an overestimation of the relative safety of SZs. However, this disagreement enabled us to implement a bias correction option in the tool to enable users to perform raw and lidar-adjusted SZ evaluations. Even with the bias correction, it is of the utmost importance that all results emerging from the use of the SSDE be validated on the ground by professionals prior to the formal designation of SZs. In the future, as airborne lidar becomes more widely available and frequently updated, efforts should be made towards basing SSD on lidar, rather than modeled geospatial vegetation height products. As fire science continues to advance our understanding of the complex nature between fuel, topography, and weather, the analysis of SSD should reflect the most updated understanding of how radiative and convective heat transfer affects SZ suitability. In addition, the practical application of SSD and SZs needs to be developed, to determine how firefighters utilize the information provided by guidelines and mapping tools.

**Author Contributions:** Conceptualization, M.J.C., P.E.D. and B.W.B.; Methodology, M.J.C., P.E.D. and B.W.B.; Software, M.J.C.; Formal Analysis, M.J.C.; Resources, M.J.C. and P.E.D.; Data Curation, M.J.C.; Writing—Original Draft Preparation, M.J.C. and P.E.D.; Writing—Review & Editing, M.J.C., P.E.D., M.P.T. and B.W.B.; Visualization, M.J.C.; Supervision, M.J.C., P.E.D., M.P.T. and B.W.B.; Project Administration, M.J.C. and P.E.D.; Funding Acquisition, M.J.C., P.E.D. and M.P.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by National Science Foundation grant number BCS-2117433 and USDA Forest Service grant number 19-JV-11221637-142.

**Data Availability Statement:** The Wildland Fire Safe Separation Distance Evaluator online tool can be accessed here: https://firesafetygis.users.earthengine.app/view/ssde (accessed on 3 January 2022). All data in this study come from freely- and publicly-available US Federal Government repositories.

**Acknowledgments:** We would like to thank Wesley Page and Alex Heeren for their feedback on the development of SSDE.

**Conflicts of Interest:** The authors declare no conflict of interest.

**Disclaimer:** The findings and conclusions in this report are those of the author(s) and should not be construed to represent any official USDA or U.S. Government determination or policy. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. government.

#### **Appendix A**

There are several factors that drive the suitability of potential SZ according to the SSDE. These include geometry of the SZ polygon (area, shape), height of the surrounding vegetation, slope of the surrounding terrain, the wind speed, and burn condition. However,

these factors are not equal in terms of their influence on SZ suitability. Quantifying the relative influence of suitability controls can provide valuable insight into the relative stability (or inversely, variability) of SZ suitability under different conditions. For example, if SZ polygon size explained the vast majority of variance in maximum within-SZ pSSD, then that would suggest SZs are fairly robust to changes in surrounding landscape conditions, and weather/fire conditions. To quantify the relative influence of different drivers of pSSD, we used random forests [58]. Random forests are an ensemble machine learning method that uses iterative subsampling of data to generate a series of decision trees to make predictions and assess variable importance. For each of the 75 random SZ polygons generated for the vegetation height analysis (Section 2.2; Figure 5) we had nine different pSSD values, resulting from unique combinations of wind speed and burn condition classes, totaling 675 data points. Using the randomForest library in *R* [59,60], we built a model to predict maximum within-SZ polygon pSSD, based on the independent variables described in Table A1. These models were run using the default parameters and were not tuned for performance, as the goal was merely to assess variable importance and not maximize predictive accuracy. The equations for calculating shape metrics were gleaned from Lindsay [61] and computed using Esri ArcGIS. Variable importance was assessed in terms of the relative proportion of mean squared predictor error that would result from the removal of a particular variable from the model.

**Table A1.** Predictor variables used in the analysis of the SZ suitability driver analysis.


The results of the SZ polygon suitability driver analysis can be seen in Table A2. Wind speed emerged as the most important determinant of maximum within-SZ pSSD by a large margin. This stands in stark contrast to the relatively minimal effect of burn condition, which was determined to be the least important predictor. These differences can be explained through the examination of ∆ values in Table 1, where increases in wind speed always result in a higher ∆, whereas increases in burn condition do not. Slope was the second most important predictor which, when taken together with the importance of wind speed, highlights the impact that convective heating, as captured by ∆, has on pSSD. SZ polygon area is the third most important predictor; however, its relatively small influence in comparison to wind speed suggests that the suitability of potential SZ polygons is not very robust to changes in wind conditions. Vegetation height is the fourth most important determinant of maximum within-SZ pSSD, meaning that SZs must take into account local vegetation conditions in order to ensure safety of wildland firefighters. The four SZ polygon shape metrics all had similar and relatively low importance, but their inclusion does have a significant effect on the amount of variance in maximum pSSD explained. For example, a random forest model built without these four metrics only explains 83% of variance in pSSD, whereas their inclusion increases that number to 94%. The ideal shape for a SZ would be approximately a circle, where SSD can be maintained from surrounding vegetation on all sides. So, a large circular SZ will possess very different suitability than an equal area linearly-shaped SZ (e.g., a road). Thus, the importance of SZ shape should not be underestimated.


**Table A2.** Variable importance for predicting maximum within-SZ pSSD, measured in terms of the percent increase in mean squared error that would result from the removal of that variable (%IncMSE).

#### **References**


## *Article* **Predicting Wildfire Fuels and Hazard in a Central European Temperate Forest Using Active and Passive Remote Sensing**

**Johannes Heisig 1, \* , Edward Olson <sup>2</sup> and Edzer Pebesma 1**


**Abstract:** Climate change causes more extreme droughts and heat waves in Central Europe, affecting vegetative fuels and altering the local fire regime. Wildfire is projected to expand into the temperate zone, a region traditionally not concerned by fire. To mitigate this new threat, local forest management will require spatial fire hazard information. We present a holistic and comprehensible workflow for quantifying fuels and wildfire hazard through fire spread simulations. Surface and canopy fuels characteristics were sampled in a small managed temperate forest in Northern Germany. Custom fuel models were created for each dominant species (*Pinus sylvestris*, *Fagus sylvatica*, and *Quercus rubra*). Canopy cover, canopy height, and crown base height were directly derived from airborne LiDAR point clouds. Surface fuel types and crown bulk density (CBD) were predicted using random forest and ridge regression, respectively. Modeling was supported by 119 predictors extracted from LiDAR, Sentinel-1, and Sentinel-2 data. We simulated fire spread from random ignitions, considering eight environmental scenarios to calculate fire behavior and hazard. Fuel type classification scored an overall accuracy of 0.971 (Kappa = 0.967), whereas CBD regression performed notably weaker (RMSE = 0.069; R <sup>2</sup> = 0.73). Higher fire hazard was identified for strong winds, low fuel moisture, and on slopes. Fires burned fastest and most frequently on slopes in large homogeneous pine stands. These should be the focus of preventive management actions.

**Keywords:** fuels; wildfire; fire behavior; fire hazard; remote sensing; LiDAR; Sentinel; modeling; simulation

#### **1. Introduction**

Climate change is expected to cause more extreme drought periods and heat waves throughout Central Europe [1]. Temperate regions, which are traditionally not prone to fire, may experience more and larger wildfires [2]. In recent years, Central Europe has already been exposed to multiple consecutive extreme drought events. Forest health in Germany has declined, increasing susceptibility to biotic and abiotic disturbances [3]. A significant rise in number and extent of wildfires has been reported for some parts of the country [4].

The vast majority of the wildfires in Germany are man-made (~95%) and tend to be much smaller than fires in typical fire-prone regions, such as North America [5]. High population density and related risks require instant fire suppression measures. In particular, the wildland–urban interface bears great risk potential [6]. However, in contrast to, for example, Mediterranean countries, in Central Europe, both the vegetation firefighting capacities and society's awareness of fire hazard are in an early stage of development. This gap between the increased susceptibility to fire through rapid environmental change and lagging fire suppression capabilities may hold challenges in the near future.

Forest fire behavior is driven by the interplay of three components: topography, weather, and fuels [7]. These components need to be quantified precisely to identify locations with elevated fire hazard potential. Topography and weather data are available

**Citation:** Heisig, J.; Olson, E.; Pebesma, E. Predicting Wildfire Fuels and Hazard in a Central European Temperate Forest Using Active and Passive Remote Sensing. *Fire* **2022**, *5*, 29. https://doi.org/10.3390/ fire5010029

Academic Editor: James R. Meldrum

Received: 29 January 2022 Accepted: 18 February 2022 Published: 20 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

and easily accessible in most cases. As fuels strongly vary in space and time, their detailed characterization often demands resource-intensive field observation [8]. Remote sensing data can help in modeling fuel data and predict continuous surfaces of fuel from punctual fuel samples [9]. In particular, dense airborne light detection and ranging (LiDAR) point clouds allow for three-dimensional assessments and help describe the vertical canopy structure [10].

Flame length (FL) and the rate of spread (ROS) are the most meaningful variables originating from fire behavior modeling [2]. Considered together, they can already support decision making. Spatial fire spread simulations may additionally provide insight into the probability of an area burning, given that a fire has started [11]. By integrating these spatial modeling outputs, one can draw conclusions about wildfire hazard [6]. Forest and fire management can benefit from this knowledge and initiate preventive measures according to local priorities [12]. In light of a prevailing need for action to adapt Central European forests to a changing climate, the efficient allocation of limited resources is crucial.

Fire behavior and hazard have been assessed in various case studies, predominantly in regions with substantial wildfire history. Botequim et al. [13] derived fuel characteristics from inventory data for even-aged management units in a Maritime pine forest in central Portugal. With the resulting fire hazard calculations, they established simple discrimination rules to implement fuel treatments. Taccaliti et al. [14] calculated basic fire behavior for black pine forests in Northeastern Italy. They used low-density LiDAR and field data to model the required fuel characteristics, yet omitted the estimation of crown bulk eensity (CBD). Stockdale et al. [15] simulated the effects of different mitigation scenarios on fire hazard in wildlife conservation areas in Alberta, Canada. The study could benefit from existing fuel maps (100 m spatial resolution) and a historic wildfire occurrence database.

This research demonstrates an end-to-end fire hazard assessment in a small study area with high-quality data from an extensive field survey and dense airborne LiDAR. It involves (i) the collection of surface and canopy fuel parameters in the field, (ii) spatial predictive modeling of fuels supported by remote sensing data, and (iii) the identification of locations vulnerable to wildfire through fire behavior and hazard modeling. Thus, we address local fuel conditions as precisely as possible, as no reliable data exists for this region. The study intends to provide a complete and structured workflow to produce spatial predictions of fire hazard in an area lacking historic fire records. Data and workflow details are shared in a way that is easy to reproduce, increasing accessibility to fire hazard information. Our study thus reacts to the current and future needs in Central European forest management and planning by supporting capacity building in the wildfire hazard domain.

#### **2. Study Area**

The Haard is a managed temperate forest located at 51.7° N, 7.2° E in North Rhine-Westphalia, Germany, with an approximate extent of 60 km<sup>2</sup> (Figure 1). Surrounded by small cities, infrastructure, industry and agriculture, it is highly frequented for recreational activities. Other land use types comprised by the Haard include a medical facility, a holiday resort, and an abandoned mine. Major roads and water bodies confine the forest in the north and west, while one public road crosses it from north to south. The study area is mostly flat, while moderate slopes can be found in northern and southwestern parts. Three dominant species form the majority of trees in the Haard, namely *Pinus sylvestris* (Scots pine), *Fagus sylvatica* (European beech) and *Quercus rubra* (red oak). All three species are economically relevant for timber production.

Considered individually, each of the three species has a different role in the wildfire hazard context. Often found on dry and sandy soils, Scots pine is adapted to limited water availability. Its crown structure allows light to reach the forest floor, which promotes the accumulation of understory vegetation and surface fuels. Scots pine is highly abundant in Northeastern Germany, the part of the country historically most affected by wildfire. European beech is the most common deciduous tree species in Germany. As a climax community species, it naturally outcompetes other species by shading the forest floor while

being shade-tolerant itself. Old beech stands produce only small surface fuel loadings dominated by litter with few woody components. Limited light availability entails higher soil moisture and sparse understory. Both conditions hinder fire spread. Young planted stands in the Haard on the other hand are very dense, row-wise plantations. Their crown base height is low, suggesting an elevated potential for flames to reach the crown via natural ladders. Red oak is native to eastern North America and became a solid silvicultural alternative to local oak species in Central Europe due to its superior growth rate. Its ability to cope with higher temperatures further makes it a popular choice for adapting forests to climate change. Regarding shading, understory, and live surface fuel loads, red oak lies in between Scots pine and European beech. Among others, these three species were reported to be strongly impacted by climatic extremes during the summers of 2003 and 2018, making them vulnerable to secondary drought effects in the following years [3]. It is conceivable that this pattern may reoccur in the future, raising their susceptibility to fire.

**Figure 1.** Our study area is a managed temperate forest dominated by Scots pine, read oak and European beech. Located at 51.7° N, 7.2° E (red dot) in densely populated North Rhine-Westphalia, Germany, it is surrounded by agriculture and industry, serving as a local recreation area. The map displays 215 surface fuel sampling locations, of which 30 were used to additionally sample canopy fuel characteristics. Wind properties were assessed at a nearby weather station.

#### **3. Materials and Methods**

#### *3.1. Field Data*

Surface and canopy fuel characteristics were sampled at points throughout the study area. The collected data enable the generation of fire behavior fuel models (FBFM) and serve as reference for predictive machine learning models.

#### 3.1.1. Surface Fuels

Surface fuel characteristics were assessed during a field campaign in the summer of 2019. A total of 215 plots were split among the three dominant tree species found in the Haard forest (beech: 73, pine: 73, oak: 69). Locations were distributed within the study area by generating random coordinates on a 25 m grid. The process was spatially stratified by a forest inventory map containing polygons of homogeneous stands dominated by either one of the three species. Each sample plot was visited physically, and centers were mapped via GPS.

Surface fuels were sampled using a standardized field method developed by Lutes and Keane [16]. Downed woody debris, humus, leaf litter and ground vegetation were sampled, using a compound cluster-plot of intercepts, point measurements and fixed-area plots. Downed woody debris was estimated based on the planar intercept principle of Brown [17], wherein three time-lag classes (1, 10 and 100 h) are defined according to diameter. Their respective particle counts along transects are translated to loadings expressed as mass per area. Humus and leaf litter depth were measured at multiple locations within each plot, averaged and converted into loadings by their respective bulk densities [16]. Similarly, vegetation loadings were estimated based on fixed-area plot observations of percent ground cover and average height of live and dead shrubby and herbaceous vegetation within the fuel bed (i.e., below two meters above ground). Field data were aggregated for each dominant species group using the median, resulting in three custom FBFMs.

#### 3.1.2. Canopy Fuels

Canopy fuel characteristics were assessed during a follow-up field survey in the summer of 2021. The main objective here was to quantify CBD, defined as the mass per unit volume of canopy biomass. It describes foliage and twigs with diameters less than 3 mm that would burn in a crown fire [18] and represents one of four canopy fuel parameters required for fire behavior modeling. Recorded field measurements include diameter at breast height (DBH), canopy height (CH), crown base height (CBH) and crown class. During this campaign, only 30 locations (10 per dominant species) were revisited. Plot locations were randomly selected from the existing set of surface fuels plots, stratified by dominant species. Each live individual within a 10 m radius around the plot center and a DBH greater 7 cm was considered. This led to recording a total of 695 individual trees. DBH was derived from circumference measurements, assuming tree stems to be perfectly cylindrical. Heights were approximated with a trigonometric method. It makes use of the fixed distance between the observer and stem and the inclination when targeting the top or bottom of the live canopy. Crowns were subjectively classified into categories by extent, height and growth form, as described by Lutes [19]. An individual was considered *dominant* if its height and extent were unmatched within the plot. Multiple trees forming the main canopy layer were considered *co-dominant*. The categories *intermediate* and *suppressed* were assigned to individuals that reached into the lower portion of the canopy or that were overtopped by surrounding crowns, e.g., due to competition.

Tree lists were input to FuelCalc [20] to derive plot-level CBD in m3/kg. Canopy characteristics described above served as input to species-specific allometric equations. This workflow for deriving CBD has been applied by other studies, such as that of Erdody and Moskal [21], as direct observations require a costly destructive sampling method [22]. FuelCalc produces CBD estimates in 1-foot (~30 cm) vertical bins and summarizes the plot-level result as the maximum 5-foot (~1.5 m) running mean [19].

#### *3.2. Remote Sensing Data*

#### 3.2.1. Sentinel-1 and -2

All Sentinel-1 and -2 tiles from 2019 covering the study area were processed to annual composites (Table 1), using Google Earth Engine [23]. The Sentinel-1 SAR GRD image collection (n = 120) was pre-processed following a framework by Mullissa et al. [24]. It applies border noise correction, speckle filtering and radiometric terrain normalization. Resulting bands include polarizations VV, VH and their ratio. Optical data from Sentinel-2 Level-2A were filtered by scene cloud coverage (<5%). Clouds in the remaining tiles (n = 24) were masked using the quality assessment (QA) band, which flags pixels that may be affected by normal or cirrus clouds. Both collections were reduced in the temporal domain by calculating three percentiles: 10th, 50th (=median), and 90th. The approach attempts to better capture the intra-annual variation in backscatter or reflectance. In particular, vegetation mapping may benefit from taking seasonal shifts of land surface properties into account.

#### 3.2.2. LiDAR

Open-access airborne laser scanning (ALS) data are available for the entire German state of North Rhine-Westphalia [25]. The latest acquisition covering the study area is from February 2020. Laser returns are declared to have a mean location error of 30 cm and a mean vertical error of 15 cm. Point cloud density is specified with 4–10 points per square meter. However, this value was found to be considerably higher in parts of the study area with a complex vertical vegetation structure. Among other meta data, each point carries its return number (first, last, nth) and classification codes for ground, non-ground and other categories, which may disqualify them for further analysis (e.g., noise).

ALS data are well suitable to extract forest structure metrics. The 3D point cloud is reduced to multiple 2D representations, describing canopy cover, height, density, or terrain. Thanks to the high point density, a partition into multiple height bins is feasible. These may have fixed widths or follow percentiles.

In total 74 metrics were derived from ALS (Table 1). Processing of the LiDAR data was facilitated by the R programming language [26] and the lidR package [27]. In total, 90 coherent 1 × 1 km tiles were retrieved from the administrations' web service and processed to 10 m resolution metrics. First, a digital elevation model (DEM) was computed from the raw point cloud. The simple triangulation algorithm only considers formerly classified ground points. Subsequently, slope and aspect were derived from the DEM using the *terra* package [28]. For further processing, the point cloud was normalized to eliminate the effect of varying elevation on vegetation structure calculations.

LiDAR-derived metrics produced for this analysis include several basic descriptive statistics that have been proven useful when modeling vertical structure for forestry applications [21,29–31]. A large portion of the metrics (>60%) is related to the height values of laser returns (Z-dimension). Some describe the point cloud as a whole (e.g., maximum, mean, standard deviation, or skewness). Others aim at revealing structural differences using (cumulative) height percentiles or the percentage of returns exceeding a certain threshold. Metrics describing the mean height of grasses, shrubs or trees as well as the vertical height gap follow the definitions introduced by the PROMETHEUS fuel type classification system [32].

Two slightly more complex metrics are the canopy height model (CHM) and CBH. Minor measurement errors in ALS allow both to be directly derived from LiDAR point clouds. The CHM was calculated using the pit-free algorithm [33]. It facilitates the triangulation of first returns in different height bins before assembling them. Further, it enriches each point with eight sub-circles that approximate the LiDAR footprint (~20 cm diameter) more realistically and create smoother results. CBH was calculated following a quantile-based approach proposed by Chamberlain et al. [31]. In theory, this method detects the height at which the first live branches are present, which intersect considerably more pulses than the stem below CBH.

Roughly one third of the metrics is devoted to canopy cover (CC). Their calculations follow the assumption made by Riaño [34] that vegetation cover is represented by the fraction of canopy (i.e., non-ground) returns from the total. Reversely, the fraction of ground returns can be obtained by subtracting CC from 1. Next to overall vegetation cover, a cumulative vertical profile at heights ranging from 0 to 30 m was computed. Again, the cover of grasses, shrubs or trees follow the PROMETHEUS fuel type classification

system and were introduced by Novo et al. [35]. As a proxy for canopy surface roughness, the Rumple index was included. It is defined as the ratio of canopy surface and ground area [36].

Two metrics are related to LiDAR return density. Besides the total number of returns per area, the density of first returns within the canopy was computed following the method by Andersen et al. [29]. Three metrics with zero variance (*Zp*[1,5,10] ) were excluded prior to the following analysis.

**Table 1.** Remotely sensed predictor variables included in this study (n = 119). The majority of datasets are derived from airborne LiDAR point clouds (n = 74) and describe vertical forest structure and terrain. Temporal composites (t. c.) were computed from all Sentinel-1 and -2 acquisitions using the 10th, 50th and 90th percentiles.


#### *3.3. Wind*

Fire behavior models require information on wind speed and direction. Both may influence fire behavior significantly. Wind characteristics were assessed using the rdwd package [39] and hourly data from a close by weather station (Figure 1). Only observations from the months most relevant to wildfire (June, July, and August) between the years 2000 and 2021 were considered.

During the summer months, the dominant wind direction is southwest (240° N). Maximum wind speeds of up to 10.7 m/s were reached, with an average of 2.1 m/s (Figure 2).

**Figure 2.** Wind direction and speed aggregated from hourly data in close proximity to the study area during months June, July, and August and from 2000 to 2021. The dominant wind direction was identified as southwest (240° N). Average wind speed is 2.1 m/s, whereas a maximum speed of 10.7 m/s was recorded.

#### *3.4. Fuels Prediction*

#### 3.4.1. Surface Fuels

To relate fuel models to spatially continuous surfaces, we used remote sensing data and machine learning frameworks. Field plot locations and their respective dominant tree species label served as training points for image classification. To delimit non-relevant (non-burnable) land cover types (water, agricultural area, urban area and bare ground), 79 training points were established manually with the support of high-resolution truecolor imagery.

Predictor variables comprised 119 layers at 10 m resolution from active and passive air- and space-borne sensors (Table 1). After pairing training locations with predictor variables, a random forest classification model was trained [40]. The number of trees was held constant at the default of 500. Hyperparameter *mtry* (number of candidates per split) was tuned with odd-numbered values ranging between 2 and 11. Model training was integrated in a forward feature selection (FFS) process as proposed by Meyer et al. [41]. The model is built successively, starting with the best-performing pair of predictors and adding further ones as long as the performance measure improves. This way, the resulting

model is as simple as possible, which is more desirable than having tens or hundreds of terms. The FFS model was optimized using the Kappa index. Overall accuracy (OA) and the class confusion matrix were further consulted for model evaluation.

If training points are clustered in space, a traditional random k-fold cross validation (CV) comes with caveats. Spatial auto-correlation may cause an overly optimistic view on model performance. In this case, however, field sampling locations were randomly distributed over a large part of the study area. This justifies the use of a simple 5-fold random CV, which was implemented in the FFS process.

As an additional criterion, the area of applicability (AOA) was computed for the final model. It aims at estimating the extent to which a prediction model and its CV error can be applied. The underlying dissimilarity index is calculated for each pixel based on its distance to the closest training point in the multidimensional predictor space [42]. Both FFS and AOA were computed using the CAST package [43], while individual model training was facilitated by caret [44].

Finally, the model prediction yielded a four-class raster dataset. A simple majority filter (3 × 3 moving window) was applied to smooth out minor classification errors and reduce their effect on subsequent fire spread simulations. Tree species classes were assigned corresponding fuel model parameters for fire behavior analysis. The remaining class represents all non-burnable land cover types and was thus not assigned any values relevant to fire behavior.

#### 3.4.2. Crown Bulk Density

Plot-level CBD values computed via FuelCalc [20] were subsequently used as target variable in a regression analysis. The same set of 119 predictors was used. Circular field plots (10 m radius) partially covered 4–9 pixels (10 × 10 m) of the predictor data. Therefore, all concerned values were extracted and weighted by their share of the plot area covered. Data were split into training and validation sets in a 70 to 30 ratio. Three predictor variables derived from Sentinel-2 (*B*01*p*10; *B*01*p*50) and ALS data (*Zp*15) were excluded, due to their near-zero variance. The training data used in this analysis had 30 observations and 119 predictors, disqualifying them for ordinary least squares regression. Many predictors are strongly correlated. This constellation calls for a regularized regression model, which aims at making a biased estimate of regression parameters. Bias is thus traded for variance [45].

We selected a ridge regression model for this task. Its regularization parameter *lambda* is dependent on the dataset. *Lambda* was tuned using a grid of 100 values ranging from 0.1 to 100. CV yielded the optimal value by considering the root mean squared error (RMSE) as a criterion. Ridge regression analysis was facilitated by R-packages glmnet [46] and caret [44]. The training data were centered and scaled prior to model training. A 5-fold random CV was implemented for the same reason as described in the previous section. Test runs revealed a non-normal distribution of model residuals. Consequently, the target variable was transformed, using the logarithmic function. In an attempt to expand the number of training samples, 15 plots of the German National Forest Inventory were added. Their tree lists were supplemented with LiDAR-derived CH and CBH and processed in FuelCalc. Predictor variables were extracted in the same way as for the remaining training points.

#### *3.5. Fire Behavior and Hazard Modeling*

Next to fuel loadings and fuelbed depth, FBFMs further consist of several constants, such as the surface-area-to-volume ratio (SAV), moisture of extinction (MOE), and heat content [47]. Constants were selected in accordance with existing standard FBFMs [47,48]. All fuel models received 1, 10, and 100 h SAV constants of 5906, 4249, and 4593, respectively. MOE was set to 30% and heat content to 18,622 Kj/kg.

FlamMap is a quasi-empirical fire spread model, widely used among forest and fire managers in operational settings [2]. Its main requirement is the fire landscape data, which consist of three terrain variables (elevation, aspect, and slope), four forest structure variables (CH, CBH, CC, and CBD) and a surface fuel classification [49]. Additional controls, e.g., fuel moisture, wind, or temperature, can be defined. FlamMap's objective is to calculate spatiotemporal fire behavior under a unique set of environmental conditions. This is especially interesting when assessing the effects of varying weather conditions, fuel properties, or fire-related forest management actions.

To evaluate fire hazard in the Haard under a variety of environmental conditions, we set up a grid of realistic scenarios. We selected two options each for the components of (i) fuel moisture, (ii) wind speed, and (iii) air temperature. This resulted in a total of eight individual scenarios (Table 2).

Like many other regions, Central Europe recently experienced extremely hot and dry summers. Extended heat waves and the absence of precipitation fostered an increased number of ignitions. Surface fuel moisture content dropped and caused higher rates of fire spread. To reflect variable moisture conditions and their effects on fire spread, Scott and Burgan [48] proposed a set of scenarios along with their standard FBFMs. Dead (D) and live (L) fuels can each be assigned four moisture levels, ranging from very low to low, moderate and high (1–4), which are expressed as representative moisture content values in percentage values. Dead fuel moisture scenarios are characterized for time-lag categories 1, 10 and 100 h and range from 3 to 14 percent. Live fuel moisture scenarios include values for herbaceous and woody fuels and range from 30 (fully cured) to 150 (fully green) percent. We selected fuel moisture scenarios D1L1 and D3L1, which should correspond to the effect of a longer and a shorter pre-fire drought period. Both are characterized by very low moisture levels for live herbaceous (30%) and live woody (60%) vegetation. However, they differ in dead fuel moisture levels for time-lag categories 1, 10, and 100 h with 3%, 4%, and 5%, and 9%, 10%, and 11%, respectively.

FlamMap uses either weather time series data or single value inputs for individual weather components. Wind speed and direction strongly affect fire initiation and spread [50]. Both can either be static inputs or may be adjusted dynamically using the WindNinja module [51]. This numerical micro-scale wind flow model is designed for fire modeling applications and simulates mechanical and thermal effects of terrain, vegetation, and air temperature on the flow. It scales down wind speed and direction to finer spatial resolutions.

To investigate the impact of different weather conditions on fire hazard, we considered two values each for wind speed and air temperature. We compare peak (10 m/s) and average (2 m/s) wind speed as well as extreme (35 °C) and moderate (25 °C) summer temperatures. Figure 3 displays exemplary WindNinja outputs for the study area from constant wind speed (10 m/s) and direction (240° N). Despite revealing local anomalies, down-scaled wind characteristics remained similar to their original input constants for the most part.

**Table 2.** Fire hazard modeling scenarios representing a range of environmental conditions. Eight combinations result from two fuel moisture scenarios (D1L1 and D3L1), two wind speed values (10 and 2 m/s), and two air temperature values (35° and 25 °C).


The two most critical and valuable landscape fire behavior outputs for decision makers are fire intensity (expressed as FL) and the ROS [2]. Both are calculated on a pixel basis. Considered individually, both characteristics bear a limited hazard potential. Stationary high-intensity fires or fast-spreading low-intensity fires only pose minor challenges for containment efforts. However, quickly moving high-intensity fires are difficult to manage and contain, resulting in an elevated fire hazard potential.

**Figure 3.** Wind direction (WD; (**A**)) in degrees north and wind speed (WS; (**B**)) in meters per second down-scaled and adjusted to the study area's terrain using WindNinja software. Constant input values were assessed in an exploratory analysis of weather station data. WD was fixed at 240° N, WS at 10 m/s.

In a second step, FlamMap employs the minimum travel time (MTT) algorithm to run fire spread simulations [52]. It requires a set of random or predefined ignition locations. Each fire burns under the same environmental conditions and stops after a specified duration. That makes MTT particularly useful for the analysis of effects of spatial patterns in fuels and topography [49]. Final perimeters are recorded for each simulated fire. Locations (pixels) that have burned more frequently throughout the simulation will receive higher conditional burn probabilities (CBP). This means they are more likely to burn in the case of a fire, but gives no indication of the probability of a fire starting. CBP values are scaled by their maximum and broken down into five ordinal classes, ranging from lowest to highest.

Next to CBP, MTT also produces a conditional flame length (CFL). It represents a weighted version of FL previously introduced in the context of landscape fire behavior. MTT distinguishes heading, flanking, and backing fire spread. Therefore, heading fire burns at a higher intensity than the other two. Throughout repeated fire spread simulations, every pixel can burn multiple times and with different intensities. MTT computes probabilities for 20 FL classes that sum up to 1. These are aggregated to a single value per pixel by applying probabilities as weights for the mid-point FL of each class (0.5, 1,..., 9.5, >9.5 m). The resulting CFL is considerably lower than FL from a landscape fire behavior calculation. Similar to CBP, CFL is subdivided into six standardized classes.

Integrated fire hazard (FH) is finally calculated from CFL and CBP, following the classification scheme proposed by the Interagency Fuels Treatment Decision Support System (IFTDSS) [53] (Figure 4). The scheme divides continuous fire behavior values into hazard categories, which are more intuitive to read. By combining both measures to a single metric, it offers a more holistic view on the present situation than other fire behavior variables.

A combination of high values in both categories leads to high FH, whereas high values in either one lead to medium FH at most. Equally to CBP, the five FH classes range from lowest to highest.

For each scenario, we ran fire spread simulations using the same set of 10,000 randomly distributed ignition locations. The maximum simulation time for MTT was set to 60 min. This constraint was selected to reflect a sensible reaction time for nearby fire departments. Local forest and fire management confirm that this value approximately matches reality. In addition to the non-burnable land cover class fire spread was restricted by barrier structures, such as major roads or water bodies surrounding the study area. These are represented by vector geometries, which serve as input data for the model.


**Figure 4.** Integrated fire hazard (FH) classification scheme proposed by the Interagency Fuels Treatment Decision Support System (IFTDSS) [53]. It categorizes continuous fire behavior variables (CFL and CBP) to present a single intuitive metric. The highest FH is only reached in places where both CFL and CBP are very high.

#### **4. Results**

#### *4.1. Surface Fuels*

Field sampling yielded surface fuel loadings for three custom FBFMs. Distributions of sampling data are displayed in Figure 5. Custom FBFMs are created from the median of each parameter (see Table 3). Significant differences between all fuel models can be observed for 1 h fuels. They are generally greater than 1 h loadings of similar standard FBFMs suggested by [48]. In turn, 10 h fuels are very similar among all groups. Red oak shows significantly more 100 h fuels with approximately twice the loading of pine or beech. Pine has the largest quantities of live fuels, especially for herbaceous vegetation. Beech, however, forming dense top canopies and shading the forest ground well, only shows small live shrub loadings and almost no live herbs. This pattern is reflected in the fuelbed depth. Pine fuelbeds are about twice as deep as the others.

**Table 3.** Fuel loading parameters from field sampling. Data were aggregated using the median to create fuel models for fire behavior calculation. Dead and live fuel loadings are measured in [kg/m<sup>2</sup> ], and height of the fuelbed is in [m].


**Figure 5.** Boxplot of sampled surface fuel parameters for each fuel model. Dead and live fuel loadings are measured in [kg/m<sup>2</sup> ], height of the fuelbed is in [m]. Beech and red oak produce slightly more 1 and 100 h fuels, while more live herbaceous biomass is found in pine stands. Higher-growing live fuels cause pine's fuelbed depth to exceed that of the others.

The FFS process found the best performing solution among 13,924 possible combinations of predictor variables. CV selected an optimal mtry of 2 and reported an OA of 0.971 with a Kappa of 0.967. The final model consisted of five variables, with *B*05*p*<sup>90</sup> being the most important one. *B*06*p*<sup>90</sup> and *Rumple* also scored high variable importance, followed by *B*05*p*10. Finally, *Zcov<sup>g</sup>* was only able to add very little improvement to the model. FFS successfully reduced the number of predictors significantly (5 vs. 119) while keeping model performance high. For comparison, a classic random forest model was built, using all available predictors. Performance was still high, but both OA and Kappa decreased by 2.8% and 3.8%, respectively.

The model confusion matrix revealed a perfect classification of non-burnable areas. This may be explained by noticeable differences in the spectral signatures and vertical structure of non-burnable land cover (e.g., water bodies, urban area) compared to forest vegetation. Minor confusion was found between tree species classes. Errors ranged between 4% and 8%, with beech having slightly larger errors than pine and red oak. Previously separated validation samples (n = 90) were tested on the model prediction. All samples were classified correctly, yielding an OA and Kappa of 1.0. Overall, the fuel model classification results were very good. However, this could be expected, as it was a simple modeling task combined with a large range of powerful predictors.

The final prediction is shown in Figure 6 alongside with the AOA. The dissimilarity of predictors from the training samples was only critical in areas less relevant for fire spread. Fewer than 1% of burnable pixels fell outside the AOA and represented mainly roads, bare ground or vegetation along water bodies. These areas are often located in close proximity to the forest edge or to existing fire barriers, which makes them less relevant for fire spread modeling.

**Figure 6.** Surface fuel model prediction (**A**) and respective area of applicability (AOA; (**B**)). Spatial patterns originating from forest management are clearly visible. Pine is the most abundant species in the study area (51%) followed by beech (32%) and red oak (17%). Burnable pixels falling outside the AOA sum up to less than 1%.

#### *4.2. Crown Bulk Density*

Field measurements of canopy fuel characteristics show clear differences between dominant species (see Figure 7). Many beech stands in the study area are young and have small DBH. They often occur in row-wise plantations with high stem density and nearly identical CBH and CH. Red oak and pine stands are older and less dense as a result of selective logging. Their DBH and CH are two to three times greater than those of beech. CBH, the height of the lowest live branch, is generally much lower for red oak compared to pine, while the opposite is the case for CH. As a consequence, the CBD estimates for red oak are lower.

**Figure 7.** Boxplot of canopy fuel parameters from field sampling for each fuel model. Parameters include crown base height (CBH; (**A**)), canopy height (CH; (**B**)), diameter at breast height (DBH; (**C**)), and crown bulk density (CBD; (**D**)). CBD was calculated via FuelCalc using the other measures as inputs. Beech stands sampled in this study area are mostly young with low, yet dense canopy layers. Contrarily, red oak stands are old with open crowns and low CBH, resulting in low CBD. Pine stands are in between the others, with higher CBH leading to higher CBD than red oak.

Ridge regression analysis was applied to predict CBD for the Haard forest. CV selected an optimal regularization parameter *lambda* of 10.7. Overall, model performance was poor which was expected, considering the small number of training samples. A model-R<sup>2</sup> of 0.59 with a RMSE of 0.054 was reported. Independent validation samples produced a higher R <sup>2</sup> of 0.73, while RMSE degraded to 0.069. Although R<sup>2</sup> is acceptable, RMSEs are large, considering a CBD training sample mean of 0.095. Differences in model performance and validation scores indicate the introduction of bias by ridge regression.

Variable importance scores indicated strong dependence on LiDAR-derived vertical structure metrics and optical predictor data. The most relevant predictors included *C*25, *B*05*p*10, *Zp*20, *NDV Ip*90, *B*04*p*90, and *DEM*. Further, the nine next relevant variables in the ranking, *Ziqr*, *Zpcum*80, and *Zp*65,...,95, all describe vegetation structure in the upper third of the tree. This coincides with relative heights at which CBD can be found at the maximum.

We tested adding training samples from NFI plots (n = 15) but were not able to improve the model. On average, their derived CBD values were significantly smaller than the existing values based on field sampling. NFI surveys include records of species and CH among many other observations. However, they do not include CBH. Supplementing NFI tree lists, for example, with LiDAR-derived CBH at 10 m spatial resolution, is rather inaccurate, especially when considering plots with heterogeneous species composition, age, and vertical structure.

Spatial prediction and AOA for CBD are shown in Figure 8. CBD values range from 0 to 0.3. They roughly follow the tree species classification, while higher densities can be observed for pine than for beech and red oak. Anomalies in CBD within homogeneous patches dominated by a single species are related to structural differences. Considering

only forested pixels, 20% fall outside the AOA. This may again be explained by the low number of training samples. A significant portion is located in areas with steeper slopes.

**Figure 8.** Spatial prediction of crown bulk density (CBD; (**A**)) and the respective area of applicability (AOA; (**B**)). Values strongly depend on the dominant tree species. Anomalies within stands are related to differences in forest structure. Pixels falling outside the AOA (20%) occur especially on slopes and are attributable to the low number of training samples for the predictive model.

#### *4.3. Fire Behavior and Hazard*

Fire behavior was computed for the Haard forest considering landscape properties, fuel metrics, and eight sets of weather conditions. Variations in fire behavior can be observed in the spatial domain. This is related to the heterogeneity of surface fuels, vertical forest structure and terrain. Further, wind speed and fuel moisture impact fire behavior, resulting in major differences among scenarios. Air temperature, however, had no significant effect on fire behavior. Therefore, only scenarios S1, 3, 5, and 7 are evaluated in the following.

Overall, the strongest fire behavior was found for S1 (Figure 9 row 1). Low fuel moisture and high wind speed led to FL exclusively larger than 2 m, sporadically even exceeding 10 m. Low to medium FL predominantly occurred in red oak and beech stands. On the contrary, pine stands in particular showed high fire intensities. The same pattern was observed for ROS. Even more so than FL, ROS was closely linked to the spatial patterns of surface fuel models. Spread rates ranging from 6 to beyond 10 m/min were common in pine-dominated stands, whereas broad-leaved stands produced significantly slower fires. Most areas in the Haard with steep slopes are populated by pine. The combined effects of fire-prone pine fuel beds and steep slopes result in the highest observations of FL and ROS within the study area.

**Figure 9.** Landscape fire behavior outputs including flame length (FL; (**A**)) and rate of spread (ROS; (**B**)) for scenarios S1, 3, 5, and 7 (rows 1, 2, 3, and 4). Fire behavior depends more strongly on wind speed than on dead fuel moisture. Even if wind speed is low and dead fuel moisture is high, individual locations in pine stands and on steep slopes show FL > 6 m and ROS > 8 m/min. NB = Non-Burnable.

The effect of high wind speed on fire behavior becomes visible when comparing scenarios S1 with S3 (Figure 9 rows 1 and 2). Overall, a reduction in wind speed from 10 to 2 m/s cut the mean FL in half from 5.1 to 2.5 m, while ROS dropped from 6.8 to 2.6 m/min. Lower wind speed in S3 diminished fire behavior disproportionately in pine stands, which showed extremes in S1. A weaker, yet still significant, reduction occurred in beech and red oak stands. Few locations with high FL and ROS in the central northern and southwestern parts of the study area still stand out in S3 and S7 (Figure 9 rows 2 and 4). These correlate well with increased slopes and related higher wind speeds, highlighting their promoting effect on fire behavior.

Comparison of S1 and S5 (Figure 9 rows 1 and 3) emphasizes the influence of dead fuel moisture. S5 features similar spatial patterns in both FL and ROS, with extreme values consistently occurring in steep pine stands. On average, however, both are reduced significantly when facing increased dead fuel moisture levels. Mean FL dropped to 3.2 m while ROS slowed down to 4.8 m/min. Similar, yet weaker, shifts apply to average wind scenarios S3 and S7. Hence, a more substantial reduction in fire behavior was recorded for a decrease in wind speed as opposed to a decrease in dead fuel moisture.

Fire behavior outputs from MTT fire spread simulations (Figure 10) show comparable patterns as landscape fire behavior results. While strong winds in S1 and S5 result in elevated CFL for the majority of the study area, average winds only lead to a few hotspots of small extent. CFL has substantially lower values than FL, which can be attributed to the weighted calculation method.

CBP for all scenarios was scaled, using the overall maximum of 0.013. Higher and highest CBP in S1 occurred in the center and the eastern part of the Haard. Both areas are characterized by large continuous pine stands indicating an interdependence. Within the limited simulation time, fire spread fastest and formed the largest perimeters when burning in this setting. Concerned locations consequently burned more frequently, leading to higher CBP. Increased fuel moisture content in S3 yielded similar spatial patterns, yet on a lower level without exceeding medium CBP. Disregarding few exceptions, S3 and S5 only produced CBP values belonging to the lowest category. This indicates that fire spread is very limited during the absence of strong winds. Even though spatial variations in CBP exist, they are hidden by the scaled classification scheme.

Similar to the fire behavior results, extreme FH only occurred during strong winds and with low fuel moisture present (Figure 11A). In case of a fire in these conditions, large pine stands in all parts of the study area appeared to be affected more severely than other species. Within pine stands, the highest FH category coincides with areas showing increased CBD or slopes. Regardless of fuel moisture, fire spread simulations featuring strong winds (Figure 11A,C) revealed a small number of FH hotspots.

Simulations for average wind scenarios S3 and S7 (Figure 11B,D) showed clearly reduced FH predictions. Low FH prevailed under both scenarios. While medium FH still emerged in pine-dominated areas for S3, it only sums up to a negligible extent in S7. It should be noted that in both cases, high FH can be found sporadically on steep slopes, despite the absence of strong winds.

**Figure 10.** Conditional flame length (CFL; column (**A**)) and conditional burn probability (CBP; column (**B**)) from MTT fire spread simulations for scenarios S1, 3, 5, and 7 (rows 1, 2, 3, and 4). CFL is presented as the mid-point of 20 FL classes weighted by probabilities. CBP was scaled by its maximum (0.013) and classified into five equal bins. Large continuous areas dominated by pine show the highest probability of burning if strong winds are present. NB = Non-Burnable.

**Figure 11.** Integrated fire hazard (FH) as a product of conditional flame length (CFL) and conditional burn probability (CBP) for scenarios S1, 3, 5, and 7 (**A**–**D**). High and highest hazard throughout the study area can only be expected for strong winds and low fuel moisture. Higher fuel moisture significantly reduces high FH and limits it to large homogeneous pine stands. FH is medium to low when wind speed is low, while large pine stands and slopes bear more hazard than other areas. NB = Non-Burnable.

#### **5. Discussion**

With this paper, we demonstrate a comprehensive workflow for the assessment of wildfire hazard in a small managed temperate forest. Surface and canopy fuels were sampled in field surveys and extrapolated in space, using predictive statistical modeling methods. We combined fuel information with other variables, completing the fire landscape to feed fire behavior models, simulate the growth of several thousands fires, and derive a hazard index.

Spatial predictions of fuel characteristics were performed based on field data and statistical modeling. Sufficient training data in combination with powerful predictor variables led to a near-perfect classification result for three surface fuel models. Dense ALS point clouds allowed for directly deriving canopy fuel metrics, such as CC, CH, and CBH, due to their relatively low measurement error.

CBD, on the contrary, required regression analysis using training data from field observations. This canopy fuel variable could not be measured directly but was estimated via allometric equations, using field measurements of other tree properties. Estimates then represented the reference for statistical modeling. Through this process, training data were subject to multiple sources of errors. Due to the limited number of samples, the resulting statistical model was weak. The value ranges of field and predicted CBD were similar to

previous studies by Andersen et al. [29] and Erdody and Moskal [21]. Their regression model validation produced RMSEs of 0.27 (R<sup>2</sup> = 0.84) and 0.024 (R<sup>2</sup> = 0.88) while including 101 and 57 field sampling locations, respectively. Cameron et al. [30] reported a RMSE of 0.098 (R<sup>2</sup> = 0.78) whilst using 52 training samples. Considering the low number of samples in the present study, a RMSE of 0.069 (R<sup>2</sup> = 0.73) was acceptable, especially as CBD is only one of multiple components of the FH assessment. Further, predicted CBD showed variations between stands with different dominant species, which is comprehensible from an ecological stand point. The ridge regression model introduced bias in exchange for reducing variance. This effect may have been reduced by applying a more complex tuning approach, such as nested cross validation. Alternatively, a better fit may be achieved through more flexibility, for example, by applying a penalized generalized additive model.

Errors in the prediction of all fuel variables may have been caused by the temporal gap between field (summers 2019 and 2021) and LiDAR (winter 2020) data acquisition. In addition, the collection of LiDAR data over deciduous trees during their defoliated phase is problematic. This may cause underestimation of CC and other metrics involving the number of canopy returns.

LiDAR-based predictor variables were crucial for modeling CBD and could contribute to the mapping of surface fuel models. Sentinel-2 data were able to support both efforts. In particular, bands 4, 5, and 6, and the NDVI were listed among the most important predictors. All involve wavelengths in the near-infrared, which are proven to describe vegetation characteristics. Temporal composites of the 10th and 90th percentiles were more important than the ones of the 50th percentile. They represent stages of the vegetative cycle, which reveal more differences among tree species than the annual median. The surface fuel classification model, for instance, took advantage of these differences by selecting *B*05*p*<sup>10</sup> and *B*05*p*<sup>90</sup> as two of its five predictors. Sentinel-1, on the contrary, did not play a relevant role. The SAR signal is able to partially penetrate the forest canopy and was expected to react to differences in CBD. Theses differences may have been too subtle to be captured, or were explained more precisely by other predictors.

As Papadopoulos and Pavlidou [50] stated, fire behavior and spread models require major assumptions, and their output accuracy heavily depends on the quality of input data. This study classified three FBFMs, each connected to one dominant tree species. Within stands, however, the species composition may slightly differ from observations made at sampling locations, altering surface fuel properties. Additional uncertainty in FH prediction arises from the accumulation of errors over the course of a fire growth simulation and the common problem of an unfeasible validation against real wildfire events [50]. Despite numerous sources of uncertainty, simulation-based FH is more robust and informative than simple fire behavior calculations, as it considers flanking and backing fire, not only heading fire [54]. CBP (and therein FH) gives no indication about the probability of a fire starting. As most fires are caused by humans, ignition density is likely to increase with proximity to man-made structures, such as roads, paths or recreational points of interest. Together with an accessibility assessment for fire suppression units, this information should be considered in any follow-up risk assessment.

The present analysis revealed wind speed and slope to be the most relevant factors of FH in the Haard forest. Further, low fuel moisture content can amplify FH. While high wind speed and low fuel moisture are difficult to prevent, forest management can address elevated FH on slopes. These areas could be converted to stands with low stem density, consisting of species that produce slow and weakly burning surface fuels (e.g., red oak). The same strategy could be applied to areas surrounding medical or recreational facilities, as they bear elevated risk. Other FH hotspots were identified in large pine stands. Among the three dominant species in the study area, pine's fuel bed produced the highest flame lengths and spread rates. Additionally, fires reached their largest perimeters in pine-dominated areas, resulting in the highest values for CBP. This may, in part, be related to the fact that local pine stands are larger and more contiguous than beech and oak stands, allowing fire to cover longer distances through a homogeneous and highly flammable

landscape. To mitigate fire spread in large pine stands, species composition and forest structure could be diversified, creating a more heterogeneous landscape. Fire barriers in the form of planted strips of broadleaved species could be introduced perpendicular to the dominant wind direction. The Haard is confined by roads and water bodies in the west and north. However, wind and, with it, fire predominantly move toward the east. FH calculations for high wind speeds suggest a need to secure the eastern forest edge to prevent potential wildfire from affecting adjacent campsites and agriculture.

Future research in Central Europe could aim at developing detailed fuels datasets with national or continental extents. These are needed to allow for wildfire hazard and risk assessments on various spatial scales. Previous studies as well as the present one have confirmed the importance of LiDAR data for the spatial prediction of wildfire fuels. ALS data, however, have high acquisition costs and are rarely accessible. Two current trends in LiDAR remote sensing may advance the quantification of fuel variables in the near future. The growing availability of drones with lightweight laser scanning instruments enables high temporal and spatial resolution in small areas. From a global perspective, space-borne LiDAR missions, such as GEDI, will improve not only biomass and CH mapping, but also CBH or CBD estimation.

#### **6. Conclusions**

We characterized the fire landscape in a small managed temperate forest by connecting field and remote sensing data. With this baseline, we calculated fire behavior, simulated fire spread and deduced spatial wildfire hazard information. By interpreting the results of this process, we identified critical areas that should receive special attention in management actions related to wildfire safety. With fire as an emerging threat to Central European forests, interest in such information may grow. We deliberately used openly accessible data and low-cost instruments for field sampling. Processing was facilitated through freely available software without advanced computational requirements. These attributes lower the hurdle for forest management to obtain FH information, following our approach.

**Author Contributions:** Conceptualization, methodology, data curation, formal analysis, validation, visualization, writing—original draft preparation, project administration, J.H.; investigation, J.H. and E.O.; writing—review and editing, J.H., E.O. and E.P.; supervision, E.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data, code, and instructions on how to reproduce this study can be accessed via https://github.com/joheisig/Haard\_Wildfire\_Fuels\_Hazard (27 January 2022).

**Acknowledgments:** We thank Harald Klingebiel at RWR Ruhr Grün for sharing his experience on fire behavior and history in the Haard. Further, we are grateful to Gregory Dillon, Duncan Lutes and Charles McHugh at the Missoula Fire Sciences Laboratory for their valuable advice on fire and fuels modeling. Thanks also goes to Lukas Blickensdörfer at the Thünen Institute for his help with forest inventory data. We acknowledge support from the Open Access Publication Fund of the University of Münster.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Considerations for Categorizing and Visualizing Numerical Information: A Case Study of Fire Occurrence Prediction Models in the Province of Ontario, Canada**

**Den Boychuk 1, \* , Colin B. McFayden 2 , Douglas G. Woolford 3 , Mike Wotton 4,5 , Aaron Stacey 6 , Jordan Evens <sup>4</sup> , Chelene C. Hanes 4,5 and Melanie Wheatley 5**


**Abstract:** Wildland fire management decision-makers need to quickly understand large amounts of quantitative information under stressful conditions. Categorization and visualization "schemes" have long been used to help, but how they are done affects the speed and accuracy of interpretation. Using traditional fire management schemes can unduly restrict the design of new products. Our design process for Ontario's fine-scale, spatially explicit, daily fire occurrence prediction (FOP) models led us to develop guidance for designing new schemes. We show selected historical fire management schemes and describe our method. It includes specifying goals and requirements, exploring design options and making trade-offs. The design options include gradient continuity, hue selection, range completeness and scale linearity. We apply our method to a case study on designing the scheme for Ontario's FOP models. We arrived at a smooth, nonlinear scale that accommodates data spanning many orders of magnitude. The colouring draws attention according to levels of concern, reveals meaningful spatial patterns and accommodates some colour vision deficiencies. Our method seems simple now but reconciles complex considerations and is useful for mapping many other datasets. Our method improved the clarity and ease of interpretation of several information products used by fire management decision-makers.

**Keywords:** colour coding; communication; forest fire; ordinal categorization; palette; risk; wildfire

## **1. Introduction**

Situational awareness and decision-making for operational wildland fire management is supported by a large amount of complex, numerical information, often covering large areas and sometimes spanning multi-day forecasts. Comprehending and interpreting that quantity of information under time-limited and stressful conditions is challenging. Among other ways, this task is commonly made faster and easier by categorizing and visualizing the numerical information. There are many ways to do so, but how it is done can help or hinder interpretation, highlight or obscure valuable information and accurately portray or distort the data.

**Citation:** Boychuk, D.; McFayden, C.B.; Woolford, D.G.; Wotton, M.; Stacey, A.; Evens, J.; Hanes, C.C.; Wheatley, M. Considerations for Categorizing and Visualizing Numerical Information: A Case Study of Fire Occurrence Prediction Models in the Province of Ontario, Canada. *Fire* **2021**, *4*, 50. https://doi.org/ 10.3390/fire4030050

Academic Editor: James R. Meldrum

Received: 18 June 2021 Accepted: 12 August 2021 Published: 18 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

The categorization of numeric indicators of potential fire activity has a long history in Canadian and other fire management agencies. Established categories are deeply integrated into fire operations and culture. Although categorization is useful, conforming to traditional schemes, such as a four-category blue–green–yellow–red sequence from low to extreme, may be less than ideal for new information products. This issue arose during our implementation of fine-scale, spatially explicit fire occurrence prediction (FOP) models. Our design process for a new scheme is the basis for this paper.

FOP is one of the pillars of situational awareness, a requirement for daily and multiday preparedness planning [1] and a key component for modelling risk [2]. Ontario's fire management agency has a lightning-caused FOP model [3] and, more recently, a humancaused FOP model [1] that was developed in collaboration between the agency's science specialists and decision-makers and external researchers.

Our objective is to provide guidance for the design of the categorization and visualization of complex information used in fire management. We begin with an overview of selected historical categorization and colouring schemes used in fire management. We then describe our method in steps including (1) specifying design goals and requirements, (2) exploring design options and seeing how they interact and (3) making trade-offs. We present a case study on designing the scheme for displaying the output of Ontario's FOP models. Although our method arose from this work, the considerations and general principles employed can be useful in many other situations. The Supplementary Material includes other applications of our method and extensions to some other considerations for the display of data. Note that we use the terms "category" and "categorization" synonymously with "class" and "classification", the latter set being conventions in fire management.

#### *Overview of Selected Historical Categorization and Colouring Schemes*

Knowing the historical origin of the schemes aids us to understand their limitations and engender improvements. Prime categorization examples are those for the outputs of two Canadian Forest Fire Danger Rating System (CFFDRS) [4] subsystems: the Fire Weather Index (FWI) System [5] and the Fire Behaviour Prediction (FBP) System [6].

The FWI System accounts for effects of past and present weather on fuel ignitability and fire behaviour and has six main numeric outputs. Three track moisture in different depths or sizes of fuel: Fine Fuel Moisture Code (FFMC), Duff Moisture Code (DMC) and Drought Code (DC). The other three indicate potential fire behaviour: Buildup Index (BUI), Initial Spread Index (ISI) and Fire Weather Index (FWI). FWI indicates the potential fireline intensity in a standard pine (genus *Pinus*) stand [5]. FWI is also used in Ontario as a general indicator of fire hazard or danger and is mapped into ordinal adjective classes: Low– Moderate–High–Extreme. The classes accompanied the introduction of the FWI System [5] and were implemented early on by Ontario [7]. Fire danger schemes are intended to "sound an alarm" about potential extreme behaviour and difficulty of control [8]. The class boundaries were set so that Extreme covered the worst 2% of historical days, and the remaining lower boundaries were set by a geometric progression [5]. For a more detailed history of the FWI classification methods in Ontario, see [9]. Other fire agencies in Canada calculated different boundaries according to their data [10,11]. The FWI System and component classifications have also been applied outside North America. For examples, see [12,13]. The current danger classes and colours used in Ontario for the FWI System's main outputs are outlined in Table 1. Examples of these are shown for roadside signs for public alerts about daily fire hazard (Figure 1a) and also for operational maps (Figure 1b). In mapping products such as Figure 1b, the FWI System values are calculated for weather station locations and are interpolated [14] across Ontario's fire management area, which is almost 50% larger than France.

**Table 1.** Current danger classes and colour scheme [15] used in Ontario for each of the six Fire Weather Index System components: Fine Fuel Moisture Code (FFMC), Duff Moisture Code (DMC), Drought Code (DC), Buildup Index (BUI), Initial Spread Index (ISI) and Fire Weather Index (FWI). **Table 1.** Current danger classes and colour scheme [15] used in Ontario for each of the six Fire Weather Index System components: Fine Fuel Moisture Code (FFMC), Duff Moisture Code (DMC), Drought Code (DC), Buildup Index (BUI), Initial Spread Index (ISI) and Fire Weather Index (FWI).


**Figure 1.** Examples of different displays of the Fire Weather Index (FWI) class used by the Ontario's fire agency: (**a**) a public roadside sign; (**b**) an agency operational map showing interpolated, colour-coded FWI and the raw FWI values for each weather station. **Figure 1.** Examples of different displays of the Fire Weather Index (FWI) class used by the Ontario's fire agency: (**a**) a public roadside sign; (**b**) an agency operational map showing interpolated, colour-coded FWI and the raw FWI values for each weather station.

The other major CFFDRS subsystem that has outputs commonly communicated using classes is the FBP System, which provides quantitative estimates of fire behaviour outputs [6]. A primary output is fire intensity, the rate of energy or heat release per unit time per unit length of a spreading fire front [16], which ranges from 1 to ~100,000 kW/m in Canadian conditions. Fire intensity is also commonly categorized into fire intensity classes (ICs). Rather than adjective classes (i.e., low–extreme) they were given five numeric labels (IC 1–5 [17]), and later six (IC 1–6 [18]). Higher IC numbers correspond with the higher intensity values, but the boundaries are not evenly spaced (Table 2). The most commonly used IC boundaries in Canada delineate distinct differences in fire type characteristics (for example, surface, torching or crowning) in mature jack pine (*Pinus banksiana* Lamb.) stands and the corresponding general effectiveness of different types of fire suppression activities (for example, hand tools, pumps and hose, airtankers). These ICs are described in the Field Guide to the FBP System [19] and the ICs are also used to map fire intensity by many fire management agencies, for various purposes (Table 2). The other major CFFDRS subsystem that has outputs commonly communicated using classes is the FBP System, which provides quantitative estimates of fire behaviour outputs [6]. A primary output is fire intensity, the rate of energy or heat release per unit time per unit length of a spreading fire front [16], which ranges from 1 to ~100,000 kW/m in Canadian conditions. Fire intensity is also commonly categorized into fire intensity classes (ICs). Rather than adjective classes (i.e., low–extreme) they were given five numeric labels (IC 1–5 [17]), and later six (IC 1–6 [18]). Higher IC numbers correspond with the higher intensity values, but the boundaries are not evenly spaced (Table 2). The most commonlyused IC boundaries in Canada delineate distinct differences in fire type characteristics (for example, surface, torching or crowning) in mature jack pine (*Pinus banksiana* Lamb.) stands and the corresponding general effectiveness of different types of fire suppression activities (for example, hand tools, pumps and hose, airtankers). These ICs are described in the Field Guide to the FBP System [19] and the ICs are also used to map fire intensity by many fire management agencies, for various purposes (Table 2).

Table 2 shows some of the different intensity values for higher-end class boundaries (i.e., additional thresholds beyond IC 6). The choice of colour when mapping can be operationally significant because colours covey information rapidly and have strong associations with levels of alarm—for example, red for danger [20] and green for calm [21]. Such psychological factors are not always considered in the visualization, however, which could lead to misinterpretation. For example, IC 4 in Table 2, which is associated with the upper limit of direct fire suppression effectiveness [17,22], is variously coloured a calming light green, a cautionary yellow or a warning orange. In Table 2 there are also cases where the same colour refers to different IC classes—for example, calming light green is used for IC 2 in Ontario, IC 3 in Alberta and IC 4 nationally. Table 2 shows some of the different intensity values for higher-end class boundaries (i.e., additional thresholds beyond IC 6). The choice of colour when mapping can be operationally significant because colours covey information rapidly and have strong associations with levels of alarm—for example, red for danger [20] and green for calm [21]. Such psychological factors are not always considered in the visualization, however, which could lead to misinterpretation. For example, IC 4 in Table 2, which is associated with the upper limit of direct fire suppression effectiveness [17,22], is variously coloured a calming light green, a cautionary yellow or a warning orange. In Table 2 there are also cases where the same colour refers to different IC classes—for example, calming light green is used for IC 2 in Ontario, IC 3 in Alberta and IC 4 nationally.

**Table 2.** Examples of the diverse classification of fire intensity and colouring used in Canada. The Ontario, Alberta and national schemes are used in daily maps. The British Columbia (BC) scheme is used in a static map of the 90th percentile of historical fire intensity. The Field Guide to the Fire Behaviour Prediction System (Field Guide) scheme is used in printed tables. Intensity classes IC 1–IC 6 are as defined in [19]; the higher classes are informal. The colours are approximate.

<sup>1</sup> Range applies to the full row except where overridden by the text in the cell. <sup>2</sup> The colours are limited by the printing ink used. <sup>3</sup> BC's IC 1–5 corresponds to BC's colour progression rather than the row label.

An additional caveat in the classification is that simplifying numerical information by aggregation into few classes has a cost. For FWI and fire intensity, the wide range within some classes is operationally significant for some decisions—for example, FWIs of 11–22 become "High". Furthermore, for interpolated maps, neighbouring points that are displayed as different classes will not have operationally meaningful differences. To compensate, maps often include the raw point values that were used for interpolation (Figure 1b). The numeric information cannot, however, be read and interpreted as quickly, thus reducing the benefit of categorization. The categorizing of data for spatial application such as FWI and FBP is common, and there are many recognized considerations (for example, [26]) and built-in solutions in geographic information systems. However, using a built-in classification option without a deep understanding may be unsuitable because potential distortions can lead to radically different interpretations, as others have noted [27]. Consequently, there is a strong need to use schemes that convey information accurately.

#### **2. Methods**

We propose five steps to categorize and visualize model outputs for use in fire situational awareness and decision-making:


We first describe these steps in general, below, and then with further detail on their application, in our case study. Although the method is described in a linear sequence, the work is partly concurrent and highly iterative, especially within Step 4.

#### *2.1. Step 1: Understanding and Scoping the Data*

Our method applies to data with a continuous numerical scale of measure (real numbers). With minor modifications, it can also apply to data with a discrete numerical scale of measure (integers) and to ordinal categorical data (for example, Very Low, Low, Low– Moderate, . . . ). The modification is that colouring with continuous gradients (described below) does not apply unless there are a great many discrete values or ordinal categories.

The technical details of the raw model output data may be straightforward, but unfamiliar units, scaling, storage or other conditions can lead to misinterpretation. The following need to be understood by the designers:

	- # It may or may not be useful to have more categories where there was a lot of data, and it may or may not be pointless to have multiple categories where there was little or no data (see examples in Section 2.2).
	- # The stored precision may not correspond with the accuracy. Measured data such as weather observations have low precision (for example, to 0.1 ◦C), but calculated data such as FWI System values should be calculated and may be stored with full machine precision (~16 decimal digits), which is well beyond the accuracy of typical fire management data.

Understanding data precision and accuracy is necessary for presenting information accurately and having decision-makers understand it easily and correctly. Regarding data storage and all subsequent calculations using data, full machine precision should be maintained to avoid accumulating rounding errors. Regarding the numbers displayed for decision-makers, the displayed precision should not exceed the data's accuracy, because that could be misleading. The numbers displayed for decision-makers should ideally have the lowest precision that is operationally significant to minimize unnecessary mental processing. Further discussion and examples are given in the Supplementary Material.

#### *2.2. Step 2: Understanding the Decision-Making Uses of the Information*

The purpose of the information is to support decision-making, so it is necessary to know who is using the information and how it is used. Working directly with fire management staff to understand their needs is necessary for ensuring that new model outputs are effectively integrated into the decision-making process [28,29].

There are key questions to consider. What decisions are being supported? Are certain parts of the range more important, needing higher attention? Is a higher resolution (smaller class size) needed in some parts of the range rather than others? For example, consider the categorizing and mapping of the accumulated 24-h rainfall from a precipitation radar [30]. For differentiating the degree and duration of the reduced fire behaviour potential, a high resolution is useful at the low end but not at the high end. Conversely, for differentiating the degree and duration of flood potential, a high resolution is useful at the high end but not at the low end. Moreover, a higher top category is appropriate.

#### *2.3. Step 3: Specifying the Design Goals and Requirements*

As stated in the introduction, the high-level goals of categorization and its visualization are simply to show the information completely and accurately and to have the information be understood quickly and easily. These goals are elaborated into criteria as follows.

	- # Are the magnitudes shown with the original or reduced precision?
	- # Are the magnitudes undistorted or distorted by categorization, scale nonlinearity or truncation?
	- # Are the relative magnitudes evident by the colouring without or with referral to the legend?

• Regarding the quick and easy understanding of information: Regarding the quick and easy understanding of information:

*Fire* **2021**, *4*, x FOR PEER REVIEW 6 of 20

Regarding the complete and accurate display of information:
