Forecasting Variations in Profitability and Silviculture under Climate Change of Radiata Pine Plantations through Differentiable Optimization

Article Forecasting Variations in Profitability and Silviculture under Climate Change of Radiata Pine Plantations through Differentiable Optimization Miguel A. González-Rodríguez 1,2,* , Miguel E. Vázquez-Méndez 3and Ulises Diéguez-Aranda 2   Citation: González-Rodríguez, M.A.; Vázquez-Méndez, M.E.; Diéguez-Aranda, U. Forecasting Variations in Profitability and Silviculture under Climate Change of Radiata Pine Plantations through Differentiable Optimization. Forests 2021,12, 899. https://doi.org/ 10.3390/f12070899 Academic Editor: Mathias Neumann Received: 7 June 2021 Accepted: 8 July 2021 Published: 10 July 2021 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). 1 CERNA Ingeniería y Asesoría Medioambental S.L., R/Illas Cíes n º 52-54-56, Ground Floor, 27003 Lugo, Spain 2 Unidade de Xestión Ambiental e Forestal Sostible, Departamento de Enxeñaría Agroforestal, Universidade de Santiago de Compostela, Escola Politécnica Superior de Enxeñaría, R/Benigno Ledo, Campus Terra, 27002 Lugo, Spain; [email protected] 3 Departamento de Matemática Aplicada, Instituto de Matemáticas, Universidade de Santiago de Compostela, Escola Politécnica Superior de Enxeñaría, R/Benigno Ledo, Campus Terra, 27002 Lugo, Spain; [email protected] *Correspondence: miguelangel.gonzalez.r[email protected] Abstract: Climate change might entail significant alterations in future forest productivity, profitability and management. In this work, we estimated the financial profitability (Soil Expectation Value, SEV) of a set of radiata pine plantations in the northwest of Spain under climate change. We optimized silvicultural interventions using a differentiable approach and projected future productivity using a machine learning model basing on the climatic predictions of 11 Global Climate Models (GCMs) and two Representative Concentration Pathways (RCPs). The forecasted mean SEV for future climate was lower than current SEV ( ∼ 22% lower for RCP 4.5 and ∼ 29% for RCP 6.0, with interest rate = 3% ). The dispersion of the future SEV distribution was very high, alternatively forecasting increases and decreases in profitability under climate change depending on the chosen GCM. Silvicultural optimization considering future productivity projections effectively mitigated the potential economic losses due to climate change; however, its ability to perform this mitigation was strongly dependent on interest rates. We conclude that the financial profitability of radiata pine plantations in this region might be significantly reduced under climate change, though further research is necessary for clearing the uncertainties regarding the high dispersion of profitability projections. Keywords: differentiable optimization; Pinus radiata; stand-level management; climate change; risk modelling 1. Introduction Climate change is intended to shift forest dynamics in the following decades [ 1 ]. Declines in forest productivity and fast changes in species suitability are among the potential negative consequences of global warming [ 2 – 4 ]. These consequences may compromise the ability of forest ecosystems for producing goods and services, leading to socioeconomic fallouts, such as scarcity in timber supply chains [ 5 ], turns in timberland value appreciation [ 6 ], and food and energy shortages in rural vulnerable communities [7]. In recent years, the concern for proactively adapting to shifts in forest productivity has provoked a scientific turnaround in the field of empirical growth and yield modelling [ 1 , 8 ]. The current research trend aims at developing growth-environment relationships through predictive modelling, mainly focusing on the site index ( SI ), the most frequent empirical indicator of forest productivity [ 9 ]. A variety of supervised learning techniques have been used for this purpose [10–12], yielding, overall, successful results (R2∼0.3–0.7). Even so, connecting future forest productivity predictions with its economic and silvicultural repercussions is still an uncertain task. In this regard, several recent studies have evaluated financial risks associated with uncertain future productivity basing on Forests 2021,12, 899. https://doi.org/10.3390/f12070899 https://www.mdpi.com/journal/forests Forests 2021,12, 899 2 of 14 optimization at stand-level [ 13 , 14 ] and forest-level [ 15 , 16 ]. These studies consist, in summary, on the numerical optimization of a certain financial indicator (e.g., the soil expectation value), which depends on decision variables associated with silviculture and investments, under varying economic and climatic conditions. According to Pasalodos-Tato [17] , most of the previous research on stand-level management optimization has relied either on dynamic programming methods or on direct search methods. Dynamic programming was the earliest of both techniques [ 18 ] and consisted basically of simplifying the optimization problem by dividing it into a series of simpler problems that were solved recursively. Though it had the major advantage of ensuring convergence to the global maximum, it also implied some important disadvantages, such as the need for discretizing decision and state variables [ 19 ]. Direct search methods were applied for first time in forestry by Kao and Brodie [20] and since then have been extensively used for stand-level optimization. In comparison with dynamic programming, direct search methods provide reasonably good solutions faster and can implement continuous decision and state variables. However, direct search methods do not ensure the convergence towards the global maximum [ 21 ]. Several families of direct search methods have been applied in recent decades for optimizing stand-level management, including the one solution vector methods, such as the Hooke and Jeeves algorithm [ 22 ], and the populations-based methods, such as Differential Evolution [23], Particle Swarm Optimization [24] and Evolution Strategy [25]. To cope with some of the disadvantages of these techniques, [ 26 ] proposed the use of differential optimization methods. The latter allow working with continuous decision variables, thus avoiding the information loss due to discretization, and produce good solutions in a relatively short computing time. In their comparative analysis, [ 26 ] found that a differentiable method was between ∼ 3 and ∼ 20 times more efficient, in terms of computing costs, than Hooke and Jeeves and Differential Evolution, respectively. Moreover, some of the observed computing limitations of dynamic programming and direct search seem to scale sharply as we increase the number of decision variables [ 21 , 26 ], as it can be the case when thinnings are implemented in addition to clearcutting-related decision variables. Concerning the specific problem of management optimization under uncertain future productivity, the usual approach consists on the use of risk metrics derived from a covariance analysis between risk factors (i.e., productivity) and profitability. Under this approach, especially frequent in the field of modern portfolio optimization [ 6 , 13 ], the covariance analysis is based on a simulation-based procedure in which the objective function is evaluated exhaustively, encompassing a wide range of combinations of silvicultural, economic and climatic conditions. Considering the high number of simulations that this task might imply, computational efficiency becomes an important concern that should guide the selection of the optimization method. In this context, the use of dynamic programming and direct search methods can lead to a certain level of oversimplification in the problem setup (i.e., a reduction in the number of decision variables and possible solutions via discretization) to reduce computing costs. In this article, we simulated the development of forest plantations in the northwest of Spain under different climate change scenarios. For each scenario, we optimized economic profitability using stand-level differential optimization. We focused our scope on a set of radiata pine (Pinus radiata D. Don) stands distributed mainly in the Spanish province of Lugo. After the simulations, we evaluated the changes in financial profitability and risk between current and future climatic conditions. As a side goal, we analyzed the changes in optimum silviculture variables, such as the rotation length. 2. Materials and Methods 2.1. Optimization Approach We optimized forest stand management following a similar methodology to that used by Arias-Rodil et al. [26] in the NW of Spain for pure, even-aged stands of Pinus pinaster Ait. In it, the development of a forest stand is simulated using the dynamic systems-based framework (frequently referred to as “state-space” approach) first used in Forests 2021,12, 899 3 of 14 forestry by García [27] . According to it, the forest stand is characterized at each moment by state variables whose evolution, described by time-dependent transition functions, is considered independent of previous states. These functions represent natural dynamics (growth and mortality), are affected by control variables which encapsulate the impact of silvicultural treatments on the state, and are complemented by output functions that translate state variables into outcomes (e.g., timber volume). An economic model provides the objective function value corresponding to these control variables and initial stand conditions (see Figure 1). SIMULATOR Output functions Objective function Control variables (management prescription) utiIiRitnt+1 = Outcomes Vd ( , )ut... Vi r,d ( )u Objective value SEV( )u DYNAMIC SYSTEM MODEL State variables at an initial age t0 H0N0G0 Transition functions State variables at age t H t( ) N t( , )uG t( , )u ECONOMIC MODEL Prices ( ) and costs ( )p C j Figure 1. Flow diagram of the simulator developed for computing the objective value (SEV) of a given management prescription. We used the most frequent setting under this approach, which characterizes the stand using three state variables: dominant height (mean height of dominant trees in the stand, in metres), H(t) , number of stems per hectare, N(t) , and stand basal area (total area of stem sections at 1.3 m, in m 2 /ha), G(t) . The dynamic of these variables is described by species-specific transition functions h , n , g : R+×R+×R+−→ R , so that, if for an initial age t0≥ 0 there is a known state where H(t0) = H0 , G(t0) = G0 and N(t0) = N0 , and no silvicultural treatments are applied in [t0 , t] , then it is verified that H(t) = h(t0 , H0 , t) , N(t) = n(t0,N0,t)and G(t) = g(t0,G0,t). Concerning the simulation of silvicultural treatments, each ( i -th) thinning is characterized by its intensity (proportion of stems removed), Ii , removal relation (ratio between the proportion of stand basal area removed and the proportion of stems removed), Ri , and timing, ti . Then, a management prescription is defined by the number of thinnings, nt∈N , and the vector (control variable) u= (I1 , R1 , t1 ,..., Int , Rnt , tnt , tnt+1)∈R3nt+1 , determining these thinnings and the rotation age tnt+1 . As Ri is usually kept in the interval ( 0,1 ] , the dominant height is not affected by treatments (note that a value of Ri> 1 would lead to thinning from above). The state variables H(t) , N(u , t) and G(u , t) can be predicted at any age with the transition functions h , n , g , and the control variable u (see [ 26 ], for details). Outputs are then obtained from the predicted values of the state variables. For instance, the merchantable timber volume (in m 3 /ha) for a certain product (defined by a limit diameter d , in cm) can be obtained from a known output function v(H , N , G , d) : for any t≥0, Vd(u,t) = v(H(t),N(u,t),G(u,t),d) Forests 2021,12, 899 4 of 14 gives the merchantable timber volume with diameter greater than d at age t . From this function, the removed timber volume at the i-th thinning is calculated as Vr,d i(u) = Vd(u,t− i)−Vd(u,t+ i), where t− i and t+ i denote, respectively, the instants before and after the i -th thinning. In the same way, the removed timber volume at the rotation age is given by Vr,d nt+1(u) = Vd(u,tnt+1). Considering that the purpose is to compare the profitability of management alternatives with different rotation lengths, the economic model considers as objective function the soil expectation value (SEV, [28]): SEV(u) = R(u)−C (1+r)tnt+1−1, (1) where R(u) and C are the discounted revenues and costs, and r is the interest rate. For each timber product considered, revenues were computed as the discounted product of stumpage prices and extracted volumes at each thinning and at final harvest: R(u) = nt+1 ∑ i=1 1 (1+r)ti na ∑ j=1 pj iVr,dj i(u)−Vr,dj+1 i(u)!(2) where na is the number of different timber products, Vr,dna+1 i(u) = 0, Int+1=Rnt+1= 1, and pj i is the stumpage price ( e /m 3 ) for product j in the i -th cut. If pj is the stumpage price at clearcutting, we assume a depreciation in thinning price due to its lower intensity and, maybe, lower removed relation, such that pj i=pj a(2−Ri)(1−Ii), where a>1 is a parameter that measures the stumpage price depreciation in thinnings. According to all previous considerations, a simulator for computing the SEV corresponding to each management prescription was developed (see Figure 1). In addition, economic and logistic constraints were considered, and upper and lower bounds for the decision variables were set, determining the admissible set Unt⊂R3nt+1 of possible values of u . Therefore, the forest stand management problem was formulated as the following Mixed-Integer Nonlinear Problem (MINLP): max SEV(u)(3) subject to u∈Unt, (4) 0≤nt≤ntmax, (5) where ntmax is the maximum number of thinnings allowed. 2.2. Transition Functions and Parameters As transition functions for estimating the time-dependent changes in the state variables we used the dynamic equations developed by Diéguez-Aranda et al. [29] and CastedoDorado et al. [30] for radiata pine in this region: h(t0,H0,t) = H01−exp(−0.06738t) 1−exp(−0.06738t0)1.755+12.44/A1, (6) Forests 2021,12, 899 5 of 14 with A1=0.5log H0+1.755A2+q(log H0+1.755A2)2−49.76A2, A2=log(1−exp(−0.06738t0)); n(t0,N0,t) = (N−0.3161 0+1.053t−100 −1.053t0−100)−1/0.3161; (7) g(t0,G0,t) = exp(A3)exp−(−276.1 +1391/A3)t−0.9233, (8) with A3=0.5t−0.9233 0−276.1 +t0.9233 0log(G0) + q5564t0.9233 0+276.1 −t0.9233 0log(G0)2. The initial values for the before-thinning state were set as follows: (i) for h(t0 , H0 , t) , t0= 20, which is the reference age for the species according to Diéguez-Aranda et al. [29] , and H0 is the site index ( S ); (ii) for n(t0 , N0 , t) , t0= 0 and N0 is the plantation density; and (iii) for g(t0,G0,t),t0=10 and G0is given by G0=exp(A4)exp−(−276.1 +1391/A4)10−0.9233, (9) with A4=4.331S0.03594 −114.3 n(0, N0,10). Thus, forest productivity was implemented within the simulations through S , which affects the initial basal area and dominant height. The inputs of forest productivity used in this study are explained in Section 2.3. The output function for the merchantable timber volume was [31] v(H,N,G,d) = 0.4046H1.013G0.9776 exp(−0.2933D−2.818 gd3.192), (10) where the quadratic mean diameter is given by Dg=100r4G πN. The necessary information for implementing silvicultural interventions and economic parameters within the simulations was provided by the Spanish forest consultancy company CERNA SL. The proposed management program comprised the initial plantation, scrub clearance, and low pruning. Considering that plantation densities for radiata pine use to vary from 800 to 1100 stems/ha in this region, we chose as initial value for N0 the mean of that interval, 950 stems/ha. The timber product considered were chip and pulpwood, sawlog and rotary veneer. The cost of management interventions and stumpage prices by product at clearcutting are shown in Table 1. Concerning thinnings, we executed simulations for management prescriptions of one and two thinnings ( ntmax = 2). The constraints for thinning-related decision variables mentioned in Section 2.1 were set as (minimum-maximum): 15%–45% for intensity ( I ), 0.35 (thinning for below)-1 (systematic thinning) for removal relation ( R ), and 10–60 years for timing, with a minimal interval of five years between cuttings. Moreover, for discounting the estimated revenues, we compared the results for interest rates of 1%, 3% and 5%. Forests 2021,12, 899 6 of 14 Table 1. Economic parameters used for the simulations. Description Value Costs Plantation (t=0) 1300 e/ha Scrub clearing (t=3) 450 e/ha Scrub clearing (t=10) 450 e/ha Low pruning (t=10) 750 e/ha Prices Chip and pulpwood (d1= 7 cm) 16 e/m3 Sawlog (d2= 16 cm) 24 e/m3 Rotary veneer (d3= 25 cm) 30 e/m3 Stumpage price depreciation parameter in thinnings 2 2.3. Future Forest Productivity Predictions of S for future climatic scenarios were obtained from a Support Vector Regression (SVR, [ 32 ]) model derived from a previous project [ 33 ]. The model was developed with data from the research plots network established by the Sustainable Environmental and Forest Management Unit (UXAFORES) of the University of Santiago de Compostela, Spain. As predictors of forest productivity, the model incorporates four variables from the Wordclim 2 bioclimatic dataset [ 34 ]: mean annual temperature, mean diurnal range (mean difference between maximum and minimum daily temperatures), isothermality ratio (proportion between mean diurnal range and maximum difference between mean monthly temperatures) and mean temperature of the coldest month. For predicting future S , we used the projections of these four variables developed by the Worldclim project [ 35 ] for the period 2041–2060. Specifically, we used the downscaled projections of a set of 11 Global Climate Models (GCMs) included in the Coupled Model Intercomparison Project Phase 5 [ 36 ] for the Representative Concentration Pathways (RCPs) 4.5 and 6.0. These pathways represent the forecasted climate dynamics under the assumption of reaching a radiative forcing (proportion of incident solar irradiance and radiated energy from Earth) of 4.5 W/m2 and 6 W/m2, respectively, by the year 2100. The future climatic projections and forest predictivity predictions were obtained for a set of 128 radiata pine stand locations in the northwest of Spain for which current productivity data were available. 2.4. Numerical Resolution and Analysis Taking into account that only a few thinnings are allowed ( ntmax is small), MINLP (3)–(5) can be solved by exhaustive search on the integer variable nt . Therefore, for each nt= 1, . . . , ntmax , the nonlinear problem (NLP) (3)–(4) is solved with the fixed value of nt , and the best of the ntmax obtained solutions is taken as the optimal solution of problem (3)–(5). Concerning to NLP (3)–(4), due to the smoothness of transition functions (6)–(8) and the output function (10), the differentiability of the objective function (1) can be proved (see [ 26 ]), and gradient-type methods can be used for solving the problem. In this study, we used the most widespread family: Sequential Quadratic Programming (SQP). Specifically, we used the implementation in the nloptr package [ 37 ]ofR[ 38 ], with a random multi-start and computing gradients by a finite difference method [ 39 ]. To speed up calculations, we did parallelization with the doParallel package [40]. The optimizations were executed for the 128 locations assuming ntmax = 2. For each location, 22 (11 GCMs and two RCPs) S predictions and three interest rates were considered, leading to a total of 8448 MINLPs (optimization scenarios), i.e., 16,896 NLPs. Finally, we analyzed the empirical distribution of SEV for each location and computed the expected shortfall (ES, [ 41 , 42 ]), a financial risk indicator preferred to other metrics (e.g., the value at risk) for non-normal distributions [ 43 ]. Following the definition given by Pfaff [44] , we computed ES for a confidence level = 1 −αas ESα=1 αZα 0qu(FSEV)du, (11) Forests 2021,12, 899 7 of 14 where qu(FSEV) is the quantile function of the SEV distribution. This indicator can be interpreted as the mean SEV below the quantile defined by α , which we fixed as 0.025 in this study, and that is equivalent to the mean financial loss above the 97.5% threshold, according to the nomenclature used by Pfaff [44] . For discussing the potential sensitiveness of SEV to favorable future climate scenarios, we also computed the symmetric of ES2.5, ES97.5. Finally, we compared the SEV estimated for current productivity (assuming no variations in future S , i.e., no climate change), SEV CP , with the mean SEV for climate change scenarios ( SEVCC ), ES 2.5 and ES 97.5 . In addition, we also tested the economic effect of not adapting silviculture to climate change, i.e., to apply the optimum silvicultural programs for current productivity to RCP 4.5 and RCP 6.0 scenarios (hereafter, “climateinsensitive” silviculture ). 3. Results 3.1. Productivity and Economic Indicators The considered future climate models forecasted, on average, an increase in the four S climatic predictors except for the isothermality, which experienced a slight decrease. The most notable shift in climatic variables was the mean temperature of the coldest month, which increased ∼ 40% with respect to previous conditions. However, these forecasts varied notably over climate models, being the mean temperature of the coldest month the sparsest for both RCPs (relative dispersion ∼ 15%). The forest productivity predictions derived from these climatic projections revealed a decreasing trend in mean S under climate change. The mean S reduced from 20.8 m (observed productivity) to 18.8 m (RCP 4.5) and 17.3 m (RCP 6.0). Moreover, the variability of these predictions increased notably, with S ranges (min.–max.) of 7.9–32.1 m for RCP 4.5 and 7.1–29.4 m for RCP 6.0 that contrast the observed range of 12.8–27.7 m. The resulting SEV under optimum silviculture for current productivity varied from − 1150 e/ ha (for r = 0.05) to ∼ 52,000 e/ ha (for r = 0.01), with a mean value of 12,800 e/ha. Concerning climate change scenarios, the SEVCC ranged, overall, from − 750 to ∼35,000 e/ha with a mean value of 10,500 e/ ha for RCP 4.5 and 9000 e/ ha for RCP 6.0. With regard to the tails of the SEVCC distribution, the ES 2.5 varied from − 1800 e/ ha to 27,000 e/ ha, while the ES 97.5 varied from 130 e/ ha to 56,000 e/ ha. Plots of SEV under current productivity vs. SEV under climate change for each interest rate-RCP combination are shown in Figure 2. As shown in Figure 2, in most of the locations SEVCC <SEVCP (points are mainly located to the right side of the identity line). In other words, in most locations, simulations under climate change led to a drop in profitability in comparison with the scenario of current productivity. The average relative decrease in profitability from current productivity to climate change scenarios varied in the ranges 15–64% for RCP 4.5 and 22–89% for RCP 6.0. These wide ranges of variation were mostly driven by interest rates, being the highest decreases in profitability associated with high values of r . Increases in SEV from current productivity to climate change (i.e., where SEVCC >SEVCP ) were scarce and mostly found in some locations with current low-average productivity. A high degree of correspondence was found between SEVCC and ES 2.5 , meaning that locations with higher SEVCC had also high ES 2.5 values. However, there were cases with high SEVCC and low ES 2.5 , and vice versa, which account for the varying dispersion rates of SEV among the different locations. The evolution of SEV from current productivity to SEVCC , ES 2.5 and ES 97.5 under RCPs 4.5 and 6.0 is shown in Figure 3. Altogether, a descending trend in SEV was noticed in the direction current productivity-RCP 4.5-RCP 6.0. A slightly decreasing trend was also found in ES 2.5 between RCP 4.5 and RCP 6.0. The estimated relative decreases of profitability based on ES 2.5 , from current productivity to climate change scenarios, were of 47–142% (also varying with r ) for RCP 4.5 and 55–156% for RCP 6.0. In contrast, ES 97.5 revealed an increase in SEV under climate change scenarios, with relative values of up to 40% for RCP 4.5 and 47% for RCP 6.0. Forests 2021,12, 899 8 of 14 Figure 2. Scatterplots of SEV under current productivity ( SEVCP ) vs. mean SEV under climate change for the two RPC’s and the three interest rates considered. The dashed line in each plot represents the identity. Forests 2021,12, 899 9 of 14 Figure 3. Parallel coordinates plots representing the change from SEV under current productivity ( SEVCP ) to ( A ) mean SEV under climate change (SEVCC), (B) future ES2.5 and (C) future ES97.5 for an interest rate = 0.03. 3.2. Optimum Silviculture Concerning silvicultural decision variables, the optimum number of thinnnings was one ( nt = 1) in all the optimization scenarios (MINLPs). However, the differences in SEV between optimum silvicultural programs of one and two thinnings were very slight, being the mean difference = 180 e/ ha. Climatic scenarios had a noticeable influence on the optimum rotation length, which tended to reach higher values under climate change (RCP 6.0 > RCP 4.5, in most of cases) in comparison to current productivity (Figure 4). As expected, the rotation length also showed a negative correlation with r . The thinning intensities and removal relations were scarcely influenced by climatic scenarios and, overall, experimented low variability. Concerning the first thinning (or the only thinning when nt = 1), the results of most of NLPs yielded values very close to the lower bounds of these decision variables, implying low intensities ( I1∼ 0.15) and thinning from below ( R1∼ 0.35). The optimum intensities of the second thinning (for those NLP with nt = 2) had a broader variation range, with some of the NLPs reaching I2∼ 0.45. As with the rotation length, the optimum thinning timings suffered a noticeable variation across optimization scenarios, especially influenced by r . The mean timing for the first thinning was, overall, 19 years and the mean timing of the second thinning was 31 years.