Modeling Marbled Murrelet nesting habitat: a quantitative approach using airborne laser scanning data in British Columbia, Canada

,


INTRODUCTION
The marbled murrelet (Brachyramphus marmoratus, hereafter "murrelet") is an auk found along the Pacific coast of North America.Unlike most auks in the Alcidae family, murrelets do this forest type to a small fraction of its original extent (Burger et al. 2020, ECCC 2021).In British Columbia (BC), coastal forests are now dominated by early and mixed seral stage forests (Wulder et al. 2009) with less than 3% of remaining old growth composed of high productivity large-treed forest (Price et al. 2021).The loss and fragmentation of nesting habitat represents the principal threat to murrelet populations, although marine threats (such as oil spills and changing prey distributions due to climate change) also pose additional risks (ECCC 2021).An annual decrease of 2.4% in the BC murrelet population was estimated from 1996 to 2018 with the largest declines found in heavily harvested landscapes (Drever et al. 2021).Due to concerning population trends, murrelets are classed as a "Schedule 1 species" under the federal Species at Risk Act in Canada (COSEWIC 2012) and as "Threatened" in Washington, Oregon, and California, under the United States Endangered Species Act (USFWS 1992).
In BC, recovery efforts for murrelets aim to retain over 70% of nesting habitat compared to a 2002 baseline, and compliance with this target is tracked through habitat mapping (ECCC 2021).Knowing what constitutes suitable nesting habitat and how this habitat is currently distributed across the landscape is a crucial requirement for the appropriate placement of reserves, informing forest management, and tracking potentially important habitat areas.Consequently, there has been significant interest from scientists, environmental agencies, and forest managers in mapping suitable nesting areas.Three principal methods are currently employed to identify potential murrelet nesting habitat in BC: (1) a model classifying suitable habitat from vegetation resource inventory (VRI) polygon information (Mather et al. 2010, BC Ministry of Forests 2022), (2) air photo interpretation (API) in which direct assessments of habitat quality are interpreted from aerial photographs (Burger 2004, Burger et al. 2009), and (3) targeted low-level aerial surveys (LLAS) during which biologists in helicopters visually assess forests for the presence of potential nest platforms and other suitable microhabitat characteristics (Burger et al. 2009).These methods are comprehensively evaluated in Burger et al. (2018), with highintensity LLAS providing the most reliable method for habitat identification, correctly predicting 85% of nests within mapped suitable habitat.However, this method is impractical to conduct over the large areas needed for murrelet management unless applied at low intensity, which is less accurate, predicting only 48% of nests within mapped areas of suitable habitat (Burger et al. 2018).Vegetation resource inventory models capture relevant macrohabitat features over large areas but do not capture the key microhabitat features that determine nest sites.No method is solely reliable for habitat prediction, so multiple methods are often used together to inform management (Burger et al. 2018).Therefore, to advance habitat management and recovery efforts for murrelets, there is a need to accurately predict fine-scale habitat suitability across coastal forests using a defensible methodology that is not unduly resource intensive (COSEWIC 2012, ECCC 2021).
Habitat suitability modeling (HSM) encompasses a range of quantitative techniques that estimate the probability of finding a species over different environmental conditions (Guisan et al. 2017).Rare and elusive species, such as murrelets, need reliable HSM for conservation planning, but modeling is constrained by small samples of species data and unreliable absence data due to low detection rates.Presence-only modeling approaches that accommodate low sample sizes have been developed for rare species' situations, with discriminative ensemble of small models (ESM) and generative maximum entropy modeling (MaxEnt) recommended as appropriate HSM methods (Pearce andBoyce 2006, Breiner et al. 2015).Both ESM and MaxEnt perform well with < 100 samples and have been shown to produce accurate models of the distribution of rare species compared to envelope or regression-based approaches (Lomba et al. 2010, Breiner et al. 2015).Importantly, presence-only models do not calculate the real probability of finding a species in a location, rather they estimate the relative probability of occurrence, i.e., the highest predicted habitat value corresponds to the most probable place to find a species within the training area.
Murrelets are a viable candidate for HSM because their nesting niche is thought to be specific to a narrow range of forest conditions.Murrelets typically nest on large branches at least 15 cm wide using epiphytes or duff as nest substrate and overhead vegetative cover for protection from predators (Burger 2002).Murrelets are fast and direct flyers that use a stall-landing and a drop take-off for nest access, requiring suitable openings in the canopy and a sufficiently high nest to enable these maneuvers (Manley 1999, Burger 2004).Nest trees in BC are typically > 30 m tall, with no nest tree recorded < 20 m, and are mostly found within 30 km of the ocean, with some outliers found up to 80 km away (Burger and Waterhouse 2009).Multi-layered canopies with potential nest platforms and cover have consistently been shown to be favorable for murrelet nesting, signifying the importance of forest structure for nesting (Plissner et al. 2015, Hagar et al. 2018, Hamer et al. 2021).
Given that murrelet nests are dependent on forest and tree structure, airborne laser scanning (ALS) data, which is routinely acquired for forest inventories and terrain mapping, offer an appropriate quantitative data source to describe murrelet habitat (Burger et al. 2020).Airborne laser scanning is an application of light detection and ranging (lidar) technology, which is an active form of remote sensing, capable of accurately measuring the distance to objects by emitting laser pulses and measuring the time taken for photons to reflect off a surface and return to the sensor (Wulder et al. 2012).By amalgamating laser distance measurements, 3-D datasets called "point clouds" are produced that represent the structure of a scanned area (White et al. 2016, Torresan et al. 2017).Airborne laser scanning is commonly flown by piloted and remotely piloted aircraft, with the scanner positioned downward to capture the underlying environmental structure.Onboard motion sensors and global navigation satellite systems (GNSS) enable ALS scans to be accurately geo-located, providing great utility for spatial mapping.A range of measurements can be calculated from 3-D point cloud data that can be used to produce spatially explicit maps of vegetation and topography structure over large areas (Wulder et al. 2011).
Airborne laser scanning data are increasingly used for avian HSM, capturing unique information not measurable from other environmental data (Davies andAsner 2014, Bakx et al. 2019).Airborne laser scanning captures a range of structural features important to the ecophysiological needs of birds, enabling ecologically relevant HSM for species such as the Northern Spotted Owl (Strix occidentalis caurina; Ackers et al. 2015, Hagar et al. 2020).Hagar et al. (2014) demonstrated the feasibility of ALS informed HSM for murrelets, with canopy cover, tree height, structural complexity, and distance to the coast used as predictors in a binomial regression model that successfully discriminated murrelet occupancy among forest stands in Oregon.However, this model cannot be replicated in BC due to the lack of reliable absence data and latitudinal differences in forest composition and structure.In BC, a murrelet habitat classification method was developed by Clyde (2017) that successfully linked LLAS survey data with tree-level ALS descriptions.Clyde's method extrapolated qualitative LLAS habitat rankings to individual trees and showed that tree-level nesting suitability can vary greatly within areas of suitable forest mapped by LLAS.However, this method did not have nest data available for training or testing, preventing assessment of its reliability for management.
Airborne laser scanning data have shown great promise for assessing murrelet habitat.However, to be operationally useful in BC, a quantitative nesting habitat model is needed that is trained on regional forest structure and nest data.Therefore, we aim to assess: (1) if ALS data are suitable for the quantitative modeling of murrelet nesting habitat in BC, (2) which modeling approach and predictor set give the top-performing habitat model, (3) how predicted habitat is distributed across the landscape, and (4) how the top-performing model compares with existing LLAS maps of nesting habitat.We address these queries and then discuss if ALSinformed habitat models are suitable to support forest management decisions for murrelets in BC.

