Deforestation leakage undermines conservation value of tropical and subtropical forest protected areas

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2020 The Authors. Global Ecology and Biogeography published by John Wiley & Sons Ltd 1Section for Geography, Department of Geosciences and Natural Resource Management, University of Copenhagen, Copenhagen, Denmark 2UN Environment World Conservation Monitoring Centre, Cambridge, UK 3International Union for Conservation of Nature, Gland, Switzerland 4Institute of Marine and Antarctic Studies, University of Tasmania, Hobart, Australia 5World Agroforestry Center (ICRAF), University of the Philippines Los Baños, Los Baños, Philippines 6European Environment Agency, Copenhagen, Denmark


| INTRODUC TI ON
Protected areas (PAs) are a cornerstone of efforts to conserve species and natural resources throughout the world (Chape, Spalding, & Jenkins, 2008). Accordingly, the 196 parties to the Convention on Biological Diversity (CBD) have adopted Aichi Target 11, committing governments to conserve ≥17% of terrestrial and ≥10% of marine areas through site-based conservation strategies by 2020 (CBD, 2010). Additionally, the integration of PAs into the wider landscape through broadly connected networks and ensuring that they are representative of important sites for biodiversity is a targeted outcome of the expansion of PA estates. Preparations to establish a post-2020 biodiversity framework are underway, and even more ambitious targets for protected area coverage, connectivity and representativeness might be forthcoming (CBD, 2018). To meet these objectives, many CBD member parties might look to expand their PA estates in the coming years (UNEP-WCMC, IUCN, & NGS, 2018).
Therefore, identification of the factors that lead to successful area-based conservation outcomes and understanding the interaction of PAs with the wider landscape will be crucial as more land becomes designated for conservation purposes.
Several studies have shown that PAs are an effective means of conservation of biodiversity (Bruner, Gullison, Rice, & Da Fonseca, 2001;Coetzee, Gaston, & Chown, 2014;Geldmann et al., 2013;Gray et al., 2016;Nagendra, 2008;Nelson & Chomitz, 2009). However, they can also result in unintended consequences that undermine conservation efforts (Bode, Tulloch, Mills, Venter, & Ando, 2015;Pfeifer et al., 2012;Renwick, Bode, & Venter, 2015;Visconti et al., 2019). One such consequence is leakage, which occurs when land uses harmful to conservation efforts (e.g., deforestation) are displaced to areas beyond the administrative boundaries of the PA (Ewers & Rodrigues, 2008). When this displacement occurs directly adjacent to the PA, it can alter the composition and structure of vegetation in the PA buffer zone, potentially limiting the ranges and dispersal capabilities of organisms living within the PA and disrupting other ecological functions (DeFries, Karanth, & Pareeth, 2010;Haddad et al., 2015;Hansen & DeFries, 2007). The severity of the impacts of leakage on biodiversity will depend on the biodiversity value (measured, for example, as irreplaceability; Pressey et al. 1993) of the buffer zone compared with the PA. Cases where irreplaceability is lower in the buffer zone are consistent with positive net biodiversity impacts of PAs, but such impact will be minimal where the irreplaceability in the buffer zone and the PA is the same and negative if the irreplaceability of the buffer zone is greater than that of the PA.
Leakage in situations where buffer zone irreplaceability is high can reduce species richness and ecosystem function in landscapes surrounding PAs, with feedbacks to species and ecosystem function within PAs. Furthermore, leakage patterns can undermine assessments on the effectiveness of PAs, because success is often determined by comparing ecosystem health indicators between the PA and its immediate surroundings (Ewers & Rodrigues, 2008;Joppa & Pfaff, 2010). The importance of PA buffer zones is recognized by the inclusion of Article 8e in the CBD, which states that parties to the convention should "Promote environmentally sound and sustainable development in areas adjacent to protected areas with a view to furthering protection of these areas" (CBD, 1992). Buffer zones also are recognized by Aichi Target 11, which highlights the importance of integrating site-based conservation measures into the wider landscape (CBD, 2010).
Despite the recognition that areas adjacent to protected areas affect the success of site-based conservation efforts, the global prevalence of leakage remains poorly understood, and few studies have attempted to quantify this phenomenon (Fuller, Ondei, Brook, & Buettel, 2019). Cases of deforestation leakage have been reported in forests of Peru, East Africa and Indonesia (Oliveira et al., 2007;Pfeifer et al., 2012;Poor, Frimpong, Imron, & Kelly, 2019). However, two studies conducted across pan-tropical selections of PAs did not find evidence of widespread deforestation leakage (Fuller et al., 2019;Lui & Coomes, 2016). These studies also offer little insight into the drivers of deforestation in cases where leakage is present, thereby precluding efforts at mitigation.
One explanation for the lack of evidence of deforestation leakage reported in these studies might be the method by which counterfactual deforestation scenarios were established (Fuller et al., 2019). To identify deforestation leakage, the rate of deforestation must be compared among a PA, its buffer zone and an unprotected control (i.e., counterfactual scenario; Ewers & Rodrigues, 2008).
The locations of protected areas are non-random. Therefore, any effort to assess the effect of protection on deforestation rates should account for this potential bias, because factors such as elevation and distance to population centres might affect deforestation rates (Joppa & Pfaff, 2009 Rodrigues, 2008;Joppa & Pfaff, 2010). This study advances previous efforts by using propensity score matching to quantify and analyse deforestation leakage surrounding tropical and subtropical PAs.
The objective of this study was to identify cases of deforestation leakage in a pan-tropical and subtropical selection of 120 PAs by comparing deforestation rates among PAs, unprotected PA buffer zones and unprotected control areas. We also sought to assess the impact of deforestation leakage from PAs across extensive areas and to assess the main drivers of deforestation in buffer zones of PAs with deforestation leakage.  (Olson et al., 2001), were not fully nested within other PAs and for which boundaries were available (i.e., polygons). If the edges of two PAs in the selection were within 20 km, we discarded the selection and made a new selection. This ensured that the 10 km buffer zones of PAs in the selection would not overlap, maintaining the spatial independence of the sample.
We assumed that a pattern of deforestation leakage would displace an equivalent amount of deforestation from the PA to the buffer zone. Accordingly, we sought to analyse deforestation leakage in buffer zones with areas relatively equal to those of their respective PAs. Previous studies analysed leakage within 10 km buffer zones (Fuller et al., 2019;Lui & Coomes, 2016;Poor et al., 2019). However, this buffer radius produces buffer zones with areas similar to that of the PA only when the PA is relatively large (c. 2,000 km 2 ). Given that most PAs within our sample were considerably smaller than 2,000 km 2 (Supporting Information Table S1), we chose to analyse leakage patterns at four buffer zone levels: 1, 3, 5 and 10 km. We created concentric buffer zones around PAs at these distances and clipped to the extent of tropical and subtropical forest biomes. We then erased areas of buffer zones that overlapped with unselected PAs in the WDPA. This ensured that all buffer zones occurred fully within tropical or subtropical forests and were unprotected.
To identify deforestation leakage, the deforestation rate must be compared between a PA and a control that approximates the area if not protected (i.e., counterfactual scenario; Ewers & Rodrigues, 2008). To create the counterfactual scenario, we overlaid tropical and subtropical forest areas with a 10 km × 10 km grid, with each grid cell representing a potential control area. We then erased cells that overlapped protected areas within the WDPA or 10 km buffer zones of the selected PAs to ensure that control areas were unprotected and were outside the influence of PAs in the sample. We used propensity score matching (PSM) to match PAs with similar but unprotected cells, thereby controlling for covariates in the comparison of deforestation rates (Joppa & Pfaff, 2010). We conducted PSM with the PS Matching extension in SPSS v.25 (IBM, 2019). We performed this analysis separately for Africa, America and Asia. We used distance to the nearest major road, distance to the nearest major population centre, mean elevation, country and terrestrial ecoregion as covariates for PSM, with exact matching on country and terrestrial ecoregion (EEA, 2016;Environmental Systems Research Institute, 2018;Natural Earth, 2018;Olson et al., 2001). We chose these covariates because they can affect rates of deforestation and have been used in many other statistical matching studies (Andam, Ferraro, Pfaff, Sanchez-Azofeifa, & Robalino, 2008;Mas, 2005;Poor et al., 2019). We used the PS Matching extension to establish callipers that defined a tolerance outside of which matches were excluded from the sample (D'Agostino, 1998). In this manner, we improved the covariate balance by eliminating matches that were relatively dissimilar (Supporting Information Table S3). Each PA was matched with up to five control areas, and for each region we used the minimum callipers that allowed each PA to match at least one control area. This ensured that each PA matched and that the best overall balance was achieved (Andam et al., 2008;D'Agostino, 1998). We set the callipers at 0.19, 0.25 and 0.15 standard deviations of the logit of the propensity scores for America, Africa and Asia, respectively.

