Introduction

Carbon offsets have become a crucial strategy for achieving carbon neutrality and addressing climate change1,2. Such offsets are usually included in carbon markets3. Among the offset programs, forest offsets are receiving particular attention as a cost-effective carbon sink and nature-based climate solution3,4,5. Forest-based offset projects regulate land use in a targeted region. Credits are generated from forest-related activities, including management, reduced deforestation, afforestation, and reafforestation6. Most studies of the effectiveness of offset policies have focused on the impacts within the treated area, with concerns about accurate measurement7,8, the assessment of additionality9,10,11, and the assurance of permanence12,13. However, spillover effects outside the boundaries of a protected area have received less consideration. This research gap may lead to biased or incorrect evaluation of project efforts14,15. Leakage (negative spillover), referring to shifts in forest loss to areas that are unprotected, might undermine the value of offsets14,16. By contrast, a positive spillover of forest gains outside the boundaries of the protected area would amplify the effects of offset policies17.

The direction of spillovers varies across project regions and types15. Leakage is a common concern for many kinds of land use restrictions. The result may be underperformance of the project’s objectives. For forest offsets, leakage can occur through two mechanisms18: (1) forest harvest activities may shift from within the project boundaries to outside areas, typically at a local scale; (2) market demand may shift due to changes in commodity prices, land prices, and the relative competitiveness of different regions, leading to the risk of global leakage19,20. Even though some attention has been paid to potential local leakage when evaluating forest offset projects, it has been only a secondary research objective; even that limited research has focused on leakage in tropical forests and reducing emissions from deforestation and degradation projects21,22,23,24, especially in the Amazon14,21,25,26. Research on the subtropical and temperate regions, as well as other types of forest offsets, is still limited.

In contrast to the leakage problem, forest offsets may yield extra benefits through the following channels15: (1) learning spillovers that transmit knowledge about more sustainable practices and technologies27,28,29, (2) social norms spillovers, such as comparisons of activities among neighboring communities26,30, and (3) the inherent ecological-physical links between protected areas and adjacent land27,31. Positive spillovers have been less studied compared to leakage. Moreover, the existing literature usually studies either positive or negative spillovers, lacking a comprehensive analysis and consideration of diverse channels.

Spillovers could play an important role in developing methodologies for forest offset programs, as they may impact the effectiveness and sustainability of these programs. A spillover aligned with the program’s objectives could enhance the overall project benefits, while a negative spillover undermining these goals can cause over-crediting, necessitating explicit mitigation measures in the methodologies. However, the evaluation methodologies regarding spillover of forest-based offsets vary widely across programs and regions. This is particularly evident in China, where spillovers have been neglected. In contrast to the stringent requirements for leakage accounting in California’s forest offset projects32, almost all of China’s projects report zero leakage estimates in the project design documents (PDDs), and claim to have positive spillover effects. In addition, China’s unique context may influence spillovers; this context includes its strict land use system, the diverse climate pattern, and the features of forest offset projects. However, to our best knowledge, the current literature lacks identification strategies and empirical evidence concerning spillovers in China, and the claims of these projects remain unproven. Given these research gaps, this study focuses on the spillovers of China’s forest offset projects.

As the world’s largest emitter, China faces substantial pressure for carbon offsets to achieve carbon neutrality. Its abundant forestry carbon sequestration resources are expected to play a pivotal role33,34,35,36. In 2012, China launched its voluntary emission reduction market, the Chinese Certified Emissions Reduction (CCER) program, inviting project applications upon announcement. Project developers were required to submit detailed PDDs and materials specifying eligible project boundaries, planned or ongoing activities, and demonstrating project additionality. By 2017, when authorities suspended new CCER project applications, the official website had listed 100 forestry projects. In 2023, the CCER program was re-launched with updated measures, guidelines, and methodologies. Although the program was temporarily halted, our results suggest that forest managers continued to maintain and improve their forests, anticipating future offset credits.

We analyzed 36 out of 100 projects, for which we collected detailed ___location data (Fig. 1a). We replicated the project areas and created 5-km buffer zones (Fig. 1b). Utilizing a grid-based sampling approach at a resolution of 30 m × 30 m and employing statistical matching, we constructed reference areas with characteristics similar to the project sites, representing the counterfactual outcome for forest status in the absence of a project (Fig. 1d). We define a positive program outcome as forest gain, also described as improved vegetative quality. This outcome is measured as the annual peak forest greenness, that is, the normalized difference vegetation index (NDVI), spanning 2000–2022. To assess the impact of CCER programs and the resulting local spillover effects, we employed difference-in-difference (DID) regressions to compare the average changes in NDVI associated with the buffer zones of the project and reference areas before and after the intervention (see details in Methods).

Fig. 1: Study area and research design.
figure 1

a The national distribution of forest offset projects across China overlaid on the 2005 NDVI map. Red circles represent project locations with documented vertex coordinates (our study sample), while blue circles indicate projects without specified coordinates. Note that a single county may have multiple project applications, and a project may span multiple counties. b An example of one specific project site (blue) and its 5 km buffer zones (light blue). c A 30 m × 30 m grid (gray) on the project area, precisely delineating its boundary. d The matched control counties (dark green) within the same province as the project locations (purple). e Our grid-based matching results (see Methods), with control pixels sampled from the control counties, excluding those within established protected areas. The matched pixels are shown in black.

