Soundscape mapping: understanding regional spatial and temporal patterns of soundscapes incorporating remotely-sensed predictors and wildfire disturbance

Increased environmental threats require proper monitoring of animal communities to understand where and when changes occur. Ecoacoustic tools that quantify natural acoustic environments use a combination of biophony (animal sound) and geophony (wind, rain, and other natural phenomena) to represent the natural soundscape and, in comparison to anthropophony (technological human sound) can highlight valuable landscapes to both human and animal communities. However, recording these sounds requires intensive deployment of recording devices and storage and interpretation of large amounts of data, resulting in large data gaps across the landscape and periods in which recordings are absent. Interpolating ecoacoustic metrics like biophony, geophony, anthropophony, and acoustic indices can bridge these gaps in observations and provide insight across larger spatial extents and during periods of interest. Here, we use seven ecoacoustic metrics and acoustically-derived bird species richness across a heterogeneous landscape composed of densely urbanized, suburban, rural, protected, and recently burned lands in Sonoma County, California, U.S.A., to explore spatiotemporal patterns in ecoacoustic measurements. Predictive models of ecoacoustic metrics driven by land-use/land-cover, remotely-sensed vegetation structure, anthropogenic impact, climate, geomorphology, and phenology variables capture landscape and daily differences in ecoacoustic patterns with varying performance (avg. R 2 = 0.38 ± 0.11) depending on metric and period-of-day and provide interpretable patterns in sound related to human activity, weather phenomena, and animal activity. We also offer a case study on the use of the data-driven prediction of biophony to capture changes in soniferous species activity before (1–2 years prior) and after (1–2 years post) wildfires in our study area and find that biophony may depict the reorganization of acoustic communities following wildfires. This is demonstrated by an upward trend in activity 1–2 years post-wildfire, particularly in more severely burned areas. Overall, we provide evidence of the importance of climate, spaceborne-lidar-derived forest structure, and phenological time series characteristics when modeling ecoacoustic metrics to upscale site observations and map ecoacoustic biodiversity in areas without prior acoustic data collection. Resulting maps can identify areas of attention where changes in animal communities occur at the edge of human and natural disturbances.


Introduction
Biodiversity is increasingly at risk due to the alteration of the environment by direct human impact and induced climate change (Dirzo et al 2014).Passive acoustic monitoring has emerged as a method to understand biodiversity dynamics in natural and human-impacted landscapes (Sueur and Farina 2015).Ecoacoustic methods summarize the acoustic environment, collectively called a soundscape.Recording soundscapes can capture unique signals of acoustic activity for a place and time and serve as proxies for community-level biodiversity (Sueur et al 2014, Alcocer et al 2022).Soundscapes can be separated into discrete components made up of vocalizing animals such as birds, insects, and frogs (biophony), weather phenomena (geophony), anthropogenic sounds (anthropophony), and, in some cases, periods when recordings do not capture any signals (referred to as quiet here).These soundscape components have been quantified directly (Mullet et al 2016), but more often, mathematical interpretations of acoustic features known as acoustic indices are used as proxies for understanding and summarizing patterns in acoustic biodiversity (Alcocer et al 2022).
The relatively discrete spatial footprint of acoustic recording units (ARUs) (Hao et al 2021), environmental conditions (Haupert et al 2022), and project constraints (Scarpelli et al 2022) limit the effective spatial coverage of ARUs and lead to difficulty in understanding regional or heterogeneous landscape biodiversity.Filling the gaps between ARUs with more continuous soundscape sampling (i.e. more ARUs) presents logistical challenges, such as site access and large data volumes.Alternatively, interpolation methods that link environmental characteristics and acoustic observations can be used to generate continuous gap-filled maps.These maps can be valuable for biogeographic uses like identifying areas and times of species occurrence (Holgate et al 2021, Desjonquères et al 2022), identifying sources of anthropogenic noise pollution in urban planning applications (Barber et al 2011), and tracking environmental disturbance (Sethi et al 2020, Rappaport et al 2022).
The spatial context of recording location is frequently incorporated in experimental design to characterize human impact (Fairbrass et al 2017) or variation in habitat and vegetation characteristics (Boelman et al 2007, Tucker et al 2014).Past research has strengthened the connections between environmental characteristics and acoustic indices (Rodriguez et al 2014, Machado et al 2017, Müller et al 2022) and derived acoustic features (Sethi et al 2020).Sethi et al (2020) used acoustic features from deep learning models to visualize different habitat types, aboveground biomass, seasonality, and daily patterns in acoustic space.Do Nascimento et al (2020) found habitat types had predictable soundscapes, and canopy cover was important in modeling acoustic indices.Recently, acoustic indices were shown to be positively correlated with bird species richness and tree species diversity (Beason et al 2023).Many ecoacoustic studies incorporate measurements of environmental characteristics, establishing the ecological foundation that ecosystem state is linked to both acoustic activity and biodiversity.However, spatially-continuous environmental measurements are required to utilize these relationships.Remote sensing and other environmental spatial variables (e.g.climate) can provide this link to generalize acoustic observations, such as detected species or acoustic indices, thus enabling an understanding of how soundscapes vary across heterogeneous landscapes in which ARU recording sites are embedded.
Several studies have explored continuous spatial mapping of acoustic activity by generating interpolated layers of ecoacoustic metrics with remote sensing and field-based observations of vegetation and habitat metrics (Krause et al 2011, Turner et al 2018, Hao et al 2021, Benocci et al 2022).An ambitious study focused on generating a continuous continental USA map of sound pressure levels to generalize acoustic point observations using a range of explanatory geospatial data (Mennitt et al 2014).However, this mapping study was limited to data from US national parks and urban cores, and it was generalized across regions with sparse observations.Given these limitations, Mennitt et al (2014) found informative relationships among geospatial variables, like nighttime lights and sound pressure levels.At a regional scale in south-central Alaska, Mullet et al (2016) spatially predicted anthropophony and biophony, resulting in anticipated relationships with landscape characteristics such as human infrastructure and activities or rivers and wetlands, respectively.More recent work has established the utility of ecoacoustic metrics in a species distribution modeling (SDM) framework (Desjonquères et al 2022), a novel direction for ecoacoustic work to aid in conservation and biomonitoring efforts by linking mechanistic relationships of communities with landscape characteristics.However, the success of acoustic SDMs relies on the ability to use reliable ecoacoustic metrics derived from a consistent, automated approach and environmental variables related to soundscape dynamics (Desjonquères et al 2022, Quinn et al 2022).
Interpolated ecoacoustic maps not only uncover more complete spatial patterns in soundscapes but can provide insight regarding the temporally dynamic nature of biodiversity, such as the response of animal communities following disturbances, like wildfires.Soundscapes have been shown to capture post-fire changes in soniferous communities in a mosaicked closed-canopy and agricultural frontier Amazonian forest (Rappaport et al 2022), biotic recovery one-year post-fire in Brazilian savannas (Duarte et al 2021), and phenological patterns of an acoustic index post-fire in sky island forested and shrubland landscapes (Gasc et al 2018).In the Western US, wildfires have emerged as a significant threat as increasing human development expands the wildland-urban interface (WUI; Radeloff et al 2018) and fires grow in size and intensity associated with drought conditions (Keeley and Syphard 2021).It is also evident that fires are an essential component of the dynamics of these habitats and that biodiversity is highest in landscapes with a mosaic of fire histories (Steel et al 2022).Soundscape observations or associated predictions following wildfire could be used to provide an improved understanding of animal community response and identify bioacoustic wildlife refugia.
Here, we use a suite of remotely sensed vegetation structure, phenology, nighttime lights, land use and land cover (LULC), climate, and geomorphology variables to predict the distribution of ecoacoustic metrics and bird species richness across Sonoma County, California, USA.Our goals were to: 1. Understand how spatial predictors derived primarily from remote sensing relate to ecoacoustic metrics.
We anticipate high relative importance of climate predictors in modeling soundscapes, as observed in previously published species distribution models (Burns et al 2020).We expect higher levels of biophony in more structurally-complex and productive forests due to niche availability.Importantly, we emphasize our modeling approaches are data-driven and not mechanistic.Relationships among predictors and responses are correlative.2. Compare how spatiotemporal predictions of ecoacoustic metrics vary across the county.We anticipate anthropophony predictions will correspond spatially with urban and suburban cores.At the same time, biophony will be related to areas with high bird, amphibian, and insect prevalence in the county, primarily in less-developed land uses.Acoustic indices should correspond to these dynamics with higher Acoustic Complexity Index (ACI) following biophony, Normalized Difference Soundscape Index (NDSI) increasing in high biophony regions and decreasing in anthropophonically rich areas, and higher values of H (Acoustic Entropy) in acoustically active spaces.3. Assess how predictions of biophony change due to varying levels of wildfire severity.For example, do elevated levels of biophony due to increased insect or amphibian activity following wildfires corroborate the trends observed after the 2017 Nuns wildfire in Sonoma County (Cook and Hayes 2020)?We anticipate biophony will noticeably change within the impacted fire perimeters.