| Deforestation estimates and drivers of leakage patterns
We used the Global Forest Change Dataset (GFCD) 2000-2017 (Hansen et al., 2013) to calculate the average yearly deforestation rate for PAs, buffer zones (1, 3, 5 and 10 km) and control areas.
Accordingly, we calculated the average yearly deforestation rate as the area deforested per area of PA, buffer zone or control per year from 2001 to 2017. The GFCD shows tree cover loss at a 30 m × 30 m resolution, with the pixel value corresponding to the first year of detected tree-cover loss. Given that each pixel can contain only one value, the data do not allow for the detection of multiple tree-cover loss events within the same pixel during the time period. We used the GFCD data mask (Hansen et al., 2013) to differentiate terrestrial surfaces from water bodies and considered only terrestrial areas for the calculation of average yearly deforestation rate (Lui & Coomes, 2016). We defined leakage as a higher average yearly deforestation rate in a buffer zone compared with the PA and its matched control area (Ewers & Rodrigues, 2008). To determine whether deforestation rates differed significantly (p < .05) among PAs, buffer zones and controls, we performed Kruskal-Wallis and Wilcoxon-Mann-Whitney post hoc tests. We defined significant leakage as an average yearly deforestation rate significantly (p < .05) higher in the buffer zone than the PA and control. We assessed potential drivers of leakage patterns with the Drivers of Global Forest Loss (DGFL) dataset, which classifies the main driver of forest loss as commodity production, forestry, shifting agriculture, wildfire or urbanization, at a 10 km resolution with one value per pixel (Curtis, Slay, Harris, Tyukavina, & Hansen, 2018). We determined the prevalent driver of forest loss by counting the number of deforestation pixels within each driver category within buffer zones of PAs that had leakage.