Full size image

Our research enriches the understanding of the role of forestry carbon offsets, and increases their credibility and reliability in mitigating climate change, particularly in subtropical and temperate regions. We contribute to the existing literature by evaluating the impacts of the voluntary carbon market, evaluating spillovers from multiple directions and channels, developing appropriate methods, and assessing the rationality of offset methodologies. We also fill the research gap on China’s offset projects, provide a reference for evaluating and managing potential spillovers in a broader context, and add encouraging empirical evidence on the effectiveness of forestry projects from a micro-level perspective.

Results

Before assessing the impacts of intervention, we explored whether the constructed reference area can accurately represent the trends of forestry before the launch of China’s voluntary carbon market. The collective trend of the projects, along with their buffers shown in Figs. 2 and 3, exhibits nearly identical pre-treatment patterns. This suggests that the reference area can serve as a counterfactual scenario with similar socioeconomics and geographical characteristics.

Fig. 2: Forestry activities in project sites.
figure 2

a The annual mean trends of NDVI in project areas and reference areas, with the red vertical line representing the year 2012. b The annual site-based NDVI differences between each of the 36 project sites and their corresponding reference areas, as represented by blue dots for each year. The red line across the plot represents the mean of these differences across all project. c The average treatment effect on the treated for individual projects estimated by difference-in-difference regression, where the whiskers extending above and below each point indicate the 95% confidence interval.

Full size image
Fig. 3: Forestry activities in buffers.
figure 3

a The annual mean trends of NDVI in the buffer zones of the project and reference areas, with the red vertical line representing the year 2012. b The annual site-based NDVI differences between each of the 36 project sites and their corresponding reference areas, as represented by blue dots for each year. The red line across the plot represents the mean of these differences across all projects. c The average treatment effect on the treated for individual projects’ buffer zones estimated by difference-in-difference regression, where the whiskers extending above and below each point indicate the 95% confidence interval.

Full size image

Forest gain in the offset project sites

A prerequisite for generating spillovers to nearby areas is that the projects must have initial impacts within their boundaries. Before interpreting the results regarding spillovers, we conducted average and project-specific tests to check the effectiveness of forest offset projects in achieving program goals. As Fig. 2 shows, the trends of the project and the reference area became divergent after 2012, which is the post-treatment period. This divergence continued to grow over time, revealing overall positive effects of the CCER program. The baseline results of overall project sites show a positive effect of 2.25% to 4.25% with p-value less than 0.01 (Supplementary Table S1). This provides some evidence that the policy intervention has incentivized project candidates to engage in additional forestry activities.

Furthermore, we conducted a site-based examination to analyze individual projects (Fig. 2c; Supplementary Tables S2, S4). Among all the available projects, 29 out of 36 project sites demonstrated statistically significant positive DID estimators, indicating positive program impacts in terms of forest gain. The majority of the estimated effects ranged between 1% and 5%, with the highest effect reaching 37%. One site showed statistically non-significant effects and the remaining five projects (in gray) had negative DID estimators at the 1% level, most of whose nearby areas also exhibiting a negative trend (Fig. 3c). The patterns observed in these negative projects suggest either the absence of forest gains or the influence of other systematic factors, including climate-related damage within a broader region. Therefore, these five significant negative projects are deemed as non-effective samples due to the lack of evidence of additional forestry activities inside the project boundaries.

The spillover effects in nearby sites

Generally, the forest status in all buffer zones follows a similar trend to that of the project sites (Supplementary Fig. S2), indicating regional ecological and physical links. As the buffer zones extend farther from the project areas, the proportion of forest types decreases and NDVI gradually declines. Due to spillover effects on surrounding areas, either forestry development or carbon leakage, trends within and beyond project boundaries are changing. Therefore, directly comparing trends between project sites and their buffer zones is insufficient to explore the treatment effects. To investigate the potential spillovers more thoroughly, we performed a matched DID analysis on the 5-km buffer zones using a grid-based sampling approach.

The full sample trends indicate that the pixels in buffer zones and their control groups exhibited consistent patterns prior to 2012, with notable changes emerging after 2012 (Fig. 3a). Over 23 years, the average forest gain spillover from the project sites to neighboring areas ranges from 0.91% to 1.60% (for NDVI and kNDVI respectively, p < 0.01, Supplementary Table S1). This effect becomes increasingly positive in the post-treatment period, reflecting strengthened forestry efforts and cumulative impacts as the CCER policy was implemented.

However, our project-specific analyses revealed variations in spillover effects. Individual site-based estimations (Fig. 3d; Supplementary Tables S3, S5) showed that, out of 36 projects, regardless the five non-effective projects, 10 projects induced negative impacts to the buffer zones, ranging from −10.27% to −0.06% with p < 0.01; conversely, 20 projects significantly induced positive impacts, ranging from 0.09% to 22.34%; while one project showed no significant effect. This highlights the complex nature of spillovers, which can vary across a wide range of settings, including the geographical ___location, the ecological setting, and the project features. Therefore, the existence, magnitudes, and even directions of spillovers should be assessed within specific settings.

