Examining Landscape-Scale Fuel and Terrain Controls of Wildﬁre Spread Rates Using Repetitive Airborne Thermal Infrared (ATIR) Imagery

: The objectives of this study are to evaluate landscape-scale fuel and terrain controls on ﬁre rate of spread (ROS) estimates derived from repetitive airborne thermal infrared (ATIR) imagery sequences collected during the 2017 Thomas and Detwiler extreme wildﬁre events in California. Environmental covariate data were derived from preﬁre National Agriculture Imagery Program (NAIP) orthoimagery and USGS digital elevation models (DEMs). Active fronts and spread vectors of the expanding ﬁres were delineated from ATIR imagery. Then, statistical relationships between ﬁre spread rates and landscape covariates were analyzed using bivariate and multivariate regression. Directional slope is found to be the most statistically signiﬁcant covariate with ROS for the ﬁve ﬁre imagery sequences that were analyzed and its relationship with ROS is best characterized as an exponential growth function (adj. R 2 max = 0.548, min = 0.075). Imaged-derived fuel covariates alone are statistically weak predictors of ROS (adj. R 2 max = 0.363, min = 0.002) but, when included in multivariate models, increased ROS predictability and variance explanation (+14%) compared to models with directional slope alone.


Introduction
Wildfires induce a variety of valuable ecosystem processes [1,2] but can inflict severe economic, social, and environmental losses. Of note, California's 2017 losses include USD 500 million in suppression costs, 47 lives lost, and over 9000 structures damaged or destroyed [3]. These statistics highlight the need to understand how fire behavior is related to underlying geospatially distributed environmental factors. Rate of spread (ROS), also referred to as fire spread rate, is a key focus of fire scientists, first responders, and fire management agencies and has been the subject of many studies seeking to quantify its dependence on environmental factors. Research on the mechanisms that govern wildfire spread are commonly conducted using laboratory [4][5][6][7] or outdoor fire experiments [8][9][10], where inputs may be controlled by researchers. Findings discerned from such studies have been extended to larger spatial scales using numerical models that simulate wildfire spread in various environments and conditions e.g., [11][12][13][14].
Laboratory studies and fire spread models have utility, however, advancements are stymied by the lack of quality landscape-scale data needed for adequate fire spread theory, model validation, and model calibration [15][16][17][18][19]. There is also a need for improved empirical quantification of the environmental drivers of wildfire spread and their control on ROS at landscape-scales during extreme wildfire events (EWE) [20][21][22][23][24]. Better understanding and near real-time estimation of wildfire ROS frequency distributions at landscape-scales could help fire scientists, managers, and emergency responders focus future research, evacuation, suppression, and mitigation efforts [19,[25][26][27][28]. the few studies besides our own that used relatively high spatial resolution remotely sensed data to empirically analyze ROS during an extreme wildfire event. They estimated ROS from time-sequential fire perimeters delineated from ortho-rectified aerial photographs during the Riba de Saelices fire in Spain and used ROS as a covariate (among many others) for predicting burn severity. Storey et al. [70] used ATIR imagery to statistically analyze spotting behaviors during several large fires in Australia, however their work was focused on analyzing spotting occurrences rather than full fire spread runs.
Repetitive-ATIR imagery collected over several recent California wildfires, made it possible to map active fire front locations and estimate fire rates of spread at high levels of spatial and temporal detail. We statistically examined landscape-scale controls on fire spread rates collected from several chaparral and oak woodland landscapes based on a statistical-empirical approach utilizing ATIR image sequences, high spatial resolution optical imagery, and DEMs. Landscape covariate relationships with ROS estimates were analyzed with regression. Within the context of primarily canopy fire spread in shrubland vegetation, the following research questions are addressed: 1.
How effective are spectral vegetation indices (SVIs) and fractional cover of growth form types, derived from high spatial resolution aerial orthoimages, as spatially explicit surrogates for fuel load distributions, in predicting ROS? 2.
How strongly associated are directional slope angle and slope orientation relative to fire spread direction with ROS?

Study Areas
Study areas for this research correspond to the spatial coverage of sequential ATIR imagery captured during the Detwiler and Thomas Fire events in California, as illustrated in Figure 1. The Detwiler Fire burned over 30,000 ha of Sierra Nevada foothills in Mariposa County from 16 July through 9 January 2018 [3]. The Thomas Fire burned from 4 December 2017 to 12 January 2018, in Santa Barbara and Ventura counties [3]. It is the second largest recorded fire (single ignition source) in California history, burning over 113,000 total hectares across generally rugged mountainous terrain. The study areas are composed of chaparral, coastal sage scrub, oak woodland, and mid-elevation coniferous forest ecosystems. can be used to precisely estimate fire rates of spread and directions [63,64]. Viedma et a [54] is one of the few studies besides our own that used relatively high spatial resolution remotely sensed data to empirically analyze ROS during an extreme wildfire event. The estimated ROS from time-sequential fire perimeters delineated from ortho-rectified aeria photographs during the Riba de Saelices fire in Spain and used ROS as a covariate (amon many others) for predicting burn severity. Storey et al. [70] used ATIR imagery to statisti cally analyze spotting behaviors during several large fires in Australia, however thei work was focused on analyzing spotting occurrences rather than full fire spread runs.
Repetitive-ATIR imagery collected over several recent California wildfires, made i possible to map active fire front locations and estimate fire rates of spread at high level of spatial and temporal detail. We statistically examined landscape-scale controls on fir spread rates collected from several chaparral and oak woodland landscapes based on statistical-empirical approach utilizing ATIR image sequences, high spatial resolution op tical imagery, and DEMs. Landscape covariate relationships with ROS estimates were an alyzed with regression. Within the context of primarily canopy fire spread in shrubland vegetation, the following research questions are addressed: 1. How effective are spectral vegetation indices (SVIs) and fractional cover of growt form types, derived from high spatial resolution aerial orthoimages, as spatially ex plicit surrogates for fuel load distributions, in predicting ROS? 2. How strongly associated are directional slope angle and slope orientation relative t fire spread direction with ROS?

Study Areas
Study areas for this research correspond to the spatial coverage of sequential ATIR imagery captured during the Detwiler and Thomas Fire events in California, as illustrated in Figure 1. The Detwiler Fire burned over 30,000 ha of Sierra Nevada foothills in Mari posa County from 16 July through 9 January 2018 [3]. The Thomas Fire burned from December 2017 to 12 January 2018, in Santa Barbara and Ventura counties [3]. It is th second largest recorded fire (single ignition source) in California history, burning ove 113,000 total hectares across generally rugged mountainous terrain. The study areas ar composed of chaparral, coastal sage scrub, oak woodland, and mid-elevation coniferou forest ecosystems.  Table 1).
The study areas are characterized by moderate Mediterranean climate regimes typical of much of the Southern California region. Hot and dry summers are succeeded by mild wet winter with average annual precipitation ranging from 32 to 53 cm. Both study areas were delineated and analyzed using the Universal Transverse Mercator (UTM) coordinate system-zones 10 and 11 north.