| Landscape-level effects of leakage
To assess the impact of leakage, we applied the framework of Ewers and Rodrigues (2008). We quantified the area of land deforested across the PA and its 10 km buffer zone over the given time period, producing a value of deforestation across the "conservation landscape". We compared this value with a counterfactual scenario, in which we calculated the area expected to be deforested in an unprotected landscape of the same size over the same period by extrapolating the area of deforestation in the corresponding control to the size of the combined PA and 10 km buffer zone. By comparing the actual area of deforestation in the conservation landscape with the expected area of deforestation for each case, we classified the pattern as avoided deforestation (actual deforestation < expected deforestation) or enhanced deforestation (actual deforestation > expected deforestation).
To assess further the potential impacts of leakage on landscape-level conservation efforts, we calculated and compared irreplaceability between PAs and their respective 10 km buffer zones. The irreplaceability metric represents the importance of a site for meeting conservation objectives for a given species or set of species (Pressey et al., 1993). A high value of irreplaceability indicates that if a site is lost, the likelihood that another site will replace its conservation value for the given species is low.
We calculated irreplaceability following the methods of Le Saout et al. (2013). This method compares the extent of a species' range within the site of interest with its total range, with a higher value of irreplaceability signifying that a site contains a higher percentage of the species' total range. The irreplaceability metrics corresponding to each species that occupies a given site are summed to calculate an overall irreplaceability value for the site. In this study, we considered irreplaceability for threatened species (i.e.,

IUCN Red List categories Critically Endangered, Endangered and
Vulnerable; IUCN, 2012) of birds, terrestrial mammals and amphibians. These taxonomic groups have been assessed comprehensively with IUCN Red List criteria. We acquired range maps from the IUCN Red List for terrestrial mammals and amphibians (IUCN, 2019) and from BirdLife International for birds (BirdLife, 2018). The spatial accuracy of these range maps varies with the extent of the range of each species. Therefore, false-positive errors are inevitable (Rodrigues, 2011) but tend to be smaller for the highly irreplaceable species (i.e., those with a small range), which drive the results.