The results are robust to alternative model specifications (Supplementary Table S6). First, we substituted the control group (the outcomes of 1:1 nearest neighbor matching) with the top three and top five nearest neighbors, recognizing that matching outcomes can affect the estimation. Second, to account for potential confounding factors arising from variations in area size, we modified our grid processing approach to adjust the number of observations of each buffer zone associated with their sizes (see Methods). Third, to address the concern about potential spatial autocorrelation in our geospatial data, we implemented spatial statistical methods and used spatially robust standard errors developed by Conley37,38. Fourth, to account for discrepancies between the CCER launch year and the project application submission time, we substituted the intervention time with the public comment start date for each project, conducting the DID regression within a staggered framework39. Lastly, to address potential temporal confounding during our long study period, particularly in light of China’s intensified forestation efforts post-2000 and the initiation of several large-scale tree plantation programs in the early 2000s40, we conducted additional estimations using shorter intervals 2007–2017 and 2009–2015.

The results are also robust to alternative indices and dataset. First, to address the non-linear relationship and saturation effect of NDVI in areas with high vegetation cover and during peak seasons, we employed kNDVI41 as an alternative vegetation index. Results remained consistent, with kNDVI matching results (Supplementary Fig. S3) aligning with trends in Figs. 2a and 3a. The estimated effects were slightly larger in absolute value (Supplementary Tables S1–S3), possibly due to NDVI limitations in dense forest areas. Second, we further validated our results using NDVI data estimated from different sources. Due to computational constraints, we used 1 km × 1 km resolution data from another source for matching and baseline regression. While the conclusions remained consistent (Supplementary Fig. S4; Supplementary Table S6, columns 8-9), the higher-resolution data used in our primary analysis yielded approximately 5 times larger effect. This is attributed to the fact that most of the projects consist of numerous small-scale polygons, with nearly 45% of the 1 km2 grids covering project sites contained less than 10% project area, resulting in a coarser representation of green biomass (see comparison in Supplementary Fig. S5). Third, to address the potential influence of the large sample size on the statistical significance of the results, we aggregated the treatment and control pixels in both project sites and the buffers to conduct baseline estimation on this smaller sample. The project-level results remain significantly positive and show only minor differences from our pixel-based estimates. These findings underscore both the robustness of our conclusions and the importance of data resolution. Additionally, the results are robust to the presence of other forest restoration and environmental programs.

Furthermore, we performed two types of placebo tests by randomly aligning the treatment pixels and intervention period using a bootstrapping procedure (Supplementary Fig. S6). We first created pseudo-treatment groups by randomly selecting samples from the control group, maintaining the original intervention period. Then, we randomly selected both treatment groups and intervention periods between 2000 and 2022. Both tests were repeated 1000 times. These placebo tests simulated scenarios where either no forestry carbon offset projects occurred in the project areas, or other policy interventions happened in different years. The results show that our actual DID estimates lie at the extreme right tail of both the empirical distributions of the placebo estimates and the associated theoretical normal distribution. This implies a positive effect that is statistically different from zero, signifying that our observed treatment effects are unlikely to be the results of random factors and instead represent genuine policy intervention effects.

Heterogeneity in spillovers

The spatial distribution of estimated spillovers for the 36 sample projects is presented in Fig. 4a, facilitating a comparison of spillover effects across different spatial scales. The density patterns reveal that as buffer distance increases, the distribution of estimators becomes more concentrated around zero. The average effects diminish, with a notable drop between the 1 km and 2 km buffers from 0.85% to 0.15%, followed by a gradual reduction to approximately zero at greater distances. This observed trend suggests a spatial gradient of forest project impacts on vegetation quality beyond project boundaries, with the primary impact occurring within the 1 km buffer zone. The pattern also indicates that the 5-km buffer zone serves as an appropriate analytical boundary for capturing spillover effects.

Fig. 4: Heterogeneity test.
figure 4

a The spatial pattern of estimated spillovers in the buffer zones. The main plot depicts the density distribution of DID estimators for buffer zones at distances of 1 km, 1-2 km, 2-3 km, 3-4 km, and 4-5 km from the project boundaries. The inset in the upper left corner displays the average effects across these buffer distances, with each bar representing the mean estimator value for a specific buffer interval. b A heterogeneity analysis based on three key features: economic status, primary ecosystem type, and project type; the whiskers extending above and below each point indicate the 95% confidence intervals for the respective estimates.

Full size image

We further investigated how project features impact outcomes by categorizing projects based on county economic development level (grouped by per capita GDP), project type (forest management = 6, afforestation = 25, and vegetation restoration = 5), and primary ecosystem type (forest = 26, grassland = 3, cultivated land = 6, desert = 1) (Fig. 4b). The positive spillovers are larger in less-developed regions, vegetation restoration projects, and farmland-predominant ecosystems.

First, regarding the economic status, the potential economic benefits and new technologies from the projects are appealing to economically disadvantaged communities42,43, leading to larger positive spillovers in less-developed regions.

Second, we discuss the heterogeneity with respect to project type. The ecosystem services provided by afforestation and vegetation restoration programs have been extensively studied44,45,46. These projects, including cropland-to-forest conversion, mine rehabilitation, and reduced deforestation, have been shown to influence local climate patterns47, mitigate soil erosion and land degradation48,49, and enhance ecosystem resilience. Such effects likely contribute to more substantial spillovers in surrounding areas.