Study area
This study focuses on two important areas of murrelet nesting habitat in BC (Pastran et al. 2022) that contain overlapping ALS and historical nest data (Fig. 1).
Area A (Fig. 1) is adjacent to Desolation Sound (DS) in the Coastal Western Hemlock (CWH), Mountain Hemlock (MH), and Coastal Mountain Heather Alpine Biogeoclimatic Ecosystem Classification (BEC) zones (BC Forest Analysis and Inventory Branch 2021).Desolation Sound lies on the inside passage of the Salish Sea in the rain shadow of Vancouver Island and has a maximum elevation of 1530 m.Area B (Fig. 1) is adjacent to Clayoquot Sound (CS) and is in the CWH, MH, and BEC zones (BC Forest Analysis and Inventory Branch 2021).Clayoquot Sound is on the exposed west coast of Vancouver Island and holds forests that span from sea level to 1000 m.Both areas contain moist maritime and very wet maritime biogeoclimatic subzones, but CS contains additional wetter areas of the hyper-maritime subzone while DS contains areas of dry maritime; and very dry maritime subzones.The landscape conditions in DS and CS differ greatly due to the legacies of timber harvesting (Fig. 1).The forests of DS are highly fragmented with approximately 80% of the treed land base harvested prior to the collection of our study data, as opposed to only 15-25% in CS (Zharikov et al. 2007).and were acquired by the Province of BC in 2020 (Government of British Columbia 2022).The average point density for DS and  (1984-2020;Hermosilla et al. 2015).Previous harvesting occurred prior to 1930 in both areas but these data were not available.
Nest locations were identified during radio-tagging projects from 1998-2002 and are a unique and important dataset for murrelet research (Bradley et al. 2004, Cooke et al. 2022a, b).Nest locations were recorded by radio-tracking birds captured at sea and are considered a random, unbiased, sample (Zharikov et al. 2006).The standard deviation of point heights in a 1-ha area.
This measures the spread of vertical vegetation structure in a forest.More complex vegetation may provide more cover and nesting opportunities.See label 3 in Figure 3.A measure of access to feeding sites.Larger distances will require more energy to reach foraging areas.
dist_coast Murrelets were tracked to nest sites where a GNSS position was taken, which resulted in a positional uncertainty of roughly ±50 m for each nest, making pinpointing nests to specific trees impossible (Waterhouse et al. 2002).Murrelets are known to occasionally nest on cliffs in our study area (Manley 1999, Bradley andCooke 2001).Out of 85 nests with ALS coverage, 6 were identified as cliff nests using ALS data and 3 m resolution PlanetScope imagery (Planet Team 2017).Cliff nests were not used in our study, resulting in a final sample of 58 nests in DS and 21 nests in CS.
Pre-existing LLAS maps of murrelet nesting habitat were provided by the BC Ministry of Water, Land, and Resource Stewardship but were only available for DS.Nesting habitat suitability was ranked on a scale of 1 (highly-suitable) to 6 (no habitat potential) using low-intensity, medium-scale helicopter surveys conducted from 2003-2019 (Burger et al. 2009).A standard survey methodology was conducted for all LLAS, with forests younger than 140 years old (as recorded in the VRI) prescreened from survey areas because it is generally considered unsuitable for murrelets (Burger 2004).Areas in DS not represented by LLAS data were defined as unranked.Unranked areas were either excluded from the LLAS due to pre-screening or are areas that have yet to be surveyed.
Forest age data were accessed from the BC VRI dataset, which records the estimated age of the leading species in a stand (Government of BC 2021, BC Ministry of Forests 2022).The VRI is expected to contain substantial inaccuracies in forest age estimates, particularly for older stands of trees with no history of harvesting (Mah andNigh 2003, Hilbert andWiensczyk 2007).
The National Terrestrial Ecosystem Monitoring System (NTEMS) is a 30 m resolution landcover change detection product produced from Landsat imagery using spectral trend analysis techniques (Hermosilla et al. 2015).This system records the year of landcover change, the predicted cause of change, and the magnitude of the change across Canada.The NTEMS change type and change year products were used to identify areas of harvesting and other stand replacing disturbance occurring between 1998 and the time of ALS acquisition.

