Optimal harvesting in forestry: steady-state analysis and climate change impact

Abstract

We perform the steady-state analysis of a nonlinear partial differential equation model that describes the dynamics of a managed size-structured forest. The harvesting policy is to maximize the net benefits from timber production over an infinite planning horizon. The existence and uniqueness of the steady-state trajectories are analysed. Closed-form steady states are obtained in meaningful special cases and are used to estimate how climate change affects the optimal harvesting regime, diameter of cut trees, number of logged trees, and net benefits in the long run.

Keywords:

• sustainable forestry
• optimal control of nonlinear partial differential equations
• size-structured forest
• climate change

1. Introduction

Forests present a renewable resource, which provides timber and energy, maintains biological diversity and offers recreation facilities, mitigates climate change, and improves air quality 3. Human intervention, natural disturbances, and climate change may cause irreversible and unfavorable changes in the forest dynamics. It is difficult to set up an experimental design to analyse the changes because of a long rotation period of forests, and a theoretical analysis can play an essential role in studying the forest dynamics. In this paper, we attempt to build up a contemporary mathematical technique for analysing long-term forest development. We employ a size-structured model for the description of the forest dynamics because the size of a tree depicts forest biology and economic demands much better than its age. Moreover, empirical studies of real stands have shown that the age cannot be used as a proxy for the tree size, as age and size are only weakly correlated 5. The model is a size-structured-controlled version of the age-structured Gurtin–MacCamy model of population dynamics and includes mortality, growth, heterogeneity of trees, logging, and intra-species competition. The considered optimal control problem (OCP) maximizes the net benefits from timber production. Similar size-structured models of forest dynamics are commonly used in describing optimal forest management 79101617.

The mathematical novelty of the paper is in providing the steady-state analysis of a nonlinear OCP for the partial differential equation (PDE) model. This analysis can assist forest management as it presents closed-form solutions for steady-state optimal harvesting regimes in meaningful cases, which in turn allows for making some qualitative predictions. To explore sustainable forest development, we assume that the planning horizon is infinite. This assumption simplifies our investigation as it excludes undesirable end-of-horizon effects and allows us to perform the steady-state analysis of the OCP. The steady-state analysis was previously used for non-optimization models of size-structured populations in 718.

One of the crucial environmental policy issues is to determine to what extent climate change affects the future evolution of forests. Forest scientists suggest that climate change primarily increases the growth rate of forests, whereas its effect on tree mortality cannot be determined unambiguously 112. Thus, we compare the forest dynamics and optimal harvesting regime for different growth rates keeping in mind that they represent diverse climate scenarios 19. In particular, we investigate how changes in the growth rate affect the optimal harvesting rate, the harvesting size (the diameter of the cut tree), the total number of logged trees, and the net benefits of the forest.

The paper is organized as follows: Section 2 introduces the model of timber harvesting in size-structured forest. To make the model more realistic, we assume that the trees may not only reach a maximal size but also a maximal age. Section 3 derives optimality conditions and reviews some theoretical results relevant for our study. Section 4 analyses steady-state regimes in the model of the controlled forest growth. Section 5 investigates the dual system of the OCP as well as shows the existence and structure of optimal solutions. Explicit formulas for the optimal steady-state harvesting regime and the harvesting size are found in special cases. Section 6 concludes.

2. Problem Statement

Let us introduce the following notations:

Let us consider the following OCP:

Determine the functions x(t, l), u(t, l), E(t), p(t), , l∈[l₀, l_m], that maximize subject to the following constraints: under the given c(t, l), k(t), and χ.

In 7, the functions g(E, l) and μ (E, l) are assumed to be dependent on the population size . Following our objectives, here we consider their dependence on E as in 9.

The objective function (1) describes the net benefits of harvesting, which are equal to the revenue from the timber production minus the expenses to plant young trees. It is assumed that the management costs to maintain the forest and harvesting costs are included in the revenue. The expression cux in (1) means that the timber revenue is linearly proportional to the amount of the logged timber ux. A similar objective function is widely used in applications, e.g. in 2.

