Assessing representation of remote sensing derived forest structure and land cover across a network of protected areas.

Protected areas (PA) are an effective means of conserving biodiversity and protecting suites of valuable ecosystem services. Currently, many nations and international governments use proportional area protected as a critical metric for assessing progress towards biodiversity conservation. However, the areal and other common metrics do not assess the effectiveness of PA networks, nor do they assess how representative PA are of the ecosystems they aim to protect. Topography, stand structure, and land cover are all key drivers of biodiversity within forest environments, and are well-suited as indicators to assess the representation of PA. Here, we examine the PA network in British Columbia, Canada, through drivers derived from freely-available data and remote sensing products across the provincial biogeoclimatic ecosystem classification system. We examine biases in the PA network by elevation, forest disturbances, and forest structural attributes, including height, cover, and biomass by comparing a random sample of protected and unprotected pixels. Results indicate that PA are commonly biased towards high-elevation and alpine land covers, and that forest structural attributes of the park network are often significantly different in protected versus unprotected areas (426 out of 496 forest structural attributes found to be different; p < 0.01). Analysis of forest structural attributes suggests that establishing additional PA could ensure representation of various forest structure regimes across British Columbia's ecosystems. We conclude that these approaches using free and open remote sensing data are highly transferable and can be accomplished using consistent datasets to assess PA representations globally.

INTRODUCTION

Protected areas (hereafter PA) are an integral component of biological conservation, designed to preserve ecosystem services and biodiversity both inside the PA and in some cases the surrounding regions (Chape et al., 2005; Watson et al., 2014). In recent decades, there has been a growing consensus of the need to conserve varying portions of the terrestrial area of the globe, with areal goals increasing over time (CBD, 2004, 2010). In the 2010s, the Aichi biodiversity target sought to conserve 17% of the terrestrial and 10% of the marine area of the globe (CBD, 2010). Nationwide, the Canadian government has set the goal of protecting 25% of Canada's terrestrial area by 2025 (ECCC, 2021). While increasing proportional ecosystem protection does in turn aid conservation, it does not guarantee the representativeness of the entire ecosystem, nor that all biodiversity within the PA will be effectively conserved (Hazen & Anthamatten, 2004).

Many conservation goals, both global and regional, are commonly based on the proportion of area protected, at least partly due to its ease of use and calculation (Brooks et al., 2004; CBD, 2010). However, while the area protected is a simple metric to report, other metrics can be more informative, with the potential to convey how effective a given PA is for protecting the inherent ecosystem services or biodiversity in the area (Butchart et al., 2015; Chape et al., 2005; Maxwell et al., 2020). Beyond areal extent, it is also relevant to consider the biases in PA placement, which are frequently located in fiscally cheaper, low‐productivity regions both globally (Joppa & Pfaff, 2009; Venter et al., 2014, 2018) and regionally, as is the case in British Columbia, Canada (Environmental Reporting BC, 2016; Hamann et al., 2005; Wang et al., 2020). The areal protection targets currently in place are much lower than what research would indicate is necessary to adequately protect biodiversity (Dinerstein et al., 2017, 2019).

In response to the proportion protected approach, a number of other methodologies have been developed to evaluate the effectiveness of PAs before these larger global targets have been met (Bolton et al., 2019; Gaston et al., 2006, 2008; Hansen & Phillips, 2018; Parrish et al., 2003). One recently identified concept in Canadian protected area management is ecological integrity. Ecological integrity is defined as an ecosystem having the expected “living and non‐living pieces for the region,” and where ecological processes occur at the expected frequency and intensity for the region (Parks Canada, 2019). Many potential ecological integrity indicators have been examined to capture biodiversity related processes within PA (Hansen et al., 2021; Hansen & Phillips, 2018). These indicators can then be interpreted manually or automatically, most often through examining temporal trends within the PA or by comparing the indicators to areas in known healthy reference ecosystems (Woodley, 1993).

Frequently, comparisons between PA and unprotected areas (UA) have been drawn in order to assess PA performance and health (Defries et al., 2005). This allows for the PA or PA network to be taken in context of surrounding and/or similar ecosystems (Wiens et al., 2009). There are, however, challenges associated with comparing the effectiveness of PA directly with UA. It can be difficult to identify suitable UA for comparison due to the increased prevalence of human pressure in UA (Geldmann et al., 2019), and the bias for PA to be in areas that would not have faced increased human pressure due to their remoteness (Joppa & Pfaff, 2009).

