Site description

The study was conducted in Santiago de los Caballeros Province, an effective area of 448 km² (19º 10’ - 19º 20’ N; 70º 45’ - 71º 05’ W), within the northcentral portion of Cordillera Central, Dominican Republic, in a region known as La Sierra (Fig. 1). La Sierra itself encompasses an area of 1800 km². Its climate varies depending on the altitude (from low montane to montane) and precipitation (from wet to very wet). Elevations above sea level range between 500 and 1600 m. Although irregular, rains manage to maintain a certain level of moisture in the soil for a large part of the year. Average precipitation ranges from 800 mm to 1600 mm annually; average annual temperature is around 24 ºC (^{Bueno-López, 2009}). The pine forests on these higher elevation slopes usually develop in shallow, carbonate, lateritic, low-producing soils, located primarily on land with rugged topography being this area more than 90% mountainous with slopes varying from 0% to 70%. Soils have a sandy texture although higher proportions of clay and silt can be observed in some places. They are usually well-drained, and where the surface is well protected with vegetation, fertility increases (^{Bueno, 2009}).

Source: Dirección de Información Ambiental y de Recursos Naturales, 2009.

Figure 1 Pinus occidentalis distribution within La Sierra, Dominican Republic. 

Sampling and analytical procedures

Selection of stands for plot installation involved placing a 2 km × 2 km grid on each of seven 1:20 000 scale maps, each map representing an area of 64 km² developed using SPOT panchromatic satellite images from 1991. The sampling frame for selecting stands included 120 non-overlapping pairs of grid central coordinates. From these, 24 central coordinates were randomly selected and undisturbed stands closest to each grid center that were unburned, reasonably free of wildlings, and free of apparent damage from insects and/or fungi were chosen as sampling units. Twenty-four temporary circular plots (350 m²), one in the approximate center within each stand, were established, encompassing different stand densities to represent the whole study area. In each plot, all trees were identified and measured for diameter at breast height at 1.30 m above ground (DBH, cm) and total height (H, m). Quadratic mean diameter (QMD, cm), density (trees/ha) and basal area per hectare (BA, m²/ha) were computed. Two trees per plot were selected, and a total forty-eight trees from dominant, codominant and intermediate classes were destructively sampled to measure biomass and C content of three components; stem plus bark, branch and foliage.

Trees were felled as close as possible to the ground (~10 cm). After felling, total tree height was measured to the nearest centimeter. Diameters to the nearest 0.5 cm and bark thickness to nearest cm were measured at stump and every two meters up the stem to a minimum top diameter of 4 cm. To partition each tree into component parts, branches were cut at their intersection with the main stem and foliage twigs with at most 1 cm diameter were separated from them. The main stem was cut into 2 m logs up to 4 cm minimum top diameter. Green weight of each component was determined in the field using heavy duty scales. Sawdust weight from main stem cuttings was also recorded and, added to its weight.

Subsamples from the main stem plus bark involved cutting 5 cm thick disks at stump height, 0.5 relative height, and 0.8 relative height. Subsamples from branches in the lower, middle and, upper portions of the canopy were also taken using 5 cm disks. Subsamples from foliage were taken from the lower, middle and upper portions of the canopy. All tissues from this subsampling stage were weighted in the field for green weight determination using high precision balances. For biomass determination, all sampled tissues were encoded, placed in Ziploc bags, and transported to the lab for further processing.

In the lab, samples were oven-dried at 105 ºC until constant weight was acquired. After drying, dry weight was measured, and dry weight / green weight ratio of components determined. Total dry biomass of each component (stem plus bark, branches, and foliage) was computed using

where:

B_i = total dry biomass of component i

B_Gi = total green weight of component i

S_i = sample dry weight of component i

S_Gi = sample green weight of component i

Summary statistics for variables from temporary plots and sample trees are shown in Table 1.

Allometric models

Graphic exploration of the data

As a first step in developing our models, we visually assessed the relations between our dependent variable, total above ground dry biomass (B_T); and predictor variables, DBH, H, and D²H, to be able to specify the mean relation and form of the error between these variables. We proceeded to plot coordinate points as a cluster to determine the nature of the mean relation and the form of the variance. If we found that the residuals were not constant (heteroscedasticity), we would adopt a power model to link the variance of the residuals to the effect variable included in the model.

Model fitting

Five different types of models were fitted to establish the relationship between B_T and predictor variables. For each type, two fitting methods, Weighted Least Squares (linear (WLLS) and non-linear (WNLLS)) and Maximum Likelihood, and two variance forms were assessed. Table 2 shows the combination of model types, variance forms, and fitting methods along with the model reference number in the study. Diameter was measured in centimeters and total height in meters.

