Optimal Harvesting Decision Paths When Timber and Water Have an Economic Value in Uneven Forests

: We developed an uneven-aged forest economic decision-making framework that combines: (i) a size-structured matrix model, based on growth and mortality predictions of a dynamic process-based forest landscape model, (ii) an optimal control model that determines the dynamics of control and state variables, which in turn are deﬁned by tree harvesting and forest stock, respectively, and (iii) a water yield function that depends on changes in the leaf area index (LAI), the latter being affected by forest management. This framework was used to simulate the effects of economic-driven harvesting decisions on water yields on a catchment of South-Western Swiss Alps when both timber and water beneﬁts are considered. Water beneﬁts are estimated as environmental prices considering current water demands for drinking, irrigation and hydropower production. We simulated optimal harvesting decisions given the initial forest structure at each 200 m × 200 m grid cells, a set of restrictions to harvesting, and speciﬁc species survival, recruitment and growth probabilities, all of which are affected by the stand’s LAI. We applied this model using different harvesting restriction levels over a period of 20 to 40-years, and accounting for single and joint timber and water beneﬁts. The results suggested that at the environmental prices estimated at the catchment area, water beneﬁts have a slight inﬂuence on harvesting decisions, but when water is accounted for, harvesting decisions would include more tree species and different diameter classes, which, in principle, is expected to favor more diverse forest structures.