Data
ATIR imagery collected by Kolob Canyon Air Services using a FireMapper TM 2.0 thermal infrared imaging system was used for ROS measurements for both study fires. FireMapper TM 2.0 is a noncryogenic sensor with a 320 × 240 frame array [64,71]. The imagery was captured in short-, mid-, and long-wave infrared wavelengths, yielding images with 5-15 m ground sample distance (GSD) depending on flight altitude above ground level [71]. The long-wave TIR imagery was utilized in this study. Imagery was acquired over the same fire front with a "race-track" pattern as described in Stow et al. [64]. All images include a GPS location and time stamp that provides each frame with temporal and geographic coordinate metadata. FireMapper imagery was geometrically corrected with ERDAS IMAGINE Photogrammetry Software Application [72] using onboard positional (GPS) and altitude (inertial motion unit) data to yield georeferenced image frames. Stow et al. [64] found that the coregistration of sequential pairs of geoprocessed FireMapper TM 2.0 imagery to be approximately one-pixel root mean square error. ATIR image sequences were prioritized for ROS analyses based on a high number of repetitive, short-interval (<12 min) flight-passes over active fire fronts burning within varied topography (different angles and orientations of slope facets during front progression) to sample a range of fire spread dynamics. ATIR observations analyzed in this research pertain to samples of maximum forward fire spread. A summary of ATIR image metadata for the sequences analyzed in this research is shown in Table 1. Prefire visible/near infrared aerial imagery from the National Agricultural Inventory Program (NAIP) was obtained from USGS Earth Explorer (https://earthexplorer.usgs. gov/) and Google Earth Engine (GEE) [73] to map and quantify fuel characteristics. The NAIP image sets were captured between July and August 2016, and consist of 0.6 m spatial resolution, four-band visible and near-infrared georeferenced orthoimages. Topographic data were derived from National Elevation Dataset (NED) 10 m spatial resolution digital elevation models (DEMs). These DEM data were obtained from the GEE data cloud.
Weather data were retrieved from the Remote Automated Weather Station (RAWS) online data library (https://raws.dri.edu/wraws/scaF.html) and from the FireBuster model Thomas and Detwiler fire forecasts [74]. FireBuster is a fire weather forecast system that provides spatially explicit weather estimates based on downscaled weather model and digital terrain data inputs and are available through an interactive webtool or downloadable 1 km raster grids (https://fwxfcst.us/firebuster/) [74]. The RAWS data and FireBuster estimates obtained for this study consisted of hourly wind speed average, wind speed max, wind direction, and relative humidity estimates nearest to the space/time domain of the study areas and ATIR image collection times. Distances from RAWS to study areas ranged from a minimum of 5.5 km (Thomas 1, 2, and 3) to a max of 13.3 km (Detwiler).

Fire Feature Delineation and Landscape Sampling Units
Fire fronts were delineated from ATIR imagery using methodologies developed by Stow et al. [63,64] (Figure 2a). ATIR imagery was contrast enhanced to aid in active fire front detection and delineation. Images influenced by large amounts of smoke and hot gases were processed using a Laplacian edge filter to enhance fire front locations and aid in their delineation. ArcGIS Pro ver. 2.5.0 [75] software tools were used to interactively delineate fire front line curves in the form of polyline features that contain attribute information associated to the front length (m), front ID, and time of ATIR image capture. Fire fronts were manually delineated in ATIR imagery by drawing polylines at the sharpest radiometric temperature descent (i.e., gradient between fire and non-fire (ambient) temperatures) ( Figure 2a).
Fire 2020, 3, x FOR PEER REVIEW 5 downloadable 1 km raster grids (https://fwxfcst.us/firebuster/) [74]. The RAWS dat FireBuster estimates obtained for this study consisted of hourly wind speed average, speed max, wind direction, and relative humidity estimates nearest to the space/tim main of the study areas and ATIR image collection times. Distances from RAWS to s areas ranged from a minimum of 5.5 km (Thomas 1, 2, and 3) to a max of 13. (Detwiler).

Fire Feature Delineation and Landscape Sampling Units
Fire fronts were delineated from ATIR imagery using methodologies develope Stow et al. [63,64] (Figure 2a). ATIR imagery was contrast enhanced to aid in activ front detection and delineation. Images influenced by large amounts of smoke an gases were processed using a Laplacian edge filter to enhance fire front locations an in their delineation. ArcGIS Pro ver. 2.5.0 [75] software tools were used to interact delineate fire front line curves in the form of polyline features that contain attribute mation associated to the front length (m), front ID, and time of ATIR image capture fronts were manually delineated in ATIR imagery by drawing polylines at the sha radiometric temperature descent (i.e., gradient between fire and non-fire (ambient) peratures) (Figure 2a). Fire spread-vectors were represented as linear features connecting sequentia fronts and their origin, spacing, and direction was identified on a sequence-by-sequ basis [64]. Evenly spaced (30 m) [64] points along a time = n (start) front curve were matically generated by local-normal polyline connections from each point to the inte tion of the time = n+1 fire front (Figure 2b). Spread vectors were assigned geometry a utes including line bearing (0-360°), distance (m), and ROS estimate (m min −1 ). Fire sp vectors were also used as landscape sampling units (LSUs) for extracting topographi fuel data for subsequent statistical analyses.  Fire spread-vectors were represented as linear features connecting sequential fire fronts and their origin, spacing, and direction was identified on a sequence-by-sequence basis [64]. Evenly spaced (30 m) [64] points along a time = n (start) front curve were automatically generated by local-normal polyline connections from each point to the intersection of the time = n + 1 fire front (Figure 2b). Spread vectors were assigned geometry attributes including line bearing (0-360 • ), distance (m), and ROS estimate (m min −1 ). Fire spread vectors were also used as landscape sampling units (LSUs) for extracting topographic and fuel data for subsequent statistical analyses.

Orthoimage Processing
Three spectral vegetation indices (SVIs) were generated from NAIP orthoimages and tested as fuel loading covariates with ROS. NAIP imagery utilized in this study were not calibrated to surface reflectance values. Normalized Difference Vegetation Index (NDVI-U) (Equation (1)), Green-Red Vegetation Index (GRVI-U) (Equation (2)), and Normalized Difference Red-Blue (NDRB-U) (Equation (3)) images were created for all fire sequences based on the following formulae: where in Equations (1)- (3), -U represents uncalibrated to surface reflectance, NIR, RED, GREEN, and BLUE are uncalibrated (surface reflectance) NAIP digital number values for near-infrared, red, green, and blue wavebands, respectively. These SVIs were selected because they exploit different waveband combinations associated with the NAIP data, and normalized indices tend to suppress terrain-related illumination and image brightness effects [76]. NAIP imagery was also used to classify and map vegetation GF types as surrogates for fuel load. Input data for classification started with generating a false-color composite image of NDVI-U (Equation (1)), Visible Brightness (VB) (Equation (4)), and the red/green band ratio (RG) (Equation (5)) [77].
where RED, GREEN, and BLUE are uncalibrated NAIP digital number values for red, green and blue wavebands, respectively. NDVI-U, VB, RG thresholds were established interactively for classifying the following GF and land cover classes: (1) shrub, (2) herb, (3) tree, and (4) rock/bare soil. Map accuracy was assessed by comparing classification products with reference data generated from visual interpretation of 100 randomly sampled NAIP pixels and prefire Google Earth imagery. The overall accuracy of image classification products was 91%. Misclassified GF pixels identified during accuracy assessment were manually edited and recoded. Polygons were drawn around misclassified pixels and reclassified to reflect the correct GF type visually observed in reference imagery. Reclassification was followed by reconducting accuracy assessment and the procedure was repeated until an accurate map was produced.

Topographic Data Processing
Slope angle and orientation (relative to fire spread direction) were assessed as landscape topographic covariates with ROS. Standard slope functions found in GIS software lack the ability to generate nonstatic directional slope (i.e., relative slope in a predefined direction). To overcome this, we developed a DEM processing routine with Python ver. 3.0.2 and the RichDEM terrain analysis library [78] designed to measure the slope inclination relative to the fire spread direction associated with spread vector bearings. Slope-degree, slope-aspect, and fire-direction (bearing direction of vector) rasters were calculated for each LSU [79]. The slope angle and aspect, and fire-direction rasters for each LSU were then computed and processed into the directional slope covariate using Equation (6): where DS is the output directional slope raster, S is a raster grid representing slope degree values, VD is the spread vector bearing (degree) raster, and A is the aspect in degrees.
The output represents slope angle as a signed degree slope, where negative (downhill) and positive (uphill) values are calculated relative to the spread vector directions. Rasters having 0 • values signify flat slopes.

Landscape Covariate Sampling and Stratification
Zonal sampling, validation, and organization of covariate data were conducted using R v.3.6.2 and ArcGIS scripting and software tools and organized as a relational database using PostgreSQL ver. 9.6.2. LSUs were created spatially coincident to spread vector using the "Buffer" tool in ArcGIS Pro. Like the spread vectors, LSUs were assigned spatial geometry attributes consisting of the fire spread bearing, the length of the unit, and the unit's size (m 2 ). The LSUs constituted the basis for automatically extracting topographic and fuel data for subsequent statistical analyses. LSU buffers were created to capture DEM and NAIP pixels spatially coincident to spread vector polylines. The mean, minimum, maximum, range, standard deviation, median, and variance of directional slope and SVI pixels were extracted for each LSU using the R "exactextractr" package [80].
GF fractional cover (GFFC) estimates for LSUs were calculated with a raster model that converts like-classified, contiguous GF pixels into polygons. The GFFC percentage of each LSU was quantified on a scale of 0.0-1.0. GFFC estimates were also used to assign a "LSU Fuel Class" type to each sample unit (Table 2 a), based on combinations of fuel/vegetation classification schemes derived from prior studies including Anderson [34], Sandberg et al. [81], and Blodgett et al. [53]. GFFC were stratified into one of two separate slope orientation groups: Slope angle > 0 • (upslope), Slope angle < 0 • (downslope) based on the mean directional slope ( Table 2 b). Google Earth was used to visualize 3-D topography and remove from statistical analyses LSUs associated with spread vectors where fire spread both up-and downslope between imaging passes. In total, 195 LSUs were removed from analyses for this reason. Samples of GFFC were stratified according to Table 2 b and regressed. Directional slope angles from all sample units were paired with fuel classes according to Table 2 c and regressed. Directional slope samples were also stratified and regressed by orientation (upslope and downslope).

Statistical Analyses
Several statistical models were employed to examine relationships between both stratified and unstratified landscape covariates and ROS estimates. These included: (1) bivariate linear regression, (2) multiple-stepwise linear regression, (3) power (log-log) regression, and (4) exponential (semi-log) regression. All models were fit and analyzed using R ver. 3.6.2 [82]. Exploratory data analyses of ROS estimates and landscape covariates were conducted to summarize data quality, distribution, and covariance [83]. Sample means, variances, minima, maxima, standard deviations, kurtosis, skewness, and correlation statistics were used to evaluate all landscape covariates separately, and cumulatively. Independent variables that exhibited significant control on ROS estimates were initially evaluated using parametric and nonparametric difference-of-means/median tests, including: two sample t-test, Wilcoxon rank sum, analysis of variance (ANOVA), and Kruskal-Wallis tests [83,84].
Bivariate linear regression models evaluated direct relationships between individual landscape covariates and ROS estimates, for each spread sequence and covariate (stratified and unstratified). All bivariate models were evaluated and compared using beta coefficients, standardized beta coefficients, adjusted coefficient of determination (adj. R 2 ), Akaike's Information Criterion (AIC), and p-value diagnostics [83,85].
Multiple linear regressions were run to determine the cumulative influence and interaction of covariates on ROS estimates. Forward-and backward-stepwise regression models were constructed using the lowest AIC, lowest Mallows CP, highest adj. R 2 , and significant F-statistic [83]. All ordinary least squares (OLS) models were evaluated using ANOVA tests to check whether they were statistically significant or not through corresponding Fand p-values [85]. Standardized coefficients of parameter estimates were also evaluated to determine the direct influence of individual covariate estimates. Covariates in multiple regression models were also tested for multicollinearity using the variance inflation factor (VIF) [86]. Multiple regression models were fit with all landscape covariates minus one GF type. One GF type at a time was removed from multiple regression models using hierarchical partitioning to clear the compositional nature of GFFC estimates [87]. The removal of any growth form variable was determined on a sequence-by-sequence basis using regression goodness-of-fit diagnostics derived from the heir.part R package (Adj. R 2 , Log-Likelihood, and root mean square error) [88] as precursors for candidacy in multiple regression models.
Throughout the OLS modeling process, we looked for violations of model assumptions and if found, data log-transformations and/or different techniques were used to model relationships. First, plots were used to validate assumptions of normality, linearity, and equity of variances [83]. Residual Quantile-Quantile (QQ) plots were used as an initial test for residual normality. To confirm QQ-plot visualizations the Shapiro-Wilk test [83] was performed to test residual normality. Residuals were also plotted against predicted ROS (y) values to determine any residual clustering or obvious patterns that indicated violation of model assumptions [83].

Fire Spread Behavior
Maps of fire front progressions for each fire spread sequence are shown in Figure 3. Frequency plots of ROS estimates with rose plots of fire spread direction frequency for the five sequences are presented in Figure 4. Although fire spread rates varied substantially within and between ATIR sequences, all sample distributions are comparably right-skewed with mean rates of spread greater than the median. Median rates of spread ranged from a low of 4.

Wind Conditions During Fire Sequences
RAWS data and FireBuster predictions associated with ATIR collection areas and periods are summarized in Table 3. The lowest recorded wind speeds were reported during Thomas 1, 2, and 3 (0.4 m/s NNW). FireBuster wind speed predictions for Thomas 1, 2 and 3 are on average 0.77 m s −1 faster than RAWS measurements. Wind direction (RAWS) for Thomas 1 and 2 range were S-SSW during the first 1.5 h of spread, and NNW for the remainder of the sequences. The fastest wind speed (8.9 m s −1 ) and lowest average reported relative humidity (12.7%) was reported during Thomas 4. Detwiler wind speeds during ATIR imaging ranged from 4.5 to 4.9 m s −1 (RAWS) and 3.1 m s −1 (FireBuster). Wind direction reported from RAWS nearest to Detwiler (13.3 km) exhibit NNE winds. FireBuster predictions for Detwiler record WNW winds. Relative humidity for the Detwiler sequence was 23-24% (RAWS) and 14.8-16.3% (FireBuster). Overall, wind speed, wind direction,

Wind Conditions During Fire Sequences
RAWS data and FireBuster predictions associated with ATIR collection areas and periods are summarized in Table 3. The lowest recorded wind speeds were reported during  FireBuster predictions for Detwiler record WNW winds. Relative humidity for the Detwiler sequence was 23-24% (RAWS) and 14.8-16.3% (FireBuster). Overall, wind speed, wind direction, and relative humidity data and predictions from RAWS and FireBuster forecasts varied little over the time domain of all sequences.

Fuel Covariate Relationships with ROS
Statistical models for exploring SVIs as spatially explicit surrogates for fuel load covariates for ROS yielded few significant correlations. Bivariate regression diagnostics by fire and ATIR sequence yielded very weak and mostly nonsignificant relationships. For the Detwiler and Thomas 4 sequences, all SVIs are statistically significant with ROS. Of the statistically significant SVI models, the maximum total explained variance is only 5.4% (NDVI-U, Thomas 4). Weak but significant SVI relationships for Detwiler and Thomas 4 are attributed to greater spatial variability in fuel coverage and type compared to other sequences. For the study areas (sequences) and SVIs examined in this research, SVIs alone were not effective predictors of ROS. Weak statistical relationships likely stem from limited spatial variability in fuel composition and density within the burning areas captured during ATIR imaging. GFFC distributions for the four Thomas fire sequences were largely characterized by homogenous, high density of chaparral shrubs. Although shrubs predominantly cover landscapes burned during the Thomas sequences, sparse cover of tree, rock/bare soil, and herbaceous vegetation types are observed, based on image-classification products. Notably, portions of Thomas 1, 2, and 4 sequences contained dense corridors of riparian vegetation along riverbeds cutting through the sequence extents. The diversity of GFs is greater for the Detwiler fire sequence compared to Thomas fire sequences. LSUs associated with the Detwiler sequence are characterized by similar shrub, herbaceous, and tree GFFC. The degree of model fit between GFs and ROS estimates is greater when samples containing zero percent fraction are removed from sample populations. Like SVIs, GFFC is weakly associated with fire spread rates. However, model coefficients and adj. R 2 metrics indicate they were slightly more effective proxies for describing spatial distributions of fuel loads and control on ROS. For example, tree GFFC explained 36% and 12% of the total variation in ROS at Detwiler and Thomas 2, respectively. Similarly, rock/bare soil beta coefficients exhibit a negative relationship with ROS estimates at Detwiler, Thomas 1, and Thomas 2 (−14.53 to −23.45), explaining 7-14% of the variance in ROS.
Stratifying GFFC samples by slope angle (upslope > 0 • and downslope < 0 • ) was largely ineffective at isolating contributions on ROS. Similar to SVI relationships with ROS estimates, weak association with GFFC (stratified and unstratified) may be attributed to limited spatial variation in prefire fuel composition within the burning areas captured by the ATIR image sequences. This observation is supported by the predominantly homogenous shrub cover portrayed in maps generated from prefire orthoimagery as depicted in Figure 5.

Topographic Covariate Relationships with ROS
Generalized trends of slope angle for entire image sequences of successive fronts are characterized by average upslope inclinations of 19.54°, 10.53°, 13.50°, and 11.06° for Thomas 1, 2, 3, and 4, respectively ( Figure 4). Fire progression for the Detwiler sequence primarily occurred downslope with an average decline of −10.01°. Of all covariates examined, directional slope is the most statistically significant predictor of ROS estimates for all study fires and ATIR sequences. Linear, exponential, and power regression diagnostics on directional slope are reported in Table 4. For Thomas sequences 1, 2, and 3, about 50% of the variance of ROS estimates are explained by directional slope (adj. R 2 ≥ 0.500) ( Table 4). Regression beta and standardized beta coefficients for all regression methods and study areas show positive linear to weakly nonlinear relationships between ROS and directional slope (β > 1). Scatter plots comparing linear, exponential, and power regression fits with directional slope data are shown in Figure 6.
Linear model coefficients are similar for Thomas sequences 1, 2, and 4 (β range of 0.066), while semi-log (exponential) regression coefficients are similar for Thomas 1, 2, and 3 (β range of 0.003). Exponential relationships for Detwiler and Thomas 4 exhibit similarity. Power regression models show similar degrees of variance explanation of ROS on directional slope as by exponential models (avg. adj. R 2 difference of 0.010). Plots of ROS on directional slope data exhibit abrupt increases in ROS per incremental increase in slope angle at or around 20° (Figure 6), confirming the strength of weakly nonlinear relationships exhibited by power and exponential regression metrics (Table 4). Comparable model fits between exponential and power regressions types are found for Thomas 1, 2, and 3 (Figure 6b-d). Linear models for Detwiler (adj. R 2 = 0.160) and Thomas 4 (adj. = R 2 0.194) yielded higher model fit than exponential and power regression counterparts. Stronger linear regression fit for the Detwiler and Thomas 4 data are confirmed when comparing regression lines-of-best-fit (Figure 6a,e). Overall, exponential and power functions were most effective for characterizing the relationship between directional slope and ROS (adj. R 2 ≥ 0.500), and consistently accounted for the steeper rise in spread rate commonly observed at slope angles of 20° and greater.

Topographic Covariate Relationships with ROS
Generalized trends of slope angle for entire image sequences of successive fronts are characterized by average upslope inclinations of 19.54 • , 10.53 • , 13.50 • , and 11.06 • for Thomas 1, 2, 3, and 4, respectively ( Figure 4). Fire progression for the Detwiler sequence primarily occurred downslope with an average decline of −10.01 • . Of all covariates examined, directional slope is the most statistically significant predictor of ROS estimates for all study fires and ATIR sequences. Linear, exponential, and power regression diagnostics on directional slope are reported in Table 4. For Thomas sequences 1, 2, and 3, about 50% of the variance of ROS estimates are explained by directional slope (adj. R 2 ≥ 0.500) ( Table 4). Regression beta and standardized beta coefficients for all regression methods and study areas show positive linear to weakly nonlinear relationships between ROS and directional slope (β > 1). Scatter plots comparing linear, exponential, and power regression fits with directional slope data are shown in Figure 6.
Linear model coefficients are similar for Thomas sequences 1, 2, and 4 (β range of 0.066), while semi-log (exponential) regression coefficients are similar for Thomas 1, 2, and 3 (β range of 0.003). Exponential relationships for Detwiler and Thomas 4 exhibit similarity. Power regression models show similar degrees of variance explanation of ROS on directional slope as by exponential models (avg. adj. R 2 difference of 0.010). Plots of ROS on directional slope data exhibit abrupt increases in ROS per incremental increase in slope angle at or around 20 • (Figure 6), confirming the strength of weakly nonlinear relationships exhibited by power and exponential regression metrics (Table 4). Comparable model fits between exponential and power regressions types are found for Thomas 1, 2, and 3 (Figure 6b-d). Linear models for Detwiler (adj. R 2 = 0.160) and Thomas 4 (adj. = R 2 0.194) yielded higher model fit than exponential and power regression counterparts. Stronger linear regression fit for the Detwiler and Thomas 4 data are confirmed when comparing  (Figure 6a,e). Overall, exponential and power functions were most effective for characterizing the relationship between directional slope and ROS (adj. R 2 ≥ 0.500), and consistently accounted for the steeper rise in spread rate commonly observed at slope angles of 20 • and greater.
All models stratified by upslope (>0 • ) and downslope (<0 • ) orientations are statistically significant. However, compared to nonstratified models, greater model fit only resulted for Detwiler and Thomas 4 (average adj. R 2 increase of 0.200). All upslope stratified models demonstrate markedly higher standardized and unstandardized beta coefficients. The opposite is true for downslope stratified groups. To minimize potential bias in slope findings incurred from samples having disparate GF types, slope samples were stratified by the LSU fuel classes shown in Table 2 c: (1) shrub, (2) mixed-shrub, (3) herb, (4) mixed-herb, (5) rock/bare soil, (6) mixed-rock/bare soil, (7) tree, (8) mixed-tree, and (9) full mix. All slope models stratified by the shrub and mixed-shrub fuel classes are statistically significant for all sequences. Minor increases in slope model fits occurred within the shrub and mixed-shrub classes for Thomas 1, 2, and 3, while most R 2 values and coefficient estimates are comparable to nonstratified models. Slope stratified by the mixed-shrub fuel class yielded the highest model fit for Detwiler (adj. R 2 = 0.422) and Thomas 4 (adj. R 2 = 0.396). However, when slope data are stratified by other fuel classes, few models are significant. Detwiler was the only sequence where directional slope stratified by the herb fuel class was statistically significant (adj. R 2 = 0.223, p = 0.002). Although the stratification of slope yielded few significant results for Thomas 1, 2, and 3 (compared to unstratified models), increase in variance explanation for Detwiler and Thomas 4 were common.

Multivariate Analyses
Multiple forward/backward stepwise linear regression models were run with all covariates as initial inputs for each sequence to gain insights on variable interactions and combined influences on ROS (Table 5). Directional slope remained the most significant predictor of ROS (average std. β coef. of 0.381). Except for the Detwiler sequence, NDVI-U was included as a significant covariate with ROS. Similarly, GRVI-U was a marginally significant predictor of ROS for two of five sequences. The Thomas 3 model explained the greatest variance in ROS estimates (54%) and included directional slope, shrub, and NDVI-U as significant covariates. Levels of statistical significance for GFFC varied largely by sequence. Rock/bare soil were included in multiple regression models for Detwiler and Thomas 2 and Thomas 4. Herb fraction for Detwiler and Thomas 1 models exhibited interaction with other covariates and ROS, and thus were included in multivariate models. Absolute std. β coefficient averages for fuel covariates were 0.12 and 0.14 for SVIs and GFFC, respectively. Overall, when compared to directional slope, all image-derived fuel covariates contributed only a small fraction to ROS variance explanation (12-17%).

Discussion
Weather, fuel, and topographic conditions and properties are the factors that control wildfire spread [29][30][31][32][33]. Most empirical studies investigating environmental controls on fire spread are limited to laboratory and outdoor fire experiments, or studies based on coarse-scale satellite images of natural events e.g., [4,8,9,66]. The primary objectives of this research were to evaluate landscape-scale terrain and fuel controls on fire spread rates derived from repetitive-ATIR imagery collected during portions of two extreme wildfire events. This study directly builds on Stow et al. [63,64] by linking detailed ATIR wildfire spread measurements to geospatial data representing fuel and topographic distributions derived from prefire NAIP orthoimagery and USGS DEMs.

Fuel Covariate Findings
Fuel properties which influence fire behavior are difficult to measure and map due to their high variability in time and space [37][38][39]. For this reason, image-derived spectral vegetation indices and growth form maps were tested as spatially explicit proxies of fuel load distributions. The low and mostly nonsignificant correlations between ROS and imaged-derived fuel covariates suggest that (1) neither the SVIs or GFFC were stable predictors of ROS for the particular wildfire events and landscapes studied here, (2) fuel load was not a significant control on spread rates at the study sites examined in this research, or conversely, (3) the temporal resolution of NAIP imagery produced inaccurate estimates of the spatial characteristics of vegetation at the study areas imaged within the extreme wildfire events. Similarly, wind speed, which cannot be incorporated at the space-time scales of our landscape statistical analysis, may dominate the unexplained variance and weaken any influence of fuel covariates on ROS (see discussion on weather below). However, when comparing model coefficients between separate study areas, GFFC depicts more consistent relationships with ROS estimates than SVIs. For example, a negative relationship was commonly found for study areas containing LSUs with larger tree or rock/bare soil fractions (std. β average of −0.613). Conversely, relationships between SVIs and ROS varied largely by study area, and no clear relationship was discernable from model coefficients. Some past studies examining the applicability of SVIs as surrogates of fuel load or biomass using moderate to coarse-satellite images report promising results e.g., [36,39,40], indicating that SVIs are more closely associated with ROS at landscape-scales when coupled with terrain variables [42][43][44] and derived from orthoimagery captured just prior to the burn event [36,58].
The inclusion of growth forms and SVIs in multivariate models did increase the predictability and variance explanation of ROS. Cross-validated R 2 values routinely increased over standard adj. R 2 metrics in multivariate models. Similarly, the stratification of directional slope by fuel class (derived from GFFC estimates) increased variance explanation of ROS for Detwiler and Thomas 4. This supports what is generally known about fire behavior, that variations in fire spread rates are likely the result of complex landscape-scale interactions between fuel, terrain, and weather characteristics [29,37,89]. The statistical significance of bivariate and multivariate regression models for fuel covariates at Detwiler and Thomas 4 show separate controls on ROS at specific times of burning within the extreme wildfire events. Detwiler and Thomas 4 both exhibited significantly higher FireBuster and RAWS wind speeds compared to other sequences ( Table 4). The ancillary wind information coupled with fuel covariate findings at Thomas 4 indicate that fire spread controls during extreme wildfire events occur at landscape scales [53]. Further, homogeneity of vegetation cover type has been linked to greater fire size and ROS by Viedma et al. [90] and Holsinger et al. [91]. To further investigate this point, we analyzed samples in the top 95th percentile of ROS estimates and found most of the associated LSUs had shrub or herbaceous cover fractions between 90% and 100%. LSUs linked to lower ROS percentile groups (e.g., bottom 20th percentile) contained larger variation of GF type. These findings correspond with many former studies, in that fuel heterogeneity or fragmentation can restrain fire propagation and ROS [90,91] and should be investigated further.

Topographic Covariate Findings
Slope directly regulates fire spread through energy transfer of flaming biomass along fire fronts and its relationship with ROS has been postulated as a curvilinear function [4,[49][50][51]. Many findings indicate slope steepness has a dramatic effect on fire spread [49,50], however, the effects of slope on ROS in chaparral fuels are largely limited to laboratory settings due to the difficulty of isolating its influence from other variables (wind and fuels) during outdoor experiments and extreme wildfire events [14]. To evaluate nonstatic terrain controls on wildfire spread rates derived from the high-temporal resolution ATIR imagery, we designed a customized directional slope covariate to calculate signed-slope angles in the relative direction of fire propagation. We also attempted to isolate directional slope by orientation and fuel classes derived from GFFC estimates ( Table 2 c).
Directional slope is the most explanatory covariate for all study sequences, accounting for more than 50% (adj. R 2 ≥ 0.500) of the variance in ROS estimates for three of five of the image sequences ( Table 4). All regression models were statistically significant for directional slope on ROS for all sequences. Slightly higher model fits for exponential and power regressions with directional slope suggest the relationship with ROS is best described as a positive nonlinear function. An important finding from our topographic analyses on ROS is that exponential regression results were consistent between study areas. For example, constants in functions derived from exponential regression (or a in ae bx ) suggest the y-intercepts of the separate models varied little among study areas (5.6-7.4 m min −1 ). Growth rates of ROS (per unit increase in directional slope) were also consistent between separate study areas. For instance, growth factors (or b in ae bx ) in exponential functions only varied from 3.0% for Detwiler to 4.8% for Thomas 2 (∆1.8%).
We further examined directional slope data for LSUs associated with the top 95th percentile of ROS estimates and discovered maximum rates of spread coincide with slope angles between 20 and 40 • ; none of these samples associated with rapid fire spread had negative slopes. These findings contradict those from Vega et al. [92] and Catchpole et al. [93], who have proposed that slope controls on ROS in shrubland fires is less than forests [14,52].
While the influence of directional slope on ROS is clear and consistent within the contexts of the study fires and space and time scales of analysis, its importance as a controlling variable may vary as a function of space and time scales [14,52]. We characterized the average ROS and slope trends for the four Thomas Fire sequences which were all generally upslope fire progressions (unlike the Detwiler Fire sequence), as shown in Table 6. The ROS per degree of upslope trends at characteristic length scales of several km (rather than tens to a few hundred meters for the LSU analyses) are also listed in Table 6. These trend characteristics show that the relationship between ROS and slope are consistent between Thomas 1, 2, and 3 (imaged on the same day). The larger ratio between ROS and slope at Thomas 4 indicates that wind was likely the predominant control on fire spread. This is confirmed by the higher wind speeds recorded at Thomas 4 and the sequence's stronger relationships between fuel covariates and ROS. In this context, the influence of wind on ROS could also apply to the Detwiler sequence, which also exhibited a higher statistical significance with fuel covariates, lower significance with directional slope, and higher average wind speed. An equally important inference is that the directional slope covariate used in this research may be capturing topographic influences on surface winds [54,58].  Excluding the Detwiler and Thomas 4 sequences, stratifying directional slope by orientation and fuel class showed similar findings as unstratified models. However, our qualitative observations of fire spread along ridgetops and valley bottoms (Figure 7) indicate that the stratification scheme used in this research did not isolate relationships effectively; or high homogeneity of vegetation type and cover largely negated the stratification approach. Moritz et al. [89] postulated that topographic "fences" and "corridors" may impede or facilitate wildfire spread and occurrence during extreme wildfire events. Our qualitative observations of topography influencing fire spread from geovisualizations of active fire fronts derived from ATIR time sequences (Figure 7) support this assumption and are also consistent with past studies describing how convective or radiative mechanisms of wildfire spread are (1) strongly influenced along ridgetops and valley bottoms, and/or (2) these features have a large effect on surface winds and fuel conditions [48,58,89,94]. Moreover, the restriction of fire spread along these topographic features do support past findings that fire "fences" occur on landscape scales [89][90][91].

Weather-Related Findings
Most of the evidence associated with weather effects on fire spread for this study is limited to general ROS trends relative to wind direction, velocity, and relative humidity from RAWS observations and FireBuster forecasts that correspond most closely to the space-time domains of our image-based fire spread estimates. RAWS data and FireBuster predictions associated with ATIR collection periods and locations are reported in Table 4. The spatial, and especially temporal scales of both data types are too coarse to be analyzed statistically in a direct manner with ROS. However, the FireBuster predictions and RAWS

Weather-Related Findings
Most of the evidence associated with weather effects on fire spread for this study is limited to general ROS trends relative to wind direction, velocity, and relative humidity from RAWS observations and FireBuster forecasts that correspond most closely to the space-time domains of our image-based fire spread estimates. RAWS data and FireBuster predictions associated with ATIR collection periods and locations are reported in Table 4. The spatial, and especially temporal scales of both data types are too coarse to be analyzed statistically in a direct manner with ROS. However, the FireBuster predictions and RAWS data do help shed light on regression analyses.
Besides reported winds for Thomas 4 and Detwiler, RAWS and FireBuster reports for other sequences show low wind speed. The incorporation of topographic and meteorological data in the spatially explicit FireBuster model [74] are also likely to produce estimates of weather variables that are more indicative of the site-specific fire spread conditions where ATIR image sequences were captured, though at coarser spatial and temporal scales. For example, wind direction reported from nearest RAWS to Detwiler (13.3 km North) are directionally opposite to fire spread directions documented by spread-vector bearings. Contrastingly, FireBuster predictions for the time and location of the Detwiler sequence characterize wind directions consistent with fire spread direction.
The most meaningful finding related to general weather versus ROS relationships is that wind speeds and ROS were substantially higher, and relative humidity lower during the Thomas 4 and Detwiler sequences than for the Thomas 1-3, as shown in Table 4. The data and estimates reveal that Santa Ana weather conditions prevailed during the Thomas 4 sequence on 09 December 2020, but not necessarily on the previous data when the other Thomas sequences were imaged. Although the burn area associated with the Thomas 4 sequence, just 15 km west of Thomas 1-3, the relationships between directional slope and ROS are substantially different between the sequences that exhibit relatively higher wind speeds (Detwiler and Thomas 4). Further, the mean and median ROS estimates for Thomas 4 are considerably higher than other sequences. Higher ROS for Thomas 4 is likely associated with the faster wind speeds reported by nearest RAWS data and FireBuster estimates. This suggests that fire behavior for Thomas 4 and Detwiler were more winddriven and the higher wind speeds may have substantially reduced the relative influence of slope on ROS [14].
The Thomas Fire covered a vast, topographically diverse area and occurred during a two-week period intermittent Santa Ana weather condition. Thus, while some fire growth areas and periods were wind-driven, such as the initial run on December 4-5, others were topographically sheltered from the Santa Ana winds and instead, driven by 'fire-induced winds' [95], or intermittently exposed to the strong ambient winds [96].

Conclusions
Through this study we demonstrated that repetitive ATIR imagery from the FireMapper TM 2.0 imaging system facilitates attainment of broader knowledge of the relationships between extreme wildfire behavior and controlling environmental factors. Additionally, such landscape-level estimates of ROS provide could provide an excellent source of reference data for validating fire spread theory and simulation models [18][19][20]. Continuing ATIR data collection during extreme wildfire events in disparate ecosystem, terrain, and climate types is of the utmost importance, and would encourage the establishment of a comprehensive data set for future fire spread research and validation purposes.
We also document procedures for relatively short-interval, time-sequential ATIR imagery collection with the FireMapper 2.0 system for quantifying ROS during extreme fire events and evaluating environmental controls on fire spread at landscape scales. We identified covariate relationships, interactions, and influences with ROS estimates that help validate and expand upon similar work conducted for larger and smaller scales. We provide evidence from the 2017 Thomas and Detwiler fires that directional slope can be a strong control variable on ROS in shrubland landscapes, and its influence can be effectively expressed as a curvilinear function. We postulate that calculating directional slope in advance of an active fire front, coupled with dominant wind data, could provide fire managers and modelers better ability to predict wildfire spread and focus fuel treatment projects, similar to findings from Viedma et al. [54] and Coen et al. [58].
This study is the first attempt at isolating and statistically analyzing landscape-scale covariate relationships with high spatial and temporal resolution fire spread measurements collected during extreme wildfire events. The Detwiler and Thomas fire events burned over 140,000 ha of land, and spread through diverse variations of terrain, fuel, and weather. Observations and results drawn from this study are necessarily limited to the specific fire spread behaviors and landscape characteristics associated with the study fires captured during ATIR imaging. Moreover, the influence of weather on ROS estimates is challenging to account for in the statistical analyses due primarily to the nature of capturing wind speed and direction data for individual LSUs and at the temporal resolution of ATIR imaging during such events. However, the effects of Santa Ana wind and humidity conditions can be discerned from the differences in ROS and model fits with directional slope when comparing results for Thomas sequence 4 relative to sequences 1 through 3 [96].
Follow-on research pertaining to spatial-scale and nonlinear relationships between ROS estimates and landscape covariates should be conducted. This could include buffering spread vectors to generate larger landscape sample units that may account for uncertainty in actual fire spread directions [64]. Spatial association effects on statistical analyses could be evaluated using Geographically Weighted Regression (GWR) and/or Eigenvector-based Spatial Filtering (ESF) regression techniques. Nonlinear, nonparametric machine learning approaches such as regression trees (RT) and random forest (RF) regression should also be explored. Similarly, this research focused on analyzing long runs of fire spread that were easily tractable from ATIR image sequences. However, focus should be given to the spatial relationships and impact of spotting behaviors on ROS that were omitted in this research to focus on the full fire runs captured by ATIR sequences.
Future research on imaged-derived fuel covariate relationships with ROS would benefit from studies of wildfires in different ecosystems with different and particularly more heterogenous fuel compositions. This would also help determine if the weak relationship between ROS estimates and fuel covariates were limited by homogenous fuel conditions within the study fire sequence zones. The inclusion of prefire Light Detection and Ranging (LiDAR) data for the characterization and mapping of fuels would also likely enhance assessment of fuel controls.
The predictive capability of directional slope and its statistical relationship with ROS and interaction with wind speed should be explored further. Moreover, wind direction and velocity data at spatial and temporal scales closer to those of the ATIR image sequences would allow better understanding of weather influences on ROS. Evaluating the influence of terrain features that impede or enhance wildfire spread could also be examined by exploring ROS relationships with topographic index metrics [97]. Funding: This research was funded by the National Science Foundation (NSF), Division of Social, Behavioral and Economic Research, Geography and Spatial Sciences program grant G00011220. NCAR is sponsored by NSF. Any opinions, findings, and conclusions or recommendations expressed in this material are the authors' and do not reflect the views of NSF.

Data Availability Statement:
The data that support the findings of this study are available upon request.