As discussed in the Introduction, our model of thinning-driven stand dynamics is based upon state dependent life history theory for individual trees, which are then extrapolated to the stand level. We focus on three trees states: size (S), reserves (R) and age (A) whose units are defined as kg, kg and years, respectively (in section 2.2, Equation 12, we transcribe size to volume to bring in line with traditional forestry approaches). A fourth state variable, competitive ability (ψ), enters into stand-level dynamics described in the following paragraph. How trees respond to A, S and R states is modelled via incorporation of four important life history assumptions: (i) innate tree growth follows a sigmoidal pattern (see Weiner and Thomas, 2001), (ii) based upon principles from evolutionary ecology, we assume that tree growth strategies have evolved to maximize lifetime reproductive success (Davies and Krebs, 1993), (iii) individual trees will acquire resources and modify their investment in growth in a state dependent manner according to (S), (R), (A) (see Equations 1-3 below) and (iv) there are no genetic constraints on the above i.e., we play the phenotypic gambit (Grafen, 1984), which assumes that selection on phenotypic variation translates directly into selection on heritable variation in the population. Finally, though we do not explicitly say so, we assume that there is an underlying tradeoff between growth and reproduction (e.g., Kassaby and Barclay, 1992) though we acknowledge that such a tradeoff is often not obvious (e.g., Harshman and Zera, 2007; Suzuki et al., 2019).

We assume that individual trees at a given combination of our three state values (S, R, A) will adopt a resource allocation strategy. Thus, a focal tree, n, in population of N trees with state variables S = s, R = r and A = a will grow in an adaptive manner as follows:

We expect that trees will mobilize reserves for growth and reproduction, at a rate that we define as alpha (α).

Where: R_min is the minimum reserves value for a tree of size (s), with s, r and a defined as above and W values are weight constants for the state variables (see Table 1 for a list of variables and their values). For a tree of a given size (s), we expect that α will increase with size, reserves and age up to some maximum that never causes the focal tree to exceed its critical reserves state value (R_crit), which is the minimal level to maintain metabolic function. Here our reasoning is: (i) trees cannot forecast weather conditions for the current year and (ii) trees are risk averse and will never mobilize so many reserves so to put their survival at risk should the current year’s weather turn out to be poor for energy harvesting (see Clark’s asset protection principle, Clark, 1994).

We define beta (β) as the proportion of mobilized reserves that go to reproduction versus growth:

Where π is the start-up cost for initiating reproduction and the W’s are weight constants for size and age states ( Table 1 ).

Here, we assume two effects of size. First, it is necessarily true that small trees have small reserves such that R will rarely exceed π, meeting the β = 0 condition. Second, small trees are not efficient at producing seeds (e.g., Sherman et al., 2019) thus, we have set the size weight (W₃) low ( Table 1 ). Regarding age weight (W₄), we further assume that older trees should strongly bias their mobilized reserves to reproduction because future discounting (e.g., Lee et al., 2011) will offset small gains from growth particularly if older trees are large, thus, we employed a small value for W₄.

We define phi (φ) as the allocation of harvested energy to growth versus reserves.

where: ϕ s = 1 − s − s m i n s (equation 3.1) and ϕ r = 1 − ( ( r − ( R m i n ) ) r ) (equation 3.2), based upon size and reserves, respectively. In addition, s ' = s + γ ( 1 ϕ s + ϕ r ) ϕ s + ( α ( 1 − β ) ξ ) (equation 3.3) and r ' = r + γ ( 1 ϕ s + ϕ r ) ϕ s − α (equation 3.4) are potential updated values for size and reserves, respectively and ξ is defined as the efficiency in mobilizing reserves to structure (see Equation 7).

Here, three conditions hold: Sufficient reserves remain after metabolic costs and reproduction start-up cost are paid and allocation to reserves does not exceed the maximum reserve level after adjusting for energy harvesting (γ). Otherwise, φ is adjusted to meet those minimum and maximum reserve levels, respectively with γ and π being functions of tree size. The rationale for the φ decision is: as trees increase in size, their growth rate decelerates (Weiner and Thomas, 2001) thus there is little gain from investing in growth vs building up reserves for reproduction (Lalonde and Roitberg, 1992).