Study area and data collection
This study used data from the Soundscapes to Landscapes project (Snyder et al 2022;soundscapes2landscapes.org, hereafter S2L), which focused on collecting acoustic data across the diverse habitats and landscapes of Sonoma County, California, USA (figure 1; 4152 km 2 ).Sonoma County is a Mediterranean coastal environment with a complex developed and natural landscape.The developed spaces along the county's WUI include cities, towns, residential areas, vineyards, and agriculture.Natural landscapes range from coastal shores to low mountainous terrain (up to 1369 m) with complex mixed oak-hardwood forests, riparian corridors, and Redwood conifer forests.The county has many protected spaces and a legacy of frequent fires in the eastern oak-hardwood forests.The S2L acoustic data campaigns targeted the bird-breeding season from late March through June 2017-2021 (figure 1).In addition to bird-breeding activity, frog abundance (e.g.Pseudacris regilla and Lithobates catesbeianus) peaks in late winter and early spring, coinciding with each species' breeding season (Cook and Hayes 2020), while insect stridulations peak in summer months (daytime activity primarily Cicadoidea and nighttime primarily Orthoptera).At each monitoring site, ARUs were deployed by citizen scientists and volunteer-based property owners.Sites were selected based on topographic lowland and highland stratification types, post-wildfire lands, and access to public or private lands.This methodology resulted in field sites with an approximate minimum distance of 500 m between sites deployed simultaneously, sampling across heterogeneous vegetation types, vocalizing species, and a range of human-impacted landscapes.ARUs recorded 1 min every 10 min with an average of 607 min (4.2 d deployment) of recordings per site with 144 min recorded daily, meaning seasonal variation was not captured at individual sites.Select sites had subsequent annual visits, which we treat as two separate sites in our study.There were two different ARU models: AudioMoths (AM) programmed with an upper-frequency range of 24 kHz (sampling rate = 48 000 Hz, digitization depth = 16-bit, signal-to-noise ratio = 44.2dB; Hill et al 2018) and LG phone recorders with an upper-frequency range of 22.05 kHz (sampling rate = 44 100 Hz, digitization depth = 16-bit).ARUs were affixed at approximately 1.5 m height to the nearest woody vegetation or post.LG ARUs make up 100% of the 2017 sites (n = 119) and transitioned to 100% AM ARUs in 2020 (n = 180), while 2021 had 24 LG sites and 442 AM sites (AM total n = 949; LG total n = 221).We used 1170 possible unique sites and filtered the potential 710 178 min of acoustic data to 474 865 min (7914 h) after sub-setting target periods-of-day described in section 2.2.We did not remove the LG ARU samples as the observations represented unique spatial and temporal observations in the dataset before wildfires.Yet, we acknowledge that including them may bias results and not reflect either ARU device fully.

Ecoacoustic metrics
We focus on seven ecoacoustic metrics representing acoustic activity (Quinn et al 2022(Quinn et al , 2023) ) and acoustically-derived bird species richness (Clark et al 2023;table 1).The ecoacoustic metrics include the soundscape components anthropophony, biophony, geophony, and quiet (ABGQ), and the acoustic indices include the ACI (Pieretti et al 2011), acoustic entropy (H) (Sueur et al 2008), and NDSI (Kasten et al 2012).The acoustic indices were selected from a combination of their ability to (1) summarize unique dimensions of acoustic data in the S2L dataset (Quinn et al 2023) and (2) represent the abundance of biologically relevant sounds (Alcocer et al 2022).ACI and H were computed in the R package seewave v.2.2.0 (Sueur et al 2008) using a frequency limit of 1-10 kHz (flim) for ACI and window length (wl) of 512 for H. NDSI was computed using the R package soundecology v. 1.3.3 (Villanueva-Rivera and Pijanowski 2018) with a fast Fourier transform window of 1024 (fft_w), anthro_min of 1 kHz, anthro_max of 2 kHz, bio_min of 2 kHz, and bio_max of 10 kHz.Importantly, our implementation of seewave's ACI differs from the soundecology implementation.Acoustic indices were calculated for each 1 min recording.
The presence or absence of a soundscape component was derived from a custom-built pre-trained convolutional neural network (CNN) trained on a withheld set of S2L acoustic data (Quinn et al 2022).The CNN classified ABGQ at 2 s increments in which ABGQ could each be present.Model classification accuracy (F(β = 0.75)) is 0.89 for anthropophony, 0.885 for biophony, 0.838 for geophony, and 0.864 for quiet (Quinn et al 2022).Presence/absence (1/0) was then averaged for each site-specific period-of-day resulting in a percentage of ABGQ activity.
Bird species richness was acoustically derived from three CNN bird-detection models that classified individual species from 54 commonly observed species in the study area that primarily represent non-coastal Combines an anthropophony and biophony-focused measurement of power spectral density.We expect NDSI levels to have lower values (towards −1) associated with anthropophony and higher values (towards 1) associated with biophony.
(Kasten et al 2012) a Conceptual foundation for soundscape components comes from Pijanowski et al (2011).
species (Clark et al 2023).For each of the 54 possible bird species, an optimal CNN and threshold per species were used to classify the presence or absence of the target species for each 2 s audio clip.We then used these presence and absence classifications at a site level to estimate the number of bird species at each site.If a species had more than three presences in a site deployment, it was considered present at that site.We used a threshold of three presences as a conservative cutoff to reduce false positive identifications (e.g.reduced inclusion of non-resident birds) following the approaches in Clark et al (2023).We calculated relative sun positions for every site and recording day using the R package suncalc (Thieurmel and Elmarhraoui 2022) with the America/Pacific time zone and site longitude/latitude referenced to the WGS84 datum.Ecoacoustic metrics were averaged to several daily periods delineated using the sun position for each site.The periods-of-day are: dawn chorus-one hour before and three hours after the astronomical end of the night (range for the end of astronomical night: 04:30-07:10, local time); midday-two hours before and two hours after solar maximum (range for solar maximum: 11:20-15:10, local time); dusk chorus-three hours before and one hour after the onset of the astronomical night (range for the onset of astronomical night: 17:40-22:40, local time); and night-two hours before and two hours after sun nadir, the darkest moment at night (range for sun nadir: 22:00-01:50, local time), all hereafter referred to as 'period(s)-of-day' .Because periods-of-day were not inclusive of the entire recording schedule, the 710 178 min was reduced to 474 865 min.Recordings corresponding to each period-of-day were used to generate site-mean values of 2 s soundscape components and 1 min acoustic indices, while bird species richness did not consider periods-of-day.