By the boundary condition (4), the actual density x(t, l₀) of the young trees equals . The restriction (5) states boundaries for the harvesting efforts and the planting trees and Equation (6) represents the density of trees at the initial moment. The restriction 0 is also common in harvesting population models (e.g. 518).

3. Optimality conditions and the structure of a solution

3.1. Optimality conditions

The OCP (1)–(6) has four unknown variables x, u, p, and E connected by two equations (2) and (3). We choose u and p to be unknown controls (independent variables). Then, the dependent variables x and E are determined from the non-local PDE boundary problem (2), (3), and (6).

Following 1420, we assume the following specifications for the given functions: , ], and ]. We also consider the unknown controls ] and , where is the space of all measurable functions bounded almost everywhere (a.e.) on the domain Ω. Then, analogous to 1620, it can be shown that for every given u and p from the domain (5), the nonlinear problem (2), (3), and (6) has a unique solution x, which is continuous in t a.e. on [0, ∞) and measurable in l on [l₀, l_m].

Theorem 3.1 (The maximum principle)

Let (p*, u*) be a solution of the OCP (1)–(6). Then where the functional derivatives are and is the unknown dual variable determined from the dual (adjoint) system as

Proof

Let us apply the method of Lagrange multipliers to the OCP (1)–(6) and choose to be a dual variable for the state equation (2). Then the increment of Lagrangean is Rewriting Equation (4) as and using integration by parts and the Taylor series, the second integral term in Equation (15) can be simplified as The variation in Equation (16) because x(0, l) is fixed by Equation (6). The transversality condition (12) and the initial condition (13) for the dual problem are obtained while deriving Equation (16).   ▪

Using (16) and obtained from Equation (3), the expression (15) becomes where γ (t) is defined by Equation (14). The first two integrals produce Equations (9) and (10) and the last two integral terms lead to the dual system (11).

The rest of the proof uses standard reasoning of the maximum principle theory 21.

All conditions and equations of Theorem 3.1 have applied interpretation. The economic explanation of the dual variable λ as the value of a tree (in situ value) requires its positivity. The transversality condition (12) is common in bioeconomics and reflects the behaviour of the process in the long run. It implies that the increase in the in situ value of a tree has to be smaller than the decrease in the discounting factor. The initial condition (13) shows that trees that reach their maximum size have no commercial value: they just die out. The conditions (8) and (10) state that the discounted in situ value of a young tree should be equal to its planting cost along an interior optimal regime. Otherwise, it is profitable to plant more trees or not to plant at all.

3.2. On the structure of solution

In this section, we present some known results on the dynamics of the growth rate g, optimal structure of the harvesting effort u, and the equation for the optimal logging diameter l, which will be used in this paper. Forest scientists estimate 112 that climate change significantly affects the growth rate of forest, whereas its effect on the mortality rate is insufficient. Moreover, vital parameters of the forest greatly depend on the growth rate. Following 1023, we assume that the growth rate g is given by where l̄>l_m, and β₀, β₁ are parameters chosen to reflect different common climate scenarios 19. We would like to point out that the model 1017 assumes l̄=l_m. The model (1)–(6) at l̄=l_m has a major drawback that will appear during further analysis of the dual problem (see Remark 1 in Section 5). Namely, it allows the trees to live indefinitely and, then, the maximum size l_m corresponds to the infinite age of a tree. To make the model more realistic, we assume that trees have not only the maximum size l_m but also the maximal age A_max. Then, the maximum size l_m Next, there are two major harvesting regimes: clear cutting logging and selective logging. Foresters and economists vigorously advocated clear cutting in the past. However, it seems that there has been a change of mind in the last decades and both harvesting regimes often coexist 23. Even Finland, where selective logging is prohibited by law, is questioning its previous evaluation 24. Speaking mathematically, a bang-bang structure of harvesting efforts u reflects selective logging when all trees with a certain diameter are to be harvested. Bang-bang principles for age-structured models are proven in 61318 and for size-structured models in 1417. Following 17, the strong bang-bang principle is valid for (1)–(6):