Third, we distinguish the primary ecosystem type within the buffer zones. The substantial spillovers observed in farmland ecosystems align with the mechanism tests on land-use changes (see Table 1), where cultivated land experienced an overall decrease while grassland and forest cover increased within the 5-km buffers. This suggests that the spillovers may be partly attributed to the promotion of farmland-to-forest conversion, potentially driven by farmers’ re-evaluation of revenue from forestry activities and offset projects when these are perceived as more profitable50,51.

Table 1 The channels of the forestry spillovers
Full size table

Conversely, a significant decrease in vegetation quality was found in the buffer zone primarily composed of desert ecosystem. Satellite imagery reveals predominant desert conditions both within and beyond the project sites (see the subplots in Supplementary Fig. S7). While local water resources are scarce, the PDD states that the project site is irrigated by rivers formed from mountain glacier melt and spring water. A plausible explanation is that increased water consumption by the afforestation project may have exacerbated water scarcity in downstream and surrounding areas. This is supported by further examination revealing a notable decrease in water surface area from 3.5% in 2010 to 1.9% in 2020, suggesting intensified drought conditions and consequently decreased vegetation quality in surrounding areas52. These findings highlight the need to consider the complex interactions between forestry activities and local ecological contexts when planning and implementing related projects.

The channels of the forestry spillovers

Overall, our analysis does not provide evidence of substantial leakage from forest offset projects in China. To investigate potential leakage, we examined land use dynamics in buffer zones of projects exhibiting negative spillovers, including non-effective projects. The leakage hypothesis would predict an increase in cultivated land use in buffer zones, reflecting a displacement of forest harvest activities from project sites to surrounding areas. However, our results contradict this expectation: a mean decrease of 1.8% is observed in cultivated land use in 2015 and 2020 compared to the 2000–2010 average, accompanied by increases in grassland and forest cover of 2.4% and 0.2%, respectively (Table 1, columns 1–3). These land use changes align closely with those observed in the full sample (Supplementary Table S7). Further investigation revealed that projects with negative spillovers were primarily implemented on degraded and eroded soils. This suggests an alternative explanation for the observed negative effects: inappropriate afforestation measures may have exacerbated local ecological stresses. For instance, certain afforestation practices could potentially degrade local soil conditions, leading to increased erosion and subsequent decline in vegetation quality53,54.

Furthermore, we examined the trends in forest loss within buffer zones during the period 2001–2023, utilizing data from Global Forest Change dataset55. We used the annual percentage of forest loss relative to the 2000 baseline forest cover in buffer zones as our dependent variable (mean: 0.39%). The DID estimation revealed that in the full sample, forest loss decreased by 0.1% (about 25% reduction) following policy implementation. Projects with negative spillovers showed a non-significant 0.04% decrease, while those with positive spillovers exhibited a significant 0.2% reduction in forest loss (Table 1, columns 4–6). These findings indicate an absence of additional forest loss or carbon leakage in the buffer zones, while supporting the positive spillover effects of forest offset projects in reducing forest loss beyond their boundaries.

Our study suggests that learning could be one of the drivers of positive spillovers in forest offset projects. This learning effect operates through knowledge dissemination and behavioral learning15,27. Forest-friendly technologies and conservation practices spread through neighboring communities and social networks, which could reduce forest diseases, pests, and wildfires beyond the project sites. Our analysis of wildfire frequency, using MODIS data from NASA spanning from 2001 to 2022, demonstrates a notable reduction in wildfire occurrences. Specifically, we observed an average decrease in wildfire detection frequency of 11% in the entire buffer zones and 0.2% per square kilometer, as detailed in Table 1, columns 7–8. This reinforces the effectiveness of learning in preventing and mitigating forest-related challenges.

In addition, information flowing among neighbors enhances their understanding of the costs and benefits of forest-related activities. When they witness profitable outcomes in nearby forest offset projects, they tend to adopt similar behaviors related to forest conservation and project applications. This is supported by the fact that 60 out of 100 CCER forest offset projects received project applications in the same cities during adjacent years. Over half of the studied projects, and 50% of the projects with evidence of positive spillovers, were located in cities with successive rounds of project applications (Supplementary Table S10).

Discussion

This study offers insights into the incentives provided by China’s CCER program for forestry activities and the resulting spillover effects. Using a quasi-experimental approach with a 23-year panel dataset from 2000 to 2022, we find that forest offset projects induce overall positive forestry spillovers to their buffer zones, with learning potentially being one of the contributing channels. Within the scope of our analysis, substantial leakage is not detected.

Several reasons may contribute to these results. First, land-use regulations in China are stringent56,57,58. For example, felling permits, with strict approval procedures, are required for forest harvesting. This makes it less likely that forest harvest activities will shift to unprotected areas. Second, the scattered and limited-scale nature of the offset projects makes it challenging for any single project to exert a substantial influence on the supply and demand dynamics of global timber or crop markets18,59. Forest offset projects in China are typically characterized by small and disconnected polygons, mostly due to the land use features and strict site selection criteria. These include requirements for non-forested land since 2005 and have certified ownership rights, with the latter complicated by China’s complex forest property system (see Supplementary Note 1 for details). These factors make it difficult to yield large fluctuations in the local timber market, as relatively small polygons are being regulated. Thus, leakage through shifts in market demand is less likely to occur.