Spatial predictors
We grouped spatial predictors as follows: anthropogenic, climate, geomorphology, phenology, and structure (table 2).Original spatial resolutions varied from 1 to 464 m.Coarser-resolution variables were downscaled, and finer-resolution variables were upscaled, realigned, and resampled using bilinear interpolation (Hijmans 2022) to a common projected reference grid with 90 m spatial resolution.Visible Infrared Imaging Radiometer Suite (VIIRS) Nighttime lights was the only dataset to use Google Earth Engine's default nearest neighbors with pyramids export method for resampling.We chose the 90 m scale to approximate a 45 m effective ARU detection radius within which sound detection was highly likely.The assumed 45 m detection The areal density of vegetation from 0 to 5 m characterizes the vertical structural complexity of the understory.
(Continued.) distance was based on prior work that found 95% of birds were detected within 40 m of the ARU in forested sites (Darras et al 2018), that vocalization detection for soundscape metrics with higher quality ARUs typically occurred within 50 m (Shaw et al 2022), and a study of sound attenuation that found the average detection distance for 80 dB white noise was 44-45 m for frequencies between 0-20 kHz (Haupert et al 2022).However, it is noted that there was significant variation in these distances.
We calculated mean and standard deviation focal window values to represent an approximate maximum distance at which sites would not overlap in recordings (Krause et al 2011).This consisted of a coarse (495 m radius, 990 m resolution) and medium scale (225 m radius, 450 m resolution).Both resolutions were chosen so calculations did not have to resample the 90 m target grid into partial cells.These multi-scale statistical summaries were calculated for selected anthropogenic and forest structure measurements to represent the effects of distance on noise attenuation due to human density and physical barriers and to capture coarser-scale spatial patterns (table 2).

Anthropogenic predictors
We used ten anthropogenic variables, including building density, distance to buildings, distance to streets, nighttime lights, and the proportion of area covered by agricultural land use (see Supplement S.1 for processing details).Building density was calculated as the sum of building footprints in each 90 m raster cell and is a measure of the intensity of human-developed space.Distance to buildings and distance to roads are measures of proximity to developed spaces and are associated with decreased biotic activities in some cases (Ware et al 2015, Mullet et al 2016, Quinn et al 2022).Nighttime lights data were downscaled from 464 m to 90 m and included as a source of human impact using the Suomi National Polar-orbiting Partnership VIIRS Day-Night Band sensor's Stray Light Corrected monthly composite (Mills et al 2013).We included agricultural cover as Sonoma County has a range of agricultural activities and may represent non-urban anthropogenic sources of sound.Agriculture in the county is primarily low-production extensive agriculture along with higher production intensive agriculture and diverse agriculture where smaller intensive and part-time agriculture occurs, such as vineyards.These land uses have lower population densities with machinery activity on altered managed landscapes that do not have the same vegetation structural complexity as natural landscapes.

Climate predictors
Climate variables were derived from the California Basin Characterization Model monthly climate data from 2017-2020 at a spatial resolution of 270 m (Flint et al 2013).Water year 2021 was unavailable, so water year 2020 was reused for sites in 2021, as the two years were climatically similar.We aggregated data into corresponding water year quarters: October-December (Q1), January-March (Q2), April-June (Q3), and July-September (Q4), resulting in four sets of six predictors for a total of 24 predictors per water year summarizing the included S2L campaign (Burns et al 2020).We included the entire water year, including quarterly values that precede measurements and follow measurements (e.g. if a site recorded early in the year relative to Q4), because each quarter summarizes the unique seasonal conditions in Sonoma County that correspond to vegetation patterns.For example, Q1 climate variables precede all ARU observations.Nevertheless, these variables capture the county's highest seasonal precipitation, linked to the end of the prior fire season and net primary production in the spring in California (Alexander et al 2023).

Geomorphology predictors
Topography and abiotic factors have a strong influence on sound and spatial patterns of biodiversity (Zarnetske et al 2019), so we included a set of six geomorphological variables derived from an airborne lidar digital elevation model (DEM) resampled from 1 m to 30 m spatial resolution.The DEM was used to generate slope and topographic position indexes (both continuous and classified).We also calculated the Euclidean distance to the coast since it is correlated with the spatial pattern of habitat types.Euclidean distance to streams was calculated from the Sonoma County lidar-derived stream centerlines filtered to 1st order streams (geomorphological data: sonomavegmap.org).

Phenology predictors
We included phenology to represent dynamics in vegetation productivity each year (Radeloff et al 2019).Phenometrics were derived using the Framework for Operational Radiometric Correction for Environmental Monitoring (FORCE), which aids in processing multispectral satellite imagery (Frantz 2019).We used FORCE to perform atmospheric correction of available Landsat 7 and 8 and Sentinel 2A and B imagery from 1 January 2016, to 1 February 2022.We then used FORCE time series analysis to stack Normalized Difference Vegetation Index (NDVI) data and perform radial basis function interpolation to standardize the data to an 8 d time sequence.These NDVI time series were then used to calculate seven annual phenology variables spanning the water year (see 2.3.2).Annual metrics (30 m) were resampled to 90 m and included dynamic habitat indices (DHI) cumulative, minimum, and variance (Hobi et al 2017) and four NDVI metrics (5th, 50th, and 95th percentiles, and the seasonal difference between the mean of June and the mean of November data; Burns et al 2020).

Forest structure predictors
We incorporated forest structure variables derived from the Global Ecosystem Dynamics Investigation (GEDI) lidar mission (Dubayah et al 2020) and Landsat time series.GEDI began acquiring data in April 2019 and is temporally aligned with the acoustic dataset from 2019-2021.However, a fundamental limitation of GEDI is the discrete 25 m shot sampling strategy as data are acquired along the International Space Station orbital path.Since our modeling approach relies on continuous 90 m gridded variables to generate map products, we developed a method to spatially predict continuous maps of GEDI structure metrics at 30 m resolution for each year that acoustic data were collected.First, we fit harmonic regression models to the Landsat Collection 2 surface reflectance time series using the continuous change detection and classification algorithm (CCDC; Zhu and Woodcock 2014) in Google Earth Engine.Next, we temporally matched high-quality, ALS-collocated GEDI shots (supplement S.1) with CCDC segments and extracted harmonic regression coefficients associated with select spectral bands and indices.We used these coefficients in a Random Forest classifier (RF; Breiman 2001) to predict GEDI structure metrics annually on 15 June (from 2017 to 2021) at 30 m spatial resolution.We then calculated the medium and coarse focal mean and standard deviation values to represent horizontal structural heterogeneity and finally resampled all structure layers to 90 m.We derived the multi-scale structural variables to serve as proxies for physical limitations on sound propagation because an area with a denser surrounding canopy will increase attenuation (Pekin et al 2012, Trimpop and Mann 2014, Do Nascimento et al 2020).However, contrary evidence exists that no structural characteristic can determine sound transmission among habitats (Darras et al 2016).These metrics also provide an additional measurement of horizontal heterogeneity in structure, capturing transitions from edge habitats related to acoustic biodiversity (Barbaro et al 2022).This resulted in 15 annual GEDI metrics (table 2) with 60 additional focal mean and standard deviation variants totaling 75 structure metrics per year.

