Use of Optimization Modeling to Assess the Effect of Timber and Carbon Pricing on Harvest Scheduling, Carbon Sequestration, and Net Present Value of Eucalyptus Plantations

: Quantifying the impact of carbon (C) and timber prices on harvest scheduling and economic returns is essential to deﬁne strategies for the sustainable management of short-rotation plantations so that they can provide timber products and contribute to C sequestration. In this paper, we present a mixed-integer linear programming model that optimizes harvest scheduling at the forest level, C sequestration, and Net Present Value (NPV) over a planning period of up to 15 years. The model included revenue from the sale of timber (pulplogs) and credits from the net C sequestered during the life of the stands. In addition, plantation establishment, management, harvesting, and transportation costs were included in the analysis. The study area comprised 88 Eucalyptus grandis W. Hill and Eucalyptus dunnii Maiden stands located in Uruguay, totaling a forest area of nearly 1882 ha. The study investigated the impact of C and timber prices on NPV, harvest schedules, stands’ harvest age, timber ﬂows to customers, and C sequestered per period. The maximum NPV among all the scenarios evaluated (USD 7.53 M) was calculated for a C price of 30 USD t − 1 , an interest rate of 6%, and a timber price of 75 USD m − 3 . This was USD 2.14 M higher than the scenario with the same parameters but that included only revenue from timber. C prices also impacted stands’ harvest age, C sequestration, and timber ﬂows delivered to end customers. On average, in scenarios that included C prices, timber ﬂows and C sequestration increased by 15.4 and 12.1%, respectively, when C price increased from 5 to 30 USD t − 1 . These results demonstrate that harvest scheduling, harvest age, and NPV are very sensitive to C and timber, and that the best economic returns are obtained when the stands are managed to maximize timber production and C sequestration. the average amount of C sequestered in periods 7, 8, and 9 was 74,909, 52,326, and 57,959 t, respectively. The above results show clearly that the optimization model aimed to maximize C sequestration at higher C prices and the production of timber at lower C prices.


Introduction
Forest ecosystems constitute the largest terrestrial carbon (C) pool in the global carbon cycle and represent one of the most significant sources to mitigate climate change effects [1]. The concentration of atmospheric carbon dioxide (CO 2 ) has significant effects on tree growth and carbon accumulation in forest ecosystems [2]. Whether forest ecosystems sequestrate more carbon into ecosystems (carbon sinks) or release more carbon from ecosystems (carbon sources) depends on the net ecosystem carbon exchange. Forests have a significant impact on the concentration of CO 2 in the atmosphere. As terrestrial ecosystems, forests act as C sinks, positively affecting the atmospheric CO 2 balance [2]. Carbon storage in forest ecosystems includes numerous components but particularly biomass and soil [3]. only; in the latter case, average rotations in intensive silviculture stands occur at the age of eight years. Additionally, land price has been increasingly growing in Uruguay in the last decades, which has led forest growers to explore strategies to increase forest productivity and efficiency along the value chain (including harvesting and transport operations), and in this way, improve the competitiveness of the forest sector [20]. As part of the tactical plan and based on forest inventory data, forest companies usually prescribe short forest rotations, no longer than eight years, as a way to maximize yield per ha [12]. Traditionally, methods to determine optimal rotation or harvest ages have not been frequently employed in Eucalyptus plantations that are managed to produce roundwood for the pulp industry under operational-constrained scenarios [21]. However, this has been changing ultimately. Most forest growers are now considering determining optimal rotations for their stands to maximize economic returns and meet the growing demand for wood in the coming years.
The forest planning problem at the stand and forest level, based on both timber production and C sequestration, has been addressed previously by several authors [21][22][23]. C pricing can result in changes to the optimal stand rotation or harvest age, and consequently, to C stocks. This is because forest stands with longer rotation or harvest ages build up a larger volume of biomass and soil organic C than forests managed on shorter harvest cycles [24]. In Eucalyptus plantations in Australia, Enzinger and Jeffs (2000) [25] showed that the addition of C payments resulted in longer optimal rotation ages than scenarios where stands were managed for timber production only. As indicated by Asante et al. (2011) [24], the maximization of NPV for a single harvest cycle is a standard criterion for determining the optimal stand rotation and harvest age. NPV is maximized when the capital value of standing timber (the percentage by which the value of the forest is growing) is equal to the interest rate (IR), assuming flows of revenues and costs in just one harvest cycle and not accounting for these flows in subsequent cycles [26]. From an economic point of view, the optimal rotation age (optimization at the stand level) and harvest age (optimization at the forest level) of single stands at the tactical level is the one that maximizes NPV for a flow of revenues and costs throughout the planning period, subject to supply and operational constraints. This is opposed to the biological optimal rotation that is obtained when MAI is maximized, which excludes the flow of revenues and costs involved in growing the forest stands [27].
We hypothesize that C and timber prices have an impact on harvest scheduling, C sequestration, and economic returns (NPV) in industrial Eucalyptus plantations located in Uruguay. The objectives of the study were as follows: (i) formulate a mathematical model able to optimize the harvest of stands in Eucalyptus plantations, including timber (pulplogs) and C stocks as a source of revenue, (ii) assess the impact of C and timber prices on harvest age at the forest level and NPV, and (iii) assess the impact of C and timber prices on spatial solutions and flow of timber from stands to pulp mills. We expect that the results from this study will help forest growers in Uruguay get a better understanding of the impact of C and timber prices on harvest scheduling, orienting the management of Eucalyptus plantations in the context of a growing global demand for timber and C sequestration.