Environmental predictors
Four ecologically justified ALS predictors (Table 1) were tested in modeling.Forest height (FH) and forest vertical complexity (FVC) were chosen because larger forests that are more vertically variable in their structure were expected to be more favorable for nesting, providing more nest platforms and overhead cover (Fig, 2, labels 1 and 3; Burger et al. 2010, Hamer et al. 2021).Two custom ALS predictors were created to capture untested vegetation structural features thought to be appropriate for describing murrelet microhabitat: canopy surface area (CSA) and internal forest gaps (empty_gaps; Hagar et al. 2014, 2018, Burger et al. 2020).Canopies with emergent trees and a variable canopy (captured by CSA) alongside gaps within and under the canopy (captured by empty_gaps) were expected to be more favorable for nesting by providing canopy openings for nest access and flight maneuvers (Fig. 2, labels 2 and 4).
Non-structural variables thought to be informative for nesting were also tested as model predictors (Table 1).These were the climate moisture index (CMI; Wang et al. 2016) and distance to the coast (dist_coast).Wetter and more humid areas (represented by larger CMI values) were expected to contain more epiphytes, which are thought to increase the probability of nesting by providing nesting substrate (Van Rooyen et al. 2011, Hamer et al. 2021) and areas closer to the ocean were expected to be more favorable because of the lower energetic cost required to reach foraging sites (Hagar et al. 2014).
All predictors were created at a spatial resolution of 100 m (i.e., each pixel covered 1 ha of forest) to correspond with the positional error (± 50 m) in nest locations.Although this approach decreased the spatial precision of the predictions, it increased the chance a pixel overlapped the true location of a nest (Guisan et al. 2017).Additionally, to make the ALS predictors (produced from data collected between 2015-2019) more representative of forest conditions experienced by murrelets at the time of nesting, any pixels recorded by NTEMS as experiencing stand replacing disturbance since 1998 were removed.

Defining available forest for nesting
The study area contained regions that are known to be generally unsuitable for murrelet nesting (e.g., clear-cuts and young, regenerating forest) that were removed from the environmental predictors prior to modeling to ensure only potentially available forest was considered (Acevedo et al. 2012, Guisan et al. 2017).
As opposed to the dominant approach in murrelet research that removes any forest younger than 140 years (as recorded by the VRI) because it is considered unsuitable nesting habitat, we exclusively used ALS and NTEMS data to identify unsuitable forest.This is justified because removing forest younger than 140 years would have erroneously excluded 21 nests from DS and 1 from CS due to imprecision within VRI polygons and age estimates.Instead, our method retained the entire nest sample and removed areas that held no murrelet nesting potential by excluding pixels that contained a forest height < 20 m or were clear-fell harvested since 1998 (Burger 2002, Mather et al. 2010, ECCC 2021).An additional benefit of our definition of available forest is that only one tree per pixel needed to meet the height criterion for the pixel to be retained, resulting in the inclusion of small areas of potentially suitable forest that are often overlooked by other habitat identification methods (Burger et al. 2018).

