Modeling sub-boreal forest canopy bulk density in Minnesota, USA, using synthetic aperture radar and optical satellite sensor data

National estimates of canopy bulk density (CBD; kg m−3) for fire behavior modeling are generated and supported by the LANDFIRE program. However, locally derived estimates of CBD at finer scales are preferred over national estimates if they exist, as the absolute accuracy of the LANDFIRE CBD product is low and varies regionally. Active sensors (e.g., lidar or radar) are better suited for this task, as passive sensors are ill equipped to detect differences among key vertical fuel structures, such as coniferous surface fuels (≤2 m high) and canopy fuels above this threshold—a key categorical fuel distinction in fire behavior modeling. However, previous efforts to map CBD using lidar sensor data in the Superior National Forest (SNF) of Minnesota, USA, yielded substandard results. Therefore, we use a combination of dormant-season synthetic aperture radar (SAR) and optical satellite sensor data to (1) expand detectability of coniferous fuels among mixed forest canopies to improve the accuracy of CBD modeling and (2) better understand the influence of surface fuels in this regard. Response variables included FuelCalc output and indirect estimates of maximum burnable fuel based on canopy gap fraction (CGF) measured at ground level and 2 m above ground level. SAR variables were important predictors of CBD and total fuel density (TFD) in all independent model calibrations with ground data, in which we define TFD as the sum of CBD and primarily live coniferous surface fuel density (SFD) 0 to 2 m above ground. Exploratory estimates of TFD appeared biased to the presence of sapling-stage conifer fuel on measures of CGF at the ground level. Thus, modeling efforts to calibrate SFD with satellite sensor data failed. Both CGF-based and FuelCalc-based field estimates of CBD yielded close unity with satellite-calibrated estimates, although substantial differences in data distributions existed. Estimates of CBD from the widest CGF zenith angle range (0 to 38°) correlated best with FuelCalc-based CBD estimates, while both resulted in maximum biomass values that exceeded those considered typical for the SNF. Model results from the narrowest zenith angle range (0 to 7°) produced estimates of CBD that were more in line with values considered typical. LANDFIRE’s estimates of CBD were weakly, but significantly (P = 0.05), correlated to both narrow- and wide-angle CGF-based estimates of CBD, but not with FuelCalc-based estimates. The combined use of field estimates of CBD, based on indirect measures of CGF according to Keane et al. (Canadian Journal of Forest Research 35:724–739, 2005), with SAR and optical satellite sensor data demonstrates the potential of this method for mapping CBD in the Upper Midwest, USA. Results suggested that the presence of live, coniferous surface fuels neither confounds remote detection nor precludes mapping of CBD in this region using SAR satellite sensor data, as C- and L-band idiosyncrasies likely limit the visibility of these smaller understory fuels from space. Nevertheless, research using direct measures of burnable SFD for calibrations with SAR satellite sensor data should be conducted to more definitively answer this remote detection question, as we suspect substantial bias among measures of CGF from ground level when estimating SFD as the difference between TFD and CBD.


Conclusions:
The combined use of field estimates of CBD, based on indirect measures of CGF according to Keane et al. (Canadian Journal of Forest Research 35:724-739, 2005), with SAR and optical satellite sensor data demonstrates the potential of this method for mapping CBD in the Upper Midwest, USA. Results suggested that the presence of live, coniferous surface fuels neither confounds remote detection nor precludes mapping of CBD in this region using SAR satellite sensor data, as C-and L-band idiosyncrasies likely limit the visibility of these smaller understory fuels from space. Nevertheless, research using direct measures of burnable SFD for calibrations with SAR satellite sensor data should be conducted to more definitively answer this remote detection question, as we suspect substantial bias among measures of CGF from ground level when estimating SFD as the difference between TFD and CBD.