Weighted regression models

Models 1, 2, 3, 4, and 5 were fitted by applying WLLS and WNLLS methods, aiming to stabilize the variance of the residuals assuming that each observation has its own variance, which is inversely proportional to the common residual variance. In equation form, this relationship can be depicted as follows:

where wi = positive weight associated with each observation and σi2 = residual variance for each observation. We also assumed, as suggested by ^{Picard et al., (2012)}, that the residual variance σ i is a power relation of the effect variable tree DBH, such that:

^{Picard et al., (2012)} also suggested that the exponent ξ must be fixed in advance, being this the main limitation of this method. To estimate coefficient ξ above, we divided the observations into five DBH classes and calculated the standard deviation of the total aboveground dry biomass in each of these classes and proceeded to plot it against the median DBH in each class on a log scale. We would check if the resulting plot would show us the points roughly aligned along a straight line, to confirm that the power model would be appropriate for modeling the residual variance. If that were the case, we would fit a linear regression of the log of the standard deviation of the biomass against the log of the median DBH for each class and consider the slope of the linear regression as an estimate of exponent ξ.

Linear and Non-linear models with a modeled variance

For models 6, 7, 8, 9, and 10 we explicitly modeled the variance of the residuals, assuming again that the residual standard deviation of the observations is a power function of effect variable DBH, such that:

Coefficients ζ and γ were model parameters estimated by means of Maximum Likelihood. Summarizing, models 6 through 10 were fitted by linear, multiple and non-linear regression of biomass against effect variables of DBH, DBH ² , H, and D²H by specifying a power model on the residual variance.

Checking hypotheses

We proceeded to produce quantile-quantile plots to visually verify that the residuals followed a normal distribution. The constant variance hypothesis of the residuals was assessed by visual inspection of scatter plots relating the residuals as a function of the predicted values, specifically making sure that the cluster of points would not show any particular trends in their structure.

Best model selection

Models developed were quantitatively evaluated by means of “100 Fold” cross validation (^{Efron & Tibshirani, 1993}). The distribution, bias, and precision of residuals was assessed to determine the accuracy of estimations (^{Vanclay, 1994}). The residuals were computed by subtracting the predicted from the observed aboveground total dry biomass values. Relative root mean square error (as percentage of average observed values) (RMSE%), relative bias (B%), mean absolute deviation (MAD), and coefficient of determination (R²) were calculated as follows:

where n = number of observations in the validation dataset

m = number of 𝛽 𝑖 parameters excluding β0

BTi = observed total aboveground dry biomass for tree i,

B^T = total aboveground dry biomass for tree i

B-T = mean total aboveground dry biomass for all sampled trees.

Ranking

To select the best model for total aboveground dry biomass estimation, we ranked the models based on four evaluation criteria (RMSE%, B%, MAD and R²). We assigned an independent score to each of these statistics based on its absolute magnitude. The lower the magnitude, the higher the rank. We then summed the assigned ranks and selected the one having the lowest score as the best model.

Carbon content

To estimate carbon content in the main stem, subsamples at each of three heights in three positions (heartwood, middle portion of sapwood, and outer portion of sapwood including bark) along the stem were pooled together after grounding. Subsamples from each of the three positions within the canopy of branches and foliage within the crown were pooled together after grounding for carbon content determination in each of these tissues. To obtain the carbon fraction (C), weights of approximately 50 mg from each of these components were placed on a Rapid CS cube (manufactured by Elementar Analysensysteme GmbH), which converts carbon compounds by oxidation at 950 °C into carbon dioxide (CO₂) which is measured via infrared (IR) detection.

Assessment of carbon content variation among tissues and its relationship with specific density

Analysis of Variance (ANOVA) was carried out to assess carbon content variation among stem wood, specifically, the outer portion of stem sapwood including the bark, foliage and branches. Within the stem itself carbon content variation between the pith, middle portion of sapwood and outer portion of sapwood including the bark at three different heights, was also examined by means of Two-Factor ANOVA. Also examined by Correlation Analysis was the relationship between C fraction and wood density.

Allometric models

Graphic exploration of the data

The results of the graphic exploration of biomass against the effect variables DBH, DBH ² and, DBH ² H are shown as scatter plots in Figure 2. The relationship between biomass and DBH (left panel) is not linear, and the variance of the biomass increases with DBH. Between biomass and DBH ² , the pattern is linear (central panel), but the variance of the biomass increases with DBH ² . The association between biomass and DBH ² H (right panel) is also linear, but again, the variance of the biomass increases with DBH²H.