Compared to the existing literature, our findings also show a considerable magnitude. Our study reveals an annual NDVI growth rate of 0.0035 in project areas and 0.0029 in buffer zones from 2000 to 2023, significantly higher than the range of 0.0006 to 0.0027 observed for China as a whole during different periods in previous studies60,61,62,63. Using a difference-in-differences approach, we identified additional average NDVI growth of 2.25% in project areas and 0.91% in buffer zones from 2012 to 2022, after controlling for common trends. These results are not only comparable to but in some instances surpass the effects found in other conservation programs and afforestation initiatives63,64,65, indicating the significant impact of CCER forest offset projects on vegetation quality.

The study offers insights for the ongoing implementation of China’s CCER market, particularly since the process of approving applications was re-started in 2023. The results suggest that forestry offset projects, due to their potential for additional carbon removal and emissions reduction, should receive priority consideration in CCER-supported projects, thereby enhancing their role in achieving carbon neutrality. Additionally, the variation in spillover effects that we found suggests that, while leakage may not be a primary concern, it should be addressed on a case-by-case basis, considering the geo-bio-physical conditions of the project area, especially in ecologically sensitive regions. Moreover, although not the primary focus of our study, our findings of forestry gains observed from 2012 both within and beyond the project sites lend support to the 2023 revision of the administrative measures for the CCER market. This revision shifted the project eligibility criterion, changing the initiation date for eligible project activities from 2005 to 2012, aligning it with the CCER program’s implementation date66. This adjustment reflects a more targeted approach to capturing the effects of the CCER initiative.

Although not all projects demonstrate forest gain or spillover effects, this does not imply that these projects are completely ineffective. The delay in project approval might be a key factor. Many project applications were submitted after 2016. However, as noted, the Chinese government suspended project approvals in 2017, restarting the process in 2023. Therefore, many potential projects have not been approved or have not yet received carbon offset credits. Even so, the anticipation of getting credits in the future induced some forest improvement. As a result, this study may underestimate the potential of CCER projects to provide financial incentives for forestry activities.

Our study provides avenues for future research on forestry offset projects. While our focus has been on investigating spillover effects within the forest context, future research could consider an expanded scope of co-benefits of these projects. The PDDs for many of the projects suggested co-benefits such as employment opportunities, poverty alleviation, climate adaption, biodiversity conservation, water conservation, air quality, and other ecological improvements. With the relaunch of the CCER market in China, it is expected that certain forest offset projects that we studied will gain approval and start receiving credits in the near future. Future studies can continue to track the progress of these projects to assess the post-approval performance, and thus evaluate the effectiveness of carbon credits in a broader context.

We acknowledge several limitations to this study. First, while we have some evidence regarding the impact of learning effects, a more comprehensive understanding of why and how forest managers respond to policies necessitates more granular data. This could include individual behavioral data obtained through field investigations and interviews, which would also facilitate exploration of other potential mechanisms, such as social norms. Second, although our estimation of net spillovers accounts for both aggregate and individual analysis, the spatial influences can occur bidirectionally, making it challenging to empirically separate these effects in a single direction within a particular site. We believe that these issues could be addressed with detailed data and specific research designs in the future. Additionally, considering the re-launch of CCER program in 2023, further research across a wider range of forestry projects in China is necessary to draw more comprehensive insights. Despite the above limitations, this research highlights the potential of China’s CCER program, and adds evidence on both positive spillover and limited leakage of forest offset projects in subtropical and temperate regions.

Methods

Figure 5 provides an overview of our methods and variables.

Fig. 5: Flowchart of methods.
figure 5

First, the project areas were identified and calibrated using remote sensing imagery. Thus, 5-km buffer zones were created for each project. Both the project and buffer polygons were rasterized into about 30 m × 30 m pixels. Second, the reference areas were constructed following data exclusion and a two-step statistical matching process: the control groups were restricted to the same provinces without other forest offset projects, and the protected areas were excluded; county matching was conducted using county-level socioeconomic covariates; and pixel matching was conducted using pixel-level geo-biophysical covariates. These steps constructed references areas for both project sites and buffers, generating the matching pixels. Third, difference-in-differences regressions were run for the entire offset portfolio using pooled data, as well as for each individual project, to estimate the impact of CCER programs on additional forestry activities of forest offset projects.

Full size image

Study sites and offset projects

From 2012 to March 2017, a total of 100 forest-related offset projects were announced. These 100 projects were distributed across 68 cities in 23 provinces, with some projects spanning multiple cities and several projects located within a single city. We collected primary information for all 100 projects and acquired PDD files for 96 of them. Among these, 36 projects provided vertex coordinates for each forestry sub-compartment, and became the focus of our study (see details in Supplementary Table S11). The representativeness of our research sample is further analyzed by comparing the geographical distribution, socio-economic contexts, forest conditions, and project characteristics of the sample projects against those lacking detailed ___location data (Supplementary Fig. S1). This analysis indicates that despite some missing sample, the projects we studied offer a broadly representativeness across multiple dimensions rather than being subject to sample selection bias (see details in Supplementary Note 2).