Site extraction of predictor data
We extracted spatial predictor data using the coordinates of each site and the underlying 90 m raster cell values.The annual temporal resolution of some predictors allowed us to extract predictor values that matched the site calendar year (structure, nighttime lights) or water year (climate, phenology) for each respective site.All other predictors (geomorphology and anthropogenic excluding nighttime lights) were extracted without temporal matching.

Spatial modeling
We implemented RF regression models to predict ecoacoustic metrics using geospatial predictors.RF can naturally account for complex predictor interactions and can be used to identify important predictors (figure 2).RF has also been successfully applied for SDM (e.g.Lorena et al 2011, Burns et al 2020) and geospatial sound level modeling (Mennitt et al 2014).RF models allowed us to (1) quantify prediction performance for each ecoacoustic metric, (2) determine which predictors were most important in explaining ecoacoustic metrics, (3) interpret the effects of top-performing predictors, and (4) spatially model ecoacoustic metrics across periods-of-day.All modeling analyses were performed in R version 4.2.1 (R Core Team 2022).

Variable selection
RF variable importance interpretation is sensitive to highly correlated predictors (e.g.Gregorutti et al 2017).A range of methods exist for RF variable selection, but agreement on efficacy varies depending on the analysis (Speiser et al 2019), and computational efficiency should be considered when not testing these methods explicitly.We selected a set of less-correlated predictors using a novel adaptation of the variance inflation factor (VIF) (Zuur et al 2010).Predictor selection consisted of three steps for each metric to reduce the potential set of 123 predictors: 1. We used a semi-informed selection method to build univariate RF regressors using R's caret's ranger RF implementation for each metric and predictor (n = 9840 models; Wright and Ziegler 2017, Kuhn 2022).Site average ecoacoustic metric values were used to keep by-metric predictor groups consistent, making variable importance comparable across period-of-day models.We calculated mean test R 2 after k = 10 cross-validation (CV) folds to estimate model fit and tuned the minimum nodes in each RF k-fold while the predictor was centered and scaled.Here, R 2 is the square of the Pearson's correlation between the test set response variable and predictor for each k-fold.2. The univariate RF test R 2 values were used to rank predictors within their respective groupings (anthropogenic, climate, geomorphology, phenology, and structure) to reduce within-group redundancy.
We used spatialRF's auto_vif function (Benito 2021) to select which variables should be retained based on a preferential ranking.We viewed this as advantageous because standard VIF procedures do not integrate the predictor's relationship to the response metric.This allows for a semi-informed selection of predictors, especially when two predictors have a similar VIF value.We used a cutoff threshold of 3, which is recommended when the ecological signal may be weak or highly variable (Zuur et al 2010).3. Finally, we used the within-group thinned predictors and the univariate RF test R 2 to rank the remaining predictors with a more conservative threshold cutoff of 10 in the auto_vif function.This step removed highly correlated predictors affecting variable importance interpretations, resulting in each metric's final VIF-selected predictor set.

Cross-validation (CV)
We implemented a geographically informed cross-fold validation (referred to here as geo-CV) method to evaluate final model performance and assess model stability.In k-fold CV, data are separated into k groups where one group is withheld to evaluate model performance as the test set, and the remaining k-1 groups are used for model training.However, when modeling involves spatial data, geographic context is a vital consideration when assessing model performance due to issues with spatial autocorrelation (Roberts et al 2017).Thus, we developed a custom geo-CV approach using geographic clustering.We used the ArcGIS Pro find point clusters function with a distance of 500 m from site to site to create site clusters (i.e.fold group IDs).We chose 500 m to match our coarse focal scale for predictors (495 m), a distance with significant acoustic attenuation so recordings would have a lower likelihood of overlapping (Haupert et al 2022).The clusters were then used in a grouped CV approach using caret's groupKFold function (Kuhn 2022) and k = 10 to generate the geo-CVs.We repeated the grouped k-fold splitting ten times using 10 starting seed values to randomize site assignment to geo-CV folds.The resulting folds were not even because cluster groups varied in the number of constituent sites; the training set size ranged from 932 (79.7%) to 1,129 (96.5%) sites among folds.

Model fitting
We used RF regression models implemented with caret's ranger RF (Wright andZiegler 2017, Kuhn 2022) to predict ecoacoustic metrics using VIF-selected predictor sets.Caret's training protocol allows for utilizing a tuning grid to optimize the RF hyperparameters mtry and min.node.size.The number of trees was fixed at 500, predictors were centered and scaled, and hyperparameter optimization was performed internally using a repeated-CV approach of 5 folds with 3 repeats.This modeling effort aimed to (1) examine the relationships of environmental measurements with ecoacoustic metrics and (2) generate highly accurate maps filling gaps among the sites in Sonoma County.This was not a mechanistic or process-based modeling effort.

Model evaluation and interpretation
We used average R 2 and root mean square error (RMSE) values of the withheld test sets in each geo-CV fold (n = 100) to assess model performance.
Variable importance was calculated using permutation importance for each RF model to assess the contribution of each predictor.When a predictor is permuted, the decrease in model performance (i.e. increase in error) is measured.The larger the reduction in performance, the more important the predictor is in the model.We assessed variable importance by calculating the mean and standard deviation of each predictor's raw permutation value across folds.The average global permutation values were then ranked to determine the overall importance of each metric's predictor set, not considering period-of-day.To summarize importance, we (1) identified the top five predictors that were consistently the most important and (2) aggregated predictor permutation values by group to understand the relative importance of each predictor group to ecoacoustic metrics.
We generated partial dependence plots (PDPs) for each metric's top five predictors to interpret the effects of variable importance on metrics.PDPs highlight relationships between predictors and metrics that are sometimes difficult to interpret due to RF models' non-parametric, complex nature (Maxwell et al 2018).We used the partial function from the R package pdp (Greenwell 2017) by modeling the predicted partial response for each predictor from 100 equally spaced observations along the predictor's training domain.
We used the 90 m, 2021 spatial predictors to predict the spatial patterns of period-of-day and metric using raster's predict function (Wright and Ziegler 2017).We created final predicted maps using the median predicted value for each period-of-day and metric combination.