If Equation 17 holds and then the optimal control u*(t, l) in the OCP (1)–(6) has the following structure: with, at most, one switching size l*(t), , .

If , then Equation (21) follows from Theorem 3.1, namely, from Equations (9) and (7).

Conditions (19) are realistic. Indeed, the first condition requires the timber price c(t, l) to increase with the tree diameter l, and the other two suggest the mortality μ (E, l) to be constant or to increase with intra-specific competition E, or the size l of a tree. Although a constant mortality rate is often assumed for practical purposes, we will consider a more general case where μ (E, l) depends on E and l as well as the special cases with and

4. Steady-state analysis

One of our objectives is to find how climate change affects the optimal harvesting rate u*, the corresponding forest dynamics x*, and the number of harvested trees H in the long run. To study this issue analytically, we employ a steady-state analysis. Namely, we will look for possible solutions of the OCP (1)–(6), which do not depend on the current time t:

It is obvious that the necessary condition for the existence of steady-state solutions (22) is that the given prices c(t, l) and k(t) do not depend on t, c(t, l)=c(l), k(t)=k.

As in 718 and other work, we skip the initial condition (6) during the steady-state analysis. It means that possible steady-state solutions (22) are not solutions to the OCP (1)–(6) in the general case. Indeed, if a steady-state solution x(l) exists, it does not necessarily coincide with the given initial function x₀(l) in Equation (6). Only if it coincides with x₀(l) on [l₀, l_m], then x(l) is a solution to the OCP (1)–(6). Steady-state solutions describe the long-term strategy of the population model (2)–(4) and can become attractive trajectories for the OCP solutions. Such attractive behaviour is common in economic optimization models 418. However, we do not discuss it for the OCP (1)–(6) due to space restrictions.

We will provide the steady-state analysis of the OCP (1)–(6) in two steps, considering the direct system (1)–(4) in this section and then the dual problem and the OCP (1)–(6) in Section 5.

4.1. Steady state of controlled forest growth

In Section 3, the structure (20) of the optimal u* is established and determined by the optimal cutting size l*. We will discuss how to find the optimal p* and l* in Section 5. In this section, we investigate the existence of the steady-state solutions x and E to Equations (2)–(4) at the given p and u. A steady-state analysis of size-structured population models was earlier explored in 718. Paper 7 considers a similar size-structured model with the possibility of natural reproduction and proves the existence and uniqueness of its solutions, the existence and uniqueness of a steady state (stationary size distribution), and conditions for the convergence to the steady state. Paper 18 presents the analytic and numeric study of steady states in the uncontrolled size-structured model (2)–(4).

Definition

The continuous function x(l), l∈[l₀, l_m], that satisfies Equations (2)–(4) at p(t)=p and u(t, l)=u(l), is called the steady state of the controlled forest model (2)–(4).

Possible steady-state regimes x(l) should satisfy the following non-local integral-differential equation obtained from Equations (2)–(4): with the initial condition The solution of the initial problem (23), (25) is

Theorem 4.1

If the growth function g(E, l) satisfies Equation (17) and then there is a unique positive value E* that satisfies Equations (23)–(25). The unique solution x(l) of Equations (23), (25) is expressed by Equation (26) for E=E*.

Proof

Substituting Equation (26) into Equation (24), we obtain the following nonlinear equation: for the steady-state E*. The function F(E) is continuous, F(0)<0 and F(E)>0 for large enough E under the first condition (27). Therefore, there exists at least one positive E* such that F(E)=0. The derivative is positive under constraints (27). Therefore, the value E* is unique. The theorem is proven.   ▪

The conditions (27) are natural and the last one, , is the first condition (19).

4.2. Qualitative dynamics of the steady state

The qualitative behaviour of the steady-state trajectory x(l) is not clear because both the numerator and denominator of Equation (26) approach zero at l→l̄ under the growth law (17). Let us find its first derivative: Thus, , if . If , then x(l) increases up to ∞ at l→l̄, otherwise, it tends to zero at l→l̄.

