The performance benchmark procedure, ex-post and ex-ante

Procedure to Define the Performance Benchmark (Section 10.4)

Overview of Establishing the Performance Benchmark

To establish the baseline and calculate the performance benchmark, the methodology required to select and match control and project plots as described in Section 10.4 of VM0047. Specifically, it required to extract a number of project plots, as described in Step 1: Select project plots, extract a wide range of donor plots, as described in Select donor pool area, and match a number of donor to project plots, the control plots, as defined in Step 2, point 3 of Section 10.4 of VM0047. Specifically, 4 control plots were assigned to each project plot based off k-nearest neighbor optimal matching approach, as described in Step 3: Evaluate match quality and finalize matching. Stocking index (SI) for each project and donor plot was calculated based off the Normalized Difference Fraction Index (NDFI), calculated as defined in Stocking index calculation. SI data from the control plots was used to calculate the dSI for the baseline.

Step 1: Select project plots

Random plots of 30 x 30 m were generated within the project area boundaries. As detailed in Evaluate project plots, an initial pool of at least 150 plots was generated, of which 102 were carried into the matching, and the 45 plots with the closest match to their control plots over the historical period (t-10 to t) were retained as the project plots.

Step 2: Select control plots for each project plot

Remote sensing data acquisition

Harmonized surface reflectance data from 2008 to 2023 were extracted for each plot. Landsat 7 (Collection 2, Level 2; LANDSAT/LE07/C02/T1_L2) was used for years up to and including 2014, and Landsat 8 (Collection 2, Level 2; LANDSAT/LC08/C02/T1_L2) from 2015 onwards.

Landsat data was first cloud filtered and shadows were removed based off Landsat cloud and cirri bands. This allowed us to collect the median reflectance value for each plot over a period spanning from 1 May to 31 August of each year. These dates were chosen to guarantee the reflectance was representative of the peak of greenness and hence the maximum potential carbon absorption of each specific plot. The bands were rescaled to surface reflectance (DN × 0.0000275 − 0.2). The stocking index is computed from this reflectance (Stocking index calculation); for the band covariates of the matching it was then normalised to unit length (L2 norm) per observation, to reduce the reflectance bias between images.

Aboveground biomass (AGB) using the ESA CCI data product from 2017 -2020.

This data product provides a remote sensing derived estimate of aboveground biomass. It is used only to fit the linear relationship between the SI and aboveground biomass that converts the ex-ante estimates of biomass into SI (Estimate the relationship between biomass and SI). The SI itself is computed from fixed endmembers and does not depend on the biomass data.

Land Cover and ancillary information to all plots

Ancillary data was used to select eligible plots from the randomly extracted first donor pool. Two land cover layers were used, both for 2018 (t = 0): the Copernicus Global Land Service Land Cover 100 m, collection 3, epoch 2018 (version V3.0.1), and MapBiomas Colombia, Collection 3.0, 30 m (Fundación Gaia Amazonas, 2025), as applied in Select donor pool area.

Stocking index calculation

The stocking index is the Normalized Difference Fraction Index (NDFI) of Souza et al. (2005), calculated as follows:

$$NDFI = \frac{GV_{shade} - (NPV + Soil)}{GV_{shade} + NPV + Soil} \qquad(3.1)$$

Where $GV_{shade} = GV / (1 - Shade)$ is the shade-normalised green vegetation fraction, and GV, NPV, Soil and Shade are fractions of the pixel obtained by linear spectral unmixing of the Landsat surface reflectance. Specifically:

GV (Green Vegetation Fraction): Represents the proportion of green, photosynthetically active vegetation in the pixel.

NPV (Non-Photosynthetic Vegetation Fraction): Represents the fraction of dead or senescent vegetation (e.g., dry leaves, stems) in the pixel.

Soil (Soil Fraction): Represents the bare soil component within the pixel.

Shade (Shade Fraction): Represents the proportion of shade in the pixel; its endmember has zero reflectance in every band.