Biophony following wildfires
Sonoma wildfire season (August-November) occurs primarily after S2L recording campaigns (March-July).Because pre-fire recordings were on a by-chance basis, we instead use the gap-filled, predicted biophony to understand how biophony changed following wildfires from 2017 through 2021.We targeted three wildfires during the five-year S2L recording campaign: the Kincade, Glass, and Walbridge fires (figure 1).Non-spatially overlapping S2L recording data include one (Wallbridge and Glass-2020) to two (Kincade-2019) years following wildfire disturbance.The Kincade fire began in October 2019, burning 77 758 acres in the eastern part of Sonoma County in mixed hardwood, conifer, shrubland, and grass terrain that overlapped in the south with the 2017 Tubbs fire, the north with the 2017 Pocket fire, and the Geysers fire in the north in 2004.The Glass fire, also in eastern Sonoma County, began in September 2020, burning 67 484 acres (26 876 acres in Sonoma) in complex landscapes dominated by hardwoods and conifer forests with smaller shrub and grasslands, agriculture, and urban land cover areas overlapping in the south with the 2017 Nuns fire.The Walbridge fire burned 55 209 acres in August 2020 in a complex rugged hill system in western Sonoma County and had a historically longer reburn interval, with the most recent historic fire occurring in 1944.The Walbridge fire land cover is dominated more by denser conifer forests (i.e.Douglas fir, coastal redwoods) than the Glass fire mixed with hardwood forests, particularly in the northern portion of the fire perimeter.
We used monitoring trends in burn severity classified difference normalized burn ratio fire severity maps to characterize the burn severity (Eidenshink et al 2007).We rescaled the classified severity maps from 30 to 90 m using a majority vote to identify burn severity classes considering unburned/low, low, moderate, and high severity areas.We also generated a 1 km wide buffer 1-2 km from the fire perimeter to assess year-to-year biophony patterns in unburned regions (i.e. a 'pseudo-control').We removed areas that burned in 2017 fires from the target fire areas and pseudo-control rings to avoid recent reburn signals (Kincade: Geysers, Pocket, and Tubbs; Glass: Tubbs and Nuns).Predicted biophony values were extracted for each fire severity class and non-burned buffer area to compare the average biotic acoustic response following wildfire.We did not assess other ecoacoustic metrics, such as anthropophony or quiet as these may relate more to the change in forest structure and, therefore, increased sound detection distances following severe burns, while biophony includes this effect as well as the reorganization of soniferous species such as insects and birds (Gasc et al 2018, Duarte et al 2021).

Aural verification
We used subsets of audio recordings to verify specific predicted patterns in our findings that targeted times or events.These data subsets were generated by sampling 30 1 min recordings to target the following questions that emerged from our findings and are detailed in supplement S.8: -What species are observed in nighttime recordings?A subset of recordings for each (1) high ACI at night (2) high biophony at night.-What species are observed following the targeted wildfires?
In addition to these verification subsets, authors and citizen scientists have listened extensively to the S2L dataset including an exhaustive manual classification of soundscape activities and bird species in 901 1 min recordings that represent biotically active recordings from across non-night times of day.

Predictor importance
The predictor selection method reduced our predictors from 123 variables to 25 for anthropophony, 26 for biophony, 25 for geophony, 14 for quiet, 26 for bird species richness, 21 for ACI, 26 for H, and 25 for NDSI (list of predictors and correlations in supplement S.3).Two predictors were shared among the eight ecoacoustic metrics (slope and VIIRS average nighttime radiance).In contrast, 27 predictors were unique to only one metric, and 51 predictors were not selected for any metric.Of the 51 excluded predictors, 49 were structural predictors, one was anthropogenic (90 m building density), and one was phenological (NDVI median).The significant reduction in structural predictors was due to the multiple spatial summary variants (e.g. 90 m, 225 m mean and SD and 495 m mean and SD; n = 75 originally) for each unique variable, leading to high levels of correlation.The only structural variables not represented in any final predictor set (i.e.all five summary variants were dropped) were PAI and RH98; however, RH95, RH98, and RH100 are highly correlated (90 m ρ = 0.999), and all approximate canopy height.
For predictor importance, we ranked predictors irrespective of period, meaning period-specific importance may differ.Overall, we observed climate predictors most frequently represented in the top five most important predictors across metrics (figure 4), with precipitation of quarter one (PPT Q1) in the top five predictors for five of eight ecoacoustic metrics (figures 5 and 6), while structure was the second most represented predictor group occurring in six ecoacoustic metrics.Anthropogenic predictors were in the top  five for anthropophony, geophony, and NDSI models.Phenology predictors were present in ACI, biophony, and bird species richness models, while geomorphology was present in all metric models except for biophony and H.
The average importance of predictor groups highlighted the prevalence of predictors in explaining ecoacoustic metrics (figure 4).In general, climate (for biophony, quiet, geophony, bird richness, and NDSI), structure (for anthropophony and H), and geomorphology (for ACI) groups had the highest prevalence of predictors in the top-five predictors for each metric.Anthropogenic predictors ranked in the top five less frequently than expected for all metrics except for anthropophony, while the phenology group was lower than expected in six metrics, geomorphology was lower than expected in four metrics, and structure frequency and climate were lower than expected in two metrics.

Predictor effects
We used the average top-five ranked predictors irrespective of period-of-day to generate PDPs for each metric (figures 5 and 6).PDPs revealed that most periods had similar patterns of predictor effect, yet with an offset.There were deviations from overall period-of-day patterns, particularly for biophony's DHI predictor effect.Anticipated patterns were seen for anthropophony with strong positive effects from nighttime lights and negative relationships with topographic slope (e.g.activity in flat areas) and less structurally complex areas.Nighttime lights had strong effects at lower values for anthropophony, geophony, and NDSI.Bird species richness had important variables primarily related to the spatial configuration of habitat in the county (e.g.climate, distance to coast, and phenology).
Accounting for correlated predictors removed in the VIF processing is essential for fully understanding the results.We have included spearman's rho correlation among the predictor variables (supplement S.3).For example, 10 structural variables with |ρ| > 0.95 were excluded for the biophony models, resulting in only canopy height (RH 95%) in the final model.Among these excluded, highly correlated structural variables was a positive relationship with VDR Bottom, indicating that increased structural skewness in the understory may be related to increased biophony.Comparatively, a positive correlation between RH 95% and VDR top may suggest that biophony decreases as the top of the canopy becomes more evenly distributed, while the positive correlation with canopy cover may suggest that a more open canopy potentially relates to increased biophony.

Spatiotemporal patterns
We found that the RF spatial predictions conveyed meaningful spatial patterns for all ecoacoustic metrics (figures 7-10).Comparison of ecoacoustic values among LULC types highlights differences in metrics across periods-of-day, particularly for ACI, anthropophony, biophony, H, and quiet (supplement S.4).Higher anthropophony was generally concentrated along the north-south urban corridors that bisect Sonoma County (figure 1), corresponding to high standard deviation values in vegetation structural variables (e.g.habitat edges and transitions).Spatial patterns of anthropophony follow expectations and are highest in agricultural and urban areas and along major roadways (figure 9).Comparatively, the lowest anthropophony occurred in forested LULCs.Temporally, high anthropophony occurred in agricultural and riparian areas at dawn and peaked at midday in all other LULCs.Midday and dusk anthropophony, though slightly lower in magnitude than the dawn levels, had relatively higher levels in less developed areas, particularly in areas zoned as land-intensive agriculture.
Biophony peaked at the interface of urban-to-natural spaces in areas where canopy height was approximately 10-15 m tall (figure 5) and during the dawn period for all LULCs.Specifically, canopy height was the most important predictor for biophony, with values decreasing when canopies were >10 m, typical of conifer forests.High levels of biophony appear along riparian corridors, in oak-hardwood forests, and along transitions from urban or agriculture to woodlands (i.e.WUI; figures 1 and 8).Biophony was lowest at night in conifer forests and midday for all other LULCs and had a secondary peak in agricultural, riparian, herbaceous, and shrubland areas at night.
The largest shift in countywide quiet was from low midday levels (i.e.large amounts of sound) to high levels at night (figure 9).This trend is most notable in non-developed forested areas, with the quietest regions in the northwestern, conifer-dominated forests (figure 1).Quiet decreases, and noise increases along transitions from less-developed areas to agricultural, residential, and urban areas where the lowest quiet levels are observed regardless of period-of-day.
Spatially predicted NDSI followed the general patterns found with biophony and anthropophony, with lower values found in urban cores similar to those described for anthropophony patterns (figure 10).However, differentiation in overlapping high anthropophony and biophony areas was poor.NDSI did capture anticipated positive values in less disturbed LULCs and negative values in more disturbed LULCs.
ACI had spatiotemporal patterns that diverged from other ecoacoustic metrics, being highest in the county's southern portion, the urban cores co-occurring with high anthropophony, some riparian corridors, and some hilltops to the northeast (figure 10).Spatial correlation was poor between ACI and biophony and bird species richness in the higher activity areas of the county (supplement S.8), and instead, it was highest across the county with anthropophony.Night had the lowest average values, but high ACI reflected select agricultural areas and hilltops to the northeast, and ACI had the most temporal variability among acoustic indices.High nighttime ACI values were associated with ARU interference and insect activity.
Broadly, H had the lowest values coinciding with urban areas, while the highest values were in the mosaicked oak-woodland regions in the east and more so in the conifer forests of the western part of the county (figure 10), capturing coarse level anthropophony and biophony patterns when grouped by LULC.Local regions of lower H values, comparable to urban areas, were interspersed in high-relief conifer forests, sparsely vegetated hills, and ridgelines.However, H had poor temporal differentiation.