Let u(l) have the bang-bang structure (20). Usually, the value u_max in (20) is much larger than the growth rate ĝ, so, on [l*, l_m]. If , then, considering formula (29) on two intervals [0, l*] and [l*, l_m] along the steady-state, we obtain that the density x(l) increases in l before the harvesting size l* and decreases to zero after l*.

Special case : Mortality μ depends on E and does not depend on l. In case (20), formula (26) becomes Expression (30) gives the explicit structure of the forest density x and allows us to analyse the harvesting yield in the next subsection.

4.3. Harvesting yield

The total number H of logged trees and the revenue R along the steady state are defined as Let us consider the case . Using (20), we obtain: where l* is the given harvesting age in formula (20) and x(l) is given by formula (30).

Let us estimate how the changing growth rate g affects H. For simplicity, we consider two growth rates g(E) and kg(E), . The value H for the second rate is By Equation (35), the number of logged trees H is larger for a larger growth rate: H<H_k, at least at and a fixed logging size l*.

Differentiating Equations (33) and (34) in l*, one can see that R increases in l* if H increases in l* and c′(l)≥0. Therefore, the revenue R is larger for a larger growth rate at the natural assumption that the tree price is not smaller for larger trees (the first condition (19)). Thus, at a fixed logging size l*, the harvesting yield R is larger and the number of trees H is larger when the growth is larger.

If we consider the OCP (1)–(5), then a change of the growth rate will also affect the optimal logging size l* of trees. This issue is investigated in the next section, where Equations (33)–(34) will be used to estimate the number of logged trees and the revenue along the optimal harvesting regime.

5. Steady-state analysis of the OCP

Section 3 states that the steady-state optimal management regime u*(l) in the OCP (1)–(6) under conditions (17) and (19) is where by Equation (21), the optimal switching size l*, , is determined as For a given x(l), the steady-state dual variable λ(l), l∈[l₀, l_m], is found from Equations (11)–(14) as where E depends on x(l) by Equation (24).

Remark 1 It is important to mention that, in the case l_m=l̄, the dual equation (38) is a singularly-disturbed ordinary differential equation with the singularity in the derivative term. It can be shown that Equation (38) at l_m=l̄ has no solution that satisfies the initial condition (39). So, our earlier assumption l_m<l̄ in Equation (18) is necessary for the regularity and solvability of the OCP (1)–(6).

5.1. No intra-species competition (linear case)

To clarify the qualitative behaviour of the optimal steady-state trajectories, let us start with the forest without intra-species competition: μ_E=0, g_E=0. Then, the forest dynamics is governed by the linear size-structured model (23)–(26) and the optimal steady state can be defined by explicit formulas.

If μ_E=0, g_E=0, then γ=0 by Equation (40), the explicit solution of the dual equation (38)–(39) is and the optimal switching size l* determined from Equation (37) is unique if c(l) is non-decreasing. Let us consider some special cases.

Case

1: Constant mortality and constant timber price: , and . Following Equation (41), for . By Equation (9), if the interior optimal harvesting size l* exists, it should satisfy Equation (37), that is, , or The left-hand of Equation (43) is always less than c, so, Equation (43) does not have a solution l* and the optimal regime is boundary. More specifically, by Theorem 3.1:

• the optimal , l*=l₀, and at ;

• the optimal p*=0 and u*(l)=0 at .

Indeed, in this case the price of a tree does not depend on its size, so, there is no reason to grow the trees and it is optimal to cut them immediately. If the future value of a young tree is less than the price k of planting a tree, then p*=0 by Equation (10).

To obtain more realistic optimal regimes, we shall assume that the selling price of a tree increases with its size. In order to keep the analytics reasonable, we consider a linear price c(l) in the next special case.

Case

2: Constant mortality and linear timber price: , c=c₀l. As shown in 1022, this case is quite realistic in forestry. Now, we can find an explicit formula for the interior steady-state harvesting size. By Equation (41), for At the natural assumption that the difference between the maximum size of the oldest tree and the asymptotic maximum size in Equation (17) is quite small: , we have The optimal harvesting age l* is now determined from and is given by The optimal l* is interior: , under the condition The condition (47) naturally holds since the size of the oldest tree is much larger than the size of young trees to be planted.