Ferraro (2009) prescribed the use of the counterfactual method based on comparing the outcomes following PA implementation with what would have happened if the PA was not implemented. The counterfactual method is widely used as a more accurate method for assessing PA management effectiveness (Coad et al., 2015; Eklund et al., 2019; Geldmann et al., 2019). This method is frequently employed by using matching methods, where treatment and control samples are similar with regards to topography, climate, land cover, or others, to select UA that directly correspond to the PA being analyzed (Geldmann et al., 2019; Ribas et al., 2020). Collecting field data across both the PA and its counterfactual UA is often time‐and‐cost prohibitive. Consequently, the increasing prevalence of freely‐available imagery has led to satellite remote sensing becoming an essential tool for PA monitoring (Nagendra et al., 2013).

The opening of the Landsat archive in 2008 (Wulder, Masek, et al., 2012) has played a significant role in the use of satellite imagery in conservation monitoring (Nagendra, 2008; Turner et al., 2015). The availability of 30‐m spatial resolution data since 1984 allows for assessment of temporal trends in satellite derived indicators (Bolton et al., 2019; Hansen & Phillips, 2018; Nagendra et al., 2013), while the global coverage allows for comparisons between similar and differing ecosystems (Nagendra, 2008; Wulder, Masek, et al., 2012). Leveraging free and open‐source optical remote sensing data products has allowed users to increasingly undertake comparisons across an entire jurisdiction's PA network (Bolton et al., 2019; Fraser et al., 2009; Pôças et al., 2011; Skidmore et al., 2021; Soverel et al., 2010), comparing them to ecologically similar UA (Buchanan et al., 2018; Turner et al., 2015). These comparisons allow for an assessment of the effectiveness of a given PA or the entire PA network at representing regional biodiversity trends (Bolton et al., 2019; Soverel et al., 2010; Turner et al., 2015).

Optical remote sensing technologies have offered a key approach to deriving indicators (Bolton et al., 2019; Burkhard et al., 2012; Fraser et al., 2009; Nagendra, 2008; Pereira et al., 2013; Soverel et al., 2010) and detecting key terrestrial processes (Turner et al., 2003) to assess PA effectiveness at conserving ecological integrity (Nagendra, 2001; Nagendra et al., 2013). These indicators derived from remote sensing technologies can be categorized and monitored at broad spatial extents and across temporal scales. Commonly used indicators include land cover proportion (e.g., forest type, wetland, and unvegetated; Parmenter et al., 2003, Olthof et al., 2006), tree species (Nagendra, 2001), habitat classification (Lucas et al., 2011; McDermid et al., 2005), spectral information (Feeley et al., 2005; Gillespie, 2005; Nagendra et al., 2010), spectral heterogeneity (Rocchini et al., 2010), and ecosystem structure (Cohen & Goward, 2004; Goetz et al., 2007; Pôças et al., 2011; Soverel et al., 2010) and function (Skidmore et al., 2021). Moreover, remote sensing technologies enable the monitoring of terrestrial processes, such as natural and anthropogenic disturbance regimes (Alsdorf et al., 2007; Bolton et al., 2019; Hermosilla et al., 2015b; Kerr & Ostrovsky, 2003), alongside biogeochemical cycles (Myneni et al., 2001), vegetation productivity (Running et al., 2004), and vegetation dynamics (Zhang et al., 2003). Diversity in forest structural attribute measurements, often derived from light detection and ranging (lidar) is also a strong indicator of biodiversity, providing habitat, influencing food quality, and mediating microclimates (Gao et al., 2014; Guo et al., 2017).

Lidar enables the accurate characterization of treed vegetation structure (e.g., canopy height, canopy cover, basal area) across forested areas by measuring the time it takes for an emitted pulse of light to return to the sensor (Lim et al., 2003). While the natural variation in vertical and horizontal forest structure has been extensively explored using lidar, comparisons between PA and UA have been less frequently drawn using these methods when compared to optical remote sensing (Nagendra et al., 2013). The lack of previous comparisons has likely been due to the frequently limited extents of lidar acquisitions, a problem that has recently been solved by generating wall‐to‐wall metrics. These wall‐to‐wall metrics can be created by combining lidar data with times series of Landsat data, generating forest structural attributes across large regions and even entire countries (Matasci, Hermosilla, Wulder, White, Coops, Hobart, Bolton, et al., 2018; Wulder, White, et al., 2012).