Nest site comparisons
Prior to HSM, all predictors (

Habitat suitability modeling
Ensemble of small models and MaxEnt models were first developed using all potential habitat predictors (Table 1) and then tested with only ALS predictors (removing CMI and dist_coast) to explore if these were sufficient to predict nest occurrence independent of any other data source.To explore the minimum information needed to predict nest occurrence, simpler MaxEnt models were built using bivariate combinations of all ALS predictors.Prediction values closer to one indicate a higher probability of nest occurrence.
The ensemble habitat prediction was created from the average of many small generalized linear models (GLMs) that predicted nesting probability using bivariate combinations of the predictors, nest locations, and pseudoabsences.(2015).
Maximum entropy models were tuned following the best practice presented in Cobos et al. (2019) with tuning parameters that achieved the highest AUC score used in the reported models.We used 10,000 randomly selected background samples from DS for model training.Maximum entropy models were internally evaluated in DS with bootstrapping running a 70% data split and 100 iterations.The final MaxEnt models were an average of 20 bootstrapped model replicates.

Model evaluation
All models were trained and internally evaluated in DS, then projected over CS to provide a fully independent transfer assessment.Evaluating our models on a fully independent dataset was crucial to assessing model performance and ecological fidelity, especially with our small sample where overfitting was a risk (Yates et al. 2018, Gregr et al. 2019).Additionally, a successful independent transfer would indicate that our model could be applied to other areas of coastal forest in BC for management purposes.To assess if models were being extrapolated to values outside the range of the DS training data, multivariate environmental similarity surface (MESS) analysis was conducted for each predictor set (Elith et al. 2010).
Two evaluation measures were used to assess model performance: area under the ROC curve (AUC) and the continuous Boyce index (Breiner et al. 2015, Guisan et al. 2017).Area under the ROC curve is a measure of classification accuracy between nests pixels and random background pixels where an AUC of 1 represents perfect accuracy and 0.5 represents a model that performs no better than a random classifier (Di Cola et al. 2016).The Boyce index assessed if our model could be used to infer habitat suitability, with a score closer to one indicating that murrelets preferentially selected forest with higher model prediction values (Hirzel et al. 2006).Preferentially selected areas were assumed to represent higher quality and more suitable nesting habitat.The top-performing model was chosen based on AUC and Boyce index scores of both the training area (DS) and the transfer area (CS; Breiner et al. 2015) with an AUC score > 0.7 and a Boyce index score > 0.5 used to identify a good-performing model (Raes and Steege 2007).Models with fewer predictors were preferred over more complex models with similar performance scores (Heinze et al. 2018).
Response curves for the top model were used to visualize the relationship between model predictors.Predictor importance was measured for MaxEnt models by measuring each predictor's contribution to the regularized gain in performance during training (Elith et al. 2011, Breiner et al. 2015) The spatial arrangement of predicted habitat from the top model was analyzed by reclassing predicted values into five classes each covering one-fifth of the prediction range.Connected pixels of the same class (hereafter called habitat clusters) were identified using the queens' case rule connecting neighboring and diagonally adjacent pixels (Hesselbarth et al. 2019).For each class, the average, standard deviation, and total area of habitat clusters were calculated and, to support interpretation, the average forest age in each class was recorded.

Comparisons with existing murrelet habitat mapping
Areas of DS that overlapped with LLAS habitat polygons were analyzed to explore how our predictions compared with existing maps of nesting habitat quality.Pixels with more than 10% of their area outside an LLAS polygon were not compared, and we excluded polygons of rank 6 because most of these areas did not meet our definition of available habitat.Predicted pixel values were plotted per LLAS rank and compared.Distributions of predicted habitat values within LLAS ranks were not normally distributed, so statistical differences were tested using a Dunn test adjusted for multiple comparisons using a Bonferroni correction (Dinno 2017).The null hypothesis that the median value of predicted habitat did not differ between two LLAS ranks was rejected if p < 0.025.

Software
All analysis was conducted in R version 4.2.0 (R Core Team 2021) using the packages: ecospat (Di Cola et al. 2016) for ESM, kuenm (Cobos et al. 2019) for MaxEnt, ROCR (Sing et al. 2005) for CS model evaluation, landscape metrics (Hesselbarth et al. 2019) for the spatial distribution of predicted habitat, stats (R Core Team 2021) for Shapiro Wilks and Wilcoxon tests, dunn.test(Dinno 2017) for LLAS comparisons, and ape (Paradis and Schliep 2019) for calculating Moran's I.
Airborne laser scanning predictors were produced using the R packages lidR (Roussel et al. 2020), and CMI was produced using ClimateNA at a 100-m resolution (Wang et al. 2016).

Nest site comparisons
The distribution of predictor values at nests and random background pixels were compared with unpaired Wilcoxon tests (Fig. 3).Forest height (FH) and FVC both showed a significant difference between nests and background pixels in both areas.
Nest sites showed a high variation in all predictor values, with similar distributions to the random background sample.Clayoquot Sound and DS nests showed similar distributions for FH, FVC, and CSA, but CS nests were closer to the ocean, wetter, and had less empty space within and under the canopy.

Overall model performance
All tested models performed better than a null model (AUC of 0.5) over the DS training area (Table 2).Ensemble of small models performed below our performance threshold when interpreted with AUC in both areas, and Boyce index scores were below our threshold in CS.The MaxEnt models with more predictors showed the greatest AUC and Boyce scores for DS but performed poorly in CS.Three MaxEnt models performed worse than a null model after transfer; however, the remaining models still provided reasonable classification accuracy with AUC scores > 0.6.

The top-performing habitat model
The top-performing model was identified as a MaxEnt model using two predictors: empty_gaps and FVC.The top model had an AUC of 0.77 with a small average standard deviation of 0.03 (Table 2).The variable FVC was the most important variable with a contribution to regularized training gain of 67.5%.The MaxEnt parameters used in the top-performing model were a regularization of 0.2, with only linear and product features used to fit the model.The probability of nest occurrence was found to (1) decrease as internal forest gaps increased and (2) increase as forest vertical complexity increased (Fig. 4, plots A and B, respectively).
Nest occurrence predictions from the top model (interpreted as habitat suitability) are shown in Figure 5.Both areas contained a wide distribution of suitability values, with DS containing more low-suitability nesting habitat with a mean value of 0.33 compared to CS with a mean value of 0.49 (Fig. 6).
Multivariate environmental similarity surface analysis showed that the only predictor in CS that did not fall within the range of the training data was CMI.For models built with all predictors, 60% of CS was environmentally similar to DS; but for models built with only ALS predictors, > 99% of CS was structurally similar to DS (MESS values > 0).The residuals of the topperforming model were normally distributed in both DS and CS with a low Moran's I value of 0.11, indicating a minimal risk of spatial autocorrelation.

Habitat spatial arrangement
Small clusters of high-suitability habitat were found across both areas with scarce habitat ranked in the prediction range of 1-0.8 (62 ha in DS and 18 ha in CS).In DS, the habitat range 0.4-0.2composed the largest areas, but in CS the higher habitat range 0.6-0.4composed the largest areas (Table 3 and Fig. 6).On average, most suitable habitat was found within the oldest forest, with lower suitability habitat found within younger forests.

Comparisons with existing murrelet habitat mapping
Top-model predictions occurring within previously mapped LLAS polygons were compared using a Dunn test.The null hypothesis that LLAS ranks shared the same median predicted habitat value was rejected for all comparisons, apart from ranks two and three that did not statistically differ (p-value = 0.049).
Our predictions broadly aligned with LLAS rankings of nesting suitability (Fig. 7), with higher quality LLAS ranks containing statistically higher median predicted habitat values.A broad range of predicted habitat values were found within every LLAS rank.Forest unranked by LLAS contained the lowest mean predicted habitat values but also the largest range of predicted habitat values.

Describing murrelet nest locations with airborne laser scanning (ALS) data
To explore ALS's capacity to describe nesting structure, we compared predictor values between nest pixels and pixels randomly sampled from the available forest in both DS and CS.This analysis indicated that only 2 100-m resolution ALS predictors, i.e., forest height and forest vertical complexity, were significantly different at nest locations compared to random forest locations for both DS and CS nests.This aligns with other murrelet studies that found measures of tree height and vertical complexity to be indicators of nesting suitability (Waterhouse et al. 2002, Plissner et al. 2015, Hagar et al. 2018).Nests sites can occur across a range of structural conditions within tall, old, forest (Silvergieter and Lank 2011a, Barbaree et al. 2014, Wilk et al. 2016), however, we did not expect such a large overlap between nest and random background forest values, with no single predictor appearing sufficiently distinct to identify nesting habitat from background forest in either area.It is expected that murrelets select nest sites based on multiple structural conditions, as opposed to a single dominant characteristic, which potentially explains the lack of clear distinction in Figure 3. Another likely explanation is that the 1-ha scale of our ALS predictors did not effectively capture key microhabitat structure, such as tree-level cover and the presence of suitable nest platforms.
For HSM, predictors are most appropriate if calculated at scales that best capture the direct features of selection, which for murrelets would be tree-level (e.g., height, overhead cover, platform availability) along with stand-level forest features (e.g., forest complexity).Clyde (2017) and Silvergieter and Lank (2011a) both found that nesting suitability can vary substantially between neighboring trees, meaning a pixel covering one ha could contain the structure of many unsuitable trees, thus diluting the signal of relevant tree-level nesting structure.However, we were restricted from conducting our analysis at finer scales by the ± 50m positional error in nest locations.Consequently, we only measured stand structure and did not directly capture tree-level

Model interpretation
Despite the range of predictor values at nest sites, our modeling results show that well-performing murrelet nesting suitability models can be built with ALS data.In both MaxEnt and ESM models, ALS predictors performed similarly to models built with both ALS and non-ALS predictors, suggesting that ALS predictors alone are sufficient for murrelet nest modeling.This may be because both non-ALS predictors, CMI and dist_coast, were worse at discriminating microhabitat differences among sites compared to ALS predictors (Hamer et al. 2021).Overall, MaxEnt models built with two ALS predictors emerged as the most appropriate modeling approach, although mixed model performance was observed post-transfer.A MaxEnt model created with the predictors internal forest gaps (empty_gaps) and forest vertical complexity (FVC) emerged as the top model, with good performance scores in DS and reasonable classification accuracy in CS.Importantly, we found that only two predictors were needed to model the distribution of nests contrary to Hagar et al. (2014Hagar et al. ( , 2018)), who required upward of five model predictors to successfully discriminate murrelet occupied forest stands (AUC 0.74).
Although our modeling approach was not designed to test hypotheses of why murrelets select certain nesting characteristics, our use of ecologically relevant predictors allows some inferences to be gleaned.Counter to our expectation, our top model indicated that high-suitability nesting areas did not contain higher values of empty_gaps, suggesting that murrelets are less likely to use areas with more internal gaps under and within the canopy (Fig. 4, plot A).However, after reflecting on the unexpected response of empty_gaps in our model, we now suspect that this predictor was inappropriate for describing the fine-scale structure that determines nest access.Murrelets are known to typically enter the canopy at, or just below, the mid-canopy level to make a stalllanding on the nest platform (Nelson andPeck 1995, Manley 1999), and although below-canopy movement has been observed (Manley and Kelson 1995), vegetative gaps at this level are unlikely to be used often by murrelets.If specific below-canopy gaps are regularly used for nest access (as we initially assumed), then the empty_gaps predictor was either too coarse to detect these individual flyways, or these flyways had little influence on nest selection.Instead, lower values of empty_gaps were likely associated with forest that contained more vegetative layers and horizontal cover, both of which have previously been shown to be positively associated with nesting suitability, potentially providing more nest platforms and cover from predators (Burger 2002, Hamer et al. 2021, ECCC 2021).
Additionally, we were surprised that the surface area of the canopy (CSA) did not appear in the top model, despite this predictor corresponding with greater canopy variability and the presence of emergent trees, two features thought to be key for murrelet nest access (Silvergieter and Lank 2011b).Instead, forest vertical complexity appeared as the most important variable in the top model, suggesting variation within the entire vegetation profile was more useful for describing nesting suitability compared to our measure of canopy variability.
To better unpack the fine-scale structural features that influence nest selection, we recommend that future ALS nest modeling test different ALS predictors, such as measures of internal spacing and horizontal cover found exclusively in the upper canopy, and predictors that directly assess emergent trees in the canopy.Airborne laser scanning voxel metrics and object detection methods are highly suitable ALS processing approaches to create such nesting predictors.

Transfer assessment
The performance of our top model decreased when transferred to a different coastal forest (CS), although reasonable performance was still observed.Drops in performance are common for HSM transfers and arise because of ecological differences between areas and statistical factors (Yates et al. 2018).
The largest ecological difference impacting our top model was the lower volume of internal gaps at both nest sites and surrounding forest in CS (Fig. 3, plot 4).Consequently, the detection of CS nests may have suffered because the relationship established in DS between internal gaps and nesting suitability was weaker in CS due to differences in forest structure.However, surveys of murrelet nests in DS and CS by Silvergieter and Lank (2011a, b) found that both areas had very similar nesting conditions, differing only by tree height, with CS nest trees on average 10 m shorter and more similar in height to neighboring trees compared to DS.This similarity in nesting structure is supported by plots 1 and 2 in Figure 3, which shows nests in both areas occurred in statistically taller and more complex forest, indicating similar nest selection behavior in both DS and CS.On balance, ecological differences between the study areas likely caused a small decrease in model performance.However, these differences were not large enough to compromise the reliability of the transferred model, and we have confidence that our model predicted real differences in nesting habitat across CS.
Statistically, model overfitting was observed in the most complex MaxEnt model with six predictors, giving high AUC values in DS but a prediction no better than chance in CS (Table 2).However, our top model with only two predictors should be more robust to overfitting due to its lower complexity (Gregr et al. 2019).
Additionally, AUC scores were expected to be lower in CS because there is considerably more suitable habitat (Fig. 6), which made the suitability of the CS nest sample less distinct from the surrounding forest, resulting in lower AUC classification accuracy.Considering these statistical factors, the transferred CS model maintained reasonable accuracy despite falling beneath our AUC performance threshold for a good model.

Spatial arrangement of predicted habitat
In total, more high-quality habitat was found in CS (Fig. 6) with a median habitat value of 0.49 compared to 0.33 in DS.High suitability forest was rarely found on the valley floor in DS, whereas in CS, most high-quality habitat was found in this area, reflecting the different harvesting histories of these areas (Fig. 5).
Although CS is not without anthropogenic disturbance, it provides a more reliable baseline of undisturbed habitat arrangement than DS (Zharikov et al. 2006, Bjorkman andVellend 2010).In both areas the most suitable habitat was found in small clusters with an average size of < 2 ha.On average, high-suitability habitat was clustered in small areas within older forest (as reported by the VRI) with low-suitability areas often associated with younger forest.This is consistent with murrelets' nesting association with old growth forests (Plissner et al. 2015, ECCC 2021).However, substantial variation in forest age was seen across every habitat class.This variation is likely caused by imprecision in the mapping and age estimation of VRI data, with many small remnant areas of older suitable forest occurring in stands mapped by the VRI as younger forest.Prior nest assessment (Silvergieter andLank 2011a, Burger et al. 2018) and modeling (Clyde 2017) support our conclusion that small high-suitability patches can be found within younger low-suitability LLAS polygons.

Comparison with existing maps of nesting suitability
Broad correspondence was seen between our predictions and LLAS habitat rankings.This gives confidence that accurate assessments of nesting habitat are being made, considering the disparate methodologies used (i.e., ALS vs qualitative aerial mapping).However, the large variation of our predictions within LLAS ranks shows some disagreement, with our top model predicting low-suitability habitat inside high-suitability LLAS polygons, and high-suitability habitat found outside of LLAS polygons (Fig. 7).Some of this error will come from inaccuracies in the low-intensity LLAS polygons that were previously shown to only predict 48% of nest within suitable habitat (Burger et al. 2018).Additional variation can be explained by our ALS model predicting finer-scale differences in habitat suitability within large LLAS polygons (often > 100 ha).This is a major strength of our ALS methodology and allows small areas of suitable forest to be mapped and targeted for management.High-suitability habitat found outside of LLAS polygons is likely explained by the prescreening conducted prior to surveys, which would have excluded small areas of suitable habitat in forest polygons that are otherwise considered unsuitable (< 140 years).Our results suggest that a substantial area of nesting habitat may occur within forests currently considered unsuitable for murrelets.
Notably, prior analysis of the LLAS dataset using high-intensity, fine-scale, qualitative mapping found that DS nests were mostly clustered within the top two LLAS ranks (Burger et al. 2018), whereas our model did not show this clustering, with nests occurring across a wider range of predicted values (Fig. 6B).

Management implications
Our top model provides a quantitative map of potential nesting suitability at a 100-m resolution that is appropriate to support reserve placement decisions and habitat identification.The continuous predictions of suitability produced by our model are particularly useful for habitat assessment.Murrelets typically nest in areas with the highest relative suitability (Wilk et al. 2016), meaning we can discern likely nesting sites by identifying the most suitable habitat in areas of otherwise similar habitat quality.Additionally, areas of disagreement between LLAS ranking and our predictions can be directly assessed by murrelet biologists for nest platforms and field verified, potentially revealing unknown suitable forest and areas of prior misclassification.Doing so would allow for more optimal allocation of resources, such as protected area placements and habitat monitoring efforts.Comparisons of our predictions within LLAS polygons and VRI polygons (Table 3 and Fig. 7), highlight that our model predicts substantial fine-scale variation in habitat quality that is not captured by existing mapping or forest inventory data, addressing the need for fine-scale predictions that better describe murrelet habitat in BC (Burger et al. 2020).Our model can also be used to track changes in structural habitat conditions over time, but this requires repeat ALS acquisitions such as those collected by forestry companies for inventory purposes.Therefore, a provincial database will be important for the efficient tracking of ALSinformed nesting habitat quality over time.
We are confident our model captures nesting suitability with sufficient accuracy to be useful for management (Raes and ter Steege 2007), however, some key limitations should be noted.Our method, unlike high-intensity LLAS, does not directly detect potential nest platforms and did not capture tree-level information, which is a core determinant of nest sites.If new nest data are to be collected, the accurate locations of nest trees should be recorded with, at most, 5-m positional error.This will enable substantially finer ALS predictors to be created that better capture the direct microhabitat features of selection, likely improving the accuracy of modeling.Additionally, our models did not incorporate any measures of nest success.Subsequently, our predictions may have been informed by suboptimal nesting locations where no chicks successfully fledged, potentially leading to the obfuscation of truly optimal habitat that is critical for population recovery.
Considering these uncertainties, our top-model predictions should be used alongside existing murrelet habitat maps when making spatial management decisions (van der Hoek et al. 2015).Areas of agreement between different habitat mapping products can be used to identify areas in which we are confident habitat predictions are accurate.Without further testing in other areas of BC with murrelet nesting data, such as Mussel Inlet and Haida Gwaii, we recommend that further application of our top model be limited to the south coast of BC, where it can inform the provincial land-use order for the recovery and protection of murrelet habitat (BC-FLNRORD 2021).However, our methodological approach is applicable to any area of potential murrelet habitat containing presence-only nest locations and ALS coverage.

CONCLUSION
In this study we found that: (1) ALS is an appropriate and effective tool for modeling murrelet nesting habitat.However, the positional error in nest locations prevented the testing of treelevel measurements, which limited our description of key microhabitat structures.(2) A MaxEnt model was able to successfully predict Marbled Murrelet nesting habitat using only two ecologically relevant 100-m resolution ALS predictors.This model performed well for the training area (DS) and, although still informative, worse for an independent area (CS).Our results aligned with previous studies of murrelet habitat associations with measures of forest height, forest vertical complexity, and internal forest gaps emerging as informative descriptors of nesting structure.
(3) Clayoquot Sound contained more suitable habitat than DS and, in both areas, the most suitable nesting habitat occurred as small clusters within older forest.(4) Our predictions broadly aligned with existing LLAS maps of nesting suitability, but captured considerable variation within each habitat rank, highlighting fine-scale differences in forest structure within forest polygons.Importantly, our model identified new areas of potentially suitable habitat that have not been captured by existing mapping.We recommend that future BC ALS models use MaxEnt with a restricted number of predictors and conduct a fully independent transfer to evaluate model overfitting.It is expected that more accurate and numerous murrelet location data would improve modeling accuracy and enable finer exploration of tree-level microhabitat features.With successful field verification, our model could be a highly useful tool to support the management and recovery of murrelet habitat in BC.

Fig. 1 .
Fig. 1.Maps showing the location of airborne laser scanning (ALS) coverage (our study area) in Desolation Sound (DS, area A) and Clayoquot Sound (CS, area B), BC, Canada.To show the forested areas within our study area, all forests recorded by ALS as over 20 m tall are shown on a blue scale.To highlight the different disturbance histories, areas harvested since 1930 are shown in red.These areas are a combination of vegetation resource inventory VRI) data (1930-2020; Government of British Columbia 2021) and Landsat landcover change attribution data(1984-2020;Hermosilla et al. 2015).Previous harvesting occurred prior to 1930 in both areas but these data were not available.