All three of the above state-dependent strategies can be summarized as: invest in growth when that leads to accelerated returns but only when reserves are sufficiently high to avoid starvation i.e., metabolic reserves fall below critical level. When trees are large and growth is constrained, trees are assumed to invest in reproduction when it is safe to do so because of the greater return on investment for the latter. Also, future discounting means that future reproductive returns from current investment in growth will be constrained for older and larger trees (Equation 2).

We employed Equations 1-3 in simulations of stage-dependent tree growth, which we refer to as Tree Adaptive Growth or TAG. Details of the TAG model follow.

We developed an agent-based simulation model using C++ language, which is shown in flow chart form in Figure 1 . Here, we simulated small stands of (N = 100 giving a starting density of 2500 stems/ha) trees that were subjected to silvicultural practices including pre-commercial and commercial thinning. We assumed that each stand was comprised of lodgepole pine (Pinus contorta) though we expect that our model could be employed for any tree species with known growth parameters or, potentially, mixed stands. The TAG model is similar in spirit (but not detail) to the DESPOT model of Buckley and Roberts (2006) in that individual trees are assumed to maximize some goal, in their case carbon gain and in ours, lifetime reproductive success.

Each stand was planted with identically aged individuals at year=0. The size (mass in kg) of each tree was randomly generated from a normal distribution (x=15, sd=5). Tree reserves (mass in kg) were generated from a uniform distribution (1,5). The simulation tracks size (S), reserves (R) and cone production (Λ), if any, each year for 200 years.

We employed a ‘weather’ concept to describe conditions that impact growth i.e., energy and nutrients that are harvested by trees for growth and maintenance. Our weather concept was created for characterizing the distribution/spectrum of favorable to unfavorable environmental conditions and such causative agents as variation in hydrological factors, heat stress or insolation. As many (most) of these factors show strong autocorrelation, we chose to use a single parameter in this proof-of-concept model, rather than breaking ‘weather’ into component effects.

At the beginning of each year, weather is randomly drawn from two vectors ς_s and ς_f, which refer to weather intensity and frequency, respectively (see Table 1 for details). Tree growth for any individual tree is a product of growth potential (Equation 4) and weather intensity for that year. We assumed that, in a given year, weather intensity is identical for every tree in a stand i.e., our model is aspatial (see Discussion).

At the beginning of each year, each tree is subjected to a random mortality event by individually drawing from a uniform distribution (0,1). If the random number falls above a preset default- annual-survival value, then the tree dies. For our runs, the default value for annual survival was set at 0.998. While this survival value is a constant, it is also possible for a tree to die if its reserves fall below a critical value, R_crit, which would happen most frequently when an individual’s R state is low, weather is poor (see Table 1 ) and the focal individual is a poor competitor (Equations 5 and 6), thus the assumption of risk averse investment and reproduction as above. In these simulations, we assumed that R_min = R_crit H.

where: H describes the degree of risk aversion.

Following the mortality evaluation, the time counter for the current year was compared with the pre-determined time step for PCT (e.g., year = 15). If the current year was determined to be a PCT year, then individual trees were removed according to specific rules - in default simulations, the quartile of smallest trees were removed i.e., thinning from below; other thinning rules were also applied as discussed in section 2.3.vii. Late PCT was applied in a similar manner for trees later in the simulation.

Trees that survived random mortality and were not thinned were allowed to grow according to a modified logistic model where the potential growth increment for a tree of size (s) was:

Where: ρ is the intrinsic growth rate, S_max is the maximum mass that an individual tree can realize (but see Bontemps and Duplat, 2012; Stephenson et al., 2014), θ is a shape parameter and ς is a senescence term that describes decreasing energy harvesting with age (e.g., Hubbard et al., 1999; Qiu et al., 2021 but see Lanner and Connor, 2001).