Study Location and Characteristics
The study area comprised about 1882 ha of Eucalyptus grandis and Eucalyptus dunnii plantations owned by Forestal Oriental SA in Uruguay. Figure 1 shows the five sites included in the study: sites b1 and b2 with 861 ha (43 • 13 97" S-64 • 24 78" W, hereafter Guichon); sites b3 and b4 with 870 ha (63 • 07 44" S-52 • 22 64" W, hereafter Trinidad), and site b5 with 717 ha (63 • 44 33.05" S-63 • 95 58" W, hereafter Sarandí del Yi). Climate is subtropical, with a mean annual temperature of 18 • C (12 • C in the coldest month, 24 • C in the warmest month) and a mean annual rainfall between 1300 and 1400 mm [28]. According to the National Commission for Agroeconomic Studies of the Land classification (CONEAT), the predominant soil in Guichon is planosols with a horizon A of 40-50 cm depth, characterized by low fertility, weak structure, low level of organic matter, slopes of 1-3%, sandy texture, medium to low risk of drought, imperfect drainage, moderately-slow to slow permeability, and good rooting ability. Soils in Trinidad are vertisols with a horizon A of 50-60 cm, and characterized by medium fertility, weak structure, medium level of organic matter, slopes of 2-5%, loam franc structure, low risk of drought, moderate drainage, moderately permeability, and good rooting ability. Soils in Sarandí del Yi correspond to vertisols, rupticos, and brunosols, with a horizon A of 50-60 cm depth and characterized by sandy-loam, loam franc texture, low fertility, low risk of erosion, moderate slopes (1-3%), weak structure, low organic matter, imperfect drainage, and good rooting capacity.

Stands Information
A forest inventory including a single 10-m fixed radius (314.16 m 2 ) plot per ha was established in the study area. In each plot, measurements included diameter at breast height (1.3 m above ground level -DBH-cm, measured with a calliper Haglöf Mantax, Långsele, Sweden), stand density (N, trees ha −1 ), and basal area (G, m 2 ha −1 ) of all trees with DBH ≥ 7 cm. Total height (H, m) was measured in the two central rows of the plots using a Vertex III hypsometer (Haglöf, Sweden); the height of the remaining trees was estimated using a logarithmic classic model ( Table 1). As shown in the table, the stands were inventoried in real life when their ages ranged between 4.5 and 9.5. This information was used to project their growth from year 0 (plantation establishment) to year 10. For example, if a stand was inventoried in year 5, the data were used to project the state variables backward (from year 0 to year 5) and forward (from year 5 to year 10) through differential equations [29,30]. State variables projected included number of trees per hectare (N), basal area per hectare (G), quadratic mean diameter (dg), and mean arithmetic heights (Hm).