We utilized the coordinate data from the PDDs to delineate the boundaries of forestry sub-compartments for the 36 projects. These boundaries were further refined using additional information from the PDDs, including project area maps and satellite images. Remote sensing imagery was also used to create a more accurate representation of the project boundaries, removing non-project areas such as rivers and residential lands, ultimately generating the final project polygons (Fig. 1). Based on the project polygons, we created concentric buffer zones at intervals of 0-1, 1-2, 2-3, 3-4, and 4-5 kilometers and grid each buffer into several 30 m × 30 m pixels for data calculation. It is worth noting that these projects typically consist of multiple small, scattered, and disconnected polygons, with very few large and independent forest areas. As a result, the 5-km buffer zones around each project are, on average, thirty times larger than the project area. This characteristic is illustrated in the lower part of the example project shown in Fig. 1b. Our high-resolution data accurately captures both the project boundaries and the scattered small polygons. In total, the 36 projects yielded approximately 2.39 million pixels (Supplementary Table S8).

Reference area construction

To identify control units for each project, we followed two steps to keep similar socioeconomic and geo-bio-physical conditions between treatment and control groups. Before matching, we conducted data exclusions and filtration. First, counties with forest projects supported by other offset programs were excluded, including those under the Clean Development Mechanism of the Kyoto Protocol and projects certified by the Verified Carbon Standard, a widely used greenhouse gas crediting program. Additionally, counties with other CCER projects were not considered in our analysis. Second, established protected areas, namely the area involving conservation, were excluded. Third, control counties were restricted to the same province as project sites to ensure a same policy environment and similar climatic conditions. This was done to mitigate potential confounding effects from other programs and policies. Forests within other forest offset projects and established protected areas are under some form of forest-friendly management and therefore do not represent counterfactual outcomes for the forest status in the absence of any intervention.

Firstly, we identified counties that were similar to the counties where the projects were implemented, based on selected conditions and statistical matching. We focused on socioeconomic covariates, including GDP, government revenue, industrial structure, and land use structure. The Mahalanobis distance67,68 was then calculated and compared for each available county for the year preceding the claimed activities year. For projects implemented in only one county, 1:2 nearest-neighbor matching with replacement was applied. For those implemented in more than one county, 1:1 matching was performed (Fig. 6a). We evaluated the quality of the matching, adjusted the relevant covariates, and iteratively optimized the control group counties (Fig. 6b). In the first few iterations, we found that some counties had been funded by other forestry carbon sequestration projects recently, indicating that our matching results were well representative. A total of 86 counties were selected into the matched control group after a series of filtration and matching, with 73 being unique. The repeated matched county mainly occurred in projects in Guangdong and Hubei provinces, where a number of cities were excluded before matching for the implementation of forest offset projects. After converting the unique matched counties into grid cells, we obtained about 350 million pixels to serve as control pixels.

Fig. 6: Spatial matching and balance assessment.
figure 6

a Delineates the distribution of project sites and their matched counties, connected by lines to emphasize the geographic pairing. b Compares the standardized mean difference of covariates between treated and control counties before and after matching. The open circles represent the standardized mean differences before matching, while filled circles (blue) indicate these differences after matching. The vertical dashed lines demarcate standardized mean differences of −0.1, 0, and 0.1, providing reference points for assessing the quality of covariate balance. The covariates include forest coverage (forest), industrial structure (1st Ind., 2nd Ind.), and economic indicators (revenue, per_GDP, GDP). c Provides an example of pinpointing individual project and matched pixels at 1 km × 1 km resolution (30 m × 30 m resolution is exhibited in Fig. 1c) where the blue pixels denote the project sites and the purple pixels denote the reference area. d Compares the standardized mean difference of covariates between treated and control pixels in buffer zones before and after matching. The covariates include land use of the pixel and the surrounding pixels (LUCC, LUCC Buffer), landform, where the proportion of each category is considered, as well as NDVI, precipitation, and solar radiation (SRAD).

Full size image

Secondly, we constructed a reference area, namely a counterfactual control area for each project and their corresponding buffer zones, based on biophysical and land-use variables, using a grid-based sampling strategy and statistical matching. This ensured equal coverage and similar features of project areas, which had multiple project types and initial land use types (non-forested, sparse-forested, cultivated, grasslands, unused, or degraded land). A nationwide 30 m × 30 m resolution grid (consistent with the NDVI base data grid) was overlaid on the project area and the matched county, using the Create Fishnet toll in ArcMap v.10.5. Thus, we performed 1:1 nearest-neighbor matching with replacement on each pixel plot from the pool of control pixels group, based on NDVI, land use, landform, precipitation, solar radiation, ecosystem, as well as land use type surrounding the project pixels (Fig. 6c, d; Supplementary Table S9). Similarly, a Mahalanobis distance metric was used to measure the similarity between treated and control pixels. This process yielding 0.6 million unique control pixels for project sites, indicating that a matched control pixel is used 4 times on average in the control group. An identical matching approach was applied to the buffer zones, with more than 70 million pixels within the 5-km buffer zones. These matched pixels were eventually synthesized into a reference area. As a result, the matching process has to be carried out individually for each of the 36 projects and their 5-km buffer zones, about 73 million pixels, which makes the process computationally demanding.

To demonstrate the effectiveness of our matching process, we compared the standardized mean difference of covariates before and after matching, a crucial step in ensuring the validity of our subsequent causal analysis. Figures 6b and d demonstrate a substantial reduction in the standardized mean differences across all covariates after matching, as evidenced by the clustering of filled blue circles closer to the zero line compared to the open circles. These confirm the successful alignment of covariates between treated and control areas, thereby strengthening the robustness of our methodological approach in constructing appropriate control groups.

