Introduction 

By applying conditions of stable equilibrium from a matrix model, a recent study ([15]) found dimensionless numbers λ₀, š‚ and š′ to provide reliable approximations of, respectively, the population growth rate, λ₀, the “sustainable/stable” harvest rate, s, where (eqn. 1):

and the proportion of trees that has to remain unharvested in order to retain the stable diameter distribution, s′ (where s′ = 1/λ₀ = 1-s), of a forest stand, by means of the following expressions (eqn. 2, eqn. 3, eqn. 4):

where R is the global amount of recruitment, G the stand basal area, p_i the transition probabilities between adjacent diameter classes, and n and w, respectively, the number and width of the diameter classes, as detailed below (the use of the circumflex above a variable here denotes a new object that represents an approximation to that variable, as defined in eqn. 2, eqn. 3 and eqn. 4).

Dimensionless numbers are numbers that represent a property of a system, but do not have physical units, such as length, area, mass or time, associated with them. Where such numbers exist for a specific system, they allow conclusions about the behaviour of the system to be extracted, without the need for a more detailed formulation relating all the different variables involved in the problem at hand. Hence, their importance and wide use in many fields, such as mathematics, physics, biology, engineering, and economics ([27]).

In financial sciences in general, indicators to quantify financial return have been widely applied, using the internal rate of return (IRR) and the return on assets (ROA - [1]). The IRR is defined as the discount rate that makes the net present value (NPV) equal to zero, and the ROA is expressed as the quotient between profit and the sum of invested capital. Attempts to derive similar return indicators for forests are discussed in detail in Knoke ([11]). Mills & Hoover ([20]) calculated the annual change in the market value of standing timber stock and land, divided by the market value of timber and land at the beginning of the year, to obtain an indicator for US timberland. Weber ([33]) computed a similar index for German tree species in order to define an indicator for “forest performance”, and Hyytiäinen & Penttinen ([9]) applied the same perspective in their analyses of forest assets in Finland.

In this context, the main aim of this study was to investigate whether there were relationships between those dimensionless numbers (in particular, eqn. 3) and the NPV and PV of sustainable timber harvests from a monospecific uneven-aged forest. In addition, we have introduced two new dimensionless numbers, ρ_NPV and ρ_PV, using the NPV or PV derived from the sale of timber throughout a harvest cycle, plus the final stocking value (as numerator), and the fair value of the standing timber under IAS 41 ([10]) as denominator. These can be interpreted as economic performance indicators for forest management, inspired by the ROA accounting concept, showing how profitable the forest is, relative to its total value, with sustainability and stability criteria. Moreover, we demonstrated that there are strong linear correlations between those numbers and š; to that end, it will be used the matrix model applied to uneven-aged managed Pinus nigra stands by López et al. ( [16], [17], [18]), which accurately describes the stand dynamics.

Matrix models, introduced by Leslie ([13]) and modified by Lefkovitch ([12]) by grouping organisms in terms of stage categories rather than age categories, have been applied to almost all the subject areas of forestry ([30], [3], [31], [14]). In forest management, matrix models have been applied, for example, to evaluate economic outcomes ([2], [5]) and ecological impacts ([32], [34]), or to optimise silviculture in mixed uneven-aged forests to increase the recruitment of browse-sensitive tree species without intervening in the ungulate population ([4]). A matrix population model works in discrete time, projecting a population from time t to t + 1 using a transition matrix (1 here represents the time step or projection interval). It is defined by the finite difference linear system of equations X(t+1) = AX(t), where X(t) and X(t+1) are column vectors containing the number of stems ha^-1 within each diameter class at time t and t+1, respectively, and A is the transition matrix. In these models, the population growth rate (i.e., the temporal rate of exponential change of the population number of individuals in the long term) is the dominant eigenvalue, λ₀ of A, and the equilibrium position (stable diameter distribution in this study) is determined by the right eigenvector, W₀, corresponding to λ₀. By asymptotic analysis (long term behaviour), we know that, independent of the initial conditions, when λ₀ > 1, the total number of stems ha^-1 of the tree population increases exponentially over time (unless harvests are conducted); when λ₀ < 1 the population decays until extinction; and when λ₀ = 1, a stable distribution proportional to W₀ is obtained. Gotelli ([6]) referred to the special case of the stable distribution when λ₀ = 1 as the “stationary distribution”; we considered the same case in reference to the stable diameter distribution of the stand W₀.

  Materials and methods 