Figure 2 Scatter plots showing the relationship between P. occidentalis total biomass and the effect variables DBH, DBH² and DBH²H. 

Weighted regression models

To stabilize the residual variance shown in our exploratory analysis of biomass against the effect variables assessed in the study, we employed as previously stated, weighted linear (WLLS) and non-linear (WNLLS) regression models as well as linear and non-linear models with a modeled variance. For the application of the weighed regression, we estimated coefficient ξ. Figure 3 shows a graph plot of the relationship of the standard deviation of biomass and median DBH in five DBH classes on a log scale.

Figure 3 Plot showing the relationship between the standard deviation of the biomass against the log of the median DBH for each of five DBH classes in logarithmic scale 

The points are approximately aligned along a straight line confirming that the power model is appropriate for modeling the residual variance. The resulting equation has a residual standard error = 0.233 on 4 degrees of freedom, coefficient of determination = 0.929, and P-Value = 0.008. The slope of the regression is 1.953, which is our corresponding ξ value. We assumed that 1.953 could be rounded off to two, and that the standard deviation of the biomass is proportional to DBH². Therefore, substituting on the weighting factor wi=DBHi-2ξ, we selected DBH⁴ as weight for the observations.

Best model selection

After performing “100 Fold” cross-validation, we proceeded to compute goodness-of-fit statistics from the 10 models predicting B_T. Once these statistics were computed, we assigned ranked points to each model based on absolute values of four of these statistics: RMSE% model 14, BIAS% model 15, MAD model 16 and R² model 17. Ranks of RMSE%, BIAS%, MAD and R² were summed and lowest summed rank score represented the best model (Table 3).

Based on ranked goodness-of-fit scores, the best model to predict total aboveground dry biomass of P. occidentalis is model 6 followed closely by models 1 and 4. Each of these three models incorporates DBH ² H as the predictor variable. Table 4 shows the parameter estimates and their corresponding standard errors for models 1, 4 and 6, along with t statistic and P-Values. Total aboveground dry biomass predictions for P. occidentalis by models 1, 4, 6, and 9 are depicted graphically as follows (Fig. 4).

Figure 4 Total aboveground dry biomass predictions for P. occidentalis by models 1, 4, 6, and 9. 

The predictions made by models 1, 4, and 6, which were fitted with this effect variable, are shown in Figure 5.

Figure 5 Residuals versus fitted (left column) and Quantile-Quantile plots (right column) for models 1 (first row), 4 (second row), and 6 (third row). 

Hypotheses assessment for WLLS and WNLLS regression models

Figure 5 depicts, from left to right, the graph of the weighted residuals against the fitted values above the quantile-quantile (Q-Q) plots for models 1 through 5, respectively. These residual plots do not show any trend. There is a little bit of structure in the residual and Q-Q plots but nothing of major concern.

Linear and Non-linear models with a modeled variance

To explicitly model the variance of the residuals, we assumed that the residual standard deviation is a power function of the effect variable DBH as depicted in equation 13. The coefficients ζ and γ would be estimated using the maximum likelihood method, but the multiplier ζ does not need to be estimated as the weights 𝑤 𝑖 are defined only within a multiplier factor (^{Picard et al., 2012}). Exponent γ is used to define the weighting of the observations and it is estimated as a parameter instead of being fixed in advance as the previous WLS models. Table 4 presents parameter estimates, their corresponding standard errors, exponent 𝜸 values and relevant statistics for model 6. Models 6 and 7 resulted with an intercept not statistically significant and were refitted without it.

Carbon content

We determined the carbon fraction on the stem including bark, branches, and foliage of individual P. occidentalis trees, through samples collected from 48 harvested trees. In total, 384, 48 and 48 randomly selected subsamples from each of above tissues, respectively, were analyzed. Table 5 shows summary statistics of carbon fraction (%) and content (kg) among tree components.

Significant differences in C% (by means of chemical analyses) among tree components were observed (P-value < 0.0001). Carbon concentrations were found to be lowest in the stem (42.059% ± 1.008% [95% C.I.]), followed by branches (46.355% ± 0.3385% [95% C.I.]), and foliage (49.80% ± 1.953% [95% C.I.]).