To determine our default value for ρ, we initiated simulations with ρ= 1.0 and then iterated a series (1.0, 1.20), increasing by 0.01 for each iteration and then plotting the stand volume against time from 0 to year 200. Since our goal was to produce a representative growth curve we visually confirmed that the TAG stand dynamics produced a sigmoidal plot as shown by existing growth and yield models like GYPSY (Huang et al., 2001, 2009) in Alberta (AB) and VDYP7 (British Columbia Ministry of Forests, Lands, Natural Resource Operations and Rural Development, 2019) in BC. Three features that we considered were, (i) the point of obvious acceleration in the growth curve, (ii) the inflection points (representing the change from the accelerating to the decelerating part of the growth curve) and (iii) maximum value of the growth curve for an unthinned (control) plot. In our case, the 15 m³/ha, year 50 and 600 m³/ha were used, respectively. Figure 2 shows the TAG output that that we used based upon these criteria. In the end, we chose ρ = 1.12, the value that gave good representation of stand dynamics. In the default model, we assumed that θ =1 i.e., there is a linear decline in per unit mass growth with tree size, yielding a classic symmetric growth increment curve. In some simulations, we modified θ to give right or left skew to the size performance curve.

We further modified an individual’s potential growth increment from energy harvesting by including competitive performance. We assumed maximum performance for a tree with no competitors whose value declines with increasing density of competitors over time according to the function:

Where: N₀ is the number of trees at planting (year = 0) and N_y is the number of trees at year y and τ’ is a shape parameter for the function and ψ_s,n,y is relative competitive ability of a tree, n, of size s in year y, which is defined as:

Where: s_n,y is the size of a focal individual tree n at year y and S_y is mean tree size for a stand at year y.

Once the potential growth increment was calculated for a given tree of a given size (s) and reserve (r), new size and new reserve state values were calculated according to the three allocation rules (α, β and φ) as discussed in section 2.1.

A tree can also increase its structural mass by mobilizing reserves (β< 1 Equation 2). In doing so, we assume that there is some inefficiency in such mobilization, which we define as ξ i.e., ξ< 1.0. Thus, growth from reserves in year y is:

Reserves are decremented when trees pay maintenance costs, which we assume are size dependent. Thus, for a tree of size S, maintenance costs are:

Where: ω is the per unit mass cost of maintenance I n this version of our theory, we set ω to 0 based upon the assumption that maintenance costs are already subsumed in the growth function (Equation 4). Future work, may explicitly separate the two parameters.

In years in which an individual invests in reproduction, i.e., β > 0, we assume that there is a linear size dependent cost to reserves from initiating reproduction, which we define for a tree n of size S in year y as:

Where: ϵ and ϵ’ are the initiation slope and intercept, respectively.

From Equation 1 through 9, we can calculate (Equation 10) the size and reserves (Equation 11), respectively of an individual tree n in year y+1 based upon its size in year y as:

and

The simulation produces an even-aged stand of trees whose size distribution changes from year to year. There are several factors that generate individual size differences among trees. First, trees vary in size and reserves at initiation. From there, based upon absolute and relative size, and reserves, trees vary in their capacity to harvest energy but also in their tendency to allocate their harvest to growth, reserves and reproduction.

The process was repeated until either the year counter reached 200 and the entire stand was harvested, or no longer held any viable trees.

For each series of simulations described below, we replicated the treatment in 10 different even-aged stands, 5 stands for non-thinned controls and 5 stands for PCT. Each stand was originated using a unique set of random number seeds, however, identical random number seeds were employed for pairs of control and treatment stands. As such, each pair of stands generated identical dynamics up until PCT was implemented, confirming that any differences observed were due to treatment effects and not stochasticity.

To make our model relevant to forestry practices, we converted tree size to volume (V) by the following physical equation:

where: D = wood density, which is species specific, in this case, for pine, is 455 kg/m3 (https://www.thecalculatorsite.com/conversions/weighttovolume.php).

For each simulation, we considered two different management tactics: (i) no pre-commercial thinning (from here on referred to as Control) and (ii) PCT as described in section 2.2.

To demonstrate the utility of our approach we used the simulations in the following manner:

Standardized relative growth (SRG) (Equation 13) is a proper indicator of the status of compensatory growth. We calculated RG to show the pattern of stand growth in thinned plots respect to control as follows:

Where V_thinned and V_Control are stand volumes for thinned and control plots, respectively. By definition, under compensation is denoted when SRG is less than 0, exact compensation happens when SRG equals to 0, and overcompensation is showed when SRG is greater than 0.

The primary goals of this paper are: (i) as a proof of concept, develop a biologically-based theoretical model for tree and stand growth (TAG), (ii) demonstrate that the TAG model generates stand dynamics similar to that found in nature, (iii) demonstrate that stands that experience thinning can overcompensate under some conditions and (iv) explore the sensitivity of the compensatory growth to some easy-to-implement management tactics (e.g. intensity of thinning) as well as inherent characteristics of different tree species (e.g. sensitivity to crowding).

Based upon evaluation of the TAG model, the following insights emerged: (i) compensatory growth, following thinning (natural or otherwise) occurs readily especially early on in the life of a stand, (ii) the window for overcompensation narrows inversely with the intensity of PCT simply because increased performance by individuals is negated by reduced numbers of such, (iii) early PCT is beneficial because increased investment in growth is most readily observed during the accelerating part of the growth curve, (iv) results from the various sensitivity analyses were in line with expectations and (v) relative performance of thinned stands were higher than those observed in nature (see below). Finally, an important insight from employing individual trees was that the ‘catching up’ response of surviving trees, post-PCT ( Figure 2 ), can be explained by such individuals gaining greater access to limiting resources and investing more of those resources into growth because it is safe to do so i.e., low risk of draining reserves. We would not have identified this crucial mechanism had we not employed an agent-based, life-history driven approach.

As shown in section 3ii, the overcompensation generated by the TAG model was nearly two to three times larger and two to five times faster, depending upon PCT treatments and site conditions, compared to observations from the Shawnigan Lake trial (Li et al., 2018), as has been observed in managed stands ( Figure 4 ). There are several biologically plausible explanations for this discrepancy. First, values from the simple inverse relationship between crowding and growth (Equation 5) that we employed may have been too high. It is possible that other factors in the environment (Bennett et al., 2017) may ameliorate stress from crowding. Second, we employed an aspatial model, for heuristic reasons, but in doing so, we implicitly assumed a sized-biased, scramble competition for those resources that are forfeited by the now-culled trees (Kato et al., 2007). This might be a reasonable approximation for mobile organisms but it is too simplistic for sessile ones and would lead to overestimates of resource accrual within the stand (Husband and Barrett, 1996). Third, we assumed no limiting resources (e.g. nitrogen) to growth other than those competed for (e.g. energy, water, etc.) in the aforementioned competition. Fourth, as is often the case in heuristic models, we assumed near instantaneous access to resources that were freed up via the thinning process. Fifth, while there is an implicit effect of year-after-year stress in our model via Equations 1-3, there is no explicit assumption that long-term stress has cumulative effects, which would make compensatory response more difficult (e.g., Lv et al. 22). If any or all of the above assumptions were relaxed, that would lead to slower and less uptake of resources by trees in the thinned plots and would necessarily reduce resource harvest and the LTH allocation value for the parameter φ, which, in turn, would reduce magnitude and speed of overcompensation. Further, we assumed that survivorship would be equal in both thinned and control plots however, windthrow could increase soon after thinning, further reducing the number of trees and thus initial thinned stand productivity (Bose et al., 2018). Finally, again, for heuristic reasons, we employed simple baseline growth curves whereas further modification using specific allometric descriptors (e.g., Kramer et al., 2018) may have reduced the gap even further between thinned and control stands though we would still expect to observe the signature of overcompensation. All of these simplifications could be remedied without changing the tenor of our approach but that would reduce the generality of our model, which was the primary goal of this proof-of-concept exercise.

Further to the discussion above, it is possible to increase the reality of the TAG model without losing its generality however doing so would sacrifice some of its elegance. For example, suppose we wish to move from single-species to mixed-species stands; that would require that we rewrite Equation 5 from the original τ = 1 − ( N y N 0 ) τ ' ψ (Equation 5) to τ = 1 − ( ∑ y = 1 V v y ∑ y = 1 V v 0 ) τ ' ψ (equation 5.1) and Equation 6 from the original ψ s , n , y = s n , y s ¯ y (Equation 6) to ψ s , v , n , y = s v , n , y ∑ v = 1 V s ¯ v , y λ v v N y (equation 6.1).

The interpretation is that we now have accounted for those different species (or varieties) of trees of which there are V types (v = 1, 2, 3…V) each with average size s in a given year y which must be further weighted by both their interspecific competitive ability (λ_v) against the focal individual and their relative density. To account for different size classes would require another summation term that would cycle through those size classes within each tree variety.

As can be seen from this discussion, it is certainly possible to modify the TAG model to include diverse forests, however the complexity from doing so increases dramatically particularly with regard to statistical analysis and/or visualization given that we would have now added more dimensions to the problem. The same can be said for some of the other simplifications that we have employed, for example, the assumption of a non-spatial world. Even if spatial representation is included, there are practical limits to degree of resolution. For example, in Guignabert et al. (2023) spatially-explicit, forest dynamics model, seedlings were binned into size classes within 10m² plots to reduce computational challenge.

At least two possible criticism of the TAG model need to be reconciled. First, TAG is generic and cannot be directly employed to make silvicultural decisions, for example, the exact amount of thinning that should be applied to a particular stand. We fully agree that this model is only the first step toward developing a comprehensive PCT silviculture system and as such we have described it as proof of concept. At this point, we cannot say how broadly applicable TAG will be across a range of tree species an environments however most of the parameters that we employed are measurable thus it should yield testable predictions. Second, we have made a number of reasonable assumptions while deriving the TAG model some of which are untested. We hope that development of this novel theory will generate interest in parameters that until-now have garnered little attention.

On the other hand, the major contribution of the TAG model development is to confirm that a simple theoretical model can provide meaningful explanation and predict the compensatory growth phenomenon including overcompensation in post-thinning forest stands and the possibility of stand productivity enhancement. In particular, overcompensation has been difficult to understand and explain in the past but here it is seen as the simple outcome from enhanced resource access and disproportionate investment resource accrual into growth by trees in thinned plots. As we have noted, the simulation results are consistent with existing observations in both short and long term, without digging dept into very detailed physiological processes or waiting for corresponding empirical observations become available. The results presented in this paper reinforce the conclusions from the statistical approach of the TreeCG model with a verification from running the operational TASS model (Li et al., 2022b). These duel approaches can help speed up our understanding and predictability of complex forest dynamics in both natural and managed stands.

A feature demonstrated in the TAG model development is the modeling approach of state-dependence (e.g., Li et al., 2020), which means that tree growth is determined by yearly changes in internal and external/environmental conditions/states, such as the internal states of tree size, age, and energy reserve, and the external states of available level of nutrients and fluctuations in weather variables. The advantage of this state-dependent modeling approach is the possible incorporation of adaptive growth/responses of trees during their long lifespan, with which the growth of trees is no longer a fixed function of tree age and site index as shown in many forest growth and yield models such as VDYP7 of BC and GYPSY of AB, Canada, but rather the results of a continuous trade-off process within a changing environment as indicated by the life-history theory. In this case, plant growth is closely related to resource capture and allocation over time (Thomas, 2002). In this paper, we finessed (i.e., no explicit optimization) state-dependent dynamic life history responses (Clark and Mangel, 2000) via the α, β and φ functions to describe the adaptive responses of individual trees to changing environment and this approach generates understandable, stand-level growth patterns. However, even using this finessed approach requires detailed calculations that take more time than that from conventional growth and yield models, but each simulation run can still be completed within a few seconds using a personal computer. The incorporation of these adaptive growth functions makes it suitable for exploring the effect of other time sensitive management operations such as fertilization and pruning, as well as addressing climate change related issues. (Note that our approach makes assumptions regarding the shape and direction of tree responses to their states. An assumption-free approach to calculating optimal response requires searching the fitness surface via backwards induction as explained by Mangel, Clark and others (Mangel and Clark 1986; McNamara and Houston, 1986). Another option for adding more realism and searching resultant complex fitness surfaces would be to employ genetic algorithms (e.g. Chubaty et al., 2013) again with the possible expense of losing some heuristics.

The state-dependent modeling approach is relatively new and unfamiliar to many foresters but promising in a sense of not only harmonizing the feature of conventional growth and yield modeling approach, but also opening the door for investigating how trees might respond to their changing environments, so as to be capable of addressing additional related issues. The state-dependent modelling approach has been widely employed in the field of behavioral ecology for representing how the behavior of animals could be explained as a result from response to specific environmental conditions. Many examples can be found from elucidating response to danger in tiny parasitoid wasps (Asobara tabida) (Roitberg et al., 2009) to explaining flexible torpor in insectivorous bats (Fjellddall et al., 2023) to evaluating impact of anthropogenic disturbance on adaptive foraging in beluga whales (Delphinapterus leucas) (McHuron et al., 2023). The research approach from behavioral ecology has recently been expanded to plants (Cahill, 2015; Simon et al., 2016), and this multi-disciplinary approach has demonstrated the benefits in enhancing our understanding of the dynamics of various ecosystems. Our current study is an example of this approach.

The general stand growth patterns simulated by our TAG model can be widely applied as long as the growth curves from natural stands display a sigmoid shape. Thus, the calibration of this life history theory-driven TAG model could be much simplified compared to statistics-driven models such as TreeCG and hence widely applied to many jurisdictions without detailed stand growth relationships. This feature could serve as a complement to experimental approaches in a sense of providing a theoretical prediction even before implementing these experiments. As such, simulation results from carefully designed TAG model experiment could help address silviculture-effect related issues such as optimal thinning strategy leading to maximized productivity and optimal spacing in achieving the goals of stand density management, and could provide useful insights for improved field experiment designs, as well as for supporting forestry policy development and practical forest management decisions to mitigate the challenges from wood supply shortage. This modeling approach can also be applied to investigating the effect of other disturbances such as fire and insect pests, as well as other species of plants and animals with sigmoid growth pattern, especially those with long lifespans.

As an example of future applications of the TAG model, consider climate change. As a starting point, we envision the following in an attempt to retain generality while focusing on specific problems: in our ‘weather’ parameter, we assumed that climate is stable and weather can be drawn from a likewise stable distribution. To deal with climate change, we would draw from two two-dimensional matrices (ς_s,y and ς_f,y) based upon data from Coupled Model Intercomparison Project (CMIP6) (https://cds.climate.copernicus.eu/cdsapp#!/search?type=dataset&text=CMIP6) or similar, that assumes a changing climate rather than a stable one. The challenge in doing so will be to extrapolate change in tree growth performance as a function of changing climate.

Finally, as is the case with many theoretical studies, we have generated many new questions, especially those related to PCT and compensatory growth. A few of those include: (i) Since the adaptive compensatory growth response depends greatly up release from crowding stress, can this feature be generalized (see Equations 5 and 6) or will they be unique among species and sites? (ii) Similarly, will such stress show as simple additive functions or will they display emergent properties in mixed stands? (iii) Since climate will change necessarily advance much faster than the evolution of tree life history traits unlike short-lived plants (e.g. Acoca-Pidelle et al., 2023), can we employ models such as TAG that include contemporary tree life history values to navigate the many management obstacles that await us, at least in the near term? (iv) Can the post-thinning stand dynamics be extended to post-disturbance, including fire and post-pest-driven stand dynamics, especially if trees are incapable of distinguishing the causes of partial mortality from different types of disturbances? (Note that Arimura (2021) reviewed plants’ response to herbivore elicitors that triggered plants unique defense mechanism that was identified as one of the compensatory growth mechanisms summarized by Prins and Verkaar (1992)). These are just a few of the challenges to sustainable forest management, however, the adaptive nature of tree growth presents an opportunity to develop effective dynamic strategies.