| Deforestation patterns
The random selection of PAs across the pan-tropical and subtropical region created a sample of 120 PAs from the 6,679 PAs that met our criteria. These selected PAs were located in 43 countries, encompassed all six IUCN management categories and had areas ranging from 2 to 14,034 km 2 (Figure 1). Propensity score matching techniques achieved good balance for all covariates (standard mean difference < 0.25; Figure 2), indicating that the selection of controls to establish the counterfactual scenario was robust (Thoemmes, 2012).
Of the 120 PAs analysed, there were 55 in which the average yearly rate of deforestation was higher in buffer zones than within the PA and its unprotected matched controls ( Table 1). The prevalence of leakage was similar for the three regions considered, with leakage in 19, 18 and 18 PAs in America, Africa and Asia, respectively. Of these 55 cases, 18 were significant. The leakage patterns were significant for five, six and seven PAs in America, Africa and Asia, respectively.
Leakage was not clearly correlated with IUCN management category, mean elevation, size, age, distance to nearest major road or distance to nearest major population centre (Figure 3).
Across the selected group, average (±SE) yearly rate of deforestation was significantly lower within PAs (0.25 ± 0.04%) than in all buffer-zone levels (1 km, 0.39 ± 0.05%; 3 km, 0.42 ± 0.05%; 5 km, 0.44 ± 0.06%; and 10 km, 0.43 ± 0.06%) and unprotected controls (0.44 ± 0.05%). The average yearly rate of deforestation was higher in unprotected controls than in 1, 3 and 10 km buffer zones, but equal to that in 5 km buffer zones. Protected areas in America had higher average yearly rates of deforestation at all buffer zone levels than unprotected controls, whereas PAs in Africa and Asia had lower average yearly rates of deforestation at all buffer zone levels than unprotected controls. The average yearly rate of deforestation within the PA was significantly higher than within unprotected

| Landscape-scale impacts of leakage
Of the 55 PAs with leakage, 43 (78.2%) resulted in a higher level of deforestation across the conservation landscape compared with the expected level of deforestation without protection, indicating enhanced deforestation (Figure 4). This includes 26 cases in which leakage was present but not significant and 17 cases in which leakage was significant. Protected areas without leakage resulted in less deforestation across the conservation landscape than would be expected without protection in 48 cases (73.8%). Overall, the existence of PAs avoided deforestation in 65 cases (54.2%) and enhanced levels of deforestation in 55 cases (45.8%). This pattern was consistent across all regions, with 21 (52%), 23 (57%) and 21(52%) PAs that avoided deforestation in America, Africa and Asia, respectively.
We calculated irreplaceability for threatened species of terrestrial mammals, amphibians and birds for the 55 PAs exhibiting leakage and their respective 10 km buffer zones (Supporting Information Table S2).
Five PAs, three in Africa and two in America, had a higher value of irreplaceability than their 10 km buffer zone (Figure 4). The difference in irreplaceability between buffer zone and PA was >.01 in four of these five cases. Of the 50 cases in which irreplaceability was greater in the buffer zone, the difference in irreplaceability was ≥.01 in six cases.
The dominant driver of buffer zone forest loss in all 18 African PAs was shifting agriculture. In America, the dominant drivers were shifting agriculture (n = 13), commodity production (n = 4) and forestry (n = 2). In Asia, forestry was the most prevalent driver of buffer zone forest loss (n = 7), followed by commodity production (n = 6) and shifting agriculture (n = 5).