Introduction
Optimization of forest resources has received considerable attention for many decades, in particular since the reintroduction of the Faustmann [1] forest optimal rotation model by Paul Samuelson [2]. In its more generic form, the Faustmann model aims to determine the optimal length of time between tree clearcuts, considering a chain of identical even-aged tree cohorts [3]. The Faustmann model has been extended in numerous papers, serving as the base in many research and practical forestry applications, see [4] for a review.
Demands for a more sustainable forestry is increasing the interest in alternative forest management practices. This is encouraging a shift back from even-aged forest management practices, based on clearcuts, to uneven-aged forest management strategies, based on selection cutting of individual trees or small groups of trees, and natural regeneration [5,6].
Uneven-aged forest management models are in general more complex than the even-aged ones. The multiple species and age classes involved in uneven-aged mixed forest stands expand the optimization model dimensions, and demand for alternative computational methods for solving forest management problems. Compared to the Faustmann rotation model, there is no simple generic and analytically solvable uneven-aged model available [3,[7][8][9]. As a result, the number of economic studies that analyse uneven-aged forest management practices is relatively small, and the comparative knowledge of their economically efficient performance still in development [10].
Different mathematical models have been used to analyze the economics of uneven-aged forest stands management. Various authors have extended the Faustmann approach to determine optimal forest cutting cycles assuming the maximization of the present value of returns from a present state to infinity, e.g., Ref. [11,12]. The work of Usher [13], representing uneven-aged forest stand dynamics by a diameter class dependent matrix model has been seminal for the development of economic models looking at the optimal stand structure. Usher's model consists of transition matrices describing how trees in a given diameter or age class move to the subsequent class assuming underlying growth functions that account for how different classes interact in the stand. Adam and Ek [14] were among the first authors building forest structure optimization models upon the transitional dynamics of uneven-aged forest stands [15]. After the work of Haight, and coauthors [16] the optimization problem has been extended to include an optimal rotation age, and since then these transition matrix models have been extensively used to estimate both the (optimal) stand structure and the cutting schedule that maximize economic returns, e.g., Ref. [17][18][19]. Those latter models are often solved by means of numerical methods [10].
Uneven-aged forest management decision-making models have mainly focused on evaluating the economic performance of harvesting schedules and strategies, for maximizing either timber returns [10] or physical yields [8]. Those models usually involve uneven size structures for single forest species, though this type of models have been also applied to determine optimal management models in mixed forest stands [18]. Multiple-benefit forestry economic models are less frequent for uneven-aged stands, even though uneven forest structures are recognized for their biological, hydrological and aesthetic benefits [19]. Some exceptions include the use of uneven-aged forest management models to analyze the economic implications of maintaining certain levels of tree biodiversity [9,19,20], or the integration of carbon benefits into forest management decisions [21,22].
The role that forests play in the provision of water services is also a research topic that is gaining attention in scientific and policy arenas [23][24][25]. The function of forests in water purification is widely acknowledged [26], along with their role in regulating the streamflow regime (e.g., peak flood reduction). Although less acknowledged, there is evidence suggesting that forest cover reduces water supply at the local small catchment level [27], however the associated processes are complex and change over time, and the magnitude and even the sign of the effect of forest cover varies across catchments [28]. Forest and water economic interactions and their role in decision-making frameworks have received relatively little attention in environmental economic literature [29], and to the best of our knowledge, these economic interactions have not been evaluated in uneven-aged and multiple-species stands.
Uneven-aged forestry strategies have a long tradition in Central Europe, particularly in France, Switzerland, Germany and Slovenia [30]. Some of those forests are located in mountainous areas of the European Alps, where seasonal drought is projected to increase under climate change, with potential negative consequences to forests in some regions [31], and where thinning (forest management) can be particularly useful to reduce impact on forest mortality [32], and potentially have an effect on water benefits. This paper develops an uneven-aged forest economic decision-making framework. Our modeling framework combines size-structured transition matrices describing the dynamics of uneven-aged and multiple-species forest stands, with an optimal control model that determines the harvesting schedules and forest structure that maximizes landowner revenues in a spatially explicit way, when both timber and water have an economic value. A novel feature of our modeling framework is that it integrates a water yield function that depends on precipitation levels, potential evapotranspiration, and marginal changes in the stand leaf area index (LAI). LAI is a dimensionless indicator of the tree crown cover that affects rainfall interception and transpiration rates [33]. Marginal changes in LAI are influenced by changes in the number of tress of specific species and size classes due to forest management. We use LAI as a proxy metric of the forest stand density, as changes in LAI due to forest management affect water yield, and potentially water provisioning value.
The optimization model is solved numerically to define optimal harvesting/thinning schedules, given initial forest structures (i.e., number of trees by species and size class in a land unit) and the total stand LAI that maximize both timber and/or water benefits. The model is solved for individual grid cells with a spatial resolution of 200 m × 200 m. Those results were then used to simulate the effects of efficient harvesting decisions on water yields on a catchment in the Central Swiss Alps that is currently affected by Summer droughts [34].

Materials and Methods
Uneven-aged forest management models rely on two main theoretical approaches: the forest growth and optimization models. The growth model simulates the stand dynamics, and the effect of management decisions on the state of the stand (i.e., defined by the stock of trees by species and size class), while the optimization model seeks the optimal management strategies (i.e., harvesting or thinning) to meet specific objectives given a set of constraints.
The stand dynamics is represented by a forest matrix growth model describing tree ingrowth (regeneration/recruitment), mortality and up-growth (transition of living trees between size classes) probabilities. This kind of models commonly assume that these probabilities depend on the basal area of the forest and the stock of trees [19,20], but also on the gap left by trees harvested from different age classes, e.g., Ref. [13]. In our model those probabilities depend on the leaf area index, as an alternative variable to the basal area for evaluating the effect of tree density on ingrowth, up-growth and mortality probabilities.
The reason of using LAI is that this is a forest management dependent variable (i.e., LAI is affected by management decisions) that has an impact on water yield [35]. The causal link between LAI and water yield is theoretically well-established and has been confirmed by experiment and observations in hydro-climatic conditions similar to that of the study area [33]. In addition, an advantage is that the relationship between LAI and water yield can be easily parameterized based on available data (see Section 2.6). We note, however, that this choice ignores some other hydrologically relevant forest properties such as canopy conductance and rooting depth. However, observations on these variables are very scarce and difficult to generalize in space and between different species.

Forest Growth Model for Uneven-Aged Stands
The forest dynamics were described using a size-structured forest matrix growth model. This type of models are an extension of population growth models [17] applied to forest stands to analyze forest management strategies that maximize economic returns, e.g., Ref. [20,21,36].
The state of a forest stand at time t (before harvest) is represented by a column vector y t = [y ijt ], where y ijt is the number of trees per unit of land that are alive in species group i = {1, 2, 3, · · · , m}, and size class j = {1, 2, 3, · · · , n} (i.e., diameter at breast height: DBH) at time t. The entire harvest at time t is represented by the column vector h t = [h ijt ].
The growth model for an uneven-aged forest stand is defined by a set of equations that predict the evolution of the forest stand state in a growing period θ. During the θ time period (5 years in this paper), a tree in species group i and diameter class j may remain into the same class, move up to a larger size class, die or being harvested. The probability a tree of species i is alive and remains in size class j is denoted by a ij , while the probabilities the same tree dies or moves up (up-growth) to the next size class (j + 1) over the θ growing period are represented by w ij and b ij , respectively. In addition the forest dynamics is affected by the expected ingrowth (i.e., the number of trees entering the smallest size class of species in group i) during the θ period.
The state of the stand at t + θ is obtained by matrix operations that depend upon the initial state of the stand after harvest and the up-growth and survival probabilities for each species group and size class. The smallest size class dynamics is further defined by the ingrowth functions (I i,t ), which considers the aggregated natural recruitment and external provision (i.e., thought tree planting) of individuals of the lowest diametric class (see Section 2.2). The state of each species in the forest stand at a κ (where κ ∈ {1, 2, 3, · · · , K}) forest unit and at the time t + θ is represented by the following equations: The state of the forest stand at a κ forest unit after harvest takes the form: where G is the matrix of coefficients defined by G = A + R. A represents the up-growth matrix, which consists of a diagonal of matrices A i , one for each species in group i, describing the probabilities of live trees to remain into the same size class or move to a larger size class over the time interval θ. R is the ingrowth matrix, which consists of sub-matrices R ik that describe the effect of species k on the recruitment of new individuals to the smallest class, and c is a vector containing ingrowth constants that accounts for the number of trees entering to the smallest diameter class for each species that do not depend on the state of the stand (e.g., planted trees): where ι ij defines the relative tree recruitment, or the share recruited trees over the total number of trees of the species in group i with a size class higher than the smallest one (j > 1) over the time period θ.
Here it is assumed that no tree planting takes place, hence c ik = 0.

Estimation of the Ingrowth and Up-Growth Matrix Parameters
Tree up-growth, mortality, and ingrowth functions were estimated based on the hypothesis that their respective probabilities are influenced by the state of the stand and density (determined using total LAI as a proxy). In the literature, the fraction of trees that stay alive and move to the next size class (b ij ) over a determined period, is frequently estimated as a linear or quadratic function of the stand density (usually defined by the stand basal area), and tree size (DBH) [20].
In this paper, both up-growth and mortality were assessed as exponential function of the combined effect of total tree density (total LAI here) and size class (d j ). Total LAI estimation accounts for specific species and size class leaf area index parameters (LAI ij ) and the stock of live trees in a species group i and size class j at the time t (see the Supplementary Material). d j represents the central diameter at the breast height of each size class.
The survival probabilities a ij were estimated as follows: In transition matrix models, ingrowth is often estimated considering both a positive linear relationship between tree recruitment and the specific species abundance, and a negative linear relationship between tree recruitment and the total stand basal area, as proxy of tree density, e.g., Ref. [19,20]. Some authors suggest, however, that a linear ingrowth relationship might be unsatisfactory [7].
In this paper, we explored non-linear functions, such as exponential, and we used total (stand) leaf area index as an alternative proxy of tree density. It is expected that higher densities would affect tree recruitment negatively, with b < 0: Up-growth, mortality and ingrowth probabilities were estimated using outputs from the simulation model FORHYCS [37], which is a coupled model combining the forest landscape model TreeMig [38] and the hydrological model PREVAH [39]. The model simulates the dynamics of forest stands (growth, mortality, regeneration and seed dispersal), which influences hydrologically relevant canopy characteristics such as LAI or canopy conductance. These characteristics affect the simulation of local water balance in the hydrological model. This impacts simulated streamflow, but also feeds back into forest dynamics, as a physiological drought index calculated from local water balance impacts forest growth and mortality in the following year. Inputs to the model are spatio-temporal meteorological data (daily), as well as spatially distributed physiographic and land-cover data.
Following the logic of the forest landscape model TreeMig [38], the forest is represented in each κ cell grid (200 × 200 m) through the number of trees of each species (30 species are simulated in the model, of which 23 are represented in substantial numbers in this simulation). For any year of the simulation period, it is thus possible to retrieve the current state of the simulated forest from the model. In this study, synthetic data (tree species, diameter class and year (y ijt )) was generated for the study region for the period 1971-2015.
The latter tree stock data was generated from the model runs used in [37], where the model was initialized by simulating the succession process for 500 years with bootstrapped meteorological forcing from the period 1980-2000. At the end of this simulation period, forest biomass and species composition reached a state of equilibrium and showed a reasonable match with observations. For more information on the simulations and the model itself, we refer to [37].
On this basis, average up-growth, survival and mortality rates for species in each group i were reckoned (see the Supplementary Material). The parameters of the ingrowth, mortality and up-growth exponential functions were further estimated using Stata 15.

Forest Management Optimization Model
The matrix growth model can be integrated into an optimization model consisting of a function that represents the forest management objective(s) given a set of constraints. For example, the objective of the landowner can be maximizing the financial returns from timber harvesting, while maintaining a sustained supply of timber [19,40]. Those management objectives can also imply optimizing simultaneously the provision of timber and other goods and services, such as carbon sequestration [21,22], tree diversity [20], or the simultaneous provision of water and carbon services, e.g., for even-aged forest [41].
In this paper, we explored the individual and joint optimization of water and timber yields, in uneven-aged forest stands, while ensuring a sustainable forest management and conservation (that is land use change, and total clear-cutting is not allowed). We considered the direct impact of different restriction levels to timber harvesting, and the effect of these restrictions on harvesting decisions, and consequently on forest state (structure) and total LAI, and indirectly on water provisioning service (hereinafter water yield).
We determined the optimal path of harvesting decisions to get a forest structure that maximizes timber or joint timber and water benefits, while satisfying a set of constrains on harvesting and forest growth. Optimal harvesting decisions and resulting forest structure depend on the initial conditions at each forest grid cell (based on virtual forest inventories). As indicated before (Section 2.2), we used the output of the ecohydrological model FORHYCS [37]), which simulates stand dynamics in a spatially distributed way), the forest growth model (depending on total stand LAI at each grid cell and growing period θ), and restrictions imposed to the maximum share of the trees that can be harvested).