Cloud (Cloud Fraction): Represents residual cloud. It is unmixed with the other four but does not enter the NDFI.

The fractions were obtained by fully constrained linear spectral unmixing (fractions non-negative and summing to one) of the six Landsat surface reflectance bands (blue, green, red, near infrared, shortwave infrared 1 and 2) into the five fixed endmembers of CODED (Bullock et al., 2020, Table 1), which adopts the GV, NPV, Soil and Shade endmembers of Souza et al. (2005) and adds a cloud endmember (Table 1). The endmembers are fixed: they are not fitted to the data of the project region, and the SI does not depend on the biomass data. The cloud fraction is computed but, unlike in Bullock et al. (2020), it is not used to exclude observations; the annual composites are cloud- and shadow-masked beforehand (Remote sensing data acquisition).

Table 1: Spectral endmembers applied, in surface reflectance (Bullock et al., 2020, Table 1; Landsat TM and ETM+ bands 1 to 5 and 7, Landsat 8 OLI bands 2 to 7).

Endmember	Blue	Green	Red	NIR	SWIR 1	SWIR 2
GV	0.05	0.09	0.04	0.61	0.30	0.10
NPV	0.14	0.17	0.22	0.30	0.55	0.30
Soil	0.20	0.30	0.34	0.58	0.60	0.58
Shade	0.00	0.00	0.00	0.00	0.00	0.00
Cloud	0.90	0.96	0.80	0.78	0.72	0.65

Data quality. Observations outside the valid range of the Landsat Collection 2 Level 2 surface reflectance product (DN 7,273 to 43,636; USGS) in any band, or with an undefined NDFI, were treated as missing data and, as for any missing value (Evaluate project plots), a plot with such an observation was removed in full: 3 project-plot candidates and 36 donor plots (14 observations outside the range, 33 with an undefined NDFI).

References cited in this procedure

• Bullock, E. L., Woodcock, C. E., & Olofsson, P. (2020). Monitoring tropical forest degradation using spectral unmixing and Landsat time series analysis. Remote Sensing of Environment, 238, 110968. https://doi.org/10.1016/j.rse.2018.11.011

• Fundación Gaia Amazonas (2025). Documento de Bases Teóricas de Algoritmo (ATBD) MapBiomas Colombia Colección 3, Apéndice Colombia – Colección 3 de Mapas Anuales de Cobertura y Uso del Suelo de Colombia, version 1.0. Map: MapBiomas Colombia, Collection 3.0 of annual land cover and land use maps, 2018, https://storage.googleapis.com/mapbiomas-public/initiatives/colombia/collection_3/coverage/colombia_coverage_2018.tif (CC BY 4.0).

• Souza Jr., C. M., Roberts, D. A., & Cochrane, M. A. (2005). Combining spectral and spatial information to map canopy damage from selective logging and forest fires. Remote Sensing of Environment, 98(2–3), 329–343. https://doi.org/10.1016/j.rse.2005.07.013

• USGS. How do I use a scale factor with Landsat Level-2 science products? https://www.usgs.gov/faqs/how-do-i-use-scale-factor-landsat-level-2-science-products

Select donor pool area

15000 random plots were generated within a 50 km radius from the project area centroid, inside the 100 km limit of Table A1 of VM0047, and outside the boundaries of the project area. Plots that did not share the landcover with any of the project plots were removed from the collection of donor plots. The remaining donor plots were used to generate a geodatabase (.gpkg) that included both donor and project plots with the Landsat information for the years from 2008 to 2018 (data from t-10 to t0), and from 2019 to 2023 (t1-t5).

The original dataset was then reduced to 2,885 plots after removing plots that did not pass the donor-pool criteria of Table A1 of VM0047 (Appendix 1, Section 10.4, Step 2.1), as applied in Table 2.