| Global prevalence of deforestation leakage
The results of this study suggest a higher prevalence (46% of PAs

| Landscape-level effects of leakage
In 12 of 55 documented cases of leakage, PAs curb deforestation within their borders to avoid deforestation across the conservation landscape.
Although PAs generally are not established with the goal of reducing harmful land uses in areas outside their jurisdiction, overall conservation success might still depend on activities beyond PA borders (DeFries et al., 2010;Hansen & DeFries, 2007;Laurance et al., 2012). This might be the case especially where PAs are small, isolated or inhabited by species with area and resource requirements that are not met fully by the PA itself (Maiorano, Falcucci, & Boitani, 2008;Woodroffe & Ginsberg, 1998). It might be necessary to give particular attention to PAs with these char-

| Drivers of deforestation leakage
Factors we hypothesized to predict deforestation leakage, such as not correlated with leakage ( Figure 3). This finding is consistent with the patterns reported by Fuller et al. (2019), further suggesting that the factors driving deforestation leakage are case specific and might be controlled by local political, social or economic factors (Lui & Coomes, 2016). Shifting agriculture was the most prevalent driver of forest loss in buffer zones of PAs with deforestation leakage. However, the DGFL data have relatively low model accuracy in differentiating between shifting agriculture and commodity production in Africa (Curtis et al., 2018). Accordingly, the drivers of leakage for several African PAs might be misclassified. Fifteen of the 18 African PAs with leakage are within F I G U R E 3 Proportional distribution of leakage and non-leakage cases across the following variables: International Union for Conservation of Nature management category, elevation (in metres above sea level), area of protected area (in square kilometres), year of establishment of protected area, distance to major population centre (in kilometres) and distance to major road (in kilometres). Numbers inside bars represent the number of cases per category landscapes in which shifting agriculture in 2010 was "low" or "very low", whereas the remaining three are within landscapes where shifting agriculture never existed or disappeared long before 2010 (Heinimann et al., 2017). Notwithstanding, commodity production (i.e., permanent conversion of forest to uses such as agriculture, mining or energy infrastructure) was responsible for eight of the nine highest rates of buffer zone deforestation among leakage cases. Although our sample size per country was relatively low, levels of deforestation leakage driven by commodity production were higher in Brazil, Indonesia and Malaysia than in other countries. This corresponds to the general deforestation trends within these countries given that their production economies have largely been driven by deforestation for the expansion of cattle and soy production in Brazil (Alix- Garcia & Gibbs, 2017;Ometto, Aguiar, & Martinelli, 2011) and oil palm in Indonesia and Malaysia (Wicke, Sikkema, Dornburg, & Faaij, 2011

| Implications for site-based conservation
Our results suggest a higher prevalence of deforestation leakage in tropical forest PAs compared with previous studies. However, it is also apparent that the majority of PAs curb deforestation within their jurisdiction, because 58% of the cases indicated that PAs significantly reduce deforestation within their boundaries compared with what would be expected without protection (Supporting Information Table S1). Additionally, in 40% of the PAs, deforestation was less than would be expected without protection across the conservation landscape, suggesting a blocking effect whereby the buffer zone receives a positive spillover of protection from the PA (Garcia, 2015). This phenomenon has been documented previously (Fuller et al., 2019;Gaveau et al., 2009;Lui & Coomes, 2016).
Additionally, PAs considered "ineffective" within our sample are not necessarily ineffective at achieving their overarching conservation goal. Protected areas are designated for myriad reasons, and some PA designations, such as sustainable forest reserves, can have a rate of deforestation that is higher than the surrounding landscape, but sustainable for the given objective.

F I G U R E 4
Classification of 55 protected areas exhibiting leakage on the basis of avoided deforestation across the conservation landscape and comparison of irreplaceability of protected areas and respective 10 km buffer zones for threatened species of amphibians, birds and terrestrial mammals. A positive avoided deforestation value indicates actual deforestation < expected deforestation across the protected area and 10 km buffer zone from 2001 to 2017. A positive irreplaceability Δ indicates irreplaceability of protected area for threatened species > irreplaceability of 10 km buffer zone for threatened species Given that PAs are likely to remain an integral conservation strategy, we recommend policy measures for buffer zone protection where leakage is deemed both high and likely to cause a considerable impact on the conservation outcome. Legal mechanisms to limit development in PA buffer zones have been established in some contexts, such as the Natural Protected Areas network in Peru (Solano, 2010)

and the Baekdu Daegan Mountain System in South
Korea (Miller & Kim, 2010). In other cases, voluntary agreements between PA management authorities and private stakeholders promote activities conducive to conservation within buffer zones (Dudley, 2008). Assessing the efficacy of existing buffer zone protection strategies to provide better guidance for the creation of effective legal frameworks that maximize compliance for buffer zone protection in various contexts might be a point of future study.

ACK N OWLED G M ENTS
This research is part of the Inspire4Nature Innovative Training Network, funded by the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement no. 766417.