As Canada progresses towards the national goal of 25% of terrestrial area protected by 2025, there is a growing need to better understand how PA compare to UA with respect to location, ecological classifications, elevations, productivity, and forest structure. In this study, we (1) examine the hypothesis that British Columbia's PA network is biased towards high‐elevation, low‐productivity regions of the province using free and open remote sensing data products, and (2) identify underrepresented forest structures in PA in the province. To accomplish this, we examined the bias in PA placement by comparing ecoregional PA coverage and land cover classes by elevation, and disturbances by latitude across PA and UA in British Columbia. We examine representative forest structural attributes by comparing the distribution of key indicators by ecological zone to determine the differences between PA and UA to find the most and least similar represented forest structures throughout the network. We conclude by highlighting the usefulness of these globally available, high quality, consistent, and transferable datasets and methods for assessing PA effectiveness.

METHODS

Study area

The province of British Columbia, Canada, covers 94.4 million ha, of which ~64% is forested (BC Ministry of Forests, 2003), and encapsulates a wide variety of biomes and ecosystems. This diversity of ecosystems is in part due to the large area as well as variations in topography and climate. The existing Biogeoclimatic Ecosystem Classification (BEC) system disaggregates British Columbia's ecosystems into zones (Pojar et al., 1987). The broadest classification delineates 16 zones, which are further broken down into subzones, variants, and phases based on microclimate, precipitation, and topography (Meidinger & Pojar, 1991; Pojar et al., 1987). As a result, BEC zones vary widely in size (ranging from 0.25 to 17.5 million ha), and in number of subzones (from 1 to 43; see Table 1).

TABLE 1

Number of subzones, total area, and percent protected by Biogeoclimatic Ecosystem Classification (BEC) zone

Zone	Zone name	No. subzones	Area (ha)	% Protected
BAFA	Boreal Altai Fescue Alpine	2	6,286,778	30.1
BG	Bunchgrass	2	257,072	11.8
BWBS	Boreal White and Black Spruce	5	16,404,142	8.6
CDF	Coastal Douglas‐fir	1	251,232	4.8
CMA	Coastal Mountain‐heather Alpine	3	3,574,039	17.9
CWH	Coastal Western Hemlock	10	10,795,067	19.5
ESSF	Engelmann Spruce–Subalpine Fir	43	17,465,113	17.8
ICH	Interior Cedar–Hemlock	12	5,538,842	10.2
IDF	Interior Douglas‐fir	12	4,488,085	5.9
IMA	Interior Mountain‐heather Alpine	2	1,257,949	29.2
MH	Mountain Hemlock	6	4,059,301	19.8
MS	Montane Spruce	8	2,863,394	9.4
PP	Ponderosa Pine	1	294,985	7.1
SBPS	Sub‐Boreal Pine–Spruce	4	2,265,365	9.5
SBS	Sub‐Boreal Spruce	11	10,337,497	6.7
SWB	Spruce–Willow–Birch	6	8,655,855	23.3

Both the British Columbian (BC Parks, 2012) and Canada‐wide (Government of Canada, 2019) PA mandates commit to conserving ecological integrity across the network. The PA network in British Columbia is designed to serve both ecological conservation and human recreation aims (BC Parks, 2012) and consists of a network of PA and PA complexes (multiple nearby PA that share the same conservation goals), with large variations in size, ranging from 0.02 to 987,899 ha (Figure 1).

FIGURE 1

Terrestrial British Columbia including Biogeoclimatic Ecosystem Classification (BEC) zones and the location of protected areas (PA) selected in this study

DATA

Biogeoclimatic ecosystem classification and PA

Boundaries for BEC zones and subzones were acquired using the bcmaps R package (Teucher et al., 2021). Two BEC subzones were entirely subsumed by PA (Boreal White and Black Spruce—Very Wet Cool and Spruce–Willow–Birch—Very Wet Cool Shrub), whereas the Sub‐Boreal Pine–Spruce–Moist Cool subzone has no PA representation.