The Harvesting/Thinning Decision Model
The optimal harvesting/thinning decision model is implemented as a finite-horizon discrete-time optimal-control problem, where the state vector reports the number of trees (y ij,t ) of species in group i and size class j at each period of time, which in turn is discrete and indexed by t ∈ {0,1,2,...,T}. T represents the evaluated time horizon (20 years and 40 years here). Our model further considers a 10-year harvesting frequency or harvesting cycle (S).
The objective function of the decision problem is framed as the maximization of the finite time horizon net present discounted value (NPV) of the net benefits delivered by uneven-aged forest stands, when both the economic value of the (standing) timber stock at the beginning and the end of the decision period were taken into account [42]. The decision problem is framed from a landowners perspective assuming that water benefits are internalised though payments to the landowner.
In the simplest case, the objective function is to maximize the NPV from timber revenues over a finite time horizon. When dealing with infinite time horizon problems, analysts usually maximize the land expectation value (LEV), which represents the perpetual income stream produced by periodic crops, such as tress, starting with bare land. This latter implies that LEV is estimated by subtracting the initial timber stock value from the NPV of an infinite sequence of forest rotations [43].
In order to deal with a finite time horizon forest optimization problem that is consistent with the optimal solutions for an infinite chain of forest harvesting cycles, we followed Samuelson's [2] recommendation to correct Fisher's [44] false solution based on a single forest cycle. Fisher argued that providing positive real interest rates, it would be optimal to cut a forest stand when its growth rate equals the rate of interest. Samuelson [2] claimed that Fisher's solution was wrong, as this solution is only correct when a forest owner does not continue with the forest activity after harvesting. Frequently, a forest owner maybe interested in replanting the forest stand, and then cutting and replanting it again in an infinite sequence of rotations. The correct solution of the optimal forest rotation problem considering this infinite sequence of rotations was proposed by [1]. Maximizing the net present value (NPV) for just one planting cycle will give a longer rotation period that when an infinite sequence of rotations were considered, implying is that case that cutting more rapidly growing young trees would increase the NPV of the forest compared to a rotation period that follows Fisher recommendations.
Our forest benefit function includes a proxy of the market land rental value. We approximated this market land rental, by quantifying the vector of the expected timber revenues (ev b h = [ev b h,ij ]) for each tree of species i and size class j, given the survival and growth probabilities defined by the uneven forest growth model presented in Section 2.2, and a timber price vector ij ]) that represents the standing value of trees (i.e., timber prices minus variable harvesting costs for each species group and size class). We used the survival and up-growth probabilities to estimate the maximum harvesting feasible, which corresponds to the harvesting rates that equal net timber growth over a growing cycle. The harvesting probabilities were estimated for each one of the size classes, and those indicate the probability ( ij ) an individual tree of a species in group i to be felled as it reaches a size class j.
For the estimation of the ev b h,ij values, we assumed that all living trees will be eventually felled, and therefore the aggregated harvesting probabilities over a T growing cycle totals one: The economic problem is to find the optimal harvesting path (h t ) and tree stock (y t ) in a harvesting cycle T that maximize the NPV of the net forest benefits over a finite time horizon, given constrains on: (i) forest (biological) growth (Equation (12)), the feasible harvest (Equation (13)), (iii) non-negativity constrains (Equations (14) and (15), and the initial forest state, defined as the initial distribution of trees by species and size class (y 0,κ ) at each forest grid cell (indexed by κ): subject to: where FC is the fixed harvesting costs, v b · y 0,κ is the stock value of the standing trees at the beginning of the evaluated period, y T,κ a vector representing the stock of trees at the end of the evaluated period in the κ forest unit, δ is the discount function (δ = 1/(1 + r), being r the discount rate).
To estimate the market land rental value after the cutting cycle T, we considered a subsequent growing cycle of a τ-years duration. It is assumed that identical chains of the growing cycle are periodically repeated. The term (1 − δ τ ) is used to approximate the stock value of the expected timber revenues as the sum of an infinite geometric series (where δ τ < 1).
The coefficient η stands for the maximum harvesting intensity allowed, (where η ∈ {0, 1}), which is used to avoid complete deforestation solutions. We considered, in that case that harvesting/thinning interventions cannot be higher than a η share of the tree stock, at every harvesting period evaluated.
For the joint maximization of timber and water benefits, we added a second term to estimate the value of felled trees and final tree stock, defined by the price vector v w = [v w ij ], and the expected benefits associated with water yield. Those price vectors account for the potential water benefits due to a decrease in the number of trees and their associated LAI, which result in an increase in water yield, as shown later.
The objective function in the presence of payments for increasing water yields can be written as follows: the expected value of each individual tree of species i and size class j in terms of their contribution to water benefits. As a decrease in tree density, hence in total LAI, favors water yields in this model, no additional water provisioning services is associated with to the stock of trees.

Forest State Dynamics
The optimal stock (forest structure) at each time period is defined by the state vector (y t ) at each period of time t. Changes in the state variable along the evaluated time periods rely upon the ingrowth, mortality, and up-growth functions, which in turn depend on the total LAI at the beginning of each θ growing period, and on the optimal harvesting decisions. Forest growth for each species in group i is estimated using Equation (1) for each one of the θ growing periods comprised in one harvesting cycle T. LAI is recalculated after each θ growing period, and up-growth, mortality and ingrowth probabilities reappraised in accordance, considering Equations (6), (7) and (9), and integrated again into Equation (1).
The decision maker must decide upon the optimal number of trees to be felled by species and diametric class at each harvesting cycle S. Harvesting is the control variable that models the decision mechanism. The control variable enters in the state motion equation through a single control vector (h t ) that influences the forest state.
The model is solved in Matlab R2016a for different finite time horizon periods (T ∈ {20, 40}), using linear programming techniques, and assuming that harvesting takes place at the beginning of each harvesting period S, every 10-years, which in turn comprise two consecutive θ growing periods of 5 years.

Case Study for the Application of the Modeling Framework
The combined forest growth and optimal-control models were applied to the Navisence River catchment area. This catchment is in the Central Swiss Alps, in canton Valais (see Figure 1), covering about 250 km 2 , which includes close to 50 km 2 of glaciers and 46.84 km 2 of forest. This catchment comprises three Municipalities: Anniviers, Chalais and Chippis; and it is divided into 5 subcatchments, defined by the location of hydropower plants, where streamflow is gauged, see [37] for details. These subcatchments differ markedly in their climatic and land-cover characteristics.
We applied the forest growth and management optimization models to grid cells of sub-catchment 3, which covers 42% of total forest area of the Navisence catchment, with altitudes ranging from 1000 m to 2400 m asl. This provides a case study with a gradient of temperatures and precipitation, allowing for spatial heterogeneity in vegetation cover and land uses.
As mentioned before, the forest growth model is parameterized with virtual forest inventory data, i.e., the output of the dynamic forest ecohydrological model FORHYCS [37]. The reason for this choice is that FORHYCS outputs have a higher temporal (annual) and spatial granularity (grid cells of a 200 m × 200 m resolution) than observations: there have been four forest inventories in Switzerland since 1981, on a grid with spacing between 1 km (first inventory) and 1.4 km (second to fourth inventories). To generate robust and reliable estimates of growth and species composition from inventory data, it is necessary to average data from a sufficiently large number of inventory sites, which is problematic given the small extent of the study area. In addition, this area is characterized by a sharp elevational and hydro-climatic gradient, so that inventory data should additionally be stratified by elevation, which further increases the data requirements. For this reason, synthetic data from a simulation model was preferred to inventory data. Results from the FORHYCS model showed a reasonable agreement with observed species distribution, biomass and canopy structure [37]. The model accounts for 23 different tree species, whose average distribution by size class is shown in Table 1. According to FORHYCS model outputs at the begin of the study period, conifers contribute 75% of total forest LAI at the sub-catchment level, while broad-leaf species cover the remaining 25%. Larch (Larix decidua), Norway spruce (Picea abies) and Scots pine (Pinus sylvestris) are the main tree species in the case study area.
We grouped the forest species into five LAI classes and seven timber price categories (see Table 1). Each LAI class is defined by specific leaf area index and allometric parameters [45] that relate DBH to the foliage weight and area (see Supplementary Material). Up-growth, ingrowth and mortality probabilities were estimated by each timber price group (see Section 2.2 and Table 2). Timber price groups were defined considering akin species with similar prices for the same DBH (see Section 2.5). Akin species share genus or belong to the same group of species such as broadleaves or conifers (Groups 2-6), or to a group of species (i.e., Group 1 or Group 7) with no commercial interest as timber or biomass, or with interest only for biomass. Optimal harvesting and forest stock were defined at each cell for each one of the 23 tree possible species, and eight different size classes. To facilitate the presentation of results, the outcomes of the model are further regrouped into the 7 timber price groups.  (1) Size classes are defined by the diameter at breast height (in cm).

Valuation of Timber and Water-Related Ecosystem Services
Only two types of forest benefits were considered: timber and water. Timber price category groups of akin species with similar prices for same DBH and timber quality. The classification of species into a timber price groups considers species which prices for different tree diameters do not depict significant differences. We included a category for those species without market interest to which we assigned a zero timber price (Group 1). We further considered a category for species only used as biomass (Group 7), three groups of conifers with different price levels (Group 2, 4 and 6), and two groups of broadleaves with different price levels (Groups 3 and 5) (see Table 2). Timber prices were estimated for different forest species and size classes in Switzerland from years 2000 to 2014. Timber price estimations were based on price data provided by the Federal Statistical Office used for the economic valuation of standing timber in Switzerland [46]), and properly updated to year 2014 using the Swiss consumer price index [47] (see Supplementary Material). Available timber prices correspond to prices at the forest gate, thus we estimated timber standing values by removing average timber harvesting cost. Average timber harvesting cost in Switzerland fluctuates between CHF 50 m −3 to CHF 55 m −3 [47]. We further considered an average harvesting cost of CHF 50 m −3 , and assume no fixed harvesting costs per hectare.
As mentioned before, forest-hydrology literature suggests that water yield decreases as the leaf area index increases [33,35]. Leaf area, on the other side, can be estimated by individual tree in relation to its DBH, and for the forest stand considering its state and dynamics. The stand LAI is, in this model, the only forest management element affecting water yield, and as tree harvesting reduces total LAI, we estimated the marginal effect of changes of LAI by species and age class on water yields, as detailed in Section 2.6.
Changes in water yield were priced considering the demands of three water-related provisioning services: irrigation, drinking water, and hydropower water, and their estimated environmental prices in the Navisence catchment area (Table 3). Environmental prices were here defined as the unit contribution of natural capital to water supply, and they were estimated as a residual value of unit revenues from each water-related service, after covering its production costs, and the returns to labor and manufactured capital as production factors, in line with a [48] recommendations. Forest water is valued then in view of those environmental prices and corrected to account for the share of forest water with an economic use. The estimated environmental price of water represents the maximum amount that could, in theory, be transferred to the owner of the natural asset (forest) that deliver the water supply services.
The irrigation water demand in the Navisence catchment is estimated considering the irrigation water use in the neighboring valley to our case study (Crans-Montana-Sierre). In that case, we calculated the unit irrigation water use by hectare of agricultural land at the 11 Municipalities included in [49] study, and assuming the same per hectare water demand for the agricultural area of the three Municipalities inside the Navisence catchment [50]. Drinking water accounts for both residential and touristic sector water demands. Those demands were estimated in view of the average consumption by inhabitant and tourist night in Switzerland and European countries, respectively [51,52]. Hydropower water demand and its environmental prices were gauged using data recorded by the hydroelectric company operating in this catchment [53], and further economic data for the hydropower sector in Switzerland [54] (see Supplementary Material). In parenthesis the environmental prices for forest water with an economic demand; Own estimations based on: (a) [50][51][52][55][56][57]for drinking water; (b) [49,58] for agricultural water, and (c) [53,54] for hydropower.