Rewriting Equation (46) as one can see that it is optimal to cut the trees of a larger size, when the growth rate ĝ is larger. Depending on the magnitude of the change in the growth rate, the rotation may be shortened or extended compared to the ex ante situation.

Let us estimate how a higher growth rate affects the harvesting revenue R and the number of logged trees H defined by Equations (33) and (34). Their dynamics is not straightforward because, when the growth rate is higher, the trees grow faster but are cut at a larger size and, as a result, we can cut less trees and have smaller revenue.

At , Equation (33) leads to For analytic simplicity, we take and l₀=0. Then, by Equation (46), and substituting it into Equation (48), the number of logged trees H is Differentiating H in g, we obtain The sign of Equation (50) depends on the relations among the given parameters μ, r, and g. We can immediately see from Equation (50) that dH/dg<0 at small μ such that μ≪r+g. On the other side, let us consider the case of small g, g≪r+μ. Then, expanding ln(x) in Equation (50) using the Taylor series, we obtain that dH/dg>0 if g<μ+r, that is, the number of logged trees H is larger for a larger growth rate only if the growth rate g<r+μ.

Analogously, the harvesting revenue is As before, we assume and l₀=0. Then we find that The first factor in (51) is the above function (49) and the second factor decreases in g, that means the larger harvesting revenue R for a larger growth rate if g<r+μ.