Background
Canopy bulk density (CBD) is routinely cited as a salient biophysical forest canopy fuel variable-along with canopy base height, stand height, and percent canopy cover-required for modeling crown fire behavior and risk (Finney 1998;Scott 1999;Keane et al. 2005;Keane et al. 2006;Krasnow et al. 2009;Engelstad et al. 2019). In this study, CBD represents burnable mass per unit canopy volume (kg m −3 ) available to crown fires, primarily live foliage (needles) and branches less than 3 mm in diameter and dead branches less than 6 mm in diameter (Van Wagner 1977;Alexander 1988;Finney 1998;Keane et al. 2005). Live broadleaf tree foliage is not considered critical canopy fuel due to low overall flammability (Van Wagner 1977). Thus, exclusion of broadleaf fuels from CBD calculations is standard practice (Cruz et al. 2003). Spatially explicit estimates of CBD, often used in modeling wildfire behavior (e.g., FARSITE ;Finney 1998), are provided by the Landscape Fire and Resource Management Planning Tools (LANDFIRE) program (Reeves et al. 2009;Rollins 2009;www.landfire.gov). These CBD estimates are national in scope, calibrated using a combination of optical Landsat sensor data, US Forest Service Forest Inventory and Analysis (FIA) ground data, and the LANDFIRE reference database (Rollins 2009), which are then extrapolated to the landscape using sophisticated statistical and environmental modeling approaches (Keane et al. 2006). In this calibration process, certain assumptions regarding both the spatial and biophysical accuracy of FIA field plot data as well as the allometry used to determine CBD were adopted (Keane et al. 2006). With respect to FIA, recent work in Minnesota, USA, provides evidence to question some of the common assumptions regarding the combined use of Landsat and FIA for mapping biophysical forest structure information, especially spatial accuracy and suboptimal subplot design with regard to 30-m satellite sensor data (Thapa et al. 2020). Regardless of the sources of mapping error, it is understood that the absolute accuracy of the LANDFIRE CBD data product varies regionally (Keane et al. 2006;Scott 2008;Krasnow et al. 2009) and should not replace locally derived fuel products if they exist at finer spatial scales (Keane et al. 2006), which is the pretext for this research.
While CBD is a salient variable for modeling crown fire behavior, it is also one of the most difficult parameters to measure accurately in the field (Alexander 1988;Keane et al. 2005). This is because the vertical distribution of biomass varies by conifer species, crown position, and shade tolerance and from stand to stand (Brown 1978;Keane et al. 2002;Keane et al. 2006). Moreover, according to Keane et al. (2005), various vertical average measures of CBD likely underestimate effective CBD fuel conditions, noting that a few patches of fuel having substantially higher biomass values than the vertical canopy average will sustain fire spread. Therefore, they included the maximum CBD among all 1-m vertical layers as an important fuel variable to model, which we follow in this study. Equally vexing is the fact that key forest fuel properties (e.g., surface fuel, CBD, and canopy base height) are practically impossible to detect or distinguish via optical satellite sensors (e.g., Landsat-8 and Sentinel-2). Logical options suggested for remote detection of vertical canopy structures include lidar and radar (Keane et al. 2001;Keane et al. 2006).
Efficient methods for the calculation of CBD from indirect measurement of canopy gap fraction (CGF) and leaf area index (LAI) have been tested under Western forest conditions. Keane et al. (2005) described and tested six indirect ground-based methods to measure CGF (AccuPAR ceptometer, CID digital plant canopy imager, hemispherical photography, LI-COR LAI-2000, point sampling, and spherical densiometer). Of these methods, the LAI-2000 instrument (LI-COR Biosciences Inc., Lincoln, NE, USA) measures transmittance of diffuse sky radiation through a vegetation canopy via separate optical sensors oriented at five zenith angles (7, 23, 38, 58, and 68°). With this instrument, CGF and LAI are modeled in three dimensions using radiation attenuation rates at the different zenith angles (Norman and Welles 1983;Perry et al. 1988). In the Keane et al. (2005) study, iterative, incremental, destructive sampling of the coniferous forest canopy (25% basal area removal per iteration) was performed to produce field estimates of CBD across a gradient of crown biomass conditions. Concurrently, they collected six intervening measurements of CGF (listed above) before and after each destructive removal of the forest basal area. Results were used to build a suite of regression models to relate CGF estimates to their direct field measurements of burnable CBD. One of their models (CID digital plant canopy imager; R 2 = 0.94, root mean squared error [RMSE] = 0.01 kg m −3 ) required the use of a tree height parameter, while the remaining models did not. Of the remaining methods tested, the summed weighted average of LAI-2000 sensor data using the top three zenith angles (7, 23, and 38°) showed the strongest relationship to the maximum burnable CBD (foliage and small branches <3 mm diameter) in the vertical profile for Western coniferous forests (R 2 = 0.70, RMSE = 0.03 kg m −3 ). Equation 1 shows the Keane et al. (2005) weighted transformation of CGF for estimating maximum burnable CBD (CBD, kg m −3 ) across LAI-2000 zenith angles Ɵ i : Field estimates of CBD may also be modeled using specific desktop software, such as FuelCalc (Reinhardt et al. 2006). Current versions of FuelCalc (e.g., v1.6, released February 2019) use lists of inventoried tree species, with the associated bole diameter at 1.37 m height and canopy position information, to model available canopy fuel (i.e., foliage plus half the fine branch material ≤6.3 mm diameter) using Brown's (1978) allometric biomass equations. In instances for which species-specific allometry for fine, burnable branch material is not available, less specific allometry is substituted. The resulting estimates of crown biomass are then scaled to fit Brown's (1978) equations (Reinhardt et al. 2006). In other cases, where no published allometry exists for a tree species, allometry for similar, substitute species is used. In such cases, canopy fuel estimates are adjusted using specific factors to account for crown position (Gray and Reinhardt 2003). It is important to note that the FuelCalc routine excludes all trees less than 1.83 m (6 ft.) in height from CBD calculations, as biomass below this threshold is not canopy fuel, but rather surface fuel (Saatchi et al. 2007). As such, a substantial quantity of burnable, coniferous fuel in the Superior National Forest (SNF) in northern Minnesota, USA (Fig. 1), straddles this categorical boundary, especially balsam fir (Abies balsamea [L.] Mill.) (Frelich and Reich 1995).
In 2015, the US Forest Service funded a pilot study (see Engelstad et al. 2019) to investigate the potential improvement in methodologies for local modeling and mapping of forest structural parameters critical for modeling wildfire behavior within the Gunflint Ranger District of the SNF in northern Minnesota, USA (Fig. 1). Study-specific field data were combined with lowdensity lidar data (0.44 points m −2 ) to model stand age, crown fuel height, live crown base height, and CBD. For CBD, field estimates were derived using three methods: US Forest Service plot-based line-intercept method (USDI National Park Service 2003), allometric modeling of component fuel biomass from biophysical tree data according to Perala and Alban (1994), and transformation of indirect measures of CGF according to Keane et al. (2005). Of the fuel variables modeled using ran-domForest (Breiman 2001) k-nearest neighbor imputation, CBD models were the weakest (adjusted coefficient of determination [Adj R 2 ] range = 0.32 to 0.48; Engelstad et al. 2019) and did not exceed the accuracy of CBD estimates provided by the LANDFIRE program (R 2 = 0.58) for the Lake States region, Zone 41 (LANDFIRE 2011). The low reported accuracy among modeled estimates of CBD is discouraging, especially since low-density lidar has been used successfully (R 2 > 0.83) for these purposes in other regions (Andersen et al. 2005;Erdody and Moskal 2010). Thus, a need exists for alternative methodologies for mapping CBD within the SNF.
For passive satellite sensors (e.g., Landsat), the approximate 2-m boundary above ground level (AGL) separating canopy fuel from coniferous surface fuel is problematic for detection efforts, as such sensors are ill equipped to directly distinguish differences in fuel height (Keane et al. 2006;Saatchi et al. 2007;Jakubowksi et al. 2013). Data from active sensors (i.e., lidar and radar) represent a potential remedy for such limitations (see Keane et al. 2001;Keane et al. 2005). However, considering the marginal lidar-based CBD mapping results reported by Engelstad et al. (2019), coupled with decadelong revisit times proposed for higher-density lidar coverage (MNGAC [Minnesota Geospatial Advisory Council] 2020), an examination of the capabilities of synthetic aperture radar (SAR) satellite sensor data for quantifying CBD within this region is warranted due to proven canopy penetration characteristics (Huang et al. 2018). Encouraging CBD mapping results have been achieved (R 2 = 0.85, RMSE = 0.67 kg m −3 ) using highresolution (i.e., 5 to 10 m, L-band and P-band), polarimetric AIRSAR sensor data within Yellowstone National Park, USA (Saatchi et al. 2007). While the authors used allometric equations (Brown 1978;Van Hooser 1983) to derive field estimates of CBD from inventoried tree data, it is was unclear whether their estimates were equivalent to burnable CBD defined by Keane et al. (2005). In any event, it has yet to be determined whether data from space-based SAR sensors (e.g., Sentinel-1 and Palsar-1) can afford similar CBD modeling results under Minnesota forest conditions. Moreover, it is important to determine whether the presence of coniferous surface fuel may influence CBD characterization, especially given the factors that affect canopy penetration depth (Ustin 2004;Schlund et al. 2019).
Hence, we combined study-specific biophysical field data with both SAR and optical satellite sensor data to calibrate models for estimation and mapping of CBD to support fire behavior modeling efforts in the SNF. Field estimation of CBD response variables included both (1) FuelCalc output and (2) weighted transformation of indirect optical field measurements of canopy gap fraction (CGF) according to Keane et al. (2005) using a LAI-2200C instrument (LI-COR Biosciences Inc.). The fusion of optical and SAR satellite sensor data has been shown to augment remote characterization of forest structure in this region (Wolter and Townsend 2011). As such, we hypothesized that inclusion of SAR backscatter amplitude imagery with optical satellite sensor data would enable increased characterization of forest CBD over that possible via low-density lidar in this region (Engelstad et al. 2019). We also investigated the potential effects of coniferous surface fuel density on our ability to distinguish and model CBD using SAR satellite sensor data, given the greater canopy penetration afforded by radar. Finally, we compared and discussed the results of CBD model calibration strategies and provided recommendations for their operational use.

Study area
The 356-km 2 pilot study area falls within the Gunflint Ranger District of the Superior National Forest (SNF), Minnesota, USA, and is a mix of both managed and wilderness forest (Fig. 1). Forest cover is diverse (e.g., five conifer genera and seven broadleaf genera) and is considered transitional between the sub-boreal Great Lakes-St. Lawrence forests and boreal forest (Heinselman 1973;Baker 1989). Non-wilderness forest areas are intensively managed for wood fiber, which has resulted in a dominance of quaking aspen (Populus tremuloides Michx.), paper birch (Betula papyrifera Marsh), white spruce (Picea glauca [Moench] Voss), and balsam fir forest associations (Frelich and Reich 1995;Wolter and White 2002). Wilderness areas have an extensive fire history that supports vast stands of pioneer forest dominated by jack pine (Pinus banksiana Lamb.), as well as remnants of old-growth white pine and red pine (Pinus strobus L. and P. resinosa Ait., respectively) (Heinselman Fig. 1 Location of Minnesota, USA, pilot study area (black box) proximal to the Superior National Forest (SNF) and the Boundary Waters Canoe Area Wilderness (BWCAW) and Lake Superior. Here, we investigated the potential of mapping coniferous forest fuel density using satellite sensor data and ground data collected between 2015 and 2016 1973; Frelich and Reich 1995). However, early twentieth century fire suppression policies led to an increase in the dominance of shade-tolerant, fire-sensitive balsam fir on this landscape (Frelich and Reich 1995;Corace et al. 2012). Several wilderness and non-wilderness areas within the northwestern portion of the study perimeter experienced the effects of a severe downburst wind event in 1999 that caused substantial wind throw damage (Rich et al. 2007). Much of the downed course woody material from this wind event persists today amidst the regenerating forest. Throughout much of the study area, balsam fir exists largely in the understory below dominant and co-dominant canopy associates (Frelich and Reich 1995;Wolter and Townsend 2011). High flammability and understory canopy position make balsam fir an effective ladder fuel for crown fire propagation (Abbas et al. 2011). Other conifer species in this study area include eastern larch (Larix laricina [Du Roi] K. Koch), northern white cedar (Thuja occidentalis L.), and black spruce (Picea mariana Mill.). The coniferous forest species within this landscape are the focus of this research, which we discuss in the methods below.