Ecohydrological Functions
The annual average water streamflow is estimated in 191.6 million m 3 (±190.3 Mm 3 ), which is used for hydropower production [33]. On the other hand, we estimated that the average annual demand of water for irrigation and drinking water represents a 3% of this streamflow We assumed only this latter to be the share of water with economic value as fresh and irrigation water. In view of the residual values presented in Table 3, we estimate a water environmental price (p w ) of CHF 0.10 m −3 ± 0.02 m −3 .
To operationalise the causal relationship between LAI and water balance, we used MODIS leaf area index/fPAR data, being fPAR the fraction of photo-synthetically active radiation (400-700 nm) absorbed by green vegetation. LAI and fPAR data were used for calculating surface photosynthesis, evapotranspiration and water yield to estimate total evapotranspiration.
Total evapotranspiration is estimated as function of LAI, total precipitation (P) and potential evapotranspiration (PET) in the catchment area. P and PET values refer to total annual precipitation and potential evapotranspiration, and treated as constant values: AET = P · (−0.014 · (LAI) + (0.526 + 0.024 · (LAI)) · (P/PET)) (17) We further assessed the contribution of marginal changes in LAI on the total evapotranspiration, and water yield (W), which is in turn estimated as the residual value of total precipitation and total evapotranspiration. The marginal contribution of changes in the leaf area index on the water yield (W ij ) is estimated considering the change in LAI every time a single tree of a species i and size class j is felled: The elements of vector of water provision prices (v w ), were estimated then as: v w ij = p w · W ij , which gives a different water yield value associated with harvesting trees of different species and size classes. Those prices were then considered in the optimization problem, to analyze if the internalization of marginal water values has an effect on the optimal distribution of species and age classes.