Wood C content within the stem varied significantly radially (P-value < 0.0001) with C concentrations from the pith (45.355% ± 2.157% [95% C.I.]) being statistically different to C content in the sapwood (40.988% ± 0.823% [95% C.I.]) and bark (39.659% ± 1.387% [95% C.I.]) at same stem heights. Results showed weak or no correlation between C fraction in any of the tissues and wood density, which was additionally estimated for each of the subsamples. Overall, if assuming general C fractions of 0.50, C stocks of this species would be overestimated by, on average, 13.25% in individual trees.

Aerial biomass in tree components

Mean dry-weight of biomass in main stem including bark, branches and foliage represent 74.5%, 13.3% and 12.1% of the mean total aboveground biomass, respectively. Individual trees had an average aboveground dry weight of 377 kg. With an estimated average 657 trees per ha at an average age of approximately 30 years (^{Bueno-López, 2009}), total aboveground dry biomass is approximately 247.69 Mg ha^-1. In a similar site for another close relative, P. caribaea (695 trees per hectare, 25 years), ^{Ambagahaduwa, Prasad, Gunatilleke, Seneviratne, and Gunatilleke, (2009)} estimated mean aboveground biomass at 194 Mg ha^-1. These authors also report that biomass distribution among main stem including bark, branches, and foliage was 60%, 17%, and 13%, respectively. ^{Subasinghe and Munasinghe (2011)} found that the average aboveground biomass on individual trees of a 27-year-old P. caribaea plantation was 369 kg.

Allometric models

P. occidentalis stem allometry for constructing volume and taper equations have been studied before (^{Bueno-López and Bevilacqua, 2010}; ^{Bueno-López and Bevilacqua, 2012}) but not to establish biomass relations. In contrast to the approach used in this study, most reported biomass equations for Caribbean region pine included bias in biomass estimation introduced using log transformation. Goodness-of-fit statistics for our models (n = 48, RMSE% = 6.17-9.96, BIAS% = 0.024-0.726, MAD = 0.05-0.07, R² = 0.86-0.95) show that total aboveground dry biomass predictions are generally good. Overall, models using DBH²H as the only predictor variable performed best based on goodness-of-fit statistics from cross-validation data.

Weighted and variance component models

Based on our ranking, model 6, fitted by modeling the variance of the residuals by Maximum Likelihood performed best, followed closely by models 1 fitted by WLLS and 4 fitted by WNLLS. Variance component models differed from the weighted regression models because the coefficients for weighing were parameters that needed to be estimated, and the least squares method could no longer be used to estimate model coefficients.

Model similarities

If we were to compare the predictions of total aboveground biomass made by these two techniques, WLS and ML, few differences between them could be discerned (Fig. 5). In fact, except for the intercept, the coefficients are very similar. Parameter coefficients for matching pairs of models - i.e., 1 and 6, 2 and 7, 3 and 8, 4 and 9, and 5 and 10 - are very close in their magnitude.

Carbon fraction among tree components

Carbon concentrations for P. occidentalis were found to be lowest in the stem and highest in needles (Table 5). Within the stem, wood C content varied significantly radially (P_Value < 0.0001) with C concentrations from the pith being greater and statistically different to C content in sapwood and bark at same stem height. We explored the possibility of developing separate equations to predict total carbon in individual trees and determined that it would be possible to accomplish these predictions from explanatory variables such as DBH and total tree height (R² = 0.83 and 𝜎=0.049). Overall, assuming that the generic C fraction is 0.50, it would lead to an overestimation of C stocks by 13.25% for this species.

Model predictions are very similar for WLS and ML models, especially for those having as DBH²H as the predictor variable (Fig. 5). Therefore, we conclude that regardless of the fitting method, aboveground dry biomass allometric equations using DBH²H as the predictor variable are the best for this species. Among the three statistical techniques used to fit B_T prediction models - namely weighted linear (WLLS) and non-linear (WNLLS) least squares as well as maximum likelihood with a modeled variance (ML), model 6 was the best, which was fitted by modeling the variance component and applying the Maximum Likelihood method. This selection was based on a ranking of four goodness-of-fit statistics, evaluated after performing 100-fold cross validation.

On average, each individual tree is storing 175 kg of carbon in the aboveground biomass components. Significant differences in the three tissues were obtained from an analysis of variance (P-value < 0.0001). Overall, assuming generic C fractions of 0.50, C stocks of this species would be overestimated on average by 13.25% in each individual tree. If we apply the conversion factor of C to CO₂ equivalent units (3.67), at this point of time, the amount of CO₂ equivalent per hectare stored in these P. occidentalis stands is 421.96 Mg. In our country, we do not have information regarding specific carbon concentration in forest trees. We consider that reporting C fraction for this tropical species, will significantly contribute to improve local estimates of terrestrial C stocks.