The donor plots were then screened on MapBiomas Colombia, Collection 3.0, for 2018 (t = 0): a donor plot was kept only if its class is held by at least 5% of the project-plot candidates. Of the 105 candidates, 85 lie in class 21, Mosaico de agricultura o pasto (mosaic of agriculture or pasture), 18 in class 12, Formación herbácea (herbaceous formation), one in class 25 and one in class 11. Classes 21 and 12 were kept, and 2,264 donor plots remained. After the data-quality screen of Stocking index calculation, which removed 36 more, 2,228 donor plots were carried into the matching.

From the time series of 2008 to 2018, years for which Landsat data was missing from a large amount of project plots were removed (Table 3). This was necessary to ensure that the matching was correctly done when selecting the control plots as described in Select control plots.

Table 2: Required factors and source data to delineate donor pool area, evaluated for time t = 0. Time variant geospatial layers used must be current as of t = 0, ±5 years.

Factor	Data Source	Rule
Jurisdictional Boundary	legal boundaries of Colombia as the jurisdictional layer	Exclude plots outside subnational jurisdiction
Ecoregion	Olson, 2001 as the Ecoregion/Biomes layer	Exclude plots that are not in the same biome as the project area
Policy Environment	public land protected GIS data	Exclude plots in protected areas
Outside registered AFLOU projects	VERRA project registry	Exclude plots falling within other project areas
Land Tenure	Resguardos & Consejos, https://data-agenciadetierras.opendata.arcgis.com/search?collection=Dataset	Removed areas that are either concessions or reserves areas from the donor area
Land Cover	Copernicus Global Land Service: Land Cover 100m: collection 3: epoch 2018: Globe (Version V3.0.1)	Exclude areas with Land Cover classification label not in any of the labels found in the project area
Land Cover (second screen)	MapBiomas Colombia, Collection 3.0, 2018, 30 m (Fundación Gaia Amazonas, 2025)	Exclude plots whose class is not held by at least 5% of the project plots
Distance from project plots	Distance from centroid (50 km radius; Table A1 of VM0047 sets a 100 km maximum)	Exclude plots outside

Evaluate project plots

Project plots were generated by extracting a pool of random plots of 30x30m from a random selection within the project area boundaries to extract at least 150 project plots. Plots were extracted so that they had 0% overlap. No overlap was achieved by assuming that each plot was represented by the specific Landsat pixel in which the random coordinates were included. To make that sure that no plot was assigned the same pixel, a QA test was performed to make sure the same pixel was not used twice. To quantify historic and initial conditions plots an originally larger pool was extracted to make sure that Landsat reflectance was always available for at least 45 plots, for at least 3 years in the t-10 to t period, specifically:

• At least one year between t-10 and t-8

• At least one year between t-7 and t-1

• At t=0

Likewise, the plots needed to have a valid measurement of Landsat reflectance for all years of the first monitoring period (t1 to t5). A total of 105 plots met these conditions and, after the data-quality screen of Stocking index calculation, 102 were carried into the matching. Since more than 45 project plots were available, the 45 with the closest match to their Control plots in the historical period (t-10 to t) were chosen. The final number of candidate project plots and the selection of the best 45 plots is now documented in the Appendix titled “Dynamic Baseline and Performance Benchmark”. The choice of selecting 45 project plots was justified by the attempt to minimize the risk of having less than 30 project plots any time in the future where the Stocking Index will be needed to be calculated, because of faulty remote sensing data.

Table 3: Summary of the amount of plots missing from the Landsat dataset for each year, used to identify which years should not be included in the analysis. Years in bold were removed because of the large amount of Landsat pixels missing in Project plots locations.

Year	Plot type	Missing Landsat
2008	Donor	617
2008	Project	15
2009	Donor	1665
2009	Project	56
2010	Donor	134
2010	Project	0
2011	Donor	27
2011	Project	0
2012	Donor	1150
2012	Project	52
2013	Donor	5819
2013	Project	112
2014	Donor	344
2014	Project	3
2015	Donor	389
2015	Project	1
2016	Donor	797
2016	Project	1
2017	Donor	993
2017	Project	11
2018	Donor	316
2018	Project	6
2019	Donor	567
2019	Project	4
2020	Donor	186
2020	Project	2
2021	Donor	1233
2021	Project	0
2022	Donor	788
2022	Project	14
2023	Donor	215
2023	Project	0