Results and Discussion
This section presents and discusses the results of the optimizations modeling framework. Those results include the estimation of the parameters of the forest growth matrix models (Tables 4-6, and Figure 2). This also includes an interpretation of the results comparing optimal harvesting decisions and the resulting forest structures with different restriction levels to harvesting over time, and when both single or joint timber and water benefits are considered (Figures 3 and 4 ). This section also addresses the effect of optimal harvesting decision on forest tree diversity ( Figures 5 and 6), and the effect of optimal harvesting decisions on water yields over time. Results also show the spatially distributed effect of the optimal harvesting decision path on water yield at the study case area (Figures 7 and 8).

Forest Growth Functions
Annual up-growth, mortality and ingrowth annual probabilities were estimated by group of species, classified according to the timber price category those species belong to (Table 1). Up-growth and mortality probabilities are estimated as exponential functions of the combined effect of total stand LAI and DBH (Tables 4 and 5). These functions allowed the estimation of specific probabilities associated with species group i and size class j.
A larger stand LAI combined with a larger diameter class is expected to increase the time a tree remains in a j size class, and hence, to reduce its up-growth probability. Likewise, it is expected that the mortality rates would decrease within tree diametric (size) class, an increase within stand density (LAI). Figure 2 illustrates the effect of LAI and tree DBH in the estimated annual up-growth and mortality probabilities for a range of diametric classes and stand LAI values for species in groups of 2 (i.e., lower price conifers) and 6 (i.e., high price conifers). In this figure it could be observed that tree mortality probabilities decrease within the tree size, while those probabilities would slightly increase within the leaf area index. In our model, both mortality and up-growth probabilities change every θ growing period depending on tree growth and harvesting decisions, which in turn affects total stand LAI and the distribution of size classes. This dynamic is considered into the optimal decision problem (Equations (11) and (16)). Table 4. Parameters of the up-growth probability functions by group of species (1) . Notes: * p ≤ 0.1, ** p ≤ 0.05, *** p ≤ 0.01; (1) Exponential function of the type: y = A · exp b·(LAI·d j ) , where y is the share of trees that move up to the next diametric class in a year, total stand LAI was considered; d j is the diameter at breast height (DBH); (2) N is the number of observations (FORHYCS model outcomes).