FVC 4 .
Internal forest gaps The volume of all unoccupied 1 m 3 voxels (3-D pixels) found below and among the canopy in a 1-ha area.Volume of empty space under and within the canopy.More gaps potentially make nest access easier.See label 4 in Figure2.empty_gaps 5. Climate moisture index CMI was calculated for the year 2000 and is the difference between annual precipitation and potential evapotranspiration.A potential proxy for epiphyte abundance and cover.Epiphytes can provide nesting substrate on large branches.CMI 6. Distance to the coast Euclidian distance to the coast for each pixel.

Fig. 5 .
Fig. 5. Maps of the continuous habitat suitability predictions made by the top-performing model for Desolation Sound (DS, area A) and Clayoquot Sound (CS, area B).Higher nesting suitability is represented by larger values (in purple), and lower suitability by smaller values (in yellow).

Fig. 6 .
Fig. 6.Histograms of the top-model predictions are shown in plot A, with Desolation Sound (DS) in purple and Clayoquot Sound (CS) in green.The area covered by each habitat prediction bin is represented on the Y axis in hectares.Plot B shows quartile boxplots of the predicted habitat recorded at nest sites in DS and CS.Exact nest values are represented as gray dots.Note that plots A and B share the same x-axis.

Fig. 7 .
Fig. 7. Quartile box plots showing the distribution of our topmodel predicted habitat values that lay within ranked low-level aerial survey (LLAS) areas.Black dots beyond the boxplot whiskers show outliers (defined as values outside of 1.5 times the interquartile range).Rank One represents the most suitable habitat.The number of pixels within each LLAS rank is represented by n.