Select control plots

Control plots were selected from the donor plots as the plots whose historical series best matched that of their assigned project plot over the t-10 to t period (Table 4). The covariates of the matching are the SI of each retained historical year (2008, 2011 and 2014 to 2018) and the principal components of the six Landsat bands of 2018 that together explain at least 80% of their variance, which are two; all nine are standardised over project and donor plots together. 4 donor plots were assigned to each project-plot candidate by optimal matching without replacement, the Manhattan distance being used as the metric of the distance between plots. To ensure that the matches reflected global rather than local optimisations, the Hungarian algorithm (linear sum assignment) was applied four times in succession over the 102 project-plot candidates, each time minimising the overall distance and excluding the donor plots already assigned. The 45 candidates with the smallest mean distance to their four control plots were retained as the project plots (Step 1: Select project plots).

The Manhattan distance was also used to calculate the weights of each control plot to their respective Project plot. Specifically:

$$w_{control,i,j} = \frac{e^{- MD_{i,j}}}{\sum_{j = 1}^{n_{i,j}\ }e^{- MD_{i,j}}} \qquad(3.2)$$

$w_{control,i,j}$= Weight of control plot j matched to project plot i (value between 0 and 1; dimensionless)

$MD_{i,j}$= Manhattan distance of control plot j to project plot i over the standardised covariates (dimensionless)

$n_{i,j}$ = Number of control plots matched to project plot i (equal to k at project start date)

Table 4: Summary of Project to Control plot match with distance (Manhattan distance) and ordinal order of control plot to project plot, defined as the closest to farthest Control Plot to its corresponding Project plot (workbook sheet T3_matching).

Step 3: Evaluate match quality and finalize matching

For the population of the 45 project plots and their matched control plots the match quality was finally estimated by calculating, for each covariate x, the standard difference of means (SDM), as follows:

$$SDM = \frac{\left| ({\overline{x}}_{wp,x} - \ {\overline{x}}_{bsl,x}) \right|}{\sqrt{\frac{\sigma_{wp,x}^{2} + \sigma_{bsl,x}^{2}}{2}}} \qquad(3.3)$$

Where:

$SDM$ = Standardized difference of means

${\overline{x}}_{wp,x}$= Mean value of covariate x in the population of project plots

${\overline{x}}_{bsl,x}$= Mean value of weighted sums of covariate x in the population of matched sets of control plots

$\sigma_{wp,x}^{2}$ = Standard deviation of covariate x in the population of project plots

$\sigma_{bsl,x}^{2}$ = Standard deviation of covariate x in the population of matched control plots, weighted by the weights of Eq. 2

Overall match results were deemed valid since the SDM of every covariate is below 0.25; the largest is 0.220, for the SI of 2008 (Table 5).

Table 5: Summary of SDM for each covariate used to calculate the matching of Control to Project plots (workbook sheet T4_sdm).

Step 4: Monitor control and project plots

The same control to project plots combination was used for monitoring SI and dSI after project start date. In case of one or more control plots will be deemed unsuitable either because not conform to the requirements of Table 2 or because missing Landsat reflectance in any year of monitoring, that plot will be substituted by the next best matching plot from the donor pool, excluding any control plot already matching with a project plot. In that case, the baseline will be recalculated for the 10 years temporal range used for the rest of the baseline and updated with the most recent 4 control to 1 project plot combinations and weights, as defined in Step 4 of the methodology (Section 10.4 of VM0047).

Step 5. Derive and evaluate slopes for time series of stocking indices