Up-Growth Functions
The parameters of the exponential ingrowth probability functions presented in Table 6 allow estimating annual tree recruitment probabilities by species group. Tree recruitment is estimated as a solely function of the stand density (represented here by total LAI). The effect of total (stand) LAI on ingrowth probability is expected to be negative and significant, as higher tree densities are likely to reduce tree recruitment. In all cases (functions by group of species), the LAI is a significant explanatory variable of ingrowth probability functions. In addition, as expected, the LAI coefficient is negative for the exponential ingrowth function estimated for species of the groups 1,2,4 and 5, and also for the ingrowth function estimated for all the species grouped together (see ALL in Table 6). This expected relationship does not hold for species of the groups 3, 6 and 7. In these three latter cases, we opted for using the exponential ingrowth probability function estimated for all the species grouped together (ALL). Though, we acknowledge that more research is needed to improve our understanding of tree recruitment by specific species. Table 5. Parameters of the mortality probability functions by group of species (1) .

Statistic
Coefficient Robust SE Coefficient Robust SE Notes: * p ≤ 0.1, ** p ≤ 0.05, *** p ≤ 0.01; (1) Exponential function of the type: y = A · exp b·(LAI·d j ) , where y is the share of trees that die in a year; total LAI is considered; d j is the DBH; (2) N is the number of observations (FORHYCS model outcomes). Table 6. Parameters of the ingrowth probability functions by group of species (1) . Notes: * p ≤ 0.1, ** p ≤ 0.05, *** p ≤ 0.01; (1) Exponential function of the type: y = A· exp(b · LAI), where y is the relative ingrowth, which is estimated as the ratio between total trees of the species i of the lowest diameter class (DBH ≤ 2.5 cm) and the number of trees of the remainder diametric classes; total LAI is considered; (2) N is the number of observations.