Carbon Sequestration and Release Assumptions
The process of photosynthesis permits green plants to uptake CO 2 from the atmosphere and convert it into organic carbon as they grow, and in turn, organic carbon is converted back to CO 2 when it is eaten or decomposed, which is known as the process of respiration [31]. In our study, the total C sequestered annually and stored in trees as organic C was computed following the procedure presented by Díaz and Rodríguez (2006) [21], and Toochi (2018) [31]. For this purpose, a basic density of 0.5 t m −3 was considered for both E. grandis and E. dunnii., as well as a C content coefficient of 0.5 (C content equal to 50% of the tree's total volume). Thus, to calculate the C content of the species, a 0.25 conversion factor (wood density x proportion of C in the tree) was used. Only above ground carbon (AGC) was considered, assuming that in short-rotation plantations of Eucalyptus, the soil organic C does not have significant variations [30]. The aboveground biomass (AGB) required to calculate AGC was obtained using a compatible system of equations developed by Hirigoyen et al. (2020) [32] as functions of dg and Ho ( Table 1). The rates of change in total AGB stocks were estimated by calculating the expected annual growth of the stand. Two cases were assessed considering the proportion of C released at harvest time and during the periods after the harvest. In the first case, it was assumed that 100% of the C stored is released at harvest time. In the second case, it was assumed that 50% of the C stored is released at harvest time, while the other 50% is stored after harvest as residues and paper products, and released linearly over the next 5 years after the final harvest [21].

Model Description and Mathematical Formulation
To optimize harvest scheduling and stand harvest age, a mixed-integer linear programming model was developed to maximize the discounted net revenue (NPV) over a planning period of up to 15 years. Eucalyptus stands in Uruguay are usually managed for periods of up to 10 years, with a minimum harvest age of 7 years (periods); this was taken into consideration in our study by adding a constraint to the optimization model. The objective function included revenue from timber (pulp logs) and C sequestered in the stand each year until the harvest age. It also included costs associated with plantation establishment, weed control, fertilization, and pest control, as well as timber harvest and transport to two alternative pulp mills: Paso de Los Toros (hereafter P. Toros) and Fray Bentos (hereafter F. Bentos). In each period t, the revenue generated from timber sales as well as plantation, harvest (clearfelling), and transport costs were discounted at an interest rate IR. A discount factor per period was calculated assuming the flow of revenues and costs occurred at the midpoint (t-0.5) of the period. The yearly C sequestered by each stand, from establishment through harvest age, was accounted as a revenue source. Additionally, using similar approaches presented in previous studies [21,33], we calculated the value of the C released back to the atmosphere at harvest time and during the years after the harvest using the same C price as that used to quantify the revenue from C sequestration. Both C sequestered and released were discounted in the same way as the other revenues and costs Sets, parameters, and variables of the optimization model are presented in Table 2. All costs, prices, and revenues in the manuscript are expressed in US dollars. The optimization, mixed-integer programming model can be expressed mathematically as follows: Harvesting cost in stand i ∈ I in period h ∈ H (USD) SPPC i Site and plantation costs in stand I ∈ I (USD) WCFA it Weed control, fertilization, and ant control costs in stand i ∈ I in period t ∈T (USD) Carbon sequestered in stand i ∈ I in period t ∈ T (t) AccC ih Carbon accumulated in stand i ∈ I that is harvested in period h ∈ H (t) ir Interest rate (decimal) Decision variables X ih Binary variable. 1, if stand i ∈ I is harvested in period h ∈ H, 0 otherwise PLD ijh Continuous variable. Pulp logs delivered from stand i ∈ I to pulp mill j∈ J in period h ∈ H (t)

Accounting variables
RevenuePLStand i Revenue from the sale of pulp logs generated in stand i ∈ I (USD) ValueCSStand i Value of carbon sequestered by stand i ∈ I (USD) ValueCS iht Value of carbon sequestered by stand i ∈ I in period t ∈ T, if the stand is harvested in period h ∈ H (USD) ECostStand i Establishment cost in stand i ∈ I (USD) HTCostStand i Harvesting and transport cost in stand i ∈ I (USD) ValueCRStand i Value of carbon released from stand i ∈ I (USD) ValueCR iht Value of carbon released from stand i ∈ I and period t ∈ T, if the stand is harvested in period h ∈ H (USD). ValueCRStandA i Value of carbon released from stand i ∈ I (USD) at the harvest time ValueCRStandB i Value of carbon released from stand i ∈ I (USD) over the 5 years after the harvest

Objective Function
The objective function maximizes the NPV over the planning period (Equation (1)): where, To compute the discounted revenues and costs by period, a discount factor δ t (Equation (4)) was used: Revenue and cost equations Equations (5) and (6) were used to compute the discounted revenue from the sale of pulplogs in the stands.
Equations (7)-(9) were used to compute the discounted revenue from the carbon sequestered in the stands.
Equations (10) and (11) were used to compute the discounted costs associated with plantation establishment in the stands.
Equations (12) and (13) were used to compute the discounted harvesting and transportation costs in the stands.
Several equations were used to compute the value of the C released at the harvest time. Equations (14) and (15) were used for the scenario where it was assumed that 100% of the C accumulated was released at the harvest time.
Likewise, Equations (16)- (19) were used for the scenario that assumed that 50% of the C accumulated was released at the harvest time, and the other 50% was released over the next 5 years after the harvest.

Constraints
In addition, three groups of constraints were included in the model to ensure that (a) every stand was cut once within periods 7 to 10 (Equation (20)), (b) no stand could be harvested before period 7 (minimum age to harvest a stand) [12,18] (Equation (21)), and (c) the volume of timber delivered to the pulp mills did not exceed the volume available of these products in the stands (Equation (22)). A constraint in Equation (20) was included since trees harvested at an early age (≤ 6 years) produce wood with a very low pulp yield, and there is a decline in marginal growth when trees are harvested after period 10.

Data Input
The area of the 88 forest stands used for the analysis ranged between 5.5 and 54.3 ha (average = 21.9 ha). Site and stand variables for the stands at the beginning of the planning period are presented in Table 1. Parameters associated with C stocks were calculated as per the description provided in Section 2.3, while a range of C prices (from 0 to 30 USD t −1 ) was used for the analyses. These C prices correspond to the price per ton of organic C stored in trees resulting from CO 2 sequestration. The carbon market of the European Union Emissions Trading System (https://sandbag.org.uk/carbon-price-viewer/ (accessed on 15 December 2020) was consulted to define the C prices used in the study.
For the analysis, the base scenario included an IR of 6%, a C price of 0 USD t −1 , and a delivered price for pulp logs of 55 USD m −3 . Unit harvesting costs ranged between 17 and 21 USD m −3 as a function of tree volume and harvesting age (periods 7, 8, 9, or 10), assuming a typical cut-to-length (CTL) system for Eucalyptus plantations in Uruguay [34].
The CTL system comprises one or more wheeled harvesters/processors that perform tree felling and stem processing tasks, followed by forwarders that extract the processed logs to the roadside. Based on previous studies [35], a single transport cost of 0.1 USD t-km −1 was used. Distances from stands to pulp mills were calculated with ArcGis 10.1 (ESRI, 2011); they ranged from 152 to 337 km (average of km) to F. Bentos pulp mill, and from 107 to 191 km (average of 157 km) to P. Toros pulp mill. In Uruguay, road transport is usually performed by 6-axis trucks with a gross vehicle mass of 48 tons and an average tare of 13 tons (net payload of approximately 33 tons).
In addition to harvesting and transport costs, the NPV calculation included costs associated with site preparation and plantation (year 0 at a cost of 800 USD ha −1 ) as well as weed control (periods 1 to 3 at a cost of 100 USD ha −1 , and fertilization and ant control between periods 1 and 3 at a cost of 95 USD ha −1 ) [34].

Model Implementation
The optimization model was implemented using the IBM DocPlex library v. 2.19.202 ( https://pypi.org/project/docplex/ (accessed on 15 December 2020), licensed under the Apache License v2.0 within the IBM Decision Optimization CPLEX Modeling for Python program. The model was run on a computer Dell Inspiron 5530 (Dell Inc., Texas, USA) equipped with an Intel®Core™ i7-8850H, 2.60GHz (6 cores) processing unit, and 32 Mb of RAM. The model script was coded using the PyCharm editor v.2020.3.3 (Professional Edi-tion), and to facilitate the runs of the optimization model, a simple GUI was programmed with the PyQt5 Library (https://pypi.org/project/PyQt5/ (accessed on 15 December 2020). All the input data were imported from Excel and converted into Pandas data frames for later use in the optimization model. After each run, all the model outputs were exported back to Excel for further analysis.

Scenarios and Sensitivity Analysis
In total, 54 scenarios were assessed and their NVPs calculated and compared. These scenarios included the combinations of three IRs (6, 8, and 10%), three timber prices (55, 65, and 75 USD m −3 ), and six C prices (0, 5, 10, 15, 20, and 30 USD t −1 ), which were used to conduct the sensitivity analyses. The above scenarios were grouped into two classes to facilitate the analyses and comparisons: 1.
Scenarios that included revenue only from the sale of timber (pulplogs) generated during clear-cut (SC1); 2.
Scenarios that included revenue from both sources, timber generated during the clear-cut, and C sequestered over the rotation length of the stands (SC2).
Sets of outputs combining harvest periods (periods 7 to 10) and C prices included timber revenue and sequestered C and its revenue, along with harvest and transport costs to both pulp mills. The outputs also included the optimal harvest age (final cut) of each stand and the corresponding value of the objective function (NPV). All statistical analyses and summaries were performed in Python, using the NumPy and Pandas libraries. All statistically significant differences between scenarios were computed using a paired t-test (p ≤ 0.05). In addition, maps with the optimal spatial solutions were created using ArcGIS®software by Esri. Figure 2 shows the relationship between NPV and C price for the range of IRs, timber prices, and C release approaches (50 or 100% of C released at the harvest time) included in the study. As expected, in SC2 scenarios, NPV was sensitive to C prices and IR. There was a clear relationship between NPV and C price, with NPVs that, on average, increased at a rate of 0.55 (IR = 6%), 1.38 (IR = 8%), and 1.61 M USD (IR = 10%) for every 5 USD t −1 rise in C price. The maximum NPV among all the scenarios evaluated (2.52 M USD) was calculated for a C price of 30 USD t −1 , an IR of 6% and a timber price of 75 USD m −3 . In addition, there was a substantial reduction in NVP when the IR increased from 6 to 10%. For example, for a C price of 15 USD t −1 , a timber price of 65 USD m −3 , and a 100% C release approach, the NVP was 3.75, 2.73, and 1.99 M USD for an IR of 6, 8, and 10%, respectively. The differences in NPV between the three IRs were statistically significant at p ≤ 0.05. For SC2 scenarios that included an IR of 8 and 10%, and a timber price of 55 USD t −1 , NPVs became positive only after a minimum C price. For an IR of 8%, this minimum C price was 5 USD m −3 irrespective of the C release approach being used. For an IR of 10%, the minimum price was 10 USD m −3 for a 50% C release approach, and 15 USD m −3 for a 100% C release approach.

Effect of C and Timber Price, IR, and C Release Approach on NPV
On average, using a 50% C release approach resulted in NPVs that were 8 (IR = 6%), 28 (IR = 8%), and 33% (IR = 10%) higher than using a 100% C release approach. These differences were statistically significant at p ≤ 0.05, excepting for IR between 8 and 10%. It is worth noting that differences in NPV between the 50 and the 100% C release approaches occurred in SC2 scenarios only, since in SC1 scenarios all the revenue was obtained solely from the sale of timber. The largest difference in NPV between the two C release approaches was obtained for a C price of 30 USD t −1 . At this C price and a timber price of 65 USD m −3 , the differences were 0.39, 0.43, and 0.43 M USD, for an IR of 6, 8, and 10%, respectively.

Effect of C Price and Timber Price, IR, and C Release Approach on Revenue
The impact of C price on revenue is presented in Table 3, considering a timber price of 65 USD m −3 , and a 50% C release approach. The results are broken down by IR, total revenue, revenue from C, and revenue from timber. As expected, total revenues increased with higher C prices and lower IRs. The minimum revenue (11.3 M USD) was obtained in an SC1 scenario with an IR of 10%, where 100% of the revenue resulted from the sale of timber. On the contrary, the maximum revenue (23.1 M USD) was obtained in an SC2 scenario with an IR of 6%, where 70% of the revenue resulted from the sale of timber and 30% from C credits.
Since prices were considerably higher for timber than for C credits, regardless of C price and IR, the revenue from timber in SC2 was always higher than the revenue from C credits. The largest gap (88% difference of contribution to total revenue or 15.1 M USD) between revenue from timber and revenue from C credits was obtained for a C price of 5 USD t −1 and an IR of 6%. In contrast, the smallest gap (28% difference of contribution to total revenue or 7.7 M USD) was obtained for a C price of 30 USD t −1 and an IR of 10%. The highest proportional contribution to revenue from C credits was 34% for an IR of 10% and a C price of 30 USD t −1 , whereas the highest monetary contribution was 6.8 M USD and occurred for the same C price and an IR of 6%. Irrespective of the C price and IR, the contribution of C credits to total revenue ranged between 6-7% and 30-34%. In addition, there was a notable variation in the revenue attributed to C credits and timber. Thus, the coefficient of variation for C revenue was 74.7 (IR = 6%), 76.7 (IR = 8%), and 77.1% (IR = 10%), whereas the variation for timber revenue was only 0.5 (IR = 6%), 1.5 (IR = 8%), and 0.6% (IR = 10%), which reflects the impact of C prices on this variation. Table 3. Effect of C and timber price, and IR on revenue. Values (in millions of dollars) are calculated for a timber price of 65 USD m −3 and considering 50% of C released at the harvest time. In parenthesis, it is the percentage of total revenue associated with timber and carbon.   Table 4 presents the volume harvested (m 3 ) by harvest age, C price, and IR. Irrespective of the IR, in SC2 scenarios, the number of stands being harvested in period 10 grew gradually as C prices passed from 5 to 30 USD t −1 , which resulted in higher volumes of timber being harvested and delivered to customers during this period. For example, for a C price of 30 USD t −1 and an IR of 6, 8, and 10%, the proportion of stands being harvested in period 10 was 100 (88 stands), 99 (87 stands), and 75% (66 stands), respectively. From Table 4, it is observed that harvest ages commenced declining when IR passed from 6 to 8 and 10%, which resulted in more stands being harvested (and a higher volume of timber delivered) in periods 7, 8, and 9. On the contrary, for a C price of 0 USD t −1 (SC1 scenarios) and an IR of 8 or 10%, the largest proportion of the stands were harvested in period 7 (minimum harvest age). In SC1 scenarios, due to the absence of revenue from C credits, a larger proportion of the stands was harvested in periods 7 and 8, particularly when the IR was 8 or 10%. The most notable case was that with an IR of 10%, where 44% of the stands (39 out of 88) were harvested in period 7. Harvest ages also had an impact on the total volume being harvested across the planning periods. Thus, for example, in SC2 scenarios with a C price of 30 USD t −1 , a timber price of 65 USD m −3 , and an IR of 6, 8, and 10%, the total volume harvested in all periods was 435,830, 435,372, and 417,993 m 3 , respectively. As a comparison, in SC1 scenarios (C price = 0 USD t −1 ), the volume harvested in all periods was only 411,592, 367,740, and 356,882 m 3 for an IR of 6, 8, and 10%, respectively.

Effect of C Price on Harvest Scheduling (Harvest Age) and Volume Harvested (Supplied) Per Period
The spatial distribution of the stands being harvested per period (Figure 3) was impacted by harvest age, C and timber prices, and IR. Spatially arranged harvest of stands aimed to minimize harvest and transport costs, and maximize returns from timber and C, thereby maximizing NPV. This also had an impact on the flow of timber being delivered from stands to customers (Figure 4). In the case of SC2 scenarios, the peak flow of timber from stands to P. Toros pulp mill occurred at longer harvest ages when C price increased from 5 to 30 USD t −1 , reaching a maximum volume (174,289 m 3 ) in period 10 when C price equaled 30 USD t −1 . On the contrary, in SC1 scenarios, the maximum flow of timber to P. Toros pulp mill occurred during periods 7 and 8, reaching a maximum of 73,807 m 3 in year 8. Regarding timber flows to F. Bentos pulp mill in SC2 scenarios, maximum volumes of 96,529, 129,289, 138,632, 186,163, and 255,142 m 3 , were obtained in year 10 for C prices of 5, 10, 15, 20, and 30 USD t −1 , respectively. In SC1 scenarios, the maximum volume delivered to F. Bentos pulp mill (78,889 m 3 ) was obtained in period 9. In all scenarios that included a C price lower than or equal to 15 USD t −1 , there was a positive volume of timber delivered to F. Bentos pulp mill. For scenarios that included a C price of 20 USD t −1 , there was no timber flow to F. Bento from coupes being harvested in period 8. The same occurred in periods 7, 8, and 9 for a C price of 30 USD t −1 . At this C price, no timber was delivered to any pulp mills since no forest stands were harvested in these periods.

Effect of Carbon Price and IR on Carbon Sequestered
As observed in Figure 5, both the total and per period C sequestered were shown to be very sensitive to C prices, and highly affected by the harvest age. In all scenarios, the maximum amount of C sequestered occurred in period 10 for a C price of 30 USD t −1 . At this C price, most stands were harvested in period 10 to maximize revenue from C credits; amounts of C sequestered totaled 293,196 t for an IR of 6%, 289,196 t for an IR of 8%, and 203,586 t for an IR of 10%. The amount of C sequestered in period 10 declined with lower C prices, which was more evident in scenarios with an IR of 8 and 10%. In the latter case, the C sequestered dropped by 53.9% when C price was reduced from 30 to 20 USD t −1 . With an IR of 6%, the C sequestered declined gradually when C price was reduced from 30 to 0 USD t −1 . Likewise, when passing from 30 to 10 USD t −1 , the drop in C sequestered occurred at an average rate of 9295 t per every 5 USD t −1 reduction in C price, while when passing from 10 to 0 USD t −1 , the drop was about 46,145 t per every 5 USD t −1 reduction in C sequestered.
For C prices lower than 30 USD t −1 , the amount of C sequestered in periods 7, 8, and 9 showed a gradual increase when IR passed from 6 to 10%. Thus, for an IR of 6%, the average amount of C sequestered in periods 7, 8, and 9 was 15,590, 3892, and 25,310 t, respectively. Likewise, for an IR of 8%, the average amount of C sequestered in periods 7, 8, and 9 was 41,785, 38,159, and 66,703 t, respectively. Finally, for an IR of 10%, the average amount of C sequestered in periods 7, 8, and 9 was 74,909, 52,326, and 57,959 t, respectively. The above results show clearly that the optimization model aimed to maximize C sequestration at higher C prices and the production of timber at lower C prices.

Discussion
The results obtained in this study show that paying for C sequestered can increase economic returns in short-rotation plantations depending on C and timber, as well as logging and market conditions. However, the implementation of forest planning that includes C sequestration as a source of revenue in Eucalyptus plantations requires regulatory changes at the international and national level [36]. Optimizing harvest schedules and harvest age of stands considering C and timber prices is key to maximize the use of the forestry resources that are subject to intensive silviculture schemes [8]. The results show that NPVs were quite sensitive to the range of C prices evaluated in the study. For any C price, adding up the revenue from C credits and the sale of timber gave the best economic returns. High C prices increased the proportional contribution of C credits to NPV and reduced the contribution of timber, as found in previous studies [8][9][10]37]. Growth rates were used to compute C sequestered and stored in aboveground biomass, and a sensitivity analysis was conducted to quantify the effect of IR as well as C and timber prices on NVP and optimal harvest scheduling (including stands' harvest age). Following similar approaches to those from previous studies [21,33], gross revenue in each stand was computed from the sale of timber as well as from credits associated with C sequestration due to forest growth.
In addition, C released at the harvest time and costs associated with the management of the stands were used to quantify NPV. Taking all these considerations into account, our optimization modeling approach proved to be a practical and effective method to determine the impact of C and timber prices on harvest scheduling and economic returns in Eucalyptus short-rotation plantations.
The harvest schedule provided by our solution approach maximized the revenue from timber and C sequestration in each stand. This plan was generated by using a mathematical programming model that optimized harvest scheduling between periods 7 and 10, and which served as the basis to program new plantations. As confirmed in our study, optimal planning models are critical to maximize the economic and silvicultural outcomes of forest management plans [33,38]. Additionally, determining optimal harvesting scheduling and harvest age, including multiple products and services from forest ecosystems, is essential for implementing sustainable forest management strategies [39]. The production of Eucalyptus in Uruguay, however, has a clear orientation for pulp, so other products were not contemplated in our study (e.g., sawn wood), since no markets exist for these products in the country. Optimal biological rotations are dictated by the annual rate of C accumulation in tree biomass. For example, for short-rotation Eucalyptus plantations in subtropical China, Zhou et al. [40], showed that the marginal annual rate of C accumulation increased from year 4 to year 10, and sharply decreased from year 13 to year 21, suggesting optimal rotations between 12 and 15 years. However, economic and biological optimal rotations determined at the stand level usually differ from harvest ages determined at the forest level, since the former is driven by several operational constraints and affected by management costs and the prices of multiple products, such as C and timber.
Regarding economic optimal rotations, Diaz-Balteiro and Rodríguez (2006) [21], developed several models using dynamic programming to calculate the optimal rotation for intensive Eucalyptus plantations in Spain and Brazil. Their solution approach did not include operational constraints as we did in our study, but it focused primarily on determining the economic rotation age at the stand level using the Land Expectation Value. Regardless of the solution approach, when including the combined revenue from timber and C sequestration, the rotation or harvest age tends to increase when C prices are high, as reported by various authors [33,37,38]; this, however, will depend on the annual rate of C sequestration. In our study, the annual rate of C sequestered increased continuously until period 10, which was set as the maximum harvest age for any stand in the plantation estate. In addition to annual rates of C sequestered, harvest ages were determined by the ratio between timber (pulplogs) and C prices. Thus, the impact of C revenue on NPV was more pronounced when C prices were high and exceeded that of pulplogs; this resulted in longer optimal harvest ages when C prices were higher and shorter optimal harvest ages when C prices were lower. IR also had an impact on the harvest age of the stands. Thus, there was a larger proportion of stands with longer harvest ages when the IR was 6% and a lower proportion when the IR was 8 or 10%. Historically, a harvest age of 8 years has been adopted by most forester growing pulp-oriented Eucalyptus plantations in Uruguay, based primarily on MAIs figures reported in previous studies [41]. The results presented in our study show that harvest ages can range between 7 and 10 years without impacting financial returns and C sequestration. However, the impact of shorter rotations (<7 years) must be further investigated since it might have negative implications on timber supply and economic returns. Additionally, in a context where forests are one of the most important carbon sinks to mitigate climate change, a substantial reduction in the optimal harvest age may impact their C sequestration capacity [42,43]. Hence, it is essential to monitor and accurately predict and quantify changes of C stocks with the assistance of local models [44] and remote sensing technologies [38].
C pricing can be used as a political and planning tool. For example, in a study that involved forests being managed to maximize discounted revenue from both timber production and C, Pukkala (2020) [9] reported that if 100 EUR t −1 of C were paid to forest land owners every 100 years, the C sequestration of Finnish forests would increase by 70% (nearly the maximum increase obtained by subsides). Although these results are consistent with our findings, some slight differences were observed. These can be attributed to the range of forest products (multiple timber products in the case of Pukkala's study, and a single timber product in the case of our study) and C prices. In the case of our study, discounted returns obtained primarily from the sale of timber when C prices were lower than or equal to 15 USD t −1 , and from C sequestration when C prices were higher than 15 USD t −1 . The high volatility of prices in the C market might be associated with uncertain C revenues, so this should be an important variable to be factored into the economic analysis of intensive silviculture strategies [39,40].
Our results reveal that the spatial arrangement and number of stands to be harvested in each period were sensitive to C prices. With low C prices (0 and 5 USD t −1 ), the model aimed to maximize NPV and returns from the sale of timber, which resulted in the harvesting of a large proportion of stands in periods 7 and 8; these stands were associated with lower harvesting and transport costs due to their yield per ha and distance to the pulp mills. On the contrary, with high C prices (20 and 30 USD t −1 ), the model aimed to maximize NPV and returns from C sequestration, which resulted in longer harvest ages (9 and 10 years) and in a large proportion of stands that reached high yields per ha in these periods.
Globally, forest planners are concerned with the sustainability of planted forests. Strategies to increase C sequestration are under consideration by many forest owners and government agencies in their forest management and logging plans [8]. Linear programming economic models of forest management that include C sequestration in intensive Eucalyptus plantations are a powerful tool to optimize their role for climate change mitigation [45]. These optimization-based solution techniques are essential to quantify the economic value and potential contribution of natural forest and plantations to climate change based on integrated returns obtained from multiple ecosystem services and forest products [46]. Silviculture geared to maximize C sequestration promotes the growth and conservation of standing biomass stocks, either by modifying the standard forest management strategies or by optimizing the rotation age [47]. In addition, these goals are achieved through the assessment and control of different silvicultural schemes, including optimal thinning and harvesting strategies, promotion of agroforestry systems, and diversification of forest uses that provide multiple ecosystems services and increase the value of the forests [48].
As far as we know, our optimal modeling proposal is the first effort in assessing the incorporation of C sequestration into forest management plans of Eucalyptus plantations in Uruguay. However, we are aware of some of the limitations of our study. Firstly, the data were collected from clonal trees in three regions in Uruguay for the two most common Eucalyptus species, plus geographical conditions, proximity to pulp mills, harvesting technology, and market conditions may affect the impact of C prices on harvesting plans and NPV. It is expected that expanding the number of forest sites and other genetic sources might overcome these limitations. Secondly, our analysis did not include decisions on initial stocking or the possibility of thinning. Of course, these decisions can be included in future studies to be able to assess their impact on C storage strategies in low-yield plantations. Additionally, our research focused primarily on investigating how C and timber prices impact harvest scheduling (harvest age) and NPV; the assessment quantified C sequestration and storage in each period, and C released at the end of the planning period. It is worth pointing out that since there is no formal C market for plantations in Uruguay, we used a rather conservative range of C prices (from 0-30 USD t −1 ) to analyze the potential impact that this parameter could have on the harvest age, harvest scheduling, and NPV, based on similar articles found in the literature. Additionally, due to the absence of applicable legal mechanisms for C crediting and payment in Uruguay, these models have limited use in the current scenario. On the other hand, our analysis focused more on highlighting the impact of C pricing on forestry planning (space and time) and its importance to Eucalyptus plantation afforestation rather than on comparing payment mechanisms (subsidies and taxes). The optimization model developed and presented in our study is deterministic and generates the same outcomes under a set of initial conditions and value of the parameters (e.g., timber and C prices, C accounting). We are aware of the limitations of this approach, and given the role attributed to forests in climate change policy and in view of potential future policy discussions, it is important to verify and quantify the potential of forest carbon sequestration in the presence of uncertainty [49]. Future research should test and implement alternative solution approaches (e.g., stochastic programming, robust optimization, etc.) that can deal with uncertainty in C sequestration projects, since forest management decisions under uncertainty could deviate from those made with the assistance of deterministic models. In addition, our analysis did not consider the impact of logging and transportation operations on C emissions. However, previous studies have shown that their contribution is small compared to the C contained in wood and biomass products, with ratios of about 40:1 and 45:1 for roundwood and biomass, respectively [50]. According to previous reports, other variables included in the study will also influence forest management and C storage strategies, including IRs, growth rates, and prices of sawn wood and biomass. Our analysis included a limited range of prices and costs, and their fluctuation may affect decisions related to C storage.
Management alternatives, such as silvopastoralism or forest plantations destined for solid wood, were not included. The accumulation and flux of C from these forest management alternatives may be different from those computed in our study. Biomass was not included in our analyses, since in Uruguay, a common practice is to leave the logging debris on site to mitigate the impact of the harvest equipment on the soil.

Conclusions
Sustainable forest management strategies that allow the provision of multiple timber products and enhance C sequestration require the assessment of scenarios to quantify the impact of carbon (C) and timber prices on harvest scheduling and optimal harvest age. Our results show that mixed-integer linear programming is a practical and effective tool for developing optimal forest management plans in short-rotation Eucalyptus plantations. These management plans must be designed and implemented with the assistance of accurate field data, remote sensing tools (e.g., LiDAR) and growth estimation models. The results of our study in Uruguay suggest that high C prices (>10 USD t −1 ) can result in longer harvest ages, as well as increased timber flows to customers, C sequestration, and economic returns in Eucalyptus stands. Our results also demonstrate that harvest scheduling, harvest age, and NPV are very sensitive to C and timber prices, and that the best economic returns are obtained when the stands are managed to maximize timber production and C sequestration.