Biophony following wildfires
Wildfire soundscape analyses revealed measurable deviations in biophony levels following wildfires (figure 11).We show dusk patterns in figure 11 to illustrate a general trend in the effect of wildfire disturbance found for all periods-of-day.Following fires, there was an increase in biophony relative to pre-fire years, particularly for moderate-and high-severity areas in the Glass and Wallbridge fires and all severities for the Kincade fire.This increase in post-fire predicted biophony is most notable for dusk and night periods and, to a lesser degree, for dawn except in the Kincade fire that shows a notable increase at dawn (period-of-day in supplement S.5).The midday patterns show an inverse pattern, with lower predicted post-fire biophony than other periods as severity increases except for the Kincade fire where levels stayed consistent.

Discussion
Maps of ecoacoustic metrics serve as valuable resources for understanding biodiversity and human impacts and may improve biodiversity monitoring efforts and decision-making.We provide an analysis that uses multiple landscape characteristics to model ecoacoustic metrics among discrete site observations to facilitate ecological interpretation of landscape characteristics on soundscape diversity and increased predicted biophony following wildfire disturbance.

Spatial patterns
Relationships in this analysis are among the observed ecoacoustic values and environmental characteristics of Sonoma County and are not necessarily mechanistic.For example, water year Q1 precipitation is the fifth most important anthropophony predictor.This importance is due to the strong spatial pattern of heavy winter precipitation in the northwestern region of Sonoma, which is sparsely populated and has low levels of anthropophony.In general, the importance of climate predictors is related to their spatial distribution in the county and not the mechanistic relationship with the ecoacoustic metric.Climate predictors were useful for building relatively accurate models, with the tradeoff of decreased interpretability when certain variables were selected and found important.Therefore we emphasize the correlative nature of machine learning modeling efforts, which may not lead to a mechanistic understanding between a predictor and response.
In general, climate predictors were the most important predictors across the various ecoacoustic metric models, a result that agrees with bird SDM work in Sonoma (Burns et al 2020).Climate is also strongly correlated with animal vocalizing behavior (Krause and Farina 2016) and encompasses various weather patterns that affect sound levels (Mennitt et al 2014).Climate predictors were three of the top five most important for predicting bird species richness.
Similarly, several SDM studies have found forest structure variables to be vital descriptors of bird species habitat (Vierling et al 2008, Goetz et al 2010) and national-scale richness patterns (Goetz et al 2014).In this study, structural predictors were the second most important predictor group in general.Interestingly, we found a negative relationship between biophony and canopy height above a certain threshold height (10-15 m).We anticipated that increased forest complexity represented by canopy height in the biophony model would be positively related to biophony based on the increased availability of resources and habitat space (e.g.Goetz et al 2007) and interactions between canopy height and other structural metrics may be informative (e.g. is canopy openness favorable for biophony).Investigation of the spatial distribution of canopy height reveals a structural divide between the northwestern, less biodiversity-rich conifer forests and the mixed LULC in the center of the county, where shorter-stature oak-hardwood forests are located.This pattern, we believe, is driving the negative relationship between some structural characteristics like canopy height and biophony.The importance of structural variables for predicting ecoacoustic metrics aligns with findings from prior ecoacoustic work (Boelman et al 2007, Do Nascimento et al 2020, Dröge et al 2021, He et al 2022).We expanded upon previous efforts by using novel structural metrics from earth observation sensors (GEDI, Landsat) that quantify vertical and horizontal heterogeneity, and we found that biophony decreased in taller forests with less complex understories relative to forests with a more open canopy cover and more complex understories (Supplement S.9).These transition zones in biophony support findings that a gradient in soundscapes exists across urban-wild gradients (Carruthers-Jones et al 2019).
We found ecosystem productivity important for predicting biophony, but the negative relationship between cumulative DHI and biophony for all periods except for midday is counter to prior findings (e.g.Radeloff et al 2019), though this may be related to the difference in extent of the studies.This may be driven by high cumulative DHI associated with coastal redwood forests and low biophony.However, the positive trend between DHI variance (i.e.variation in seasonal productivity) and biophony follows previous work demonstrating that breeding birds in temperate ecosystems utilize abundant resources only available in some seasons (Hurlbert andHaskell 2003, Radeloff et al 2019).Further, we also note that decreased biophony and increased quiet in more productive and structurally complex habitats may also be due to physical effects on sound propagation (e.g.Haupert et al 2022).One supporting line of evidence is the strong negative relationship between quiet and negative TPI values, representing valleys, where sound transmission would Among relationships not in the top five predictor results, we found spatial patterns of biophony increased with the distance from building structures and roads (Supplement S.7), corroborating previous work (Mullet et al 2016, Buxton et al 2019).Notably, this finding differs from our previous study in this area (Quinn et al 2022), which found that average site-based biophony values were not strongly modulated by distance from the road.This difference may be due to the use of the non-linear RF models here or increased biophony beyond approximately 1400 m from roads (supplement S.7), a distance we noted may significantly affect the modeling trend (Quinn et al 2022).
For predicting anthropophony, building density was more important than the distance to the nearest building, while the distance to roads was even less important.These findings are intriguing as anthropophony is primarily car-related noise in our dataset, but the lower importance of road proximity may be due to a lack of sites close to active roads.Our results also support previous sound-level mapping efforts that demonstrated the importance of including nighttime lights in models (Mennitt and Fristrup 2016), particularly, in our case, for anthropophony and bird species richness.The nighttime lights data highlight Sonoma's cities, where a quick rise in anthropogenic noise occurs and where we observe lower bird species richness.
We recommend future ecoacoustic mapping efforts validate mapping patterns with aural verification as hypothesized patterns were able to be investigated through post-modeling recording.This demonstrates specific investigation of soundscapes using emergent patterns from continuous mapping.Many of the most important predictors in our analyses are available at near-global extents, thereby facilitating the transferability of our methods to other regions where acoustic data are present, though may require investigation into the data's coarser spatial resolution.While we used climate datasets specifically generated for California, USA, several similar options exist globally, such as monthly climate variables from WorldClim (worldclim.org) or CHELSA (Karger et al 2017).We also used predictors from global datasets like VIIRS nighttime lights and high-resolution DEMs.Forest structure can be derived using our continuous GEDI modeling process, but coarser (1 km) aggregated global canopy height maps are available (Dubayah et al 2021).