Statistic
The parameters of the exponential up-growth and mortality, probability functions are used to estimate the annual probabilities that an individual tree of species in group i and size class j moves up to the next diameter class (up-growth), dies over a period of a year, respectively. The ingrowth exponential functions, provide in turn the annual probability a tree in group species i enters into the lowest size class (j = 1) over a period of a year. The forest growth model parameters b ij , w ij , and ι ij used in Equations (3) and (4) are therefore estimated by multiplying the above referred annual probabilities by θ (which as indicated before takes a value of 5 years here).

Optimal Forest Structures and Management When Timber and Water Benefits Are Maximized
We explored the optimal harvesting paths considering four consecutive harvesting cycles (S), where S ∈ {20, 30, 40} years, to move from current situation (not optimal), which is in turn defined by the initial forest states, to the optimal distribution of tree species and size classes at each harvesting/thinning period. This later defining the transition to the forest stock that maximizes the NPV of future forest benefits for both timber and water yield.
For the application of the model, we assumed that harvesting/thinning operations take place every ten years, embracing two consecutive θ growing periods (where θ equals to 5 years). Figure 3 shows the main results of the optimization framework in terms of the number of total trees (including all size classes and species) and the number of trees felled for each one of the harvesting cycles considered and each one of the grid cells at the case study area (K = 520). This latter figure reflects the results for a moderate harvesting restriction level (η = 0.5), and for the join timber and water benefits maximization scenario. In general terms, the results suggest that reducing the density of the stand (in view of the number of trees) is optimal from a landowner's benefits maximization standpoint. When the consecutive harvesting cycles are compared, we can observe that harvesting/thinning patterns are similar for the 20, 30 and 40 harvesting cycles, while for S = 10 they imply a higher harvesting intensity. In other words, the results suggest that increasing harvesting/thinning efforts to reduce current tree density would maximize the NPV of join timber and water benefits.  The initial forest state (y 0 ) at each one of the cell grid has an important influence in the net present value of timber benefits (Figure 4a). The value of those benefit assuming a T = 40 time horizon and a 4% discount rate accrues 2443.4 Swiss Francs (CHF) per hectare (± 1220.4 CHF/ha) when only timber is considered, and 2442.5 (± 1220.8) CHF/ha when both timber and water benefits are jointly maximized. The effect of water benefits are relatively small when compared to timber ones, and in an hypothetical scenario where only water has an economic value the NPV of water benefits mounts up 456.2(± 1005.2) CHF/ha (Figure 4b).
Despite the relatively small importance of water benefits compared to the timber ones, our results show that integrating water values (scenario W and TW, for only water and for join water and timber benefits, respectively), encourages more diverse harvesting decisions, involving all groups of forest species ( Figure 5). For instance, if water were the major (or only) benefit from forest management, species belonging to the group 1 (with no economic interest for timber) will be also managed. This is expected to enhance a higher diversity in tree species. Harvesting restrictions (defined by η) do not seem to have a significant effect on harvesting decisions in terms of the distribution of species harvested ( Figure 5). In all cases, species of the group 2 (i.e., conifers with the lowest timber price) which are the dominant species in the study case area (mainly Pinus sylvestris and Picea abies) are the ones with a higher harvesting/thinning rates. Harvesting intensities are, nonetheless, relatively low, considering that the optimal decisions imply that no more than 15% of existing stock by group of species would be harvested. The results also suggest a different story when harvesting decisions are evaluated in terms of size class. The results show that when higher harvesting rates are allowed (i.e., η = 0.8), it would be optimal to harvest trees of relatively small diameters at breast height (i.e., 10-20 cm).
Harvesting decisions will affect the long-term dynamics of the species considered ( Figure 6). In general, the overall tendency on stand dynamics is to reduce current density (defined by the number of trees), as Figures 3 and 6 show. Optimal harvesting decisions depend on different factors, such as timber prices and growing rates, while the species dynamics would be affected also by recruitment and growing rates, describing growing cycles driven by survival and up-growth probabilities, and harvesting decisions. Given the harvesting restriction it is expected that the diversity of tree species will be maintained over time ( Figure 6). Nonetheless, both most abundant (e.g., Pinus sylvestris or Sorbus aria), and most profitable species (i.e., Larix decidua with a higher timber price) are favored, in detriment of broadleaved species with a lower commercial interest.