Boundaries for all PA in British Columbia were obtained from the Canadian Protected and Conserved Areas Database (available from https://www.canada.ca/en/environment‐climate‐change/services/national‐wildlife‐areas/protected‐conserved‐areas‐database.html), current as of December 2020, and includes the International Union for Conservation of Nature (IUCN) classification for each PA. PA were selected for analysis following the criteria outlined in Bolton et al. (2019). Only parks that belonged to IUCN classes Ia, Ib, II, and IV were selected, as these categories are considered strictly protected for conservation purposes, excluding monuments and areas that can include anthropogenic disturbances caused by natural resource development. PA smaller than 100 ha in size were also excluded from the analysis, as these mainly occurred in urbanized areas. After selection, 745 suitable parks managed under various jurisdictions (provincial, federal, NGOs), comprising 15.4% of the total terrestrial area of British Columbia, were studied. An equal sample of pixels equal to the area of PA or UA—whichever was lower—was randomly selected from both PA and UA for each BEC subzone. This sampling regime follows the counterfactual approach, accounting for bias in topography, climate, and climax species due to the methods used to delineate BEC zones and subzones (Pojar et al., 1987).

Digital elevation model

The Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) digital elevation model (GDEM V2, 30 m) was used to examine biases in PA land cover and ecological classification by elevation (Tachikawa et al., 2011).

Landsat derived datasets

Land cover, forest disturbances, and forest structural attributes for British Columbia were derived from the annual Landsat best‐available‐pixel (BAP) composites from 1984 to 2019 at 30‐m spatial resolution generated using the Composite2Change (C2C) approach (Hermosilla et al., 2016). These composites are generated by annually selecting the optimal observations, free from atmospheric effects (haze, clouds, cloud shadows), for each pixel from the catalog of available Landsat‐5 Thematic Mapper (TM), Landsat‐7 Enhanced Thematic Mapper Plus (ETM+), and Landsat‐8 Operational Land Imager (OLI) imagery acquired during Canada's growing season using the scoring functions defined in White et al. (2014). The annual BAP composites are further refined by applying a spectral trend analysis over the Normalized Burn Ratio at pixel level in order to remove unscreened noise, detect changes and fill data gaps with temporally‐interpolated values, resulting in annual, gap‐free, surface‐reflectance image composites (Hermosilla et al., 2015b). During this process, forest disturbances are detected, characterized and attributed to a disturbance agent (i.e., wildfire, harvest, non‐stand replacing disturbances) using a Random Forests classification model via the object‐based analysis approach (Hermosilla et al., 2015a) with an overall accuracy of 92% ± 2% (Hermosilla et al., 2016).

Annual land cover information for Canada was produced using the BAP composites following the Virtual Land Cover Engine framework (Hermosilla et al., 2018). This framework integrates post‐classification probabilities, forest disturbance information and forest successional knowledge with a Hidden Markov Model to ensure logical land cover transitions between years. The classification comprises 12 land cover classes organized as either non‐vegetated or vegetated. Non‐vegetated classes included water, snow/ice, rock/rubble, and exposed/barren land. Vegetated land cover classes discriminated among non‐treed and treed vegetation (land‐cover level). Vegetated non‐treed classes comprised bryoids, herbs, wetland, and shrubs. Vegetated treed land cover classes included wetland‐treed, coniferous, broadleaf, and mixed wood. Independent validation of the land cover maps indicated an overall accuracy of 70.3% ± 2.5%.