Temporal patterns in predictive maps
The temporal patterns in ecoacoustic metrics are linked to changes in animal and human activities across the diurnal cycle (Fairbrass et al 2019).However, we observed H, NDSI, and, to a lesser degree, ACI poorly differentiated daily acoustic activity while soundscape components captured anticipated trends.We provide predictive maps highlighting the changes in soundscape structure with period-of-day (e.g.highlighting acoustic activity hotspots) (Holgate et al 2021).For example, across period-of-day, we see that biophony remains high in urban-adjacent residential and agricultural lands and eastern mixed oak/hardwood forests.At the same time, similar hotspots of bird species richness co-occur with biophony in the center of the county, outside of the highest urban density in Santa Rosa (figure 8).Anthropophony and quiet follow anticipated spatial patterns related to commuter traffic and have opposing trends throughout the day.Quiet is affected by some high levels of biophony at night in the northeast and along riparian corridors when significant insect activity was observed to be the dominant sources of biophony following aural verification on a subset of recordings (supplement S.8).The decreased quiet, during increased anthropophony and at times of increased biophony, describes silences and not a naturally quiet landscape.Comparatively, naturally-quiet landscapes contain biophony and geophony and are valued for human enjoyment (Farina et al 2021).Quiet, here, helps identify acoustically active landscapes and avoid less ecologically relevant spaces, especially for future soundscape analysis.
Our audio sampling scheme of 1 min every 10 min effectively captured detailed diel patterns in soundscapes (Quinn et al 2023).Our ability to discern patterns could be due to the volume of data combined with our use of average period-of-day values for each location or solar positioning to determine periods that better summarize daily trends.We initially attempted to model hourly patterns and found this potentially introduced unintended temporal smoothing effects on biotic soundscape patterns.Light intensity at sunrise is one of many factors hypothesized to drive dawn bird chorus patterns, and using more informed, non-hourly temporal windows can aid in capturing relevant biotic soundscape dynamics (Farina et al 2015).Comparatively, our coarser temporal averaging may have reduced the effects of brief sounds on metric values (Bradfer-Lawrence et al 2023).Though we did not model seasonal changes in soundscapes here, modeling approaches like these could highlight seasonal shifts in biophony based on seasonal predictors like phenology and quarterly climate.Furthermore, some predictors, like humidity or temperature, could be modeled at period-specific levels to capture their effects on acoustic patterns and animal behavior.

Acoustic indices
Among the three acoustic indices included here, NDSI had the best model performance and most consistent spatial patterns with anthropophony and biophony, the two soundscape components the index was designed to track.NDSI values captured anticipated trends across forested (higher values, more biophony) and developed areas (lower values, more anthropophony).H and NDSI showed similar patterns differentiating urban versus natural areas, but this differentiation did not follow anticipated trends relative to biophony or bird species richness.For example, H and NDSI both trended lower, indicating higher anthropophony in urban areas, which was anticipated and captured some of the patterns of biophony activity but struggled to differentiate areas where both biophony and anthropophony were present.Other applications of NDSI have suggested that the presence of both human and animal sounds relative to lower animal sounds can result in similar values (Benocci et al 2022), which may explain why NDSI was lower in riparian corridors and the urban-adjacent high biophony regions.Riparian corridors also had higher levels of predicted geophony which may be further lowering NDSI values due to the overlap in frequency between geophony and anthropophony.Comparatively, ACI poorly captured the transition to urban areas where ACI was positively related to anthropophony and did not capture the patterns in biophony or bird species richness, indicating a lack of robustness against human noise (Pieretti et al 2011) and low reliability in capturing biotic activity.Interestingly, geophony, NDSI, and H have minimal variation across periods when compared to the range of values for anthropophony, biophony, quiet, and bird species richness, reinforcing prior findings that acoustic indices are sometimes not robust against anthropophony and uncoupled from the biotic-signal in urban settings (Fairbrass et al 2017) and heterogeneous landscapes.The only exception is that ACI has observably lower levels at night and a daily trend in more developed LULCs, but the spatial patterns are challenging to interpret.
Reliable biodiversity monitoring tools require a consistent ability to interpret patterns that relate to meaningful biodiversity patterns.Our results help emphasize the complex nature of acoustic indices, mainly when used to interpret patterns across heterogenous LULC and different periods-of-day in the presence of confounding sounds (Alcocer et al 2022).In contrast, we demonstrate that maps of soundscape components (e.g.biophony), which are detected directly from the presence of acoustic events by our CNN, can capture anticipated spatial patterns and discern temporal patterns in acoustic records.

Biophony following wildfires
The predicted increase in biophony 1-2 years post-wildfire is similar to results from prior work that found rise in biophonic activity following wildfires (Gasc et al 2018, Duarte et al 2021).Duarte et al (2021) found daytime soundscape differences in burned and unburned Brazilian savannas were minimal one year after fire, but differences remained at night, a pattern also observed in post-fire areas of the arid sky-island ecosystems of the southwestern United States (Gasc et al 2018).We found similar outcomes: biophony was higher 1 year post-fire for dawn, dusk, and especially at night, while midday biophony was lower in more severely burned areas than in pre-burn and unburned areas, which follows general temporal patterns in the county.Additionally, bird species richness was generally lowest at night across the county, with higher richness concentrated outside of these fire perimeters, implying the increased post-fire nighttime biophony is due to increased activity in insect and possibly amphibian activity, considering the low prevalence of mammal calls in the S2L dataset (Quinn et al 2022) and high prevalence of insects in non-dawn recordings with high biophony (supplement S.8).This hypothesis is aligned with prior work in Sonoma County that found an increase in wetland amphibian species one breeding season after the 2017 Nuns fire (Cook and Hayes 2020) and may be investigated in future work by comparing wetland fractions relative to ARU locations.
Based on 2021 prediction maps (figure 10), the only acoustic index that captured a change in acoustic activity in wildfire-impacted areas was ACI, where values across all periods were higher than surrounding regions (most notable for the Kincade fire) due to insect activity.We focused analyses on bird breeding season and biophony patterns in post-wildfire areas, which may have different patterns when explored across seasons (Gasc et al 2018).Pre-and post-fire differences in animal communities are complex, and changes could influence quantification in biotic sound propagation due to vegetation change, which we do not explicitly isolate in our analyses.Future studies could investigate this by modeling structure-only relationships.