Table 1 .
The six predictors tested in the murrelet nesting models and the habitat characteristic they are designed to capture.Note: ALS = airborne laser scanning.

Table 1
Comparisons between distributions were made using a nonparametric unpaired Wilcoxon test due to the non-normal distribution of pixel samples.
(Barbet-Massin et al. 2012, Breiner et al. 2015)012), background samples (pseudoabsences) were taken from 10,000 random points across the available forest and were weighted equally with the nest sample in model training.Generalized linear models were fit using quadratic relationships with no assumed interaction for each combination of predictors in DS.These models were then averaged (weighted by each models' performance) to give a final predictive model.Bivariate model performance was assessed through internal testing with a 70% training split and 20 replicates using the Somers' D score, which gave more weight to GLMs with good classification accuracy(Barbet-Massin et al. 2012, Breiner et al. 2015).Only models with a Somers' D score > 0.6 were retained to build the ensemble model and model averaging was calculated following the methods presented inBreiner et al.
. Spatial autocorrelation was assessed by examining the model residuals (calculated for each nest as one minus the predicted occurrence) for normality using a Shapiro Wilks test and computing Moran's I value, (values closer to zero indicate less risk of autocorrelation; Halvorsen 2013, Daniel 2014).

Table 2 .
Habitat modeling results showing the modeling approach, predictor set, and model performance for the training area (DS, Desolation Sound) and transfer area (CS, Clayoquot Sound).The top-performing model is highlighted in bold.Note: AUC = area under the ROC curve; ESM = ensemble of small models; ALS = airborne laser scanning; MaxEnt = maximum entropy; FH = forest height; CSA = canopy surface area; FVC = forest vertical complexity; CMI = climate moisture index; empty_gaps = internal forest gaps; dist_coast = distance to the coast.

Table 3 .
The size and variation of habitat clusters across Desolation Sound (DS) and Clayoquot Sound (CS).Note: Avg = average, SD = standard deviation.