For robustness checks, we expanded the matching process for project sites to include 1:3 and 1:5 nearest neighbors. Given the computational costs and the need for rigorous validation stemming from the abundance of observations in buffer zones, we applied this replacement approach to all 36 project sites. Our findings indicate that this replacement yielded nearly identical estimates.

Identification of additional forestry activities and spillovers

Quasi-experimental methods are widely used for evaluation of forest conservation. These methods can be summarized as matching, before-after control-impact (difference-in-difference), regression discontinuity, and synthetic control69,70. Because this study focuses on spillover effects, and the areas inside and outside the project boundary are correlated to some extent, regression discontinuity design is not applicable in this case71. Therefore, we use matching and difference-in-difference (DID) models to evaluate the spillover effects of forest offset projects. We also need to mitigate confounding effects that are contemporaneous with the intervention and could affect the outcome. To do so, we pooled the data and used two-way fixed effects panel regression to eliminate historical trends and unobserved environmental and social characteristics:

$${{NDVI}}_{{it}} = \alpha +{\beta \cdot {Treat}}_{i}+{\gamma \cdot {Post}}_{t}+{\delta \cdot }({Treat}_{i}\times {{Post}}_{t})+{\mu }_{i}+{\theta }_{t}+{\varepsilon }_{{it}}\,$$
(1)

where \({{NDVI}}_{{it}}\) represents vegetation status for project i in year t and is logarithm-transformed; \({{Post}}_{t}\) defines the pre- and post-intervention period; and \(\gamma\) captures temporal trends; \({{Treat}}_{i}\) defines treatment and control groups; \(\beta\) indicates baseline differences between these two groups; \({{Treat}}_{i}\times {{Post}}_{t}\) is the DID estimate term; \(\delta\) quantifies the causal impact of the intervention; \({\mu }_{i}\) is the individual fixed effect of project i; \({\theta }_{t}\) is the fixed effect of year t; \(\alpha\) is the constant term; and \({\varepsilon }_{{it}}\) is the error term. The analysis initially follows this baseline econometric specification.

The post-intervention period is defined as 2012-2022, since the year 2012 marked the launch of China’s voluntary emission reduction market, signaling potential financial support and incentivizing early action. Furthermore, project activities often commence before formal submission, and emission reduction estimations consider pre-application activities. Notably, Figs. 2a and 3a demonstrate substantial vegetation index growth starting in 2013, indicating a rapid response to the policy announcement. These factors collectively justify 2012 as the appropriate start year for analyzing the CCER program’s impact.

The DID models are also employed to investigate the channels of forestry spillovers, including changes in land use, forest loss, and wildfire activity:

$${{lucc}}_{{it}}=\alpha +{\beta \cdot {Treat}}_{i}+{\gamma \cdot {Post}}_{t}+{\delta \cdot }({Treat}_{i}\times {{Post}}_{t})+{\mu }_{i}+{\theta }_{t}+{\varepsilon }_{{it}}$$
(2)
$${{loss}}_{{it}}=\alpha +{\beta \cdot {Treat}}_{i}+{\gamma \cdot {Post}}_{t}+{\delta \cdot }({Treat}_{i}\times {{Post}}_{t})+{\mu }_{i}+{\theta }_{t}+{\varepsilon }_{{it}}$$
(3)
$${{fire}}_{{it}}=\alpha +{\beta \cdot {Treat}}_{i}+{\gamma \cdot {Post}}_{t}+{\delta \cdot }({Treat}_{i}\times {{Post}}_{t})++{\mu }_{i}+{\theta }_{t}+{\varepsilon }_{{it}}$$
(4)

where \({{lucc}}_{{it}}\) represents the proportion of five types of land use in the buffer zones of the project i in year t, including forest land, residential areas, unused land, grassland, and water areas; \({{loss}}_{{it}}\) represents the forest loss detected in the buffer zones of the project i in year t; \({{fire}}_{{it}}\) represents the wildfire frequency in the entire area or per square kilometer in the buffer zones of the project i in year t.

To manage computational complexity while preserving the integrity of pixel-level data and maintaining robust statistical inference, we employed several strategic approaches. The baseline regression for project areas utilized all available pixels, totaling approximately 107 million observations; however, applying the same regression method to buffer zones would have increased the number of observations by more than 30 times. To address this, we implemented a clustering approach for buffer zones. Pixels within each project’s buffer zone were grouped into a fixed number of parcels based on proximity, with the mean value of each parcel serving as a single observation. We set the number of parcels at 10,000 (variations in this value were tested and found to have negligible impact on estimates), resulting in approximately 16.5 million observations for the baseline regression. In our robustness checks, we employed an alternative clustering method to account for the size of buffer areas across projects, where every 25 pixels (150 m²) were grouping into a single parcel based on proximity27. This adjustment enhances the comparability of results across projects with differing spatial extents, thereby improving the robustness of our analysis against scale-dependent effects. Furthermore, for placebo tests on both project sites and buffer zones, we utilized a subsampling technique23. Specifically, for each iteration, we randomly selected a number of pixels from the control pixels equivalent to 10% of the treated pixels to construct pseudo-treatment, and this process was repeated 1000 times in total for each of the two types of placebo tests, totaling 4000 times, to ensure statistical reliability.

Research data and variables

The study period is 2000–2022. This provides sufficient temporal coverage to evaluate the changes in outcomes: 12 years before and 10 years after the implementation of CCER. Data used in this study are the county-level socioeconomic variables and 1 km × 1 km resolution geo-bio-physical variables.