Future research directions
A limitation of our study was the lack of taxonomic granularity in automated biophony classification that required aural verification.Automated classification of taxonomic groupings (e.g.insects, amphibians, and birds) will allow for a more informed interpretation of diel patterns and post-wildfire changes in biotic activity, such as by revealing changes in taxa between burned and unburned areas with lower effort.Our bird species CNN was also not focused on classifying coastal bird species, resulting in model-driven lower species richness values along the coast.We also discuss biophony as a valuable class, assuming it reflects ecologically valuable patterns (e.g.Pijanowski et al 2011); however, removing generalist species from biophony may allow for more informative interpretations.Including finer detailed taxonomic data would require extensive improvements in labeling CNN training data and proper representation within the recording archive for balancing training classes.Our study follows prior work analyzing coarse-soundscape classifications, and over the time of development, semi-or unsupervised learning techniques (e.g.Sethi et al 2020) have enabled more accessible analysis of sub-genus taxonomic groups.These techniques still require extensive technical expertise and resources, and further work is needed to balance the effectiveness of such models in different acoustic recording domains to ensure that work is transferable.
Sound attenuation has been studied extensively in select environments for ARUs (Dixon et al 2020), but to balance complexity and interpretability, we assumed a single effective detection distance here across heterogeneous habitats and conditions though investigating the use of a structure-related coefficient of sound attenuation across multiple habitats is needed.We also attempted to account for variation in sound attenuation by including predictors such as multi-scale horizontal heterogeneity of forest structure and topographic position.Moreover, we focused analyses on frequencies below 10 kHz for ACI, H, and NDSI and 11 kHz for ABGQ.This upper-frequency limit could affect our interpretation of low biophony at night, missing active nocturnal animals such as bats (20-120 kHz), and is a crucial area for additional exploration in ecoacoustic studies.Another effect of the 11 kHz limit is that spreading attenuation may drive lower frequency sound attenuation across habitats (up to 80% of attenuation below 10 kHz) instead of habitat differences which account for a smaller portion of attenuation at these frequencies (Haupert et al 2022).Darras et al (2018) found that 95% of bird detections were included within 40 m of an ARU in agroforestry tropical forest plots, while Dixon et al (2020) found Audiomoth ARUs were capable of discerning recorded bird calls up to 100 m.However, LULC types in these studies were less structurally complex than those in Sonoma County, which may have allowed for increased detection distances.
Future remote-sensing ecoacoustic mapping will benefit from the availability of global datasets like GEDI canopy height or biomass to capture landscape characteristics with improved accuracy.Other available or near-future remote sensing datasets such as the surface water and ocean topography could provide spatial distributions of wetlands where wetland maps are unavailable related to amphibian distribution, or NASA's Earth Surface Mineral Dust Source Investigation could be used to tie disturbance events (e.g.environmental remediation, forest beetle outbreak extents) to soundscape patterns.

Conclusions
Researching the relationships among landscape characteristics and ecoacoustic measurements is valuable in understanding how patterns of ecoacoustically-derived biodiversity and human activities vary over time and across diverse landscapes.Such research can enable more effective interpretation and communication of results to wildlife managers, urban planners, or future ecoacoustic scientists.The S2L study deployed over 1,200 ARUs from 2017-2021 recording more than 700 000 min of acoustic data from across Sonoma County (4152 km 2 ) and work here utilized extensive high-performance computing resources.We aim to provide insight that allows researchers to better target areas and times of interest where audio collection may be most impactful for their needs using ecoacoustic metrics that reflect meaningful ecological signals.
Our efforts highlight the ability to spatially interpolate acoustic metrics from site-level recordings using predictive models informed by remote sensing.A particularly novel method used in this study is the incorporation of continuous forest structure metrics derived from spaceborne lidar and multispectral time series data for spatially-predicting ecoacoustic metrics.While much work remains in analyzing the influence of structure and phenology on ecoacoustic metrics, our approach enabled a practical analysis focused on quantifying changes in biophony associated with wildfire burn severity.The importance of specific structure and phenological metrics may vary in other studies, but, in our case, using the relative prevalence of predictor groups facilitates comparisons with important environmental characteristics identified in future studies.Ecoacoustic metric maps can also aid in identifying areas where animal communities may be affected by proximity to development or intense disturbances like wildfire.
Deploying a high-density network of ARUs currently is cost-and time-expensive; thus, continuous ecoacoustic spatiotemporal mapping can help to fill gaps among ARUs and provide a more complete picture of the spatial patterns of biodiversity.Furthermore, continuous acoustic maps can be used for urban planning efforts, such as estimating the influence of development projects, identifying areas to prioritize for sound mitigation efforts, and comparing anthropophony among different zoning and land use areas.Continuous maps of ecoacoustic diversity can also provide information on wildlife refugia in changing landscapes, quantify the value of naturally quiet or biodiverse landscapes, prioritize areas for restoration or conservation, help measure the effectiveness of management decisions, and aid scientists and wildlife managers in targeting areas and times of interest for more intensive field data collection.The prediction maps we have developed can help establish a baseline acoustic distribution in Sonoma County that will be valuable as landscapes change due climate and associated disturbances like wildfire, as well as anthropogenically driven land use change.

Figure 1 .
Figure 1.Sonoma County S2L site density for the 2017-2021 seasons.Three wildfire boundaries used in the wildfire case study are shown.

Figure 2 .
Figure 2. Spatial modeling overview.Italics reflect a method and product (red) or method detail (black) and arrows indicate the flow of processing.RF = random forest; CV = cross-validation; PDP = partial dependence plot.

Figure 3 .
Figure 3. RF model performance by ecoacoustic metric summarized as the test R 2 for the 100 final models for each period-of-day.Center points represent the mean R 2 with one standard deviation whiskers.

Figure 4 .
Figure 4. Prevalence of predictor groups for all models in top-five ranked predictors per model and metric (n = 400 model folds per metric).

Figure 5 .
Figure 5. Partial dependency plots (PDP) for the top five predictors, irrespective of period-of-day, for soundscape metrics, with units if applicable.PDPs show the sensitivity of the ecoacoustic metric to the predictor.The larger the range the more influential the predictor on the response.For example, increased Nighttime Lights, up to approximately 5 nW cm −2 sr −1 , increases predicted anthropophony.Predictors shown represent the top-five predictors for each metric and each line represents each unique geo-CV fold (n = 100).Color coding of predictor groups matches figure 4.

Figure 6 .
Figure 6.Partial dependency plots (PDP) for the top five predictors, irrespective of period-of-day, for soundscape metrics, with units if applicable.PDPs show the sensitivity of the ecoacoustic metric to the predictor.The larger the range the more influential the predictor on the response.Predictors shown represent the top-five predictors for each metric and each line represents each unique geo-CV fold (n = 100).Bird species richness reflects non-period specific prediction.Color coding of predictor groups matches figure 4.

Figure 7 .
Figure 7. Bird species richness patterns across Sonoma County, showing highway, freeway, and primary arterial roads in light gray/blue.

Figure 10 .
Figure 10.ACI, H, and NDSI patterns for dawn and night (the two periods with the largest differences in these metrics) across Sonoma County, showing highway, freeway, and primary arterial roads in light gray/blue.Scales are normalized within each metric where minimum = highest minimum predicted value (ACI = 57, H = 0.77, NDSI = − 0.4) and maximum = lowest maximum predicted value across periods (ACI = 77, H = 0.97, NDSI = 0.6).

Figure 11 .
Figure 11.Wildfire predicted annual biophony for the dusk period within the burned perimeter of three historic wildfires stratified by MTBS burn severity and a non-burned, pseudo-control sample 1-2 km away from the fire perimeter.Post-burn years are outlined in black dashed boxes.

Table 1 .
Descriptions of the ecoacoustic metrics used as response variables in this analysis.

Table 2 .
Predictor variables used for geospatial modeling, with predictor group, units, description, and source.

Table 2 .
(Continued.)PAVD 5-10 ah m 2 /m 3 The areal density of vegetation from 5 to 10 m characterizes the vertical structural complexity of the middle or upper canopy depending on canopy height.VDR bottom ah -Vertical distribution ratio (VDR; Goetz et al 2007) of the bottom third of the REP.VDR relates to the vertical skewness of the understory relative to the top of each vertical profile (e.g. higher VDR is associated with a a Calculated predictor focal mean and standard deviation for medium (radius = 225 m) and coarse (radius = 495 m) scales.This resulted in four additional variants of each predictor.b Buildings: hub.arcgis.com/datasets/sonomacounty::buildings/about.c Streets: Merrick & Company Project Number 50 013 734.d LULC, Coast, Streams, DEM: sonomavegmap.org.e VIIRS: (Mills et al 2013).f BCM: (Flint et al 2013).g Phenology: (Frantz 2019, Radeloff et al 2019).h GEDI: (Dubayah et al 2020).