To be included in the dataset, SI values needed to be available at time t for the project plot i and all of its matched control plots i,j. Only plots for which Landsat data was available at each time t used for the analysis was in the derivation of 𝛥𝑆𝐼𝑐𝑜𝑛𝑡𝑟𝑜𝑙,𝑡 and 𝛥𝑆𝐼𝑤𝑝,𝑡. The rate of increase in stocking index in the control and project plots, ${\mathrm{\Delta}SI}_{control,t}$ and ${\mathrm{\Delta}SI}_{wp,t}$, was calculated as the slope of the weighted linear regression of the accumulated time series of SI values. For each project plot, the weight o each control plot was determined by how close the control plot SI matched their respective project plot SI in the t-10 to t period. Specifically, weights of SI values for control plots in the time series are calculated as:

$$w_{control,i,j,t} = w_{control},i,j\frac{1}{\sum_{t = 0}^{t}{n\_ rs_{t}}} \qquad(3.4)$$

Where:

$w_{control,i,j,t}$ = Weight of control plot j matched to project plot i at time t (dimensionless)

$w_{control},i,j$ = Weight of control plot j matched to project plot i (value between 0 and 1; dimensionless)

$n\_ rs_{t}$ = Number of project plots and matched control plots (i,j) with values assessed at time t

A representation of the weights and assignment of control plots for each project plot is provided in Stocking index per plot and year.

The SI was calculated as defined in Eq. 1 for each project plot and each time stamp. The SI for each batch of 4 Control plots linked to 1 Project plot was calculated as the weighted average of their NDFI, using the weights as defined in Eq. 4.

The values of SI were used to calculate the dSI, that is, the slope of the linear relationship between time and SI. This slope represents the rate of increase in stocking index in either the control (${\mathrm{\Delta}SI}_{control,t}$) or project plots (${\mathrm{\Delta}SI}_{wp,t}$).

The Summary of Project and Control plots SI for each Project to Control plot combination is given in Stocking index per plot and year. The summary includes the weight of each control plot and its SI, used to calculate the average Si control based off Eq. 4 and Eq. 1.

Step 6: Derive performance benchmark

To derive the Performance Benchmark (PB) two tests needed to be applied first. First, the ${\mathrm{\Delta}SI}_{control,t}$ and ${\mathrm{\Delta}SI}_{wp,t}$ needed to be deemed statistically different. This was needed to avoid accounting for additionality in any case project activities did not result in significant deviations form the baseline. To do so, a Z test was performed as follows:

$$Z = \frac{\left| ({\mathrm{\Delta}SI}_{wp,t} - \ {\mathrm{\Delta}SI}_{control,t}) \right|}{\sqrt{{SE}_{{\mathrm{\Delta}SI}_{wp,t}}^{2} + \ {SE}_{{\mathrm{\Delta}SI}_{control,t}}^{2}}} \qquad(3.5)$$

Where:

${\mathrm{\Delta}SI}_{control,t}$ = Average annual increase (slope) in stocking index SI in control plots through time t

${\mathrm{\Delta}SI}_{wp,t}$ = Average annual increase (slope) in stocking index SI in project plots through time t

${SE}_{{\mathrm{\Delta}SI}_{wp,t}}^{2}$ = Squared standard error of the average annual increase (slope) in stocking index SI in project plots through time t

${SE}_{{\mathrm{\Delta}SI}_{control,t}}^{2}$ = Squared standard error of the average annual increase (slope) in stocking index, SI, in control plots through time t

t = 1, 2, 3, …, t years elapsed since the project start date

The test was deemed successful (dSI significantly different) if the value of Z was equal to or greater than 1.96.

The second test required evaluate whether ${\mathrm{\Delta}SI}_{wp,t}$ or ${\mathrm{\Delta}SI}_{control,t}$ were negative. A negative dSI signified a biomass loss, therefore if project activities produced negative no additionality was to be accounted.

Where parameters ${\mathrm{\Delta}SI}_{control,t}$ and ${\mathrm{\Delta}SI}_{wp,t}$ were significantly different, ${\mathrm{\Delta}SI}_{wp,t}$ was positive, and the ${\mathrm{\Delta}SI}_{control,t}$ was significantly different from 0 (p value < 0.05), the PB was calculated as the ratio of the average ${\mathrm{\Delta}SI}_{control,t}$ to the average ${\mathrm{\Delta}SI}_{wp,t}$, as defined in Eq. 6.