NDVI (Normalized Difference Vegetation Index) is the key dependent variable of this study, which is widely used for expressing forest status including density and health72, forest greening73, forest growth74, and forest productivity75. It is measured by the difference between the near-infrared (NIR) and the visible red light reflectance (RED) from satellite imagery. Mathematically, NDVI is calculated using the formula:

$${NDVI}=({NIR}-{RED})/({NIR}+{RED})$$
(5)

NDVI values range from −1 to +1. Negative values typically indicate non-vegetative cover such as water, snow, or clouds. Values near zero suggest bare soil or rock, while positive values indicate vegetation cover, with values closer to 1 representing denser and healthier vegetation. As plant chlorophyll absorbs red light for photosynthesis, greater plant growth leads to increased red light absorption and NIR reflection. Consequently, higher NDVI values generally indicate better plant growth and higher biomass68. In this study, we utilize two sources of NDVI data. Our primary source is 30-m resolution NDVI data derived from all available Landsat 5/7/8/9 remote sensing imagery. This dataset is processed to remove clouds and shadows, and the annual NDVI is calculated as the maximum value from all valid Landsat observations throughout the year. Additionally, we employ 1-km spatial resolution NDVI data derived from the SPOT/VEGETATION PROBA-V 1 KM PRODUCTS, which provides an annual perspective on NDVI across China from 1998 to 2019. This moderate-resolution data is generated using a maximum value synthesis method, capturing the peak monthly NDVI values for each year, with values ranging from 0 to 1. Notably, NDVI can exhibit a saturation effect in dense forest areas, where the index becomes less sensitive to vegetation changes. This could lead to underestimation of effects and increased uncertainty in our analysis. To ensure a more accurate and robust analysis, we also employed kNDVI, an index proposed by Camps-Valls et al. that corrects for saturation effects.

The high-resolution 30-m NDVI dataset is sourced from the National Ecosystem Science Data Center (nesdc.org.cn)76. The 1-km resolution NDVI data, land use and the landform dataset come from the Resource and Environmental Science Data Platform of the Chinese Academy of Sciences (resdc.cn). The China Multi-temporal Land Use Remote Sensing Monitoring Dataset is based on Landsat satellite imagery, where the available periods of 2000, 2005, 2010, 2015, and 2020 are used. This dataset employs a two-tiered classification system. The primary level encompasses categories including cultivated land, forest land, residential and industrial land, grassland, water area, and unused land. The secondary level classifies these into 23 distinct types. The landform dataset, derived from Landsat satellite imagery, is constructed through fusion of remote sensing products with classifications based on slope, aspect, and elevation.

The precipitation variable comes from the 1-km monthly precipitation dataset for China (1901–2022)77. We aggregated the monthly data to obtain the annual precipitation. The downward shortwave radiation is obtained from the China Meteorological Forcing Dataset (CMFD), which is estimated directly from 2-meter temperature, surface pressure, and specific humidity78. The nationwide established protected area and ecosystem dataset is sourced from the Geographic Remote Sensing Ecological Network Platform (gisrs.cn), which is compiled based on information provided by relevant authorities. It is worth noting that a minimal number of missing values may occur, primarily in peripheral areas of the map. To address this issue, we employed a nearest-neighbor imputation method, replacing any missing values with data from the closest available pixel.

The forest loss variable is derived from the Global Forest Change 2000–2023 dataset55. In this dataset, forest loss is defined as a stand-replacement disturbance or a change from forest to non-forest state. Using 2000 as the baseline year, the data detects loss primarily in the year 2001–2023 at a high spatial resolution of 1 arc-second per pixel (approximately 30 meters per pixel at the equator).

The wildfire data is sourced from the Moderate Resolution Imaging Spectroradiometer (MODIS) Burned Area Product provided by the U.S. National Aeronautics and Space Administration (NASA). MODIS, situated on the Terra and Aqua satellites, captures data across 36 spectral bands or groups of wavelengths, allowing for daily earth surface observations. The MODIS Burned Area Product provides daily information on global burning events at a 1-km spatial resolution from 2001 to 2022, which is extensively utilized as a data source for analyzing wildfire activities79,80. We organized the data to pinpoint China’s observations using latitude and longitude coordinates. Subsequently, we calculated the annual wildfire frequency within project and reference sites, as well as their associated 5-km buffer zones, considering the quality of individual hotspot pixels. In our study sample, an average of 1.12 active fires is detected per project annually, with an average of 30.38 detections within the buffer zones. The results of log-transformed total frequencies and frequencies per square kilometer were utilized to analyze the spillover effect of wildfire prevention.

The socioeconomic variables were obtained from the China Statistical Yearbook (County-Level) and Statistical Bulletin on National Economic and Social Development of provinces. We gathered and organized data on GDP, population, general budget revenue of local governments, and added value of secondary and tertiary industry. We transformed these data into log-transformed GDP, per capita GDP, and general budget revenue, and calculated proportions of secondary and tertiary industry, representing the industrial structure. We standardized the county using coding provided by the Ministry of Civil Affairs, unified different units, and corrected data with noticeable anomalies. For a few missing values, linear trends were used to impute them. Additionally, the proxy variable for county-level land use was the proportion of land used for forestry, extracted for each county from the aforementioned land use dataset.