Stand population dynamics model

The matrix model uses the transition probabilities obtained in López et al. ([16], [17]) for uneven-aged managed Pinus nigra stands, considering three levels of tree diameter growth (Quality I, faster; Quality II, medium; and Quality III, slower diameter growth - Tab. 1, Fig. 1). The stands are located in the Spanish Iberian System, a mountain range extending about 400 km along the north-eastern edge of the central plateau, concentrating 60% of the area occupied by Pinus nigra in Spain ([7]). Based on data from those stands, the model simulates outcomes for 594 different forest management planning scenarios: three diameter growth levels (Qualities I-III) × six basal area levels (G = 21 to G = 26 m² ha^-1) × 33 recruitment levels (R = 200 to R = 840 stems ha^-1). Group selection cutting was the silvicultural system used because Pinus nigra is not a shade-tolerant species, and the gaps opened in the canopy through single-tree selection cutting would not be large enough to ensure regeneration.

Fig. 1 - Diameter growth models. I: site index 20, D = 51.68(1 - e^-0.015259t)^1.255111; II: site index 17, D = 46.645633(1 - e^-0.014318t)^1.337062; and III: site index 14, D = 40.644134(1 - e^-0.013838t)^1.456382. (D): dbh in cm; (t): time in years. The coefficient of determination was always greater than 0.999.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Since the harvesting operations in the study area generally took place every 10 years, this was the time step adopted in the model. Considering this time step, trees were grouped into n diameter classes of equal width, w = 6 cm: (0, 6), (6, 12), (12, 18), …, (30, 36), …, with the last class being (48,→; more than 48 cm) for Quality I, (42,→; more than 42 cm) for Quality II and (36,→; more than 36 cm) for Quality III. Therefore, an individual tree in class k could remain in class k or progress to class k+1 during the projection interval (t, t+1). The number of trees in each class changed in each projection interval, because some were harvested, some remained in the same diameter class and others grew past the boundary to the next diameter class. In such conditions, p_k is the probability that an individual tree in class k at time t (initial) will appear in class k+1 at time t+1 (final time of projection); the recruitment coefficient, r_k, is the number of offspring (stems ha^-1) living at time t+1 produced in the interval of projection (t, t+1) by an average tree in class k at time t (similar to many standard size classified matrix models, r₁ = 0 - [3]); h_k(t) defines the proportion of harvested trees in class k, natural mortality included; x_k(t) and x_k(t+1) describe the stem densities in class k at the initial and final times of projection. Finally, R = r₂ x₂ + r₃ x₃ + … + r_n x_n denotes the global amount of recruitment at each time step, defined as the total number of trees entering into the first diameter class from t to t+1, generally above a certain threshold diameter. Notice that 0 < p_k < 1 for k = 1, 2, …, n-1; r_k ≥ 0 for k = 2, 3, …, n, with R = 0, if and only if r_k = 0 for all k; and x_k > 0 for k = 1, 2, …, n.

By analysing the dynamics of the projections, we found that the model is described accurately by the matrix model (eqn. 5):

where (eqn. 6):

I is the identity matrix of order n, H(t) = diag[h₁(t), h₂(t), …, h_n(t)] is a diagonal matrix with the harvest rates h_k(t), including natural mortalities, and X(t) and X(t + 1) are column vectors indicating the stem densities at the initial and final times of projection, respectively.

Sustainable/stable harvesting strategy and stable diameter distribution

The “sustainable/stable” harvesting strategy is aimed at reaching the stable diameter distribution of the stand for each harvest, yielding the harvest rate, s, shown in eqn. 1 for all the diameter classes and time steps; the population growth rate, λ₀, is the dominant eigenvalue of the transition matrix without harvests (i.e., matrix A), and s also includes the natural mortalities. This is a sustainable harvesting strategy too (in the sense of population persistence over time) because s corresponds to the case in which the right eigenvector associated with the dominant eigenvalue equal to 1 of the transition matrix with harvests, i.e., matrix A(I-H), is W₀.

On the other hand, the equilibrium is defined by the stable diameter distribution of the stand W₀, which is simultaneously the right eigenvector of the matrices A and A(I-H) corresponding, respectively, to the dominant eigenvalues λ₀ and 1. Thus the components of the stable diameter distribution of the stand N₁, N₂, …, N_n (stems ha^-1) were obtained by solving the linear system AW₀ = λ₀ W₀, yielding (proportional to) the vector W₀ = (N₁, N₂, N₃, …, N_n-1, N_n), where (eqn. 7):