Effect of Optimal Harvesting Decisions on Water Yield
As expected, integrating marginal water benefits into the decision making problem in general reduces the forest density (in number of trees and in terms of LAI), in the short to medium terms (>20 years), with a marked increase in water yields respect to the initial situation. Increased thinning affects forest growth dynamics, fostering recruitment which tend to increase LAI in the long term, with the subsequent effect on water yields (see Figure 7).
Water and timber relative prices are assumed to be similar in the future than they actually are. In this context, the effect of increasing water prices as result of future water scarcity due to climate change [59] needs to be analyzed within the optimal harvesting decision model. To use this modeling framework as a more functional forest and water management decision tool, beside potential water benefits based on current demands, we need to better integrate potential changes on downstream water demand and welfare changes due to the forest interventions analyzed.
Expected changes in water yield vary within the harvesting period and the maximum harvesting rate allowed. Figure 8 depicts the spatial distribution of changes in water yields, after the first and second harvesting respectively, and considering higher (η = 0.30) to lower (η = 0.80) harvesting restrictions. The results show marked differences between the west and east sub-catchment areas, specially for the second harvesting period. These latter results seem to be connected to the orientation of the grid cells, which affect the species original composition, according to the FORHYCS model [37]. This model contributes to the yet scant literature on forest hydro-economic models, and to the best of our knowledge this is the first model that considers simultaneously timber and water benefits to find the optimal harvesting decision paths in uneven forest structures. The model is defined in a more deterministic fashion. Thus, future developments need to improve the integration of multidimensional risk and uncertainties in the forest management model, including the effect of uncertain precipitation levels on the marginal benefit of water yields due to tree harvesting.
The results must be interpreted considering the method chosen to parameterise the forest growth model: as described in Section 2.2, the output from a process-oriented forest landscape model (number of trees in each year and grid cell, clustered by species and size class) was used as virtual inventory data. Thus, effects of assumptions and simplifications in the forest landscape model are propagated to the growth model used here. An important aspect is that the allometric parameters relating leaf area to tree size are considered time-invariant. In reality, the specific leaf area of a tree can be expected to vary based on tree age climatic conditions and stand structure [60].
Another important aspect to consider is that the simulations with the forest landscape were carried out using meteorological forcing corresponding to present-day conditions . Therefore, the forest growth model parameterization used here might not be representative of future conditions. Another simplification of this study is that LAI is assumed to be the only variable affecting water yield. In reality, other stand variables, such as rooting depth and canopy resistance, may have a strong effect on local water balance [33]. However, simulating these variables and their interactions with local water balance is a complex task [61,62]. For example, Chen and coauthors [62] report on a case where thinning led to increased transpiration through complex biophysical interactions. In this study, as the model only simulated moderate thinning, we do not expect that such effects would have a large impact on the results.

Conclusions
Forest management decisions necessarily deal with complexity and trade-offs, in particular when the heterogeneity of forest attributes, growth and yield potential, and multi-benefit forestry are considered. Traditional forest management decision models are based on relatively simplistic forest representations, namely mono-specific stands with trees of the same age. Ongoing switches in forest management policies and strategies towards enhancing ecosystem services provision, demand for more complex models that account, for example, spatial heterogeneity and uneven forest structures. This paper develops an optimal decision model that considers uneven forest structures and dynamics, and spatial heterogeneity, while proposes an approach to integrate timber and water benefits into decision-making. This basic model can be extended to the integration of other ecosystems services that depend on forest structure dynamics, such as carbon sequestration, and water quality protection.
From a more general modeling perspective, water scarcity does not seem to be a priority problem in temperate forest areas of Europe, though climate change may revert this situation in the future, demanding for water-sensitive forestry [31]. The model developed in this paper is better suited for forests of sub-arid to arid areas (e.g., Mediterranean forest) that demand for water-sensitive land use and management strategies. Water scarcity is not yet an issue at the case study area where the model was applied. This latter is reflected in the relatively low environmental price of water estimated. Nonetheless and despite the low value associated with water provisioning service, this paper shows that the integration of water values into the decision problem can potentially lead to a different consideration of what would be an efficient decision. The results of this paper suggest that optimizing forest management for timber and/or water yield can drive to changes on the stand structure (reducing or increasing biodiversity levels over time). The integration of water values can favor the maintenance of more diverse forest stands, at the time it may enhance forestry revenues. This latter specific result highlights the need for developing models that can discriminate the potential effect of single species or specific group of species in the provision of ecosystem services. Timber values for instance are normally differentiated by single or a homogeneous group of species.
The model developed in this paper and its application to the case study area consider also that single or homogeneous groups of species can deliver different levels of water benefits, and henceforth the integration of water provisioning service values can promote different forest structures and species distribution.
More research is needed to integrate the potential impacts of climate change on forest and water yields and forest growth, and their potential economic effects.