Wall‐to‐wall, 30‐m forest structure metrics (i.e., Lorey's height, total aboveground biomass, elevation covariance, and canopy cover) were also annually derived from the BAP composites using the imputation method described in Matasci, Hermosilla, Wulder, White, Coops, Hobart, Bolton, et al. (2018) and Matasci, Hermosilla, Wulder, White, Coops, Hobart, and Zald (2018). This method uses lidar and field plot data to estimate forest structure metrics from topographic and Landsat spectral predictors, using a k‐nearest neighbor approach. Reported accuracy for the structure metrics indicated a root‐mean‐square error ranging from 24.5% to 65.8% and an R 2 ranging from 0.125 to 0.699 (Matasci, Hermosilla, Wulder, White, Coops, Hobart, Bolton, et al., 2018).

Forest cover classes (deciduous, broadleaf, mixed‐wood, and wetland‐treed) were used to generate land cover masks to restrict the comparison of forest structural attributes to treed pixels. Pixels with harvest activity disturbances detected post‐1985 were also removed from forest structural attribute rasters in both PA and UA, in order to restrict analysis to non‐anthropogenically disturbed areas. All data sets are displayed in Figure 2.

FIGURE 2

Visualizations for all layers included in the analysis for Garibaldi Park and surrounding region (red outline) in British Columbia for 2015

Analysis

To determine bias in ecosystem representation in British Columbia's PA network, we compared BEC zone, land cover, and disturbance proportions within and outside the PA network. We employ the counterfactual approach by examining BEC zone and land cover as a function of elevation, and secondly compiled disturbance rates on a latitudinal gradient across the province. Forest structural attributes were then examined at a finer ecosystem classification level, statistically comparing PA versus UA across similar ecosystems. Forest structural means across BEC zones were calculated to determine which forest structures need additional representation in British Columbia's current PA network. All data manipulation and analysis was conducted in R version 4.1.1 (R Core Team, 2020) or Python version 3.8.10 programming languages.

Ecosystems, land cover, and disturbances

Biogeoclimatic Ecosystem Classification zones and land cover classifications were aggregated to both PA and UA in order to determine the proportion of each zone under the protected classifications, to examine progress towards the Aichi biodiversity targets. In this analysis, zones were used to examine categorical data (land cover and disturbance) for the period of 1984–2019. Land cover and BEC zones were further examined along an elevation gradient, at 50 m increments. Histograms of area by elevation were generated in order to examine the areal magnitude alongside the proportional coverage of land cover and BEC zones. This allows us to examine the amount of area protected at each elevation, as well as the differences between PA and UA. Forest disturbances (including harvesting) were aggregated along a latitudinal gradient at increments of 0.5°.

Forest structural attributes

The t‐tests for PA versus UA were conducted on all pixels selected for analysis by BEC subzone and forest structural attribute for 2015, and the Bonferroni correction was applied. The Bonferroni correction avoids spuriously significant results in multiple comparison tests by dividing the significant p value (0.01) by the number of tests (Bonferroni, 1936). Within each BEC zone, higher proportions of significant tests will indicate dissimilar subzones in each forest structural attribute. The mean values for PA and UA forest structural attributes were calculated, in order to examine the differences in their distribution and determine which structures and zones differ between PA and UA. Values were also converted into z scores to determine the greatest standardized vector magnitude when comparing canopy cover, elevation covariance, and forest height between PA and UA.

RESULTS

British Columbia's ecosystems are protected at varying rates across the province (Figure 3). Of the 16 ecosystems present in British Columbia, seven are protected at rates above the Aichi biodiversity target (17%). Only two zones (Boreal Altai Fescue Alpine and Interior Mountain‐heather Alpine) are currently protected at rates above the Canadian 2025 protection targets (25%). Zones with Douglas‐fir (Pseudotsuga menziesii) as dominant old‐growth components (Coastal Douglas‐fir and Interior Douglas‐fir) are the least proportionally represented zones in British Columbia, with 4.9% and 6.4% protected, respectively (Figure 3).

FIGURE 3

Areal proportion of Biogeoclimatic Ecosystem Classification (BEC) zones protected in British Columbia

As elevation increases, increasing terrestrial area is protected within the PA network until ~4000 m, upon which all terrestrial area is protected (Figure 4). When comparing between PA and UA, zones are protected at differing proportions. For both Figures 4 and 6, values across a given elevation show the proportion of that elevation represented by each BEC zone or land cover class in both (b) PA and (c) UA. Zones commonly found at high elevations, such as the Boreal Altai Fescue Alpine, are predominantly located in PA, however, little terrestrial area is found at these elevations. In low elevations, proportions of area protected also differ, with Coastal Western Hemlock having a large proportion of coverage in PA, while in UA, boreal black and white spruce are underrepresented. Generally, the remaining ecosystems are found at similar rates in both PA and UA (Figure 4).

FIGURE 4

(a) Histogram of area protected in British Columbia by elevation. Proportion of Biogeoclimatic Ecosystem Classification (BEC) zone by elevation for both (b) protected areas and (c) unprotected areas aggregated to a bin width of 50 m. (d) Histogram of area unprotected in British Columbia by elevation

FIGURE 6

(a) Histogram of area protected in British Columbia by elevation. Proportion of land cover by elevation for both (b) protected areas and (c) unprotected areas aggregated to a bin width of 50 m. (d) Histogram of area unprotected in British Columbia by elevation

Protected land cover also varies by proportion (Figure 5). Non‐vegetated classes of snow/ice, exposed/barren land, and rock/rubble have higher than average proportions protected, while mixed wood and broadleaf land cover classes are underrepresented. All other classes are found at rates similar to the overall proportion of the province protected (~15%; Figure 5).

FIGURE 5

Areal proportion of land cover class protected in British Columbia

Similar to BEC zones (Figure 4), land cover also varies with elevation (Figure 6). Expectedly, snow/ice make up a large proportion of PA at higher elevations. At lower elevations in UA, mixed wood forest is a more common forest type than in PA, while wetland classes (wetland, wetland‐treed) are less frequent in the 400–900 m elevation range in UA compared to PA.

Figure 7 shows the elevation distributions of BEC zones and land cover classes. Generally, BEC zones are found at similar elevation profiles in both PA and UA. Alpine BEC zones (Interior Mountain‐heather Alpine, Boreal Altai Fescue Alpine, and Coastal Mountain‐heather Alpine) are found at similar elevations across PA and UA, while other zones such as Sub‐Boreal Pine–Spruce, Ponderosa Pine, and Bunchgrass vary in their elevation profiles. Land cover classes show differences in the wetland, wetland‐treed, and mixed wood classes. The wetland classes are found at lower elevations in PA than UA, while the mixed wood class has more variation in PA.

FIGURE 7

Elevation boxplots for Biogeoclimatic Ecosystem Classification (BEC) zones (a), and land cover classes (b). Whiskers indicate first quartile minus the interquartile range and third quartile to the interquartile range. Box and interior vertical line indicate first quartile, median, and third quartile, respectively

Overall, the burned area of forested cells is similar between PA (2.5% overall) and UA (2.3%), while harvesting is much higher in UA (7.2%) than in PA (0.33%). Harvesting is more common at lower latitudes in UA than at higher latitudes. Fire shows similar, but not identical patterns across varying latitudes, with higher wildfire proportions at high latitudes and between 51–53° N (Figure 8).

FIGURE 8

Proportion of area disturbed by latitude from 1984 to 2019 in protected areas (a), and unprotected areas (b). Proportion of terrestrial area that is protected at each latitude (c)

Figure 9 shows the subzonal proportional significance (p < 0.01) grouped by ecosystem for the 496 comparisons of forest structural variables. Higher percentages confirm ecosystems that had increased number of dissimilar subzones for the specific indicator and shows that at least half of all subzones in each ecosystem are significantly different (the exception being Ponderosa pine, which consists of a single subzone that is not significantly different in canopy structure). Median proportional significance values for canopy height, canopy cover, and aboveground biomass are universally significantly different between PA and UA within the same ecosystem.

FIGURE 9

Summary of percentage of ecosystem subzones that have significant p values from a two‐tailed t test with the Bonferroni correction (n = 496) applied at a significance level of 0.05. Boxplots represent median, interquartile range (IQR), and extreme values (1.5 × IQR)

Forest structural attributes vary between PA and UA in British Columbia (Figure 10). The largest differences between PA and UA are found in canopy structure in the Coastal Douglas‐fir BEC zone, with the PA having much higher canopy structure values. As shown in Figure 9, forests are commonly significantly different when comparing PA versus UA across all attributes. When examining the forests on a BEC zone level, only one BEC zone has a >5% difference in vertical forest structure (coefficient of variation in vegetation returns), six BEC zones have >5% difference in canopy cover, and five BEC zones have a >5% difference in canopy height. Ponderosa pine has large differences in canopy cover and canopy height (>5%), but minor differences in elevation covariance (only 0.25%; Table 2). PA in the Ponderosa Pine, Interior Mountain Heather Alpine, and Coastal Douglas‐fir have more aboveground biomass than in UA in corresponding areas (Figure 10).

FIGURE 10

Z scores of forest structural attributes in protected areas (PA), unprotected areas (UA), and their differences across Biogeoclimatic Ecosystem Classification (BEC) zones in British Columbia

TABLE 2

Mean values of forest structural attributes in protected areas (PA), unprotected areas (UA), as well as the percent difference between the means

Zone	Elevation covariance			Canopy cover (%)			Canopy height (m)
	PA	UA	% Difference	PA	UA	% Difference	PA	UA	% Difference
BAFA	0.39	0.39	0.01	46.94	48.47	3.17	14.03	13.11	−7
BG	0.38	0.38	−1.41	61.8	58.71	−5.26	23.00	21.58	−6.59
BWBS	0.37	0.38	2.62	67.18	64.72	−3.8	15.83	15.69	−0.9
CDF	0.31	0.33	6.35	89.38	83.69	−6.79	27.35	26.58	−2.89
CMA	0.38	0.39	1.63	62.65	64.01	2.13	19.91	20.01	0.48
CWH	0.34	0.33	−3.83	83.94	85.26	1.55	21.58	22.34	3.4
ESSF	0.37	0.37	0.34	61.7	64.97	5.04	18.90	18.91	0.07
ICH	0.36	0.36	−1.75	81.24	83.52	2.73	22.39	22.18	−0.98
IDF	0.36	0.36	−0.3	67.42	67.9	0.71	21.98	22.49	2.29
IMA	0.38	0.36	−3.72	68.17	62.07	−9.83	22.53	21.06	−6.98
MH	0.36	0.36	0.25	76.87	77.85	1.26	19.42	18.31	−6.07
MS	0.35	0.35	0.31	57.99	60.41	4.01	20.64	20.86	1.04
PP	0.36	0.37	0.25	57.92	48.93	−18.36	19.88	18.03	−10.24
SBPS	0.36	0.35	−1.97	32.98	34.63	4.76	18.00	18.70	3.75
SBS	0.37	0.37	−0.4	62.24	67.25	7.45	18.51	18.67	0.82
SWB	0.39	0.39	−1.22	56.67	57.71	1.8	13.78	13.67	−0.83

Note: Zones with more than a 5% difference are shown in boldface type.

DISCUSSION

The recent global availability of freely‐available, open‐source, consistent, and accurate remote sensing data products allow researchers to examine issues of representation of PA compared to UA, and regional ecosystems in novel ways (Bolton et al., 2019; Hansen & Phillips, 2018; Soverel et al., 2010). Additionally, the capacity to track forest structural attributes, a key indicator of forest biodiversity (Guo et al., 2017), across wide swaths allows for informed decisions on potential locations of new PA that capture previously underrepresented forest structure conditions (Noss, 1999). By applying this analysis to an entire PA network across BEC zones (or other ecological classifications), it becomes possible to determine which BEC zones need additional representation (the proportional metric), as well as the types of forest structures that should be represented to ensure adequate biodiversity protection (Lemieux & Scott, 2005).

Internationally, biodiversity preservation targets aim to protect a proportion of the total terrestrial area (CBD, 2010). Frequently, new PA are placed in high‐elevation, low‐productivity ecosystems both globally (Joppa & Pfaff, 2009; Venter et al., 2014, 2018), and in British Columbia, as confirmed by our analysis of ecosystem (Figure 3) and land cover (Figure 5) proportions. Alpine ecosystems are more commonly protected (Figure 3), as are the land covers commonly present within them (rock/rubble, snow/ice, exposed/barren land; Figure 5). As elevation increases, these ecosystems and land covers begin to dominate the proportional representation (see Figure 4 and 6). Differences between elevation profiles in land cover classes and BEC zones were also found, with the largest difference being that wetland classes were found at lower elevations in PA (Figure 7).

In high‐elevation ecosystems, Boreal Altai Fescue Alpine dominates the PA proportions above 3000 m, replacing the Coastal Mountain‐Heather Alpine ecosystem found in UA (Figure 4). These zones were both protected at rates above the average (Figure 3), and above the Aichi biodiversity targets. Interior Mountain‐heather Alpine had large differences in canopy cover and canopy height, while Boreal Altai Fescue Alpine only showed large differences in height. The Coastal Mountain‐heather Alpine did not any have large forest structural attribute differences (Table 2).

Distribution of disturbances followed a similar pattern to that reported by Bolton et al. (2019). Thus, the area affected by wildfires is comparable between PA and UA and at mid latitudes (51–53° N), while harvesting activity is more prevalent in UA and at low latitudes (Figure 8).

Our analysis shows that the majority of structural attributes were significantly different between the protected and unprotected forest stands across BEC subzones (Figure 9). In the south, Coastal Douglas‐fir, a zone with a single subzone, had large variation between PA and UA in two of four forest structural attributes examined. The unprotected forests were significantly less tall, had significantly less canopy cover, and significantly higher elevation covariance (vertical forest structure; Figure 10). In addition, it was the least protected BEC zone by area, with only 4.9% of the total terrestrial area protected. In this specific BEC zone, not only does additional area need to be protected to meet national goals, different forest structures need to be included in new PA (Paillet et al., 2010).

Utilizing this information on the proportion of BEC zones protected (Figure 3), as well as their forest structural attributes (Table 2), it is possible to identify which forest structures need to be added to the PA network in British Columbia. Those BEC zones with large differences (identified as being >5% change from PA to UA) suggest additional protection would encapsulate these underrepresented forest structures. For example: the forests in the Bunchgrass zone have large differences in both canopy cover and canopy height, with the PA having larger values in both attributes (Table 2). New PA in this BEC zone could contain forests with shorter and more open forests. A future avenue of research could be to incorporate forest structural attributes into spatially optimized PA placement approaches (Christensen et al., 2009).

The advent of free and open global data sets can allow for the monitoring of protected area health across the globe (Nagendra et al., 2013). Analyzing large amounts of free and open data using open‐source software approaches offers previously unseen perspectives into protected area representativeness. There are some challenges associated with this, namely: optical imagery archives being scarce in some regions due to imagery acquisition policies (Wulder et al., 2016), clouds and atmospheric interference (Li & Chen, 2020), lack of aerial lidar data available, and varying land cover legends used in differing regions and mapping products (Herold et al., 2008). New data and satellite missions are being introduced that can meet these challenges at a spatial resolution of 30 m or less. These new missions, including Landsat‐9 and Sentinel‐2, and spaceborne lidar such as GEDI (Dubayah et al., 2020) and ICESat‐2 (Neuenschwander et al., 2020), can provide global coverage of various new datasets: forest structural attributes (Potapov et al., 2021) through similar imputation methods to Matasci, Hermosilla, Wulder, White, Coops, Hobart, Bolton, et al. (2018); global land cover maps (Potapov et al., 2020; Zanaga et al., 2021); and forest disturbance maps (Hansen et al., 2013). These novel datasets provide clear opportunities for regional to global analyses of PA versus UA to be conducted concerning forest structure.

Future research monitoring protected area health using satellite remote sensing could focus on implementing essential biodiversity variables (Pereira et al., 2013) into their monitoring scheme. Advancing research towards these variables would not only benefit PA monitoring projects, but also biodiversity monitoring projects across the globe. Beyond this, examining the recovery of forest structural attribute following disturbances in both PA and UA could assess the effectiveness of PA for promoting regeneration.

CONCLUSION

We identified biases in British Columbia's PA network for PA to be placed in high‐elevation BEC zones, commonly dominated by low‐productivity land covers. We examined the disturbance regimes of PA versus UA by latitude, finding that wildfires are similar, while harvesting differs across the province. We then compared the forest structural attributes across all BEC subzones, finding that the majority of subzones have significantly different forest structures. Beyond this, we identified BEC zones with large variation in mean forest structural attributes. When new PA locations are decided upon in British Columbia, they could take forest structure into consideration, as wall‐to‐wall coverage of forest structural attributes becomes available. Novel datasets can allow this methodology to be applied across large regions, in order to identify PA biases and underrepresented forest structures.

AUTHOR CONTRIBUTIONS

Evan R. Muise, Nicholas C. Coops, and Stephen S. Ban conceptualized research. Evan R. Muise and Nicholas C. Coops conducted data curation. Evan R. Muise conducted investigation and formal analysis. Evan R. Muise and Txomin Hermosilla conceptualized methodology. Evan R. Muise wrote original draft and visualizations. All authors involved in review and editing.

CONFLICT OF INTEREST

The authors declare no conflict of interest.