for 1< i < n-1, and (eqn. 8):

NPV function

Considering that the main aim for the stands in the study area is the production of timber, the objective function is the NPV of all management operations over a planning horizon of T years, discounted to the beginning of the period, that is (eqn. 9):

where v_k(t) is the stumpage value that corresponds to class k at year t (€ stem^-1) and i is the discount rate ([17], [18]). The first summation represents the income derived from the sale of timber, and the second is the final stocking value. Harvesting costs are not considered because it is assumed that harvesting operations are not paid for the owner, but for the buyer of the stumpage. The costs for the owner of the land are the opportunity cost of invested capital and the opportunity cost of land. None of them are explicit costs, but they are incorporated into the NPV function. It is well known that the NPV is a justifiable management objective for a single economic goal ([26]).

The corresponding function of that management regime maintained in perpetuity is given by (eqn. 10):

In the case of the equilibrium defined by the “sustainable/stable” harvesting strategy, since h_k(t) = s, x_k(t) = N_k, for k = 1, 2, …, n and t = 0, 10, 20, …, the above equations become (eqn. 11, eqn. 12):

Since the stumpage prices of Pinus nigra in Spain remain fairly stable through time ([19]), simplifying eqn. 11 and eqn. 12 we obtain (eqn. 13, eqn. 14):

where ∑ N_k v_k represents the fair value of standing timber under IAS 41 ([10]). By dividing eqn. 13 and eqn. 14 by this value, we obtain (eqn. 15, eqn. 16):

which means that, for previously defined values of i and T, there is a positive linear correlation between the NPV weighted by the fair value of standing timber, ∑ N_k v_k, that is ρ_NPV, and s (resp. between the PV weighted by ∑ N_k v_k, that is ρ_PV, and s). Note that ρ_NPV and ρ_PV are dimensionless numbers. Furthermore, given that there is a strong positive linear correlation between s and š‚ ([15]), it is expected that, for previously defined values of i and T, there must be a strong positive linear correlation between these variables ρ_NPV and ρ_PV, and š, as we demonstrate below.

Stumpage value model

The stumpage value model applied in this study comes from López et al. ([17], [18]). In Tab. 2, we provide the stumpage values we used to obtain the economic data ([22], [29]). The regression models adjusted to these data appear in Fig. 2. The stumpage prices of Pinus nigra in Spain remained fairly stable between 1994 and 2011 ([19]).

Fig. 2 - Stumpage value models: Quality I: v = D^3.186471 exp(-7.704952 - 8.678687 · 10^-3 D); Quality II: v = D^3.114196 exp(-7.476506 - 9.903125 · 10^-3 D); and Quality III: v = D^2.987053 exp(-7.110977 - 1.078752 · 10^-2D); stumpage values (v) in € stem^-1 and diameter (D) in cm. The coefficient of determination was always greater than 0.9999.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

The last row in Tab. 2 shows the average price per cubic metre for trees belonging to the different diameter classes. This price was obtained by adding up the four rows above in Tab. 2 that show the contributions to the final value of a tree for the four industrial destinations according to the percentage of wood of each diameter class going into each product type. The prices of particle board, poles, sawlogs, and high-quality sawlogs were 6.0, 39.5, 19.6 and 33.2 € m^-3, respectively. It is necessary to point out that the demand for timber poles in Spain is very high, and therefore they reach a higher price than larger-sized sawlogs.

  Results 

By substituting transition probabilities p_i, recruitment R, basal area G, together with w and n into the eqn. 3, we obtained š, as summarised in Tab. 3, for selected scenarios.

On the other hand, by applying eqn. 7 and eqn. 8, we obtained the components N_k of the stable diameter distributions W₀ of the stand, as summarized in Tab. 4 for some scenarios. The λ₀ values were obtained for each scenario under conditions of stable equilibrium from eqn. 5, using the transition matrix without harvests, and the s values came from eqn. 1.

On the basis of the above data, by substituting N_k, s, š, i, T and v_k into eqn. 13 to eqn. 16, we obtained NPV, PV, ρ_NPV and ρ_PV, as summarised in Tab. 5 and Tab. 6, for i = 3%, T = 90 years and several scenarios. Additionally, Fig. 3 and Fig. 4 depict the points (š, NPV) for T = 90, (š, PV), (š, ρ_NPV) for T = 90, and (š, ρ_PV), respectively (198 points for each case).