Using estimates of the parameters μ and g from the biogeochemical process model GOTILWA (Growth Of Trees Is Limited by Water, http://www.creaf.uab.es/gotilwa%2B) and reasonable values 0.01–0.05 of the discount rate r, we conclude that various qualitative dynamics is possible: the harvesting revenue may increase or decrease with the growth rate depending on the chosen climate change scenario. Hence, a more detailed sensitive analysis is required for specific parameters.

5.2. Nonlinear case with intra-species competition

Let μ_E>0 and g_E<0, that is, the intra-species competition increases the mortality rate and decreases the growth rate. Then, the solution of Equations (38)–(39) can be represented as where γ depends on λ (l), l∈[l₀, l_m], via Equation (54).

Lemma 5.1

If μ_E>0 and g_E<0, then γ<0 in Equation (38).

Proof

Applying the integration by parts to Equation (40), we obtain By the mean value theorem, there is a point l˜, , such that under the lemma conditions and because x and λ are both non-negative. The lemma is proven.   ▪

Lemma 5.2 (on solvability of the dual system)

At and the initial problem (38)–(39) possesses the unique solution λ (t), for any given x.

Proof

Substituting Equation (40) into Equation (52), we obtain the operator equation , with the operator A as We can show that the operator A is contracting 8 at the lemma conditions. Estimating the last integral, we obtain where denotes the norm in the space C[l₀, l_m] of the continuous function λ (t), and means the norm in the space of continuous μ or g, and the factor Cγ≪1 under the lemma conditions. The condition is equivalent to β₁≪1 in Equation (17) which is natural 310.

Let Ω_R be the sphere of the radius R in the space C[l₀, l_m] and . Then Therefore, the operator A is contracting if C_A<1 (which holds for small enough μ_E and g_E).

Next, estimating the difference , , , we obtain that Therefore, the operator A is invariant of the sphere Ω_R if we take

So, we obtain that the equation possesses a unique solution by the contraction mapping principle 11. The lemma is proven.   ▪

5.3. The optimal steady-state regime

As it is shown in Sections 4 and 5.2, both direct and dual problems have steady-state solutions under natural conditions.

Theorem 5.3 (on the optimal steady state)

Let and . Then, the optimal steady-state regime exists and is determined by Equations (26) and (36) with the optimal logging size .

If then the optimal p*=0 and u*(l)=0 (harvesting is not profitable), otherwise,

If the timber price c(l) is constant, then l*=l₀ (i.e. it is optimal to cut young trees immediately).

If c(l) increases, then the optimal logging size l* can be interior: . In particular, l* is interior, if one of the following conditions holds:

• (i)  or (ii) and .

Proof

By Section 3.2, the optimal harvesting rate u*(l) and the steady-state density x(l), l∈[l₀, l_m], are completely determined by l*, . The optimal harvesting size l* is defined from the nonlinear equation (37). By Lemma 5.2, the steady-state dual variable λ (l), , is uniquely determined from Equation (52) at a given x. Substituting Equations (36) and (52) into Equation (37) and denoting , we obtain one nonlinear equation G(l*)=0 with respect to the scalar unknown l*, where the function is continuous and negative at l=l_m. So, G(l₀) must be positive to have an interior . Applying the mean value theorem, we can simplify Equation (54) at as (here l<j(w)<w): Let us first assume that . Since γ<0 by Lemma 5.1, then from Equation (54) where . So, G(l₀)<0, equation G(l*)=0 has no interior solution, and the optimal harvesting size is boundary: l*=l₀ (compare to Case 1 of Section 5.1).

Let us show that the value G(l₀) can be positive if c(l) increases. Indeed, formula (55) leads to Next, analogous to inequality (56) we obtain (here . The first term is smaller than the second one if one of the conditions holds: (i) or (ii) and . Then, G(l₀)>0 and a value l*, l₀<l*<l_m, exists such that G(l*)=0. The theorem is proven.   ▪

Theorem 5.1 demonstrates the existence of the interior optimal cutting size l*>l₀ when the timber price c(l) increases. The case of a constant timber price has no practical value because the price c(l) includes also planting costs and, therefore, is negative at the beginning and turns positive thereafter. In practice, the price c(l) may be constant, starting with some size l_s from which trees/timber are traded and have a commercial value. Then, the size l_s will be the optimal cutting size. Lemma 5.1 and Theorem 5.1 lead to the following qualitative result: the optimal cutting size of the trees at the presence of intra-species competition is smaller compared to the case without competition.

6. Conclusion

The steady-state analysis is a common tool in the mathematical economics (also known as the balanced growth analysis and comparative static analysis). In this paper, we demonstrate that the steady-state analysis represents a prospective mathematical technique for optimization models of physiologically structured populations. Such models usually employ complex nonlinear PDEs. It is challenging to use these equations for making reliable predictions about policies of rational forest development, and the steady-state analysis may shed light on these issues.

The suggested technique determines time-independent infinite-horizon solutions of the model under study. Such solutions can describe the sustainable development of managed forests. Using known results on a bang-bang structure of optimal harvesting, we have reduced the OCP to a functional equation with respect to the optimal logging size of a tree and investigated the structure of its solutions. The value of the presented results for a forest practitioner consists of finding closed-form solutions for particular cases and presenting some qualitative predictions.

This paper also attempts to estimate the qualitative impact of climate changes on the forest dynamics. Climate change affects the growth rate of forests and, thus, we compare the optimal forest harvesting for various growth rates and obtain nontrivial harvesting regimes. The major findings are

• The optimal cutting size of a tree is smaller in the case of intra-species competition, compared to the case without a competition.

• For a fixed harvesting regime (fixed logging size), the timber revenue and number of logged trees are larger for a larger growth rate.

• Changes in the growth rate also affect the optimal harvesting regime. It appears to be optimal to cut the trees of a larger size when the growth rate is higher. As a result, the harvesting revenue along the optimal harvesting regime may increase or decrease in the growth rate depending on the climate scenario.

Acknowledgements

The authors express their deep gratitude to Prof. Saber Elaydi for organizing the 3rd International Conference on Mathematical Modeling and Analysis of Populations in Biological Systems, San Antonio, TX, October 7–9, 2011; participants of the conference for their useful remarks on the presented topic; Editor in Chief Dr J. M. Cushing and two anonymous reviewers for their valuable comments. The work is supported by the NSF DMS grant 1009197 and the Ministerio de Ciencia e Innovación grant Econ2010-17020.