$$PB_{t} = \frac{{\mathrm{\Delta}SI}_{control,t}}{{\mathrm{\Delta}SI}_{wp,t}} \qquad(3.6)$$

Where the slope coefficient of the control plots, ${\mathrm{\Delta}SI}_{control,t}$, was insignificant (P > 0.05) or negative, ${\mathrm{\Delta}SI}_{control,t}$ was set equal to zero and the PB was assigned the value of 0. When the slope coefficient of the project plots, ${\mathrm{\Delta}SI}_{wp,t}$ was insignificant (P > 0.05) or negative, the PB was assigned the value of 1 (non-additionality).

Procedures to Determine ex-ante Performance Benchmark

Estimate the relationship between biomass and SI

Since the Remote Sensing data is missing for future timeseries, it was required to estimate the SI from ex-ante estimates of biomass. This was possible because the SI is meant to be a proxy of biomass. To achieve that in the most accurate way, a dataset of historical biomass estimates and SI was generated. The biomass data was retrieved for the plots for a time series spanning from 2017 to 2020 from the ESA-CCI Biomass dataset, described in Remote sensing data acquisition; observations below the 2.5th or above the 97.5th percentile of biomass were removed, leaving 769 plot-by-year observations (Biomass and stocking index observations). The SI was calculated with Eq. 1, the same index used for the historical period (t-10 to t) and for the ex-post estimates of ${\mathrm{\Delta}SI}_{wp,t}$ and ${\mathrm{\Delta}SI}_{control,t}$.

This data was used to fit the linear regression (Table 6; r = 0.389, p < 0.001, n = 769):

$$SI_{ex-ante}=\ 0.004026590*Biomass\ +0.39875586 \qquad(3.7)$$

Table 6: Linear regression of the stocking index on aboveground biomass (Mg ha⁻¹), fitted to the paired observations of Biomass and stocking index observations; it gives the coefficients of Eq. 7 (workbook sheet biomass_si_regression).

Calculate ex-ante SI

The baseline ${\mathrm{\Delta}SI}_{control,t}$ was assumed to be equivalent to the ${\mathrm{\Delta}SI}_{control,t}$ calculated for the timeseries from t-10 to t. Since the ${\mathrm{\Delta}SI}_{control,t} < 0$, a ${\mathrm{\Delta}SI}_{control,t} = 0$ was assumed for the first 10 years. From year 11 the ${\mathrm{\Delta}SI}_{control,t}$ was calculated as the average for the previous 10 years. Yet, since the ${\mathrm{\Delta}SI}_{control,t}$ was 0, the average was 0 as well, and so a ${\mathrm{\Delta}SI}_{control,t}$ was used for baseline throughout the ex-ante period.

For ${\mathrm{\Delta}SI}_{wp,t}$, the SI was calculated by applying the linear regression between SI and biomass calculated from the historical data, as defined above in Estimate the relationship between biomass and SI on the ex-ante estimates of biomass per hectare. The SI of a time series of 10 years prior time t was then used to calculate the ${\mathrm{\Delta}SI}_{wp,t}$. For years before t10, all years from year 1 to t10 were used (therefore a time series with less than 10 years until 10 years were available from ex-ante biomass estimations). Results are given in Ex-ante stocking index by management stream, the summary table of the ex-ante SI calculated from Eq. 7 using ex-ante estimation of AGB per hectare for the two different management typologies: Acacia to Eucalyptus (AOE) and Acacia to Pine (AOP). These estimates were used to calculate the dSI for timeseries of up to the 10 years prior, for each year t. For years 1 to 9, the time series was from t0 to t.

Calculate ex-ante Performance Benchmark

The ex-ante performance benchmark was then calculated using the same formula and rationale as previously described in Step 6: Derive performance benchmark and is summarized for each year in Ex-ante performance benchmark.