Field sampling
Personnel from the Gunflint Ranger District of the SNF generated and field-marked 110 random plot locations throughout the pilot study area (Fig. 1) as candidates for field sampling. We visited plots (n = 61) that were within 1 km of a road (Fig. 2) in 2015 and 2016 to collect biophysical tree data and conifer CGF data, each needed to calibrate models for estimating burnable fuel biomass (primarily standing conifer foliage and twigs 0 to 3 mm in diameter) using two different approaches, described below. Approximately 20% of these 61 randomly marked plot center locations fell in large canopy openings or on the boundary of such an opening. In such cases, we shifted plot centers toward adjacent intact forest by~42 m (i.e., 30-m pixel diagonal dimension) to accommodate potential pixel registration errors and the spatial breadth of our field sampling design (Fig. 3). All plot locations were recorded using a dual-frequency Trimble Geo7X (H-Star) GPS receiver (Trimble Inc., Sunnyvale, CA, USA) and later differentially corrected (carrier-phase, mean 2dRMS = 0.48 m) using the Grand Marais, Minnesota, CORS base station (https://www.ngs.noaa.gov/ CORS_Map/; site ID: GDMA).
The presence of broadleaf tree foliage during the growing season is problematic for forest fuel mapping using satellite sensor data and equally problematic for estimating conifer CGF. Broadleaf foliage partially conceals coniferous fuels from above while also causing erroneously high estimates of CGF from below canopy. In both cases, ground-to-space fuel density calibrations would be confounded. Therefore, we collected field CGF data during dormant, leaf-off conditions for broadleaved tree species (late October 2015 and early May 2016) to to investigate the potential of mapping coniferous forest fuel density using satellite sensor data and ground data. Plot locations and roads displayed on pre-leaf flush 7 May 2017 Sentinel-2 near-infrared imagery, in which lakes are black, conifer forest is dark gray, leaf-off broadleaf forest is bright gray, and mixed conifer-broadleaf forest is intermediate gray maximize the visibility of conifer foliage and branch structures from below the forest canopy. Similarly, dormant seasons are optimal for space-based visibility of these same coniferous fuel structures (Wolter et al. 2008). At each plot, we collected two types of ground data: CGF and tree species biophysical data (i.e., bole diameter at 1.37 m height, or 10 cm if height ≤ 1.37 m; tree height; height to first live branch; and canopy diameter). On the rare occasion that trees were exactly 10 cm tall, we used the base diameter instead. We measured field CGF using a dual-sensor LI-COR LAI-2200C instrument (LI-COR Biosciences Inc.) during dawn and dusk hours to avoid direct sunlight. We placed one of the two LAI-2200C sensor wands at a fixed position in an open area (less than 100 m from a given forest plot) with an unobstructed view of sky conditions to collect measurements every 5 s. Concurrently, the other wand was used below the forest canopy to collect CGF measurements at each of nine grid locations per plot (3 × 3 at 5-m spacing), where the center grid location was the plot center (Fig. 3). Hence, all below-canopy measurements were time-synchronized (within 2.5 s) to the open sky sensor. The grid was oriented north and south on flat terrain, or up and down slope, according to Keane et al. (2005). We recorded CGF for all nine grid locations at both ground level and at 2 m above ground level (AGL). Ground-level CGF measurements capture all fuel (i.e., surface and canopy fuels) within the LAI-2200C sensor's upward field of view, while 2-m AGL measurements capture just canopy fuel. In each instance, we used the top three zenith angles of the LAI-2200C sensor: 7°(1), 23°(2), and 38°(3). Time-synchronized CGF digital readings from the open area and nine grid points were stored onboard their respective wands and later downloaded to a computer for post-processing to both eliminate spurious readings (transmittance >1 for each zenith angle) and to correct sub-canopy readings for varying open sky conditions. Resulting sets of skycorrected measurements were averaged by plot at each separate level (0 m and 2 m AGL) and zenith angle range (0 to 7°[1], 0 to 23° [12], and 0 to 38°[123]) to provide more robust spatial estimates of CGF (LI-COR 2011).
We collected biophysical tree data for all trees (conifer and broadleaf) greater than 10 cm tall (4584 trees across 61 plots) within a single circular subplot concentric with the plot center (Fig. 3). The radius of each circular subplot varied by plot according to the product of the sine of the second LAI-2200C sensor zenith angle (23°) and the maximum tree height identified for each plot (Fig.  3). The resulting average ground plot circular area was 102 m 2 with a standard deviation of 56.8 m 2 . We measured tree heights using a Trimble Laser Ace 1000 rangefinder (Trimble Inc.) to within circa ±20 cm precision (Jamali et al. 2014). We collected crown spread measurements for trees taller than 5 m using a densitometer instrument, according to Wolter et al. (2009), and recorded bole diameters by species to estimate basal area (m 2 ha −1 ). Conifer field data were used in two ways: (1) as FuelCalc inputs for determination of plot-wise maximum burnable CBD and (2) as means to estimate total bole, branch, and foliar biomass according to Perala and Alban (1994) to explore relationships between overstory canopy composition (i.e., broadleaf, mixed coniferbroadleaf, and percent canopy cover) and coniferous surface and lower canopy biomass (≤3 m AGL). For the latter analyses, tree-wise biomass values were stratified using respective live canopy cone volumes to distribute biomass vertically into 1-m increments (see Engelstad et al. 2019), which we summed by increment and normalized by ground plot area.

Field estimates of coniferous fuel density
We converted indirect LAI-2200C CGF data from each field plot to maximum, burnable fuel density (kg m −3 ) using the Keane et al. (2005) transformation equation (Eq. 1). Fuel density estimates included CBD (fuel ≥ 2 m Ground plot sampling design (2015 to 2016) to quantify coniferous forest fuel density for calibration with satellite sensor to enable mapping across the Superior National Forest, Minnesota, USA. Nine-point grid (black dots) at 5-m spacing are the canopy gap fraction (CGF) sample locations used to estimate coniferous fuel density. The central gray-shaded circle is a forest height-dependent fixed-radius plot (average area = 102 m 2 ; standard deviation = 56.8 m 2 ) for measurement of tree biophysical dimensions (species, bole diameter, crown spread, and height to first live branch). Field data are used to calculate both canopy bulk density (kg m −3 ) via the FuelCalc routine and to estimate understory (≤3 m above ground) coniferous biomass (kg m −2 of ground area) AGL) and total fuel density (TFD; fuel ≥ 0 m AGL), which is primarily the sum of live, coniferous surface fuel density (SFD; 0 to 2 m AGL) plus CBD. We generated estimates of TFD and CBD across three LAI-2200C zenith angle range combinations: TFD 1 , TFD 12 , TFD 123 , CBD 1 , CBD 12 , and CBD 123 . Respective fuel density differences, TFD minus CBD, served as indirect field estimates of burnable SFD. It is important to note that transformation of CGF to fuel density (kg m −3 ), whether CBD, TFD, or SDF, includes all fuels visible to the LAI-2200C sensor, which consisted primarily of live coniferous fuel structures.
Field estimates of CBD were also determined allometrically using plot-wise biophysical tree data and the FuelCalc routine (CBD FC ; kg m −3 ). As mentioned above, in the absence of species-specific allometry, FuelCalc makes use of allometry from structurally similar forest species to determine a plot's CBD profile (Lutes 2020), from which we extracted the maximum CBD. In northern Minnesota, most of the conifer species found in the pilot study area lacked the necessary allometry to calculate CBD. Hence, FuelCalc substituted the following species: Abies lasiocarpa (Hook.) Nutt. for Abies balsamea, Pinus contorta (Dougl.) for Pinus banksiana, Pinus ponderosa (Dougl.) for Pinus resinosa, Picea engelmannii (Parry) for Picea glauca, and Pseudotsuga menziesii (Mirb.) Franco for Picea mariana.

Satellite sensor data
Optical imagery. We acquired images from three optical satellite sensors (Table 1) for this study: five Landsat-8 Operational Land Imager (OLI) scenes, one SPOT-5 multispectral (XS) image, and one Sentinel-2 Multispectral Instrument (MSI) image. All OLI images were from the worldwide reference system (WRS-2) path 26 row 27 footprint and were downloaded from the US Geological Survey's Earth Explorer portal (www.earthexplorer.org). Images came processed to surface reflectance, orthorectified, and geo-corrected as part of the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS; Vermote et al. 1997;Masek et al. 2006). A SPOT-5 winter image (with 53-cm snow ground cover) from 13 March 2013 was downloaded from the Earth Explorer portal (https://earthexplorer.usgs.gov/) under a temporary data purchase agreement between the US Geological Survey and SPOT Image Corporation (https:// www.spotimage.com/). Finally, a MSI image from 7 May 2017 (pre-leaf flush) was downloaded from the Earth Explorer portal (https://earthexplorer.usgs.gov/) as an orthometrically and radiometrically corrected top of atmosphere reflectance product. The relative humidity recorded by the Seagull remote automated weather station (48.120536°, −90.838725°) at the time of MSI image acquisition (1030 Central Daylight Time) was~40%. Hence, further corrections to surface reflectance were not necessary, as our goals did not include date-specific change detection (see Song et al. 2001).
Landsat-8 OLI images for CBD, TFD, and SFD work included three leaf-off winter dates with snow ground cover (4 March 2014, 121-cm snow depth; 1 December 2015, 42-cm snow depth; and 19 February 2015, 53-cm snow depth), for which the Minnesota State Climatology Office provided date-specific estimates of snow depth (https://climateapps.dnr.state.mn.us/doc/snowmap.htm). Table 1 Synthetic aperture radar (SAR) and optical satellite sensor characteristics and image dates used in this Minnesota, USA, study to quantify coniferous forest fuels using field data collected in 2015 and 2016. Images used to classify overstory canopy composition are in boldface. Sensors include the Phased Array L-band Synthetic Aperture Radar (PALSAR-1) aboard the Japan Aerospace Exploration Agency's Advanced Land Observing Satellite (ALOS); the European Space Agency's Sentinel-1 C-band SAR sensor (interferometric wide swath mode; IW) and their Sentinel-2 MultiSpectral Instrument (MSI) optical sensor; the Space Agency of France's fifth Satellite Pour l'Observation de la Terre (SPOT-5) multispectral (XS) optical sensor; and the US Geological Survey's optical Operational Land Imager (OLI) sensor aboard the Landsat-8 satellite. Polarization for SAR sensors includes horizontal send and receive (HH), vertical send and receive (VV), and two cross-polarizations (HV and VH). Dashes (-) denote multiple images from the same sensor No precipitation events were recorded for at least 2 weeks prior to each winter image date. Two leaf-on Landsat-8 images (8 June and 14 October 2014) were acquired for the determination of overstory canopy composition. The 14 October imagery coincided with peak foliar senescence (autumn leaf coloration) of quaking aspen and paper birch, which served to augment multispectral contrast with other forest species in the overstory (Wolter et al. 1995). The 8 June OLI image served as a growing-season contrast to October OLI sensor data. Optical remote sensing data from dormant seasons are optimal for observing conifer canopies otherwise obscured by broadleaf tree foliage (Sayn-Wittgenstein 1961; Wolter et al. 2008). Moreover, winter optical imagery capturing snow ground cover conditions is ideal because potential spectral confusion between canopy and forest floor is minimized (Wolter et al. 2012).
In addition to optical sensor bands from each date, we calculated three indices for use as predictors of burnable fuel density: normalized difference vegetation index (NDVI; Rouse et al. 1974), shortwave infrared to nearinfrared ratio (SWIR:NIR; Vogelmann and Rock 1988), and SWIR to visible ratio (SVR; Wolter et al. 2008) ( Table 2). We also included ten mapped estimates of forest basal area (total, conifer, broadleaf, and seven conifer species; Wolter and Townsend 2011) at 30-m spatial resolution to the pool of optical predictors. Thus, there were 61 initial optical sensor predictors and ten structure estimates as predictors available for calibrations with ground data.
Synthetic aperture radar imagery. We acquired SAR backscatter amplitude imagery (Table 1) from two satellite sensors (Sentinel-1 and Palsar-1) to complement optical sensor data for calibrating fuel density models (TFD, SFD, and CBD). We downloaded Sentinel-1 dualpolarity (vertical send, vertical receive [VV]; vertical send, horizontal receive [VH]) C-band SAR interferometric wide (IW) imagery from 16 May 2016 from the European Space Agency's Copernicus Open Access Hub website (https://scihub.copernicus.eu/dhus/#/home) as ground range detected images georeferenced and resampled to a common pixel spacing (10 m) in both range and azimuth. These SAR data captured pre-leaf flush conditions for deciduous tree species in the study area. Leaf-off Palsar-1 dual-polarization (HH, HV) Lband data from 6 November 2010 (no snow cover) and 22 December 2010 (38-cm snow cover) were downloaded from the Alaska Satellite Facility's NASA distributed active archive center as high-resolution (12.5 m), terrain-corrected backscatter amplitude imagery (https:// www.asf.alaska.edu/about/asf-daac/).
In addition to the original six SAR bands (Table 1), we calculated three ratios to highlight potential differences in canopy volume structures (Joshi et al. 2017). Ratios included May Sentinel-1 send-receive polarization ratio (VH:VV) and two Sentinel-1 to Palsar-1 cross-frequency, cross-polarization ratios (VH C :HV L and VV C :HH L , referred hereafter without frequency subscripts C and L). We created 24 additional SAR predictors from a subset of the original bands and ratio images described above using 3 × 3, 5 × 5, and 7 × 7 low-pass filters to spatially dampen high-frequency, signal-dependent noise (Sader 1987;Rauste 2005). Thus, there were 31 SAR variables available for fuel density calibration procedures, along with a subset of 25 optical sensor predictors from leafoff seasons (Table 1). Table 2 Landsat-8 Operational Land Imager (OLI), SPOT-5 multispectral (XS), and Sentinel-2 Multispectral Instrument (MSI) bands; wavelength intervals (λ); and indices used as predictors of forest structural attributes in this Minnesota, USA, study to quantify coniferous forest fuel using field data collected in 2015 and 2016. Indices include normalized difference vegetation index (NDVI; Rouse et al. 1974), shortwave infrared to near-infrared ratio (SWIR:NIR; Vogelmann and Rock 1988), and SWIR to visible ratio (SVR; Wolter et al. 2008). Equations for NDVI, SWIR:NIR, and SVR are shown using sensor-specific band numbers defined below each sensor column. Dashes (-) denote where SPOT-5 lacks an equivalent predictor band Prior to ground-to-satellite calibration analyses and mapping, we identified and masked forest change areas due to harvest (i.e., clear-cutting, November 2010 to May 2017) from all image data.
Ground-to-satellite model calibration Partial least squares regression. We performed model calibrations between ground data and satellite sensor data using iterative exclusion partial least squares regression (xPLS; Wolter et al. 2008, Wolter et al. 2012). The xPLS algorithm was originally adapted for use with the PLS regression routine in SAS (SAS Institute Inc., Cary, NC, USA; Wolter et al. 2008) and later modified to run in Matlab v. R2010a (Mathworks Inc., Natick, MA, USA; Wolter et al. 2012;Singh et al. 2013). The xPLS model-building routine iteratively withholds individual predictor variables and performs intermediate cross-validations (Gong 1986) to determine the optimal number of latent components to use. This is important because, as latent dimensionality increases, lower order latent components-often describing random measurement error-are more likely to retain collinearity issues (see Wolter et al. 2008 andGeladi andKowalski 1986 for a detailed discussion of PLS regression and latent variable structure). The predictor variable that, when withheld, results in the greatest improvement in RMSE is permanently removed from further analysis. The model-building process proceeds iteratively until the search for improvements in model RMSE is exhausted, at which time prediction accuracy is assessed using the results of cross-validation, including predicted residual error sum of squared (PRESS) statistic (Geladi and Kowalski 1986) and Adj R 2 .
Fuel density model calibration. Using the xPLS routine, respective field estimates of TFD and CBD were each analyzed independently against the candidate group of 77 optical and SAR image predictors to calibrate parsimonious models (fewest predictors) for estimating fuel density metrics (i.e., TFD 1 , TFD 12 , TFD 123 , CBD 1 , CBD 12 , CBD 123 , and CBD FC ). We also used the three respective differences between TFD and CBD to build models for surface fuel density (SFD 1 , SFD 12 , and SFD 123 ). Given the procedural differences in deriving fuel density estimates (i.e., CGF and FuelCalc), we chose to perform separate xPLS regression sequences for each response variable. In each instance, final model calibrations were assessed based on results of leave-one-out cross-validation (PRESS), Akaike's information criterion (AIC; Akaike 1973), Bayesian information criterion (BIC; Schwartz 1978), and other standard metrics (Adj R 2 and RMSE). Relationships between respective field estimates of fuel density were evaluated using Pearson correlation plots (Müller and Büttner 1994). For canopy fuels, we used each final set of image predictors for mapping the respective estimates of CBD across the SNF study area.

Understory conifer biomass and overstory forest structure
We performed two sets of analyses to compare observed levels of live coniferous biomass present in the forest understory (e.g., boles, branches, and foliage ≤3 m AGL) to factors known to affect forest understory light environment: (1) percent canopy cover and (2) Anderson (1976) level-II forest types (Zavitkovski 1976;Messier et al. 1998;Légaré et al. 2002). While LAI-2200C CGF data measured from the narrowest zenith angle (7°) are suitable for measuring forest canopy cover (Paletto and Tosi 2009), our CGF data were insufficient in this respect as they were intentionally collected during dormant seasons, when all broadleaf species were leafless, to optimize visibility of coniferous fuels. Hence, we downloaded the 2016 continuous tree canopy cover dataset from the Multi Resolution Land Characteristic (MRLC) Consortium website (https://www.mrlc.gov/data/type/ tree-canopy). Understory biomass was calculated using our biophysical tree data and local allometry (Perala and Alban 1994) and scaled to ground plot area, as FuelCalc excludes a substantial portion of this biomass, especially below 2 m AGL. Next, we classified overstory forest cover into two Anderson level-II types (broadleaf and mixed broadleaf-conifer) using 8 June 2014 and 14 October 2014 Landsat-8 OLI sensor data and unsupervised classification techniques, for which 2015 leaf-on color infrared aerial photography (https://www.dnr.state.mn. us/maps/landview/) served as ground truth.
Conifers were not present at nine plots, reducing sample size for understory biomass comparisons from 61 to 52. Of these 52 plots, we identified two potentially influential outliers based on results of a multiple regression procedure performed on overstory composition and canopy cover against understory coniferous biomass. Each of these plots had Cook's distance (D) metric (Cook 1979) greater than four times the mean D, which is a conservative outlier threshold (see Kutner et al. 2005). In one instance, we noted that understory biomass was substantially higher than all other plots, while canopy cover from our MRLC source was 70%. This was suspect since we only measured six trees greater than 6 m tall on that plot. On the other plot, understory biomass was also high, while field data did not show an inordinate number of trees with either heights or height to first live branch below 3 m AGL. Hence, we removed these plots from further analysis. We tested the remaining field data for mean contrast among two levels of vertical biomass (0 to 2 m AGL and 0 to 3 m AGL) using Tukey's honestly significant difference test (HSD; Tukey 1953) across two separate overstory groups: forest type and canopy cover. Of the 50 field plots, there were 19 with broadleaf overstory and 31 with a mixed broadleaf-conifer overstory, while 13 plots had overstory canopy cover in the 47 to 66% range and 37 plots in the 67 to 78% range.

Fuel density calibrations with satellite sensor data
We performed model calibrations between groundbased fuel density estimates (TFD, SFD, and CBD) and satellite sensor data in separate xPLS regression sequences (Tables 3 and 4). Among these, all xPLS regression attempts between CGF-based estimates of SFD and satellite image predictors failed to resolve viable solutions and were aborted. For TFD and CBD, all ground-to-satellite calibrations yielded strong models based on statistical metrics except TFD 1 (Table 5). Of these, the four CBD models performed similarly (Adj R 2 = 0.95 to 0.98, RMSE = 0.05 to 0.17 kg m −3 , and PRESS < 0.47) and were superior to the TFD group of fuel models according to these metrics (Table 5, Fig. 4). The TFD group of models (TFD 1 , TFD 12 , and TFD 123 ) derived from weighted groundlevel CGF data were less complex in that they used fewer of the initial 66 satellite-based predictors (11, 16, and 23 predictors used, respectively) with lower latent dimensionality (i.e., orthogonal components) with respect to the response variable (Table 5). However, all of the TFD models resulted in overprediction bias for low values of fuel density and the opposite trend among the upper range of fuel density compared to field estimates (Fig. 4).
Final reduced sets of image predictor variables retained from the initial set of 66 candidate predictors (31 SAR, 25 leaf-off optical, and 10 basal area estimates from Wolter and Townsend 2011) for the seven models varied from 11 for TFD 1 to 38 for CBD 123 (Tables 3 and  4). The two CGF-based sets of fuel models (TFD 1-123 and CBD 1-123 ) revealed some trends in predictor variable retention. Only one of the TFD models (123) retained any of the OLI predictor variables (green band from 19 February 2015; OLI 3), while the CBD models collectively retained 13 of the OLI predictors. On average, each set of three CGF-based TFD and CBD models retained similar proportions of SAR predictors (48.8 and 43.7%, respectively; Table 4). However, TFD models, on average, used a greater proportion of cross-frequency SAR ratio predictors (18.0% versus 5.8%; Table 4), while the three CBD models retained a higher proportion of Lband predictors on average (15.1% versus 4.5%). Overall, CGF-based CBD models collectively retained a more uniform mix of both SAR and optical variables (Table  3). No clear trends in predictor variable retention by zenith angle range emerged among the three LAI-2200C sensor zenith angle ranges (0 to 7°, 0 to 23°, and 0 to 38°) used for TFD and CBD model development.
However, the TFD 1 model was substantially different from the CBD 1 model in this respect, retaining 0 versus 7 C-band and 1 versus 5 L-band SAR predictors in favor of cross-frequency SAR ratios. The CBD FC model showed variable retention trends similar to that of the CGF-based CBD models (Tables 3 and 4).
Last, the TFD and CBD models retained different numbers of the basal area estimate variables from the Wolter and Townsend (2011) study (Table 3). The TFD models retained two to four of these basal area estimates, while CGF-based CBD models retained five to seven and CBD FC retained four. All seven of these fuel density models retained the white spruce basal area estimate predictor and, all four CBD models retained the red pine basal area estimate predictor. All but two of the fuel density models (TFD 12 and CBD 1 ) retained the eastern white pine basal area predictor. There were no clear predictor retention trends among fuel density models by CGF measurement elevation or LAI-2200C zenith angle range ( Table 3).

Comparison of fuel density estimates and distributions
To keep modeled fuel density results in perspective, the reader should note that plot-wise field estimates of CBD derived via FuelCalc (CBD FC ) were treated here as quasiground truth (modeled, not actual), even though Western species substitutions were necessary (R. Keane, US Forest Service, Rocky Mountain Research Station, Missoula, MT, USA; personal communication). That said, correlations between our seven estimates of fuel density and the CBD estimates from LANDFIRE (CBD LF ) are shown in Fig. 5 along with their respective data distributions. As expected, estimates of CBD and TFD derived from the weighted transformation of CGF according to Keane et al. (2005) were more strongly correlated among respective levels of CGF measurement than between them (Figs. 4 and 5). All six CGF-based estimates of fuel density show substantially weaker but significant correlation (Fig. 5) with FuelCalc-derived estimates of CBD (predicted CBD FC ), for which TFD 12 and TFD 123 showed the strongest correlations (r = 0.51, P = 0.001). Correlations between CBD LF and CGF-based estimates of fuel density were either not significant or marginal (r = 0.25 to 0.31, P = 0.1 to 0.05). The relationship between the CBD LF and CBD FC estimates was not significant (Fig. 5).
Measures of TFD and CBD distribution shape (skewness, SK; and kurtosis, KU) across our 61 ground plots vary among the different field fuel estimates (Figs. 5 and 6). According to Bulmer (1979), the TFD 12 , TFD 123 , CBD 12 , CBD 123 , and CBD LF estimates of fuel density have distributions that may be considered symmetric around the mean (i.e., SK between −0.5 and +0.5), while TFD 1 and CBD 1 show moderate, positive skew (SK = 0.56 to 0.59). The Table 3 Final sets of image predictors (left side) and model coefficients for respective total fuel density (TFD; kg m −3 ) and canopy bulk density models (CBD; kg m −3 ) calibrated using dependent ground data collected in 2015 and 2016 in the Superior National Forest, Minnesota, USA. Model subscript notations 1, 12, and 123 reference the zenith angle range used for canopy gap fraction (CGF) measurements (0 to 7°, 0 to 23°, and 0 to 38°, respectively), while FC represents FuelCalc-derived field estimates of CBD. Dashes (-) indicate image predictors removed during model calibration. Image predictors include synthetic aperture radar (SAR) backscatter amplitudes and ratios, optical sensor bands and ratios, and previously derived estimates of forest basal area (BA) by species or group: Abies balsamea (ABBA), conifer (CON), broadleaf (BRD), Larix laricina (LALA), Picea glauca (PIGL), Picea mariana (PIMA), Pinus resinosa (PIRE), Pinus strobus (PIST), and total basal area (Total). Where necessary, image month is indicated. Combinations of SAR send and receive polarization are indicated with H (horizontal) and V (vertical). Spatially filtered (low-pass) versions of original resolution SAR variables have notations 3 × 3, 5 × 5, and 7 × 7. Optical sensor indices include the normalized difference vegetation index (NDVI), shortwave infrared to near-infrared ratio (SWIR:NIR), and SWIR to visible ratio (SVR). Sensors include Phased Array L-band Synthetic Aperture Radar (PALSAR-1); Sentinel-1 C-band SAR; Sentinel-2 optical MultiSpectral Instrument (MSI); Satellite Pour l'Observation de la Terre (SPOT-5) multispectral (XS) optical sensor; and Landsat-8's optical Operational Land Imager (OLI) Palsar-1 (P1)   Table 3 Final sets of image predictors (left side) and model coefficients for respective total fuel density (TFD; kg m −3 ) and canopy bulk density models (CBD; kg m −3 ) calibrated using dependent ground data collected in 2015 and 2016 in the Superior National Forest, Minnesota, USA. Model subscript notations 1, 12, and 123 reference the zenith angle range used for canopy gap fraction (CGF) measurements (0 to 7°, 0 to 23°, and 0 to 38°, respectively), while FC represents FuelCalc-derived field estimates of CBD. Dashes (-) indicate image predictors removed during model calibration. Image predictors include synthetic aperture radar (SAR) backscatter amplitudes and ratios, optical sensor bands and ratios, and previously derived estimates of forest basal area (BA) by species or group: Abies balsamea (ABBA), conifer (  estimates show an opposite trend (KU = 3.93), which suggests longer, relatively flatter tails compared to normal (leptokurtic). The CBD LF estimates extracted for our field plot locations show a distinctly uniform distribution trend (KU = 0.26), which is considered highly leptokurtic (Fig. 6).

Understory conifer biomass and overstory forest structure
Of the 3738 conifer trees measured in this study across 61 field plots, 70% (2617) were balsam fir; of these, 50.4% were <2 m tall, with an additional 15.6% in the 2 to 3 m range and 9.7% in the 3 to 4 m height range (i.e., 75.7% of balsam fir were ≤ 4 m tall). Mean understory coniferous biomass (needles, branches, and boles) within two separate vertical aggregations (0 to 2 m and 0 to 3 m AGL) did not show any significant differences (P = 0.291 to 0.978), according to Tukey's HSD test, when grouped by two broad overstory forest types (broadleaf and mixed broadleaf-conifer) or two overstory percent canopy cover ranges (47 to 66% and 67 to 78%).

Canopy gap fraction and canopy bulk density
The ability to model and map CBD via indirect CGF measurements hinges on the relationship between such measures, leaf area index (LAI), and canopy biomass (Brown 1978;Keane et al. 2005). However, CGF measured by the LAI-2200C and similar instruments tends to underestimate LAI, especially for conifers (Bolstad and Gower 1990;Gower and Norman 1991). Thus, branches and species-specific clumping of needles around shoots blocking the direct view of some fraction Table 4 The number (n) and relative proportion (% of total) of satellite predictors (optical, estimates of basal area [BA], and synthetic aperture radar [SAR]) retained by each canopy bulk density (CBD) and total fuel density (TFD) model during calibration in this Minnesota, USA, study (2015 to 2016) to quantify coniferous forest fuels. Subscript notations 1, 12, and 123 reference the zenith angle range of canopy gap fraction measurements (0 to 7°, 0 to 23°, and 0 to 38°, respectively) used to estimate field fuel density, while FC represents FuelCalc-derived estimates of CBD. A further breakdown of the SAR predictors into C-band, L-band, and C to L ratios is shown in italics. Dashes (-) indicate groups of SAR predictors not used during calibration for two fuel density models  Table 5 Satellite image to ground total fuel density (TFD; kg m −3 ) and canopy bulk density (CBD; kg m −3 ) model calibration results using field data collected in 2015 and 2016. Results used to evaluate the potential for mapping coniferous fuel density across the Superior National Forest, Minnesota, USA. Models for TFD and CBD are denoted T and C, respectively, while numbers 1, 12, and 123 reference the zenith angle range used for canopy gap fraction measurements: 0 to 7°, 0 to 23°, and 0 to 38°, respectively.  of needles in an instrument's optical field of view have been suggested (Gower and Norman 1991;White et al. 1997). However, Fassnacht et al. (1994) found no improvements in modeling conifer LAI when specific clumping bias factors (sensu Gower and Norman 1991;Stenberg et al. 1994) were employed. Interestingly, Gower and Norman (1991) speculated that the impact of needle-bearing branches on CGF and LAI is likely small due to foliar masking effects. Keane et al. (2005) went further and suggested that skyward projection of foliage may largely block the view of small and large supportive branch tissues. Their speculation arose from the fact that only small differences in predicted CBD values (i.e., maximum CBD measured within any 1-m increment up to tree top) were obtained when using LAI-2000 CGF measurements (collected at 2 m AGL) across different aggregations of canopy biomass data (e.g., foliage alone; foliage + small, burnable branches; and foliage plus all branches). Thus, the use of CGF as a basis for the determination of CBD, while not perfect, may be robust to species-specific differences among coniferous fuel structures.

Disparity among CBD estimates: CGF, FuelCalc, and LANDFIRE
We developed field estimates of CBD in the SNF via two methods: allometric modeling via the FuelCalc routine and transformation of CGF according to Keane et al. (2005), in which the authors meticulously calibrated CGF-to-CBD models based on the iterative destructive sampling of Western conifer stands to determine the amount of biomass removed and amount of biomass that remained. We extended the use of their methodologies and equations to model field estimates of CBD using the 0 to7°(1) and 0 to 23°(12) zenith angle ranges for comparison to their 0 to 38°(123) models. Of the CBD models, CBD 123 most closely approximated field estimates of CBD obtained using FuelCalc (CBD FC ;Figs. 4 and 5). However, the CBD 123 model produced estimates that far exceed values believed to be typical for this region (see Engelstad et al. 2019). In this study, the upper range among CBD estimates increases substantially with increasing LAI-2200C sensor zenith angle (Fig. 4). In this respect, the CBD 1 model produced a more realistic range of estimates, while the CBD 12 model was intermediate between the extremes (Fig. 4, Table 5).
Fire managers on the SNF suspect that actual CBD values do exceed that which LANDFIRE portrays due to the abundance of lower canopy fuels (e.g., 2 to 4 m AGL), but lack a definitive set of sampling data to reinforce such conjecture (P. Johnson, USDA Forest Service, Grand Marais, MN, USA, unpublished data). Fig. 4 Regression results between field estimates and satellite-derived estimates of total fuel density (TFD; kg m −3 ) and canopy bulk density (CBD; kg m −3 ) used for mapping fuel within the Superior National Forest, Minnesota, USA, following field data collection in 2015 and 2016. Models of TFD (coniferous fuels above ground level) and CBD (coniferous fuels 2 m above ground level) based on weighted canopy gap fraction are denoted T and C, respectively, while numbers 1, 12, and 123 refer to zenith angle ranges (0 to 7°, 0 to 23°, and 0 to 38°, respectively) of the LAI-2200C instrument. FuelCalc-based estimates of CBD are denoted FC. The dotted black line represents unity between field-measured and satellite-calibrated estimates of fuel density Moreover, according to our findings, modeling the distribution and abundance of such fuels in the SNF based on overstory canopy characteristics may not be a viable mapping option. Current versions of FuelCalc conservatively overestimate CBD by using the maximum CBD found within any 3.5-m layer of the canopy, a value that often exceeds twice the profile average Scott 2008). This is in response to both fire spread sensitivity to maximum CBD (Keane et al. 2005) and underlying issues with fire spread equations (van Estimates of TFD (all coniferous fuel above ground level) and CBD (coniferous fuel ≥2 m above ground) based on weighted canopy gap fraction are denoted T and C, respectively, while numbers 1, 12, and 123 refer to zenith angle ranges (0 to 7°, 0 to 23°, and 0 to 38°, respectively) of the LAI-2200C instrument. FuelCalc-based field estimates of CBD are denoted FC. Pearson correlation coefficients with corresponding P-value significance codes (*** for 0 to 0.001; ** for >0.001 to 0.01; * for >0.01 to 0.05; · for >0.5 to 0.1; and no symbol for >0.1 to 1) are shown above the diagonal, while the diagonal shows distributions of fuel density estimates. Bivariate scatter plots with fitted lines are below the diagonal. Linking field fuel estimates and satellite sensor data facilitates ground-to-satellite calibration of fuel models for fuel mapping purposes Wagner 1977; Rothermel 1991). In this study, however, average estimates for CBD FC were over nine times greater than average CBD LF across our study area (Table  5). Hence, we suspect that the actual CBD values for the SNF region lie somewhere between these estimates. Irrespective of such conjecture, it is curious that Pearson correlations between the two CBD estimates are neither significant nor marginally significant (Fig. 5).
Whether such offsets are linked to the substitution of functionally similar conifer species (Western for Eastern) within the FuelCalc routine or some other factor such as inappropriate use of canopy closure rather than canopy cover during the calculation of CBD LF driving values downward (Scott 2008) is difficult to determine. Early in our investigation, we questioned FuelCalc's use of Douglas-fir (Pseudotsuga menziesii [Mirb.] Franco) allometry over Engelmann spruce as a substitution for black spruce. However, even higher estimates of CBD resulted when we substituted Engelmann spruce allometry. Similar questions arose for other species having low abundance on this landscape. However, for conifer species that do dominate this landscape (Fig. 7), we had no defensible reason to question FuelCalc's substitutions. Nevertheless, while we assumed FuelCalc-derived estimates of CBD were quasi-ground truth in this research, we suspect that some Western species substitutions for Eastern species that lacked detailed allometry may represent a potential flaw in these analyses. Hence, detailed allometric equations that afford accurate, region-specific estimation of CBD for common Eastern conifer species remain a research need.

Effects of CGF measurement level on fuel density estimation
The collection of CGF data from 2 m AGL to estimate CBD (Keane et al. 2005), LAI (Fassnacht et al. 1994), and other forest canopy parameters (Paletto and Tosi 2009) is standard practice, while the collection of ground-level CGF for these purposes is not. This region of the SNF has a high abundance of understory coniferous fuel (mainly balsam fir) that straddles the 2-m AGL fuel boundary between surface and canopy fuels. As such, we used groundlevel CGF measurements to enable overall assessment of the sensitivity of satellite sensor data (especially SAR) to canopy fuel biomass (CBD) and coniferous surface fuel density (SFD). That said, the most dramatic differences among CGF-based fuel density models were between levels of CGF measurement (0 versus 2 m AGL) rather than between the LAI-2200C sensor's fields of view (0 to 7°[1], 0 to 23° [12], and 0 to 38°[123]). In regression Fig. 6 Skewness and kurtosis measures (both unitless) of frequency distribution shape for respective estimates of fuel density (labeled dots). Shape statistics for total fuel density (TFD) and canopy bulk density (CBD) based on weighted canopy gap fraction are denoted T and C, respectively, while numbers 1, 12, and 123 refer to zenith angle ranges (0 to 7°, 0 to 23°, and 0 to 38°, respectively) of the LAI-2200C instrument. Shape statistics for FuelCalc-based estimates of CBD are denoted as FC, while LANDFIRE CBD is denoted LF. Dashed lines indicate zones of relative normality for the respective shape metrics. Field data for these analyses were collected in 2015 and 2016 and combined with satellite sensor data to calibrate models for mapping fuel density across the Superior National Forest, Minnesota, USA analyses, field-based estimates and satellite-based estimates of fuel density for all TFD models deviated from unity, showing both additive and multiplicative shifts between dependent TFD variables and independent, satellite-based estimates of TFD (Fig. 4). In contrast, all satellite-based estimates of CBD followed unity with field estimates of CBD via CGF (Fig. 4).
For TFD models, two complementary explanations for deviations from unity with ground-based estimates arise. First, based largely on failure of all SFD model calibrations, we suspect that contributions of smaller, understory conifer saplings to estimates of TFD were not effectively detected by either optical or SAR satellite sensors, in spite of the relative retention difference among SAR predictors observed for TFB and CBD calibrations ( Table 4). As such, TFD model calibrations between ground and satellite sensor data remained possible because some portion of canopy fuel structures were visible to SAR sensors and, to a lesser degree, by optical sensors, which we posit was not the case for SFD. Second, with respect to ground-level CGF measurements, we suspect that TFD is biased by live, coniferous surface and lower canopy fuels (e.g., saplings) when present. If true, this may preclude field estimation of SFD (i.e., TFD minus CBD) via CGF metrics, which comports with prior research (Gower and Norman 1991;Fassnacht et al. 1994).
With respect to CBD, we suspect that burnable fuel in lower layers of the forest canopy above 2 m AGL (saplings and lower branches of mature trees) in this sub-boreal landscape represents a substantial fuel biomass gap that remains hidden from remote detection, space-based SAR sensors notwithstanding. While generally accepted as dogma for passive optical sensors such as Landsat OLI (Keane et al. 2005), we suspect that lower canopy fuel detection via the current suite of SAR sensors is also limited due to a combination of C-and L-band wavelength idiosyncrasies with respect to forest cover (Le Toan et al. 1992;Imhoff 1995;Sadeghi et al. 2018). For instance, Lband (23.6 cm λ) backscatter from relatively dense, mature forests of typical basal area (e.g., 20 to 30 m 2 ha −1 ) corresponds, in decreasing order, to the abundance of larger canopy components (boles and larger branches); groundtrunk double bounce; and, to a lesser degree, depending on canopy cover, the ground (Wang et al. 1995;Ustin 2004). Hence, smaller canopy components, such as overstory foliage and understory conifer saplings (e.g., balsam fir) remain largely undetected at L-band wavelengths (see Imhoff 1995). For larger trees, however, known relationships between bole dimension, basal area, and CBD (Duveneck and Patterson III 2007;Fernández-Alonso et al. 2013) may partially explain the retention of L-band predictors during CBD calibrations (Table 3). At C-band Fig. 7 Distribution of forest basal area (BA; m 2 ha −1 ) showing the relative abundance of conifer and broadleaf forest associates among 61 random field plots sampled in 2015 and 2016 within the Superior National Forest, Minnesota, USA. Shown are three forest BA levels (0 to 5, 5 to 10, and ≥10 m 2 ha −1 ) for eight conifer species and two generalized forest groups: conifer (CON) and broadleaf (BRD). Sampled conifer species include balsam fir (ABBA), white spruce (PIGL), black spruce (PIMA), tamarack (LALA), jack pine (PIBA), red pine (PIRE), eastern white pine (PIST), and northern white cedar (THOC) (5.54 cm λ), for the same forest scenario, the scattering phase center is primarily a function of smaller components in the upper forest canopy, including foliage and small branches (Wang et al. 1993;Pulliainen et al. 1999;Sexton et al. 2009). This limits canopy penetration to a fraction of overall canopy depth (Kellndorfer et al. 2004;Sadeghi et al. 2018). Therefore, conjecture regarding the relative significance of SAR predictor retention differences between TFD and CBD model calibrations (Tables 3 and  4) would be unwise.

Recommendations for the fire modeling community
Given the discussion above, we believe that CBD 1 and CBD 12 , derived via transformed CGF data according to Keane et al. (2005), represent the most promising estimates of burnable canopy fuel in the SNF. In spite of potential allometric differences in burnable biomass between Eastern conifers and their Western substitutes, the majority of CBD FC estimates lie within the combined range of these two indirect estimates of CBD (Fig. 4). Moreover, the SAR sensor data used in this study to scale field estimates of CBD to the SNF landscape do not appear biased to the presence of underlying coniferous surface fuels. Hence, we recommend testing these CGF-based estimates of CBD as well as CBD FC in the Upper Midwest within the FARSITE modeling framework to evaluate and understand potential fire behavior sensitivity to these different fuel inputs. What seems to be of high importance when modeling fire behavior is having CBD values high enough to propagate fire across a landscape (Scott 2008). Critical CBD thresholds that, when exceeded, support active crown fire behavior vary according to fine fuel moisture content and 10-m wind speed (Agee 1996;Cruz et al. 2005). For Western coniferous forests, this may range from 0.05 kg m −3 at 4% fuel moisture and 25 km h −1 wind speed to as high as 0.30 kg m −3 at the same fuel moisture and~2 km h −1 wind speed (Cruz et al. 2005). Hence, with respect to fire behavior modeling in the SNF, our estimates of CBD certainly seem high enough to propagate fire at almost any wind speed. However, the potential importance and impact of various differences in fuel distribution and density among our estimates of CBD across this landscape (Figs. 5 and 6), which may affect rates of fire spread (Agee 1996), remain untested.

Ladder fuels
Remote detection of coniferous ladder fuels (surface and lower canopy) in the Upper Midwest remains a challenge using either low-density lidar (Engelstad et al. 2019) or dual-polarity C-and L-band SAR satellite sensor data. The distribution and abundance of ladder fuels within the SNF landscape (Fig. 8) represent a substantial gap in knowledge that affects the efficiency of efforts focused on fuel reduction treatments to mitigate wildfire risk. Research on linkages between the abundance of understory biomass and combinations of overstory canopy type and percent cover (Messier et al. 1998;Légaré et al. 2002) begs the question of whether modeling the presence of ladder fuels may be possible in this region. We suspect that positive understory light-to-biomass relationships likely hold true regarding broadleaf biomass in this region. However, the results of our analyses do not support the possibility of modeling balsam fir ladder fuels via this mechanism due to its strong shade tolerance (Corace et al. 2012). Low-density lidar sensor data (i.e., average of 0.44 points m −2 ) have not yet facilitated sufficient characterization of such ladder fuels in this region (Engelstad et al. 2019) to enable sound management decisions. Resolution of this pure broadleaf overstory. One goal of this research was to determine the impact of coniferous surface and lower canopy fuels on satellite-based detection and mapping of canopy bulk density (CBD; kg m −3 ) in this region using indirect measurements of CBD as ground truth gap in fuel information will likely only be resolved once higher-density lidar data (e.g., 8 points m −2 ) are available on a more frequent basis-a project that is currently underway in Minnesota (MNGAC [Minnesota Geospatial Advisory Council] 2020).

Conclusions
In an attempt to capture coniferous fuel structures beyond optical detection limits, as suggested by some researchers (Keane et al. 2005;Keane et al. 2006), we combined dual-polarity C-and L-band synthetic aperture radar (SAR) with leaf-off optical satellite sensor data to achieve deeper detection depth into forest fuel structures. Satellite data were coupled with both FuelCalcbased estimates of CBD and indirect, CGF-based, field estimates of CBD, TFD (all burnable fuels above ground level), and SFD (i.e., TFD minus CBD) to calibrate fuel models to evaluate fuel structure mapping capabilities. Among the indirect methods, the slope of fuel density estimates (observed versus predicted) was highly sensitive to the elevation of CGF measurements (0 m versus 2 m AGL), while the magnitude of fuel density estimates was sensitive to the field of view of the LAI-2200C CGF sensor. We believe that the presence of understory balsam fir (e.g., ≤3 to 4 m tall) attenuated the ground-level view of the LAI-2200C instrument, causing substantial bias in predicting TFD. This partially contributed to the failure of all SFD model calibrations with satellite sensor data, in addition to known canopy penetration limits associated with C-band SAR (Kellndorfer et al. 2004;Sadeghi et al. 2018). At the same time, it is likely that lower strata biomass structures (boles and large branches) are invisible at L-band wavelengths (see Imhoff 1995). Thus, for larger trees, we suspect that the relationship between bole dimension, basal area, and CBD (Duveneck and Patterson III 2007;Fernández-Alonso et al. 2013) may be linked to the retention of L-band predictors during CBD model calibrations (Tables 3 and 4).
Models of CBD derived via CGF data did not suffer equally from suspected lower fuel bias issues, as all showed better unity between field-derived and satellitepredicted estimates of CBD. However, the maximum range of modeled CBD remained sensitive to the zenith angle ranges of the LAI-2200C instrument. While the 0 to 38°zenith angle range produced statistically superior CBD calibrations with satellite sensor data and were most similar to CBD FC estimates, both CBD 123 and CBD FC produced suspiciously high upper ranges of CBD. Canopy gap fraction measured at 0 to 7°and 0 to 23°zenith angle ranges each produced a more reasonable range of CBD values for this region. However, the importance of spatial variability among these CGFderived CBD estimates on actual fire behavior modeling remains unknown. With respect to FuelCalc, substitution of Western for Eastern conifer species' allometry in the calculation of CBD represents an unknown effect on these analyses. Hence, a need exists for the development of region-specific CBD allometry.
A substantial amount of coniferous biomass in the SNF region resides in the understory layers below the dominant overstory canopy. Much of this biomass straddles the 2-m AGL categorical boundary between surface fuels and canopy fuels, each of which is largely hidden from the view of optical and SAR sensors. Unlike other regions of the USA, the abundance of understory coniferous biomass, especially balsam fir, in the SNF appears to be insensitive to both overstory forest type (broadleaf versus broadleaf-conifer mix) or percent canopy cover. This contradicts current dogma (Zavitkovski 1976;Messier et al. 1998;Légaré et al. 2002) and is likely due to the strong shade tolerance of balsam fir.