Fig. 3 - First row: set of points (š, NPV) for T = 90; second row: set of points (š, PV). Left: Quality I; middle: Quality II; right: Quality III (198 total points for each case, 33 points for each basal area level, six points for each recruitment level, and i = 3%).

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Fig. 4 - First row: set of points (š, ρ_NPV) for T = 90; second row: set of points (š, ρ_PV). Left: Quality I; middle: Quality II; right: Quality III (198 points for each case, and i = 3%).

  Enlarge/Shrink   Download   Full Width  Open in Viewer

The discount rates applied fell between 2% and 4%, as is typical in forest management studies, and the time horizon used to calculate the NPV ranged between 50 and 130 years. The influence of small variations in the discount rate i, or the time horizon T, on the results was minimal because they did not change the behaviour of the expressions shown in Materials and Methods (chapter “NPV functions”).

Regarding the regressions, we obtained the following results:

• As already mentioned, for previously defined values of i and T, the dimensionless numbers ρ_NPV and š have a strong positive linear correlation (see Fig. 4). For i = 3%, in particular, the total proportion of variance in ρ_NPV explained by š was r² = 0.99974 for Quality I, r² = 0.99965 for Quality II, and r² = 0.99956 for Quality III. It is important to note that such r² values do not depend on the planning horizon T.
• Similarly, for each discount rate i, the dimensionless numbers ρ_PV and š have a strong positive linear correlation, with the same r² values as described above (Fig. 4).
• For previously defined values of basal area (between G = 21 m² ha^-1 and G = 26 m² ha^-1), i and T, the variables NPV and š have a strong positive linear correlation (Fig. 3). The same result holds for the variables PV and š. The best approximations always occurred for the lowest level of diameter growth and the highest basal area scenarios, with r² being very close to 1 in all cases.
• For previously defined values of recruitment level (between R = 200 and R = 840 stems ha^-1), i and T, the variables NPV and š have a strong negative linear correlation (Fig. 3). The same result holds for the variables PV and š. The best approximations always occurred for the lowest level of diameter growth and the lowest recruitment scenarios, with r² being very close to 1 in all cases.

  Discussion and conclusions 

This paper proposes a simple and direct method based on dimensionless numbers, obtained under conditions of stable equilibrium from a matrix model, to provide reliable approximations of the NPV, PV, ρ_NPV and ρ_PV of a monospecific uneven-aged stand. In fact, these variables were related to š (see eqn. 3), which is a dimensionless number based on the following variables: (a) the transition probabilities p_k, which can be calculated from counts of observed transitions between adjacent diameter classes from t to t+1; (b) the basal area of the stand G, and the global amount of recruitment R. Although Pinus nigra regenerates well in the study area under its own canopy ([23], [28]), there is a critical basal area that allows sufficient light transmittance for regeneration. This maximum value is around 20-30 m² ha^-1 ([24], [25], [21]). On the other hand, basal area clearly below 20 m² ha^-1 (before harvesting) could compromise natural regeneration ([25]). Consequently, in addition to three levels of diameter growth (Qualities I-III), we considered six levels of stand basal area from G = 21 to G = 26 m² ha^-1, and 33 levels of global recruitment, from R = 200 stems ha^-1 (scarce recruitment) to R = 840 stems ha^-1 (abundant recruitment), creating a total of 594 planning scenarios. The corresponding problems simulated a wide variety of typical management cases in the study area; and (c) the number n and width w of the diameter classes. As shown in López & Fullana ([15]), these scaling factors n and w could have a limited influence on š and therefore on the NPV, PV, ρ_NPV and ρ_PV values.

The ρ_NPV and ρ_PV dimensionless numbers should not be interpreted as a kind of IRR, but as economic performance indicators for forest management inspired by the ROA accounting concept, using the NPV or PV derived from the sale of timber throughout a harvest cycle, plus the final stocking value (as numerator), and the fair value of standing timber (as denominator), showing how profitable the forest is relative to its total value, with sustainability and stability criteria. Application of the NPV (or PV) implies the introduction of a discount rate exogenously defined. In this regard, we applied the same criteria as in the International Financial Reporting Standards (IFRS), in particular IAS 41 ([10]), to estimate that exogenous discount rate. In fact, IAS 41 ([10]), in applying the expectation value approach to forest valuation, when referring to the process of deriving the present value of expected net cash flows for a forest stand, prescribed the use of a determined discount rate in each case, recommending that assumptions about these rates should be internally consistent with the underlying economic factors of the sector. By using a discount rate in accordance with international accounting standards, the possibility of manipulation should be avoided or limited. Thus, to consider the incorporation of eventual risks, in this study we initially chose a wide range of discount rates from 2% to 4%, typical of forest management scenarios ([9], [8], [11]). Note that (eqn. 17):

It is not surprising therefore that the highest ρ_NPV and ρ_PV values were obtained for Quality I, followed by Quality II and Quality III (see Tab. 6 and Fig. 4). In fact, since there is a strong positive linear correlation between these numbers and š, and considering that the core of this dimensionless number š is given by the quotient (eqn. 18):

under the same conditions for R and G, the faster the flow of the stems through the diameter classes (i.e., the higher product of the transition probabilities), the greater š and consequently the greater the ρ_NPV and ρ_PV values. Likewise, within each site quality (where the product of the transition probabilities may be assumed to be constant), the key factor underlying the ρ_NPV and ρ_PV equations is the quotient R/G, which means that the greater R/G, the more profitable the forest will be (i.e., the G = 21 m² ha^-1 with R = 840 stems ha^-1 is the most favourable scenario in this context). Two arrows in Fig. 4 point in the directions of the higher and lower R/G values.

Furthermore, the strong positive linear correlation between ρ_NPV and ρ_PV and š means that, stratifying the forest into uniform areas based upon tree diameter growths, the same linear regression line can be used in each stratum to estimate the ρ_NPV and ρ_PV dimensionless numbers. The same result holds for the variables NPV or PV and š for each basal area or recruitment levels, although, in the latter case, the corresponding strong linear correlation is negative.

Since the model rests on sustainability and stability criteria, it is important to investigate the sensitivity of the results to deviations from the equilibrium (stable diameter distribution); in particular, whether the strong linear relationship between the dimensionless numbers, ρ_NPV and ρ_PV, and š is sensitive to changes in the equilibrium position. In this regard, there have been two recent studies involving Pinus nigra stands in the study area ([17], [18]), both maximising the same NPV function, constrained by the same matrix model to describe the population dynamics, and with the same initial condition (the stable diameter distribution), as used in this study. Additionally, in the second one the final condition was assumed as the stable diameter distribution of the stand, and a wide variety of typical management scenarios in the study area, obtained by combining three rates of diameter growth (Qualities I-III), three levels of basal area (G = 22, 24 and 26 m² ha^-1), and three levels of recruitment (R = 200, 520 and 840 stems ha^-1), were considered. The first one (hereafter S1) showed that the economically optimal harvesting strategies underlying the solutions for that model changed the stand towards distributions substantially deviated from the stable diameter distribution; and the second (hereafter S2) revealed that the stand diameter distribution did not deviate substantially from the equilibrium position over time. The NPV values obtained for the optimal harvesting strategies from S2 were always below those of the optimal harvesting strategies from S1 (where no restrictions were introduced for the final state), and above the NPV values obtained under sustainability and stability criteria.

In each case, the ρ_NPV and ρ_PV dimensionless numbers were computed, and a regression analysis was conducted, using the same š values as in this study. The results showed that, in both S1 and S2 strategies the strong linear relationship between the dimensionless numbers ρ_NPV and ρ_PV, and š was not sensitive to changes in the equilibrium position, with r² values being very close to 1 in all cases. The lowest r² values always occurred for the lowest level of diameter growth, with the minimum value of r² being 0.9986 for the strategy S1.

In conclusion, this study showed that, under sustainability and stability criteria, there is a strong influence of the dimensionless number š on the NPV and PV of the forest harvesting strategies. Additionally, ρ_NPV and ρ_PV might be interpreted as economic performance indicators for forest management, inspired by the ROA accounting concept, showing how profitable the forest is relative to its total value, with sustainability and stability criteria. Those dimensionless numbers would also allow the extraction of conclusions about the behaviour of NPV of the forest harvesting, capturing information of all the global variables involved in the problem (e.g., recruitment, basal area and transition probabilities between diameter classes). Moreover, the study showed that there are strong linear correlations between the variables, NPV, PV, ρ_NPV , ρ_PV, and š.

  Acknowledgments 

The authors appreciate the helpful comments offered by the reviewers. The authors received no funding for this work, and have declared that no competing interests exist.

  References