# Modelling the growth stress in tree branches: eccentric growth vs. reaction wood

Arnoul VAN ROOIJ, Eric Badel, Jean-François Barczi, Yves Caraglio,
Tancrede Almeras, Joseph Gril

<figure><img src="image_1.png" /></figure>

## To cite this version:

Arnoul VAN ROOIJ, Eric Badel, Jean-François Barczi, Yves Caraglio, Tancrede Almeras, et al. Modelling the growth stress in tree branches: eccentric growth vs. reaction wood. Peer Community Journal, 2023, 3, pp.e78. 10.24072/pcjournal.308/. hal-03748026v6

HAL Id: hal-03748026

https://hal.science/hal-03748026v6

Submitted on 21 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
RESEARCH ARTICLE

Published
2023-08-31

**Cite as**
Arnoul Van Rooij, Eric Badel, Jean-François Barczi, Yves Caraglio, Tancrede Almeras and Joseph Gril (2023)
*Modelling the growth stress in tree branches: eccentric growth vs. reaction wood*, Peer Community Journal, 3: e78.

**Correspondence**
vanrooij.arnoul@gmail.com

**Peer-review**
Peer reviewed and
recommended by
PCI Forest & Wood Sciences,
https://doi.org/10.24072/pci.forestwoodsci.100108

<figure><img src="image_3.png" /></figure>

This article is licensed
under the Creative Commons
Attribution 4.0 License.

# Modelling the growth stress in tree branches: eccentric growth vs. reaction wood

Arnoul Van Rooij<fnref n="1" />,², Eric Badel<fnref n="1" />,³, Jean-François Barczi³, Yves Caraglio³, Tancrede Almeras<fnref n="1" />,⁴, and Joseph Gril<fnref n="1" />,⁵

Volume 3 (2023), article e78

https://doi.org/10.24072/pcjournal.308

## Abstract

This work aims to model the mechanical processes used by tree branches to control their posture despite their increasing weight loading. The two known options for a branch to maintain its orientation are the asymmetry of maturation stress, including reaction wood formation, and eccentric radial growth. Both options can be observed in nature and influence the stress distribution developed in the branch each year. This so-called "growth stress" reflects the mechanical state of the branch. In this work, a growth stress model was developed at the cross-section level in order to quantify and study the biomechanical impact of each process. For illustration, this model was applied to branches of two 50-year-old trees, one softwood *Pinus pinaster*, and one hardwood *Prunus avium* (wild cherry tree), both simulated with the AmapSim discrete element software. For the wild cherry tree, the computed outputs highlighted that the eccentricity of radial growth seems to be as efficient as the formation of reaction wood to maintain postural control despite the increasing gravity. For the pine tree, eccentric radial growth appears to be less efficient than the formation of reaction wood. But although it does not necessarily act as a relevant lever for postural control, it greatly modifies the profile pattern of mechanical stress and could provide mechanical safety of the branch. This work opens experimental perspectives to understand the biomechanical processes involved in the formation of branches and their mechanical safety.

<figure><img src="image_4.png" /></figure>

<footnote n="1">¹Laboratoire de Physique et Physiologie Intégratives de l'Arbre en environnement Fluctuant, Site INRAE Crouël : 5 Chemin de Beaulieu, 63039 Clermont-Ferrand, France // Site Clermont : 1 Impasse Amélie Murat, 63178 Aubière, France, ²Institut Pascal, Campus Universitaire des Cézeaux, 4 avenue Blaise Pascal, TSA 60026 / CS 60026, 63178 Aubière Cedex, France, ³Botanique et Modélisation de l'Architecture des Plantes et des Végétations, Bd de la Lironde TA A-51/ PS 2 34398 Montpellier cedex 5, France, ⁴LMGC, CNRS, Université of Montpellier, Montpellier, France</footnote>
Contents

<table><tr><td>Abbreviations and notations (in order of occurrence)</td><td>2</td></tr><tr><td>1 Introduction</td><td>3</td></tr><tr><td>2 Material and methods</td><td>4</td></tr><tr><td>2.1 Numerical model</td><td>4</td></tr><tr><td>2.2 Realistic growth data</td><td>9</td></tr><tr><td>3 Results</td><td>11</td></tr><tr><td>3.1 Prunus avium</td><td>11</td></tr><tr><td>3.2 Pinus pinaster</td><td>13</td></tr><tr><td>3.3 Influence of branch orientation: the stationarity hypothesis</td><td>14</td></tr><tr><td>4 Discussion</td><td>15</td></tr><tr><td>4.1 Prunus avium: heavily loaded hardwood</td><td>15</td></tr><tr><td>4.2 Pinus pinaster: lightly loaded softwood</td><td>16</td></tr><tr><td>4.3 Influence of branch orientation: the stationary hypothesis</td><td>17</td></tr><tr><td>4.4 Vertical bending moment vs horizontal bending and torsion moments</td><td>17</td></tr><tr><td>4.5 Limits of the model</td><td>17</td></tr><tr><td>Conclusion and perspectives</td><td>19</td></tr><tr><td>Acknowledgements</td><td>19</td></tr><tr><td>Fundings</td><td>19</td></tr><tr><td>Conflict of interest disclosure</td><td>20</td></tr><tr><td>Data, script, code, and supplementary information availability</td><td>20</td></tr><tr><td>References</td><td>20</td></tr><tr><td>A</td><td>22</td></tr><tr><td>B</td><td>23</td></tr><tr><td>C</td><td>24</td></tr></table>

Abbreviations and notations (in order of occurrence)

<table><tr><td>NW, TW, CW</td><td>Normal Wood, Tension Wood, Compression Wood</td></tr><tr><td>(x, y, z)</td><td>Local reference coordinates associated with the section</td></tr><tr><td>O</td><td>Centre of the section</td></tr><tr><td>r, R</td><td>Radial polar coordinate and Radii of the cross section (m)</td></tr><tr><td>e(R), e(R)</td><td>Eccentricity at the stem radius R, integrated eccentricity up to r = R</td></tr><tr><td>(x', y', z')</td><td>Local reference coordinates associated with the section, centred on the pith</td></tr><tr><td>σ</td><td>Stress (MPa)</td></tr><tr><td>σ₀</td><td>Induced maturation stress (Mpa)</td></tr><tr><td>S</td><td>Cross section area (m²)</td></tr><tr><td>N, M</td><td>Loads (N): normal force parallel to z' and bending moment around y'</td></tr><tr><td>E</td><td>Module of elasticity in L direction (GPa): MOE</td></tr><tr><td>μ</td><td>Induced maturation strain</td></tr><tr><td>ε, a, b</td><td>Strain, at the center, local curvature</td></tr><tr><td>Kᵢ</td><td>Structural stiffness of the cross-section</td></tr><tr><td>Fᵢ</td><td>External coefficients (maturation and load)</td></tr><tr><td>θ</td><td>Circumferential position in section (rad)</td></tr><tr><td>σ₀(θ)</td><td>Maturation strain in the new ring at circumferential position θ</td></tr><tr><td>α</td><td>Mean maturation stress in the new ring</td></tr><tr><td>β</td><td>Differential stress in the new ring</td></tr><tr><td>Rₓ'y'</td><td>Radius of the cross section at the instant of appearance of the point (x', y')</td></tr><tr><td>λN, λM, νM, νN</td><td>Load power law: allometric coefficient</td></tr><tr><td>λb, νb</td><td>Change of curvature power law: allometric coefficient</td></tr></table>
<table><tr><td>σ<sub>NW</sub>, σ<sub>TW</sub>, σ<sub>CW</sub></td><td>Maturation stress in the normal wood, tension wood and compression wood</td></tr><tr><td>μ<sub>NW</sub>, μ<sub>TW</sub>, μ<sub>CW</sub></td><td>Maturation strain in the normal wood, tension wood and compression wood</td></tr><tr><td>N&#x20D7;<sub>n</sub>, M&#x20D7;<sub>n</sub></td><td>Loads of growth unit n: normal force and bending moment around y</td></tr><tr><td>N<sub>z</sub>, M<sub>x</sub>, M<sub>y</sub>, M<sub>z</sub></td><td>Loads of growth unit n: projection of N&#x20D7;<sub>n</sub> on z&#x20D7; and bending moment M&#x20D7;<sub>n</sub> around x&#x20D7;, y&#x20D7;, z&#x20D7;</td></tr><tr><td>m<sub>n</sub></td><td>Mass of the growth unit n (kg)</td></tr><tr><td>g</td><td>Gravity constant: g = 9.8 m.s<sup>-2</sup></td></tr><tr><td>G<sub>n</sub></td><td>Centre of gravity of the growth unit n</td></tr><tr><td>E<sub>d</sub>, E<sub>g</sub></td><td>Air-dry, green MOE</td></tr><tr><td>ρ</td><td>Density</td></tr><tr><td>μ<sub>strain</sub></td><td>1/10<sup>6</sup></td></tr><tr><td>D<sub>n</sub>, D<sub>n+1</sub></td><td>First and second diameter the growth unit n</td></tr><tr><td>D<sub>f</sub></td><td>Deflection of a growth unit</td></tr><tr><td>L<sub>n</sub></td><td>Length of the growth unit n</td></tr></table>

# 1. Introduction

From a mechanical point of view, wood in trees fulfils three major functions: construction of the tree structure, postural control of trunk and branches and breaking resistance to external stimuli (Thibaut, 2019). These three functions are provided by the way wood cells differentiate and accumulate during the wood formation process. Each axis of a tree can be considered as an inclined beam, consisting of a succession of conical growth units (Barthélémy and Caraglio, 2007). It is built in two steps: first, primary growth resulting in new growth units that increase the length of the initial axis; and secondary growth resulting in thickening of already existing units by addition of annual rings. These two interactive and additional processes lead to a specific pattern of mechanical stress, called 'growth stress', which can be analysed as the superposition of support stress and maturation stress (Archer, 1976; Fournier et al., 1991a). The support stress results from the continuous increase of the weight supported by the axis over the years. It reaches maximal levels close to the stem and vanishes near stem periphery, where the recently formed wood contributes to the support of recently produced biomass only. Maturation stress is set up at the end of the cell-wall maturation process, when molecular components such as lignin polymerise, generating growth forces by small dilatation or contraction restrained by the rigidity of the previously formed wood cells (Alméras and Clair, 2016). An evaluation of the maturation stress can be obtained by measuring the strain associated to stress release at stem periphery, where no support stress is present (Nicholson, 1971; Yang et al., 2005; Yoshida and Okuyama: 2002). The circumferential heterogeneity of this peripheral stress is needed to regulate stem curvature. In most cases, a tensile maturation stress is produced in the newly formed 'normal wood'. But observations on inclined trunks (Alméras et al., 2005; Coutand et al., 2007; Thibaut and Gril, 2021), seedlings (Hung et al., 2016) and branches (Fisher and Stevenson, 1981; Huang et al., 2010; Hung et al., 2017; Tsai et al., 2012) have evidenced a clear difference between hardwood and softwood. Hardwoods produce 'tension wood', inducing a much higher tensile stress on one side, while for softwood, a compressive stress is induced in 'compression wood'. The first pulls, the second pushes. In the most usual case of inclined stems restoring their vertical orientation, tension wood is formed on the upper side while compression wood is formed on the lower side of the trunk. But other situations can be encountered depending on the biomechanical requirements of the tree (Wang et al., 2009). In addition to their participation in the postural control of tree stems, these two types of so-called 'reaction wood' are characterised by specific anatomical pattern (not discussed here) and specific physical and mechanical properties.

As an alternative to complex experimental approaches, growth stress modelling plays an important role in the understanding of the phenomena involved in the orientation process of a stem. The history of biomechanical models began with (Kübler, 1959) who proposed an analytical formulation of growth stress for a perfect cylinder made of a homogeneous and transversally
isotropic wood. Later, (Archer and Byrnes, 1974) took into account an asymmetry of the matura-
tion stress, and (Fournier et al., 1991a,b) proposed a semi-incremental version of these models,
allowing to take into account a potential gradient of mechanical parameters (stiffness, matura-
tion). By associating their previous model to the loading induced by the tree weight, (Fournier
et al., 1994) made the connection between growth stress and stem orientation. This model has
been adopted and developed by several authors in order to study the orientation process of
stems. (Yamamoto et al., 2002) added a primary shoot and returned to curvature calculations.
(Alméras and Fournier, 2009) introduced the notion of gravitropic performance (capacity of the
tree to correct the bending moment induced by its weight) and proposed criteria of long-term
stability. (Huang et al., 2005) and (Alméras et al., 2005) impoved the model by introducing a sec-
ondary growth asymmetry and its resulting pith eccentricity, as well as stiffness heterogeneity,
allowing to quantify the effectiveness of eccentricity, maturation, stiffness gradient and initial
radius in the curvature regulation process. They highlighted that the main factor in the gravit-
ropic process is the spatial distribution of the maturation stress. Still in line with Fournier's 1994
model, (Alméras et al., 2018) recently developed analytical models of longitudinal growth stress,
taking into account different configurations, like eccentricity or maturation gradient, and evolu-
tion laws, like evolution of stiffness per additional layer. Finally, based on the same philosophy as
established by Kübler, tree-scale and finite-element models have emerged (Ancelin et al., 2004;
Fourcaud et al., 2003).

(Huang et al., 2010)'s model has been used to understand how eccentric growth and reaction wood are involved in branch orientation (Huang et al., 2010; Hung et al., 2017; Tsai et al., 2012; Wang et al., 2009), but all these studies were based on the current state of the branch, without consideration of the previous history: although some of them quantified the roles of maturation and eccentricity in the regulation of curvature, none did evaluate their capacity to ensure a given growth scenario.

Unlike trunks, which usually seek verticality, after the first stages of growth, branches tend to grow in a stationary way at a fixed angle to the vertical. Therefore, in this framework, we focus on understanding how branches can control their orientation, through the study of two growth parameters: eccentric growth and compression wood. The aim is to check by calculation what option is mechanically possible and safe for the branch. For this purpose, we developed a semi-incremental biomechanical model of growth stress at the cross section level that takes into account the eccentricity and maturation gradients during the construction of branches. Using the digital models of a softwood (Pinus pinaster Aiton) and a hardwood (Prunus avium L), the impact of each of these two growth parameters on the stress state was evaluated.

2. Material and methods

**2.1. Numerical model**

2.1.1. General hypotheses. The problem was set in the framework of the beam theory. From a geometrical point of view, branches generally show profiles that suit to this type of analytical framework: a slender shape and no important diameter variations. The shape effects due to twigs and other local biological phenomena (cavity, nodes, etc.) were neglected. The same set of hypotheses as in (Alméras et al., 2018) was adopted. In this study, we focused on the behaviour in the longitudinal direction (parallel to the main axis). Horizontal bending and torsion loads were not considered. Only the vertical bending moment (caused by the weight) was considered; these hypotheses on the loading modes are discussed later.

2.1.2. Geometrical settings. The object of study was the cross-section of a branch, placed within a plane locally orthogonal to the pith. The local reference frame of the section is (x̄, ȳ, z̄), with z̄ the longitudinal direction of the axis, and x̄ placed in a vertical plane and facing upwards (Fig 1). The shape of the cross-section was assumed to be circular at any stage of development, described by the successive depositions of wood rings. The term of 'ring' refers here to the volume occupied by wood cells produced by the cambium during a certain duration of time, not necessarily annual: it must be taken in a numerical meaning. These rings possibly could present
an eccentricity resulting from asymmetry of secondary growth. Since the model only takes into account vertical bending, the eccentricity was set along the x axis, as expressed by the following equation:

$$ (1) \qquad O(t) = \int_{0}^{R(t)} e(r)dr = \bar{e}R(t) $$

with $O(t)$ the position of the geometrical centre and $R(t)$ the radius of the section at time t, $e(r)$ the eccentricity when the stem radius was $r$ and $\bar{e}$ the integrated eccentricity up to $r = R$. The eccentricity can vary in the interval $[-1, 1]$. A zero eccentricity corresponds to a centred section, while -1 or 1 corresponds to maximum eccentricity resulting from secondary growth only on the lower or the upper side of the section, respectively. In the following, the position $x'$ in the pith reference frame is needed. By calling $x$ the vertical position in the geometrical reference frame, we deduce from equation (1):

$$ (2) \qquad x = x' - \bar{e}R $$

2.1.3. Computation of the mechanical behaviour. We developed a radial incremental method. For each radial increment, the longitudinal stress was computed in order to satisfy the static equilibrium of the cross section:

$$ (3a) \qquad \left\{ \begin{array}{l} \displaystyle \int_S \delta\sigma dS + \int_{\delta S} \sigma_0 dS = \delta N \end{array} \right. $$

$$ (3b) \qquad \left\{ \begin{array}{l} \displaystyle \int_S \delta\sigma x dS + \int_{\delta S} \sigma_0 x dS = -\delta M \end{array} \right. $$

where $S$ is the cross-section area, $\delta S$ is its increment, $\delta\sigma$ is the increment of stress $\sigma$ in the already formed wood, in response to the maturation stress $\sigma_0$ generated in the new wood layer. $\delta N$ and $\delta M$ are respectively the increment of external force $N$ and bending moment $M$, that are applied on the cross-section. For illustration, the geometric situation for $K$ rings and an increment of stem radius $\delta R$ is proposed in Fig 1.

The stress $\sigma$ is linked to the strain $\epsilon$ by a classical pre-stressed Hooke's law:

$$ (4) \qquad \sigma = E (\epsilon - \mu) = E\epsilon + \sigma_0 $$

where $E$ is the longitudinal Young's modulus, $\mu$ is the maturation strain and $\sigma_0$ is the maturation stress. In the context of the beam theory, the planar cross-sections remain so (Euler-Bernoulli assumption). The strain field is then described by the deformation $a$ at the centre of the pith and the curvature $b$ relative to the y-axis, as follows:

$$ (5) \qquad \delta\epsilon = \delta a + x\delta b $$

where $\delta\epsilon$, $\delta a$, $\delta b$ are the increments of $\epsilon$, $a$, $b$, respectively. The stress increment $\delta\sigma$, in the already formed wood, where no maturation occurs anymore, can then be deduced:

$$ (6) \qquad \delta\sigma = E\delta\epsilon = E(\delta a + x\delta b) $$

From these considerations, the system (3) becomes (details of the calculation are given in Appendix A):

$$ (7a) \qquad \left\{ \begin{array}{l} K_0\delta a + K_1\delta b = \delta F_0 \end{array} \right. $$

$$ (7b) \qquad \left\{ \begin{array}{l} K_1\delta a + K_2\delta b = \delta F_1 \end{array} \right. $$

with

$$ (8) \qquad K_0 = E\pi R^2, \quad K_1 = E\pi\bar{e}R^3, \quad K_2 = E\pi R^4 \left(\bar{e}^2 + \frac{1}{4}\right) $$
<figure><img src="image_5.png" /><figcaption>Figure 1 - Geometrical representation of a section with K numerical rings and a radial increment δR between rings (k - 1) and k.</figcaption></figure>

$$\delta F_0 = - \int_{\delta S} \sigma_0 dS + \delta N, \quad \delta F_1 = - \int_{\delta S} \sigma_0 x dS - \delta M$$

The calculation of the coefficients δF₀ and δF₁ depends on the formulation of the maturation stress. The maturation stress was assumed to vary circumferentially as follows:

$$\sigma_0(\theta) = \alpha + \beta \cos \theta \quad (9)$$

where the mean stress α and differential stress β were defined differently in softwood and hard-wood species:

$$\begin{array}{ll} (10a) & \left\{ \begin{array}{l} \text{Hardwood: } \alpha = \frac{\sigma_{TW} + \sigma_{NW}}{2}; \quad \beta = \frac{\sigma_{TW} - \sigma_{NW}}{2} \\ (10b) & \text{Softwood: } \alpha = \frac{\sigma_{CW} + \sigma_{NW}}{2}; \quad \beta = \frac{\sigma_{NW} - \sigma_{CW}}{2} \end{array} \right. \end{array}$$

σTW (resp. σCW) being the maturation stress in the tension wood (resp. compression wood), and σNW that in the opposite wood. One gets :

$$(11a) \qquad \delta F_0 = -\pi R (2\alpha + e\beta) \delta R + \delta N$$

$$(11b) \qquad \delta F_1 = -\pi R^2 (3\alpha e + e^2\beta + \beta) \delta R - \delta M$$

From equations (8), (11a) and (11b), the components of the system (7) are known. By inversion, δα and δb can be obtained according to the following equations (see details in Appendix B):
$$
\begin{array}{l}
(12a) \qquad \delta a = \frac{4}{ER} \left[ \left( 3e\bar{e} - 2e^2 - \frac{1}{2} \right) \alpha + \left( \bar{e}e^2 - e\bar{e}^2 + \bar{e} - \frac{e}{4} \right) \beta \right] \delta R \\
\qquad \qquad + \frac{4}{E\pi R^3} \left[ \bar{e}\delta M + \left( \bar{e}^2 + \frac{1}{4} \right) R\delta N \right] \\
(12b) \qquad \delta b = \frac{-4}{ER^2} \left[ (3e - 2\bar{e})\alpha + (e^2 - e\bar{e} + 1)\beta \right] \delta R - \frac{4}{E\pi R^4} (\delta M + \bar{e}R\delta N)
\end{array}
$$

Once $\delta a$ and $\delta b$ are known, the stress increment $\delta\sigma$ at any position given by $(x', y')$ can be obtained from equation (6). The stress distribution at this position can be obtained as the sum of the initial maturation stress and all the stress increments undergone by the material point since its creation.

$$
(13) \qquad \sigma(x', y', R) = \sigma_0(x', y') + \sum_{k=k_{x'y'}}^{K} \delta\sigma_k
$$

where $\delta R_k = r_k - r_{k-1}$ for a succession of ring radii $0 < r_0 < \dots < r_k < \dots < r_K = R$, $\delta\sigma_k$ is the corresponding increment, and $k_{x'y'}$ designates the ring containing the point.

Analytical formulations. Using equations (12b) and dividing by $\delta R$, we get the following equations when $\delta R$ tends to zero:

$$
\begin{array}{l}
(14a) \quad \displaystyle\frac{da}{dR} = \frac{4}{ER} \left[ \left( 3e\bar{e} - 2e^2 - \frac{1}{2} \right) \alpha + \left( \bar{e}e^2 - e\bar{e}^2 + \bar{e} - \frac{e}{4} \right) \beta \right] \\
\qquad \qquad + \frac{4}{E\pi R^3} \left[ \bar{e}\frac{dM}{dR} + \left( \bar{e}^2 + \frac{1}{4} \right) R\frac{dN}{dR} \right] \\
(14b) \quad \displaystyle\frac{db}{dR} = \frac{-4}{ER^2} \left[ (3e - 2\bar{e})\alpha + (e^2 - e\bar{e} + 1)\beta + \frac{1}{\pi R^2} \left( \frac{dM}{dR} + \bar{e}R\frac{dN}{dR} \right) \right]
\end{array}
$$

Using equation (13) and dividing again by a vanishing $\delta R$, we obtain the following equation in-
volving the partial derivative $\partial\sigma/\partial R$:

$$
(15) \qquad \sigma(x', y', R) = \sigma_0(x', y') + \int_{R_{x'y'}}^{R} \frac{\partial \sigma}{\partial R'}(x', R') dR'
$$

where $R_{x'y'}$ is the radius of the section at the instant of appearance of the point with coordinates $(x', y')$.

On the other hand, the expressions of axial force $N(R)$ and bending moment $M(R)$ are re-
quired to compute the evolution of the stress distribution in the cross section. For this purpose,
we assumed that both vary as a power function of the radius of the branch. This resulted in the
following allometric laws:

$$
\begin{array}{l}
(16a) \qquad \left\{ \begin{array}{l} N = \lambda_N R^{\nu_N} \\ M = \lambda_M R^{\nu_M} \end{array} \right. \\
(16b) \qquad \left\{ \begin{array}{l} N = \lambda_N R^{\nu_N} \\ M = \lambda_M R^{\nu_M} \end{array} \right.
\end{array}
$$

where $\lambda_{N,M}$ and $\nu_{N,M}$ are allometric coefficients. The $\lambda$-coefficients are directly proportional to the weight of the branch part supported by the cross section (the branch itself and the other axes of higher orders). The $\nu$-coefficients express the kinetics of the secondary growth: a small $\nu$ refers to an early secondary growth while a higher one refers to a later diameter increase.

The calculation of $\sigma$ requires also the knowledge of the temporal variation of the curvature $b$. In order to simplify the analyses, we mainly studied stationary cases, i.e. we assumed that the branch maintains its orientation and remains straight. This assumption results in $db/dR = 0$. Physiologically, this equation expresses that the branch always compensates its weight increment at each deposition step of a new wood layer, corresponding to an additional weight. However, we can consider two cases for which the branch does not build up in a stationary way: i) the passive
bending (under its own weight) case, and ii) the up-righting case (i.e. the action of maturation is stronger than the additional weight). In both cases, the resulting change in curvature has been modelled by (Alméras and Fournier, 2009) and (Alméras et al., 2018). It can then be written as follows:

$$ (17a) \qquad \left\{ \begin{array}{ll} \frac{db}{dR} = -4 \frac{\beta}{ER^2} & \text{Up-righting:} \\ \\ \frac{db}{dR} = 4 \frac{\lambda_M \nu_M}{E \pi} R^{\nu_M - 5} & \text{Passive bending:} \end{array} \right. $$

For the next computations, we used the following general law:

$$ (18) \qquad \frac{db}{dR} = \lambda_b R^{\nu_b} $$

As a remark, even if this equation bears some resemblance to (16), it does not express any notion of allometry and is used here only for convenience. Combining (14), (15), (16) and (18), the total stress can then be computed as:

$$ (19) \qquad \sigma^i(x', y', R) = \sigma_0^i(x', y') + S_1 \ln \left( \frac{R}{R_{x'y'}} \right) + \frac{S_2}{S_3} (R^{S_3} - R_{x'y'}^{S_3}) \\ \qquad \qquad \qquad \qquad + \frac{S_4}{S_5} (R^{S_5} - R_{x'y'}^{S_5}) + \frac{S_6}{S_7} (R^{S_7} - R_{x'y'}^{S_7}) x' $$

where $S_1 = 4 \left[ (3e\bar{e} - 2e^2 - \frac{1}{2})\alpha + (\bar{e}e^2 - e\bar{e}^2 + \bar{e} - \frac{e}{4})\beta \right]$ is driven by the maturation process, $S_2 = \frac{\lambda_N \nu_N}{\pi} (\bar{e}^2 + \frac{1}{4})$, $S_3 = \nu_N - 2$, $S_4 = \frac{4}{\pi} \lambda_M \nu_M \bar{e}$ and $S_5 = \nu_M - 3$ by the branch loading (geometric evolution of the branch), $S_6 = E \lambda_b$ and $S_7 = \nu_b + 1$ by the branch orientation.

For each radius r, the remaining unknowns are the mean stress α, the differential stress β and the eccentricity e. Equation (14b) can be rewritten as:

$$ (20) \qquad (3e - 2\bar{e})\alpha + (e^2 - e\bar{e} + 1)\beta = -\frac{1}{\pi R^2} \left( \frac{dM}{dR} + \bar{e}R \frac{dN}{dR} \right) - E \frac{R^2}{4} \frac{db}{dR} $$

Thus by fixing two parameters, the third is directly determined. The maturation parameters α and β are determined by the maturation stress σ$_{NW}$ in normal wood and σ$_{TW}$ or σ$_{CW}$ in reaction wood according to equation (10).

We considered two possible configurations for the simulations in the next section:

(1) First, we applied a constant eccentricity (so that $\bar{e} = e$) and we fixed the stress level in the normal wood. In that case, the maturation stress of the reaction wood was given by equations (10):

$$ (21a) \qquad \sigma_{TW} = \frac{-2}{\pi R^2 (1 + e)} \left( \frac{dM}{dR} + eR \frac{dN}{dR} \right) + \sigma_{NW} \left( \frac{1 - e}{1 + e} \right) + \lambda_b \left( \frac{ER^2}{2(1 + e)} \right) R^{\nu_b} $$

$$ (21b) \qquad \sigma_{CW} = \frac{2}{\pi R^2 (1 - e)} \left( \frac{dM}{dR} + eR \frac{dN}{dR} \right) + \sigma_{NW} \left( \frac{1 + e}{1 - e} \right) - \lambda_b \left( \frac{ER^2}{2(1 - e)} \right) R^{\nu_b} $$

(2) Second, we fixed the maturation parameters and we observed how the eccentric growth could, or not, maintain the orientation of the branch. In this configuration, equation (20) became a two degrees equation in e that could be solved numerically.

In these two configurations, using data on the support allometries $\lambda_N, \lambda_M, \nu_M, \nu_N$, we can calculate the stress in the reaction wood and/or the eccentricity with different ($\lambda_b, \nu_b$), then deduce the growth stress profile in the section (eq. 19). In the next part, we see how the allometric coefficients can be obtained from data generated by growth model.
## 2.2. Realistic growth data

2.2.1. *Tree architecture modelling.* Numerical experiments were carried out using two reference models: a softwood, maritime pine (*P. pinaster*) and a hardwood, wild cherry (*P. avium*) (Fig 2). Their growth follows the architectural model of Rauh (Hallé et al., 1978). This implies that the branching is rhythmic, the axes are monopodial and the branches are orthotropic. These digital trees were computed with the AmapSim software (Barczi et al., 2007). The input of this software are architectural parameters which were provided by observations and field studies: (Coudurier et al., 1993) and (Heuret et al., 2006) for *P. pinaster*, (Caraglio, 1996) and (Barthélémy et al., 2009) for *P. avium*. The choice of these species was based on the availability of temperate species in the AmapSim database. The two trees were modelled over 50 years in open-growth conditions, which did not correspond to the same ontogenic stage of development, but allowed both trees to be considered mature. In the final state, the pine (resp. cherry) was 18,2 m (resp. 14,1 m) high. The diameter at the base was 40 cm for both species. The insertion height of the first branch was 14,3 m for pine and 4,6 m for the birch. The branches of interest were the main branches; those that were directly attached to the trunk. In addition, only branches that were more than 20 years old were studied, so that they had a consistent loading history. Finally, 33 branches for the pine and 45 for the birch were selected. For each of the branch groups, the distributions of length L, radius r and insertion angle with the trunk θ are shown in Table 1.

<table><caption>**Table 1** – Geometric distribution of branches of interest</caption><thead><tr><td>Species</td><td>L<sub>m</sub> (m)</td><td>r<sub>m</sub> (m)</td><td>θ<sub>m</sub> (°)</td></tr></thead><tbody><tr><td><i>Pinus pinaster</i></td><td>5.3 ± 0.4</td><td>5.2 ± 0.3</td><td>70 ± 0.01</td></tr><tr><td><i>Prunus avium</i></td><td>7.9 ± 1.4</td><td>8.1 ± 0.7</td><td>80 ± 0.05</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Figure 2** – AmapSim representation of aerial architecture of 50-years old trees. (a) *Prunus avium* and (b) *Pinus pinaster*</figcaption></figure>

2.2.2. *Loading scenarii: allometric laws.* Each tree was composed of axes organised hierarchically according to their order: 1 for the tree seed, 2 for the trunk, 3 for the main branches, 4 for those attached to them, etc.. Each axis was described as a succession of growth units, which
were sections of cones, identified by a number (in order of appearance), and defined by a parent number, an order, a start and end diameter, the coordinates of the centres of both initial and final sections as well as their length (Fig 3). Note that the description provided by AmapSim did not include the internal structure of the growth units, such as eccentricity. To avoid unnecessary complications, the coordinates of the centres were taken as those of the pith. From the model data, the moments and normal forces can be computed in each growth unit, at any time of the tree's existence. In addition to a part of its own weight, each unit is subjected to the weight of its offsprings - this term referring to any growth unit that would fall if the studied one was cut. The normal force $\vec{N}_n$ and bending moment $\vec{M}_n$ supported by the growth unit n can be written:

$$ (22) \qquad \vec{N}_n = \frac{1}{2} m_n \vec{g} + \sum_{\substack{k>n \\ k \text{ offspring}}} m_k \vec{g} $$

$$ (23) \qquad \vec{M}_n = \overline{G_n G'_n} \wedge \left( \frac{1}{2} m_n \vec{g} \right) + \sum_{\substack{k>n \\ k \text{ offspring}}} \overline{G_n G'_k} \wedge (m_k \vec{g}) $$

where $G_n$ is the centre of gravity of the current growth unit, $G'_n$ is the centre of gravity of its second half. On the downstream side of $G_n$, $G_k$ is the centre of gravity of an offspring of number $k > n$, $m_i$ is the mass of growth unit i and $\vec{g}$ is the gravity vector. Once $\vec{N}_n$ and $\vec{M}_n$ were computed in the absolute coordinates used for the description of the whole tree, they were projected in the local coordinates system ($\vec{x}', \vec{y}', \vec{z}$), with $\vec{z}$ of the chosen cross section. In the following, in accordance with the development of the previous section, $N_z$ refers to the projection of $\vec{N}$ on $\vec{z}$ and $M_y$ to the projection of $\vec{M}$ on $\vec{y}'$.

Power law regressions were performed to recover the allometric coefficients $\lambda_M, \lambda_N, \nu_N, \nu_M$. A summary of the analysis process is proposed in Fig 3.

For the selected branch groups, the distribution of all allometric coefficients are presented in Fig 4.

In *Pinus*, there was a large variation in $\nu$-coefficient, with $\nu_M$ varying by almost a factor 2 in the studied sample; indicating very variable secondary growth kinetics. In *Prunus*, the range of variation of the allometric power coefficients was smaller, which depicted a higher homogeneity of secondary growth kinetics. For both species, a great diversity in $\lambda$- coefficients was observed, which depicted a significant variability in the loading history. This is particularly interesting as the branches showed geometric determinants that did not vary over large ranges (Table 1). Also, these coefficients do not appear to vary as a function of geometric parameters. This reflects the complexity of predicting the loading of a branch from the determinants of the main axis, and shows the importance of branching. In both cases, these variations in the $\lambda$-coefficients result in a factor 4 in the bending load between the lightly loaded and the heavily loaded branches. The average values of each allometric and final geometry, indicated in table 2, will be used for the simulations.

2.2.3. *Material data.* The stress values in the normal wood were fixed according to the average maturation strains advised by (Thibaut and Gril, 2021). Similarly, the green wood MOE were given by the correlation between dry and green MOE identified by (Thibaut and Gril, 2021): $E_g = 0.89 * E_d$. Dry MOE were provided by the tropix database of CIRAD (Gérard et al., 2011). The density of green wood was approximated by the density of water $\rho = 1000 \text{ kg.m}^{-3}$. These inputs are summarised in Table 2.

In the following section, the case of stationary growth ($\nu_b = 0$) will be considered principally and analysed thoroughly. Situations of changing curvature will be then considered briefly.
<figure><img src="image_7.png" /><figcaption>**Figure 3** – Illustration of the different steps of computation of the allometric law *M*<sub>y</sub> = *f*(*r*). First the modelled branch of *P. avium* (a) was divided into a succession of conical units (b). This allowed to compute the bending moment *M*. Then, the repetition of the computation each year, provided the relationship between the branch radius *r* and the bending moment *M*<sub>y</sub>, represented in the graph (c). The fitted curve provides the final allometric law.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 4** – Statistical distribution of allometric coefficients for modelled branches: (a) *Pinus pinaster* (b) *Prunus avium*. *λ*<sub>M,N</sub> : weight factor. *ν*<sub>M,N</sub> : kinetic of secondary growth factor (See equation 16 for more details).</figcaption></figure>

3. Results

**3.1. Prunus avium**

Fig 5 shows the results of the simulations for *P. avium*, when one of the factors (eccentricity or reaction wood) is set to zero. On Fig 5.a, the stress on the whole section is represented. In
<table><caption>**Table 2** – Mean input characteristics of the branches.</caption><thead><tr><th>Species</th><th>λ<sub>M</sub></th><th>λ<sub>N</sub></th><th>ν<sub>M</sub></th><th>ν<sub>N</sub></th><th>r</th><th>μ<sub>NW</sub></th><th>E<sub>d</sub></th><th>E<sub>g</sub></th></tr></thead><tbody><tr><td>Pinus pinaster</td><td>-6.4e6</td><td>5e4</td><td>3.2</td><td>2.5</td><td>0.05</td><td>410</td><td>8.8</td><td>7.9</td></tr><tr><td>Prunus avium</td><td>-2.6e7</td><td>9.5e3</td><td>3.6</td><td>2.7</td><td>0.08</td><td>712</td><td>10.2</td><td>9.1</td></tr></tbody></table>

$λ_{N,M}$ and $ν_{N,M}$ correspond to the allometric evolution of the normal load and bending moment, $r$ (m) is the radius at the basal part of the branch, $μ_{NW}$ ($μ_{strain}$) is the maturation strain in the normal wood, and $E_{d,g}$ (GPa) is the dry and green modulus of the material.

this case, the branch maintains its orientation through the formation of reaction wood only (no eccentric growth). The area near the pith is under compression (red), while the periphery is under tension (blue), with a higher tension on the upper side, allowing to maintain the orientation. Fig 5.b shows the interpolation of the stress distribution of Fig 5.a on the main axis $y=0$. Fig 5.c represents the maturation stress in the tension wood throughout the growth of the branch. This stress becomes greater as the branch grows. The symmetric case, with no formation of reaction wood but eccentric growth, is presented in Fig 5.d-f. This example illustrates that eccentricity alone could theoretically provide the orientation control. Fig 5.f shows the evolution of the eccentricity through the radial growth of the branch. Like the reaction wood stress in the previous case, the needed eccentricity increased when the branch grew. The pattern of stress distribution of Fig 5.d is quite similar to that in Fig 5.a, with compression near the pith and tension at the periphery, but the section is off-centred and the tension at periphery is the same all around the section, confirming the absence of reaction wood.

<figure><img src="image_9.png" /><figcaption>**Figure 5** – *Prunus avium*: the orientation of the branch is maintained by the two different processes: (a-c) the maturation stress provided by the formation of reaction wood; (d-f) the eccentric growth; (a,d) 2D visualisation of the growth stress in the whole section; (b,e) Growth stress profile on diameter $y=0$; (c,f) Parametric representation of the tropic driver, maturation stress (c) and eccentricity (f).</figcaption></figure>

Fig 6 shows the combination of the two factors. For each of them, three different scenarios were proposed. In Fig 6.a, the reaction wood controlled the orientation. Different eccentricities, ranging from -0.5 to 0.5 were imposed. The resulting stress patterns are represented in Fig 6.a.i: the higher tension on the upper side maintained the posture. This tension stress becomes higher
as the eccentricity becomes hypotrophic. This is confirmed by the evolution of the maturation stress induced by reaction wood through branch growth in Fig 6.a.ii. The situations where the eccentricity controlled the posture are shown in Fig 6.b. Where uniform tension was imposed ($\sigma_{TW} = 2\sigma_{NW}$, $\sigma_{TW} = 3\sigma_{NW}$), the eccentricity pattern became particular: we observed a decrease during the first year, followed by an increase (Fig 6.b.ii). This is explained by the growth scenario: at the beginning of the development, fixing a uniform reaction wood formation tended to right-up the stem, while a stationary orientation was imposed. Therefore, the eccentricity process counteracted this righting up movement, leading to the initial decrease. As the branch grew, the effect of reaction wood decreased and the branch tended to bend forward: the eccentricity counteracted this trend, leading to the final increase. This coordination problem may probably be specific to our scenario that imposed a stationary orientation throughout the entire growth the branch, including the first stages of development.

<figure><img src="image_10.png" /><figcaption>(a) Stress in tension wood $\sigma_{TW}$ is the main driver of postural control. Different eccentricities are applied : solid line, epitrophic eccentricity $e = -0.5$, dashed line, no eccentricity, $e = 0$ and dotted line, hypotrophic eccentricity $e = 0.5$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(b) Eccentricity $e$ is the main driver of postural control. Different maturation gradients $\sigma_{TW} - \sigma_{NW}$ are applied : solid line, no maturation gradient ($\sigma_{TW} = \sigma_{NW}$), dashed line, $\sigma_{TW} = 2\sigma_{NW}$ and dotted line, $\sigma_{TW} = 3\sigma_{NW}$.</figcaption></figure>

**Figure 6** – Different possible options to maintain the orientation of *Prunus avium* branches: stress in reaction wood (a) or eccentricity (b). For each option, the subfigure (i) represents the total stress on diameter y=0, and the subfigure (ii) shows the evolution of the tropic driver (e or $\sigma_{CW}$) vs. the radius R of the branch each year.

## 3.2. *Pinus pinaster*

For *Pinus pinaster*, we used the same approach. The set of results is presented in Fig 7 and Fig 8. When no eccentricity was involved (Fig 7.a-c), a light compression stress was observed on the lower side of the section. When the branch grew, the compression stress increased (Fig 7.c).
In case of no reaction wood formation (i.e. homogeneous maturation stress), the distributions of growth stress and eccentricity (Fig 7.d-f) were quiet similar to the previous example with the birch tree: tension in periphery, compression near the pith, and an increasing eccentricity with branch growth.

<figure><img src="image_12.png" /><figcaption>**Figure 7** – *Pinus pinaster*: the orientation of the branch is maintained by the two different processes: (a-c) the maturation stress provided by the formation of reaction wood; (d-f) the eccentric growth; (a,d) 2D visualisation of the growth stress in the whole section; (b,e) Growth stress profile on diameter y=0; (c,f) Parametric representation of the tropic driver, maturation stress (c) and eccentricity (f).</figcaption></figure>

The combination of the two factors is shown in Fig 8. As in *Prunus avium*, different eccentricities were imposed (Fig 8.a): the more epitrophic the eccentricity, the higher reaction wood maturation stress. Although the different compression stress levels were close, the dynamic of this stress within the growth of the branch was different (Fig 8.a.ii). Also, the stress pattern exhibits a difference near the pith (Fig 8.a.i), with some tension in this area for eccentricity *e* = 0.5. In case of a uniform reaction wood maturation (Fig 8.b), the profile remained quite similar to birch tree. We could not impose a too low compression stress because of the above-mentioned coordination incompatibility.

**3.3. Influence of branch orientation: the stationarity hypothesis**

In order to evaluate the relevance of the stationarity hypothesis, i.e., the branch keeps the same orientation, different growth scenarios were considered. For each branch, the case of active up-righting or passive bending was modelled (using equation 17). Passive bending was driven by increasing weight, calculated on the modelled branches. Up-righting was driven by the maturation gradient, which was set at 400 µstrain (σ ≈ 3.2 MPa) for pine and 700 µstrain (σ ≈ 6.2 MPa) for cherry (the gradient was of the order of magnitude of normal wood stress). The results are shown in Fig 9. In cherry, no major change of the stress pattern was observed. In contrast, the pattern changed greatly for pine. For a passive-bending branch, a 'V' profile and the absence of compression wood were observed. For up-righting, the previously-mentioned profile with tension at the pith was observed.
<figure><img src="image_13.png" /><figcaption>(a) Stress in compression wood $\sigma_{CW}$ is the main driver of postural control. Different eccentricities are applied : solid line, epitrophic eccentricity $e = -0.5$, dashed line, no eccentricity, $e = 0$ and dotted line, hypotrophic eccentricity $e = 0.5$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(b) Eccentricity $e$ is the main driver of postural control. Different constant maturation gradients $\sigma_{NW} - \sigma_{CW}$ are applied : solid line, no maturation gradient ($\sigma_{NW} = \sigma_{CW}$), dashed line, $\sigma_{CW} = 0$. Dotted line is the maximal absolute value of $\sigma_{CW}$ before divergence of the computation.</figcaption></figure>

**Figure 8** - Different possible options to maintain the orientation of *Pinus pinaster* branches: stress in reaction wood (a) or eccentricity (b). For each option, the subfigure (i) represents the total stress on the diameter $y=0$, and the subfigure (ii) shows the evolution of the tropic driver ($e$ or $\sigma_{CW}$) vs. the radius $R$ of the branch each year.

# 4. Discussion

## 4.1. *Prunus avium*: heavily loaded hardwood

Regarding the stress distribution (Fig 5), using either eccentric growth or reaction wood formation led to realistic stress values (except near the pith, which is an intrinsic limit of our model. This specific point is discussed in section *Limits of the model*). In the case with no eccentricity, a tensile strain of $\mu_{RW} \approx 2140\mu_{strain}$ ($\sigma_{TW} \approx 19.5$ MPa) was obtained, quite similar to literature values, for much smaller branches: on 4 cm plagiotropic branches of eight tree species, (Tsai et al., 2012) reported an average strain in reaction wood of around 2100$\mu_{strain}$, with some values up to $\approx 5000\mu_{strain}$. When combined with uniform eccentricity (Fig 6.a), it seems safer to promote the growth on the upper side: it minimises both high tensile stress and area with high compression stress. Interestingly, the worst case (hypotrophic eccentricity $e = -0,5$, more, solid line in
<figure><img src="image_15.png" /><figcaption>**Figure 9** – Distribution of growth stresses on diameter y=0 for different orientation scenarii. Solid line, up-righting movement; dashed line, stationarity; dotted line, passive bending.</figcaption></figure>

Fig 6.a) led to levels approaching the limits, but previously observed (Huang et al., 2005; Tsai et al., 2012): $\mu_{RW} \approx 4970\mu$ strain ($\sigma_{TW} \approx 45.4$ MPa)).

In the absence of reaction wood (Fig 5.d-f), the eccentricity alone ensured the orientation. The maximal value was around 0.6, which corresponds to an average eccentricity of 0.35. Although the absolute value is realistic, average eccentric found in literature is on the opposite side. For example, (Hung et al., 2017) performed measurements on 10 plagiotropic branches of *Koelreuteria elegans* (Seem.) A.C.Sm. The average radius was 2.6 cm, and the average eccentricity had an average value of -0.37, with a maximum at -0.54. Unpublished data on more than 150 branches from six different temperate species showed very different patterns, depending from the species, but eccentricity was mainly hypotrophic. These results suggest that eccentricity doesn't usually counteract gravity: its action counteracts the one of tension wood, leading to non-optimal stress patterns. An extensive measurement campaign on branches would be needed to clarify this point.

Also the combination of radial growth eccentricity with uniform maturation stress showed the same tendency as the dual combination (uniform eccentricity): a higher maturation stress led to a larger eccentricity. Comparing all simulations, the most theoretical optimal case was a constant positive eccentricity (dotted line in Fig 6.a), which is again not what is showed by experimental observations. It raises interesting question on the main mechanical driver of branch construction. From a biological point a view, it could be less "costly" to produce tension wood than to produce a large amount of wood for an eccentricity result. But this hypothesis was never investigated. Also, more work is needed to understand how tension wood and eccentricity are linked in angiosperm trees: since they may have some uncoordinated action, we can wonder if they have common triggered factors. Finally, these results suggest that growth eccentricity doesn't have the same role in branches and trunks: (Alméras et al., 2005) showed that eccentricity in leaning stems explains ≈ 29% of the of the curvature!

## 4.2. *Pinus pinaster*: lightly loaded softwood

First of all, the values of the stress distribution were much lower than for *Prunus avium*. This was explained by the size of the modelled branches: the average bending moment is much higher for cherry tree than for pine, by a factor of roughly 10 (see $\lambda_M$ and $\lambda_N$ in Table 2). The effect of each factor alone (Fig 7) suggested that maturation is a much more efficient option than eccentricity. To ensure the same growth scenario, the eccentricity alone rose to about 0.8, which is close to a theoretical limit, whereas maturation alone led to low maturation strains in compression wood (<500 $\mu$strain, corresponding to 4 MPa). Besides, this eccentricity was not in the
direction of what is commonly observed. This point remains logical, because without compres-
sion wood, the epitrophic eccentricity is the only way to counteract the effect of gravity.
A uniform eccentricity combined with reaction wood formation led to quite similar patterns (Fig
8.a): for this range of loading, the eccentricity had little influence on stress distribution. Consid-
ering that the density of elastic energy is proportional to the square of the stress, the pattern
produced a low level of stored elastic energy, possibly reducing the risk of mechanical failure.
Also, although eccentricity did not bring much variations in the value of the maturation stress, it
considerably modified the shape of the resulting stress profiles (Fig 8.a.i). Indeed, these profiles
can become 'crenellated' (Fig 8.a.i, dashed curve for zero eccentricity, solid curve for e = -0.5)
or include tension at the pith (dotted line for e = 0.5). It seems that before producing tension
in the pith, an efficient configuration could be reached by generating compression below the
pith and tension above. Ideally, this may be a very relevant option for branches. These results
about the mechanical strategies of branches should be confronted to experimental measure-
ments. Otherwise, these pattern changes could also be an optimisation of the residual strength
of wood: compression wood is known to have high compressive strength conferred by its high
lignin content and cell wall structure. Generating some tension at the pith allows the branch to
create more compression wood. To answer this question it would be necessary to take into ac-
count strength parameters in our stress computation model. Adding a damage-elastoplastic law
would also allow to study the effects of stress relaxation and to observe if some profiles, that
are here not optimal for maintaining the branch orientation, could possibly become optimal for
resisting breakage.

Using eccentricity combined to formation (Fig 8.b) leads to usual patterns, with compres-
sion near the pith, tension on the upper side and compression on the lower one. Eccentricity is
epitrophic: this is the opposite to what is usually observed: unpublished data on 20 branches
(average radius of 3 cm) of *Pinus nigra* showed an average eccentricity of -0.2. This non-intuitive
result is partly explained by our hypothesis of uniform stiffness, as will be discussed later.It is
also explained here by the change of sign between normal wood and compression wood. In the
early stages of growth, as long as the stress in the compression wood is lower than in the normal
wood, the best option to maintain the orientation is to do epitrophic eccentricity. Once the stress
in the compression wood becomes higher than in the normal wood, it is more efficient to do hy-
potrophic eccentricity. Our scenarios do not allow us to reach stress levels in compression that
are higher than the stress in normal wood. This is due to the above mentioned incompatibility
of our scenario.

4.3. Influence of branch orientation : the stationary hypothesis

In both trees, the orders of magnitude are compatible with a mechanical safety margin for
the branches. Apart from modified tropisms (change of light environment, weight change by
loss of part of the branch, etc.), the maintenance of the orientation is quite common for real
branches. However our simulations suggest that if, for any reason, branches need to modify
their orientation, they can do it without taking too much mechanical risk.

4.4. Vertical bending moment vs horizontal bending and torsion moments

One of the hypotheses of our model was that the vertical bending moment (M_y) prevails
over the torsional M_z and horizontal bending M_x moments. This allowed to consider only one
direction of eccentricity and to avoid all the non-linear terms generated by the torsional com-
ponents. We evaluated the maximum values of the three moments for all modelled branches of
each species for comparison purpose. The results are presented in Fig 10. They show that for
each species, the vertical moment displays much higher values than the torsional and horizontal
bending moments and validates our initial hypothesis.

**4.5. Limits of the model**

The hypothesis of homogeneous wood stiffness in the whole branch section is questionable.
Systematic stiffness differences have been observed between wood types (tension wood or com-
pression wood vs normal wood). (Alméras et al., 2005) have studied the variation of Young's
<figure><img src="image_16.png" /><figcaption>**Figure 10** – Comparison of maximum moments for the branches modelled in this research. Mx: horizontal moment; My: vertical moment; Mz: torsional moment.</figcaption></figure>

modulus in the section of leaning stems from 14 angiosperms and 3 gymnosperms, all coming from different families. For the angiosperms, the average Young's modulus of tension wood was higher than in normal wood by 15%, while for the gymnosperms, the Young's modulus was 38% lower in compression wood than in normal wood. This heterogeneity of rigidity plays a role in the postural control of the stems (Alméras et al., 2005; Huang et al., 2010; Hung et al., 2017). In our case, either a higher rigidity in tension wood or a lower one in compression wood would make the branch bend upward. In the current formulation of the model imposing an homogeneous stiffness, an almost equivalent effect would have been obtained by an initial offset in the eccentricity. Calling this offset tentatively 'compensating eccentricity' $e_c$ (Fig 11), the model computed a total eccentricity, $e$, combining $e_c$ and the "real" eccentricity needed to maintain the orientation. Therefore, in case in formation on one side, the eccentricity displayed needs to be offset by $e_c$ to correspond to more realistic situations. This explains, for instance, why the simulations for the softwood resulted in hypertrophic eccentricity while it is well-known that inclined softwood stems usually exhibit hypotropic eccentricity. Although data are missing to approximate the value of this parameter, and further work is needed to assess theoretically the possible equivalence between rigidity variations and eccentricity, the available information on relative stiffness of normal wood and suggests a more important effect in gymnosperms than in angiosperms.

<figure><img src="image_17.png" /><figcaption>**Figure 11** – How the hypothesis of a uniform wood stiffness impacts the initial position of the pith. $e_c$: pith offset equivalent to a centered pith with heterogenous wood properties; $E_{TW}$, $E_{NW}$, $E_{CW}$: Longitudinal Young's modulus of tension, normal, compression wood, respectively.</figcaption></figure>

The evaluation of the stress during the first stages of branch development is another issue of the model. In almost every stress profile, a tension or compression peak is generated in the
pith. It generally exceeds the wood strength, which is not compatible with branch sustainability. This point could be corrected in two ways. First, the role of the bark could be taken into account. Its importance in postural maintenance was clearly highlighted (Clair et al., 2019; Ghislain et al., 2019). Our model could include the mechanical action of bark in the early stages of branch development. This improvement would require additional data about the mechanical behaviour of the bark but would bring more realistic stress predictions and limit the artefacts at the pith. A second exciting perspective would be to take into account the elastoplastic behaviour of wood. By imposing a realistic plastic strain limit, the peak at the pith would then disappear; the increments would be spread over the middle part of the section, thus modifying the non-realistic patterns previously observed.

Finally, modelling the evolution of normal force and bending moment loads by allometric laws remains questionable. Indeed, the orientation of the branch may vary with time, which implies variations of the effect of weight. For example, modelling a constant increase of the normal force is inappropriate if the inclination of the branch decreases with time. An improvement of our model could be the construction of loads based on equivalent length allometries taking into account the mass of the branch, and the computation of the loads for each position in the right reference frame.

## Conclusion and perspectives

A semi-analytical growth stress model has been developed in the context of branch development. This model was applied to test the effectiveness of two well-known biomechanical processes of woody plants to control the orientation of their axes: eccentric radial growth and formation. For the hardwood branches, the computations highlighted that the eccentricity needed to maintain orientation did not corroborate the observations reported in literature. This suggests that this parameter probably provides another function than the orientation control, like the improved bending strength of the branch that provides it a greater mechanical safety. For the softwood branches, although the model showed that eccentric radial growth did not play a major role in maintaining the branch's orientation, it does modify the shape of the stress profiles in the cross section of the branch. A few odd and critical profiles, crenellated or with tension near the pith, have been identified. Their analysis provided exciting perspectives for further experimental works to gather real data.

Now that a complete model is available, it becomes crucial to start experimental investigations on branches in order to compare the outputs with real in situ observations. Especially, we need to evaluate the relevance of the different biological processes used by branches to ensure their mechanical sustainability over the years.

From a biological point of view, a key point for understanding branch sizing is the question of biomass costs. Building additional wood on one side or forming are carbon sinks with possible trade-offs. In order to investigate this point, our model could help by affecting a cost to the production of as well as to eccentric growth. The resulting computations could then help to understand the relevance of some options and would lead to coupling the biomechanical point of view to other biological considerations.

## Acknowledgements

Preprint version 3 of this article has been peer-reviewed and recommended by Peer Community In Forest and Wood Sciences (Dreyer, 2023; https://doi.org/10.24072/pci.forestwoodsci.100108).

## Fundings

Funding was provided by CNRS, we thank them for the PhD support granted to A. Van Rooj. No external funding used for the research.
**Conflict of interest disclosure**

The authors declare that they comply with the PCI rule of having no financial conflicts of interest in relation to the content of the article.

**Data, script, code, and supplementary information availability**

Script and codes are available in (Van Rooij, 2022; https://doi.org/10.57745/RECOB0).

**References**

Alméras T, Thibaut B, Gril J (2005). Effect of circumferential heterogeneity of wood maturation strain, modulus of elasticity and radial growth on the regulation of stem orientation in trees. *Trees* **19**, 457-467. https://doi.org/10.1007/s00468-005-0407-6.

Alméras T, Fournier M (2009). Biomechanical design and long-term stability of trees: Morphological and wood traits involved in the balance between weight increase and the gravitropic reaction. *Journal of Theoretical Biology* **256**, 370-381. https://doi.org/10.1016/j.jtbi.2008.10.011.

Alméras T, Clair B (2016). Critical review on the mechanisms of maturation stress generation in trees. *Journal of The Royal Society Interface* **13**, 20160550. https://doi.org/10.1098/rsif.2016.0550.

Alméras T, Jullien D, Gril J (2018). Modelling, Evaluation and Biomechanical Consequences of Growth Stress Profiles Inside Tree Stems. In: Plant Biomechanics: From Structure to Function at Multiple Scales. Ed. by Anja Geitmann and Joseph Gril. Cham: Springer International Publishing, pp. 21-48. https://doi.org/10.1007/978-3-319-79099-2_2.

Ancelin P, Fourcaud T, Lac P (2004). Modelling the biomechanical behaviour of growing trees at the forest stand scale. Part I: Development of an incremental transfer matrix method and application to simplified tree structures. *Annals of Forest Science* **61**, 263-275. https://doi.org/10.1051/forest:2004019.

Archer RR (1976). On the distribution of tree growth stresses. II. Stresses due to asymmetric growth strains. *Wood Science and Technology* **V10**, 293-309. https://doi.org/10.1007/BF00350833.

Archer RR, Byrnes FE (1974). On the distribution of tree growth stresses - Part I: An anisotropic plane strain theory. *Wood Science and Technology* **8**, 184-196. https://doi.org/10.1007/BF00352022.

Barczi JF, Rey H, Caraglio Y, Reffye P, Barthélémy D, Dong QX, Fourcaud T (2007). AmapSim: A Structural Whole-plant Simulator Based on Botanical Knowledge and Designed to Host External Functional Models. *Annals of Botany* **101**, 1125-1138. https://doi.org/10.1093/aob/mcm194.

Barthélémy D, Caraglio Y (2007). Plant Architecture: A Dynamic, Multilevel and Comprehensive Approach to Plant Form, Structure and Ontogeny. *Annals of Botany* **99**, 375-407. https://doi.org/10.1093/aob/mcl260.

Barthélémy D, Caraglio Y, Sabatier S (2009). 4.1 Crown Architecture of Valuable Broadleaved Species. *Valuable broadleaved forests in Europe* **22**, 87.

Caraglio Y (1996). Le développement architectural du merisier. Forêt Entreprise 107, 72-80.

Clair B, Ghislain B, Prunier J, Lehnebach R, Beauchêne J, Alméras T (2019). Mechanical contribution of secondary phloem to postural control in trees: the bark side of the force. *New Phytologist* **221**, 209-217. https://doi.org/10.1111/nph.15375.

Coudurier T, Barthelemy D, Chanson B, Courdier F, Loup C (1993). Premier résultats sur la modélisation du pin maritime pinus pinaster ait.(pinecae). *Architecture des arbres fruitiers et forestiers*, 306.

Coutand C, Fournier M, Moulia B (2007). The Gravitropic Response of Poplar Trunks: Key Roles of Prestressed Wood Regulation and the Relative Kinetics of Cambial Growth versus Wood Maturation. *Plant Physiology* **144**, 1166-1180. https://doi.org/10.1104/pp.106.088153.
Dreyer E (2023). An important contribution to the description of growth stresses in branches of adult trees based on a new model and an optimisation process with digitised branches. Peer Community in Forest and Wood Sciences, 100108. https://doi.org/10.24072/pci.forestwoodsci.100108.

Fisher JB, Stevenson JW (1981). Occurrence of Reaction Wood in Branches of Dicotyledons and Its Role in Tree Architecture. Botanical Gazette **142**, 82-95. https://doi.org/10.1086/337199.

Fourcaud T, Blaise F, Lac P, Castéra P, Reffye P (2003). Numerical modelling of shape regulation and growth stresses in trees. *Trees* **17**, 31-39. https://doi.org/10.1007/s00468-002-0203-5.

Fournier M, Chanson B, Guitard D, Thibault B (1991a). Mécanique de l'arbre sur pied : modélisation d'une structure en croissance soumise à des chargementspermanents et évolutifs. 1. Analyse des contraintes de support. Ann. For. Sci. **48**, 513-525. https://doi.org/10.1051/forest:19910503.

Fournier M, Chanson B, Guitard D, Thibault B (1991b). Mécanique del'arbre sur pied :modélisation d'une structure en croissance soumise à des chargements permanents et évolutifs.2. Analyse tridimensionnelle des contraintesde maturation, cas du feuillu standard. Ann. For. Sci. **48**, 527-546. https://doi.org/10.1051/forest:19910504.

Fournier M, Baillères H, Chanson B (1994). Tree biomechanics: growth, cumulative prestresses, and reorientations. *Biomimetics* **2**, 229-251.

Gérard J, Guibal D, Paradis S, Vernay M, Beauchêne J, Brancheriau L, Châlon I, Daigremont C, Détienne P, Fouquet D, Langbour P, Lotte S, Thévenon MF, Méjean C, Thibaut A (2011). *Tropix* 7. https://doi.org/10.18167/74726F706978.

Ghislain B, Alméras T, Prunier J, Clair B (2019). Contributions of bark and tension wood and role of the G-layer lignification in the gravitropic movements of 21 tropical tree species. *Annals of Forest Science* **76**, 107. https://doi.org/10.1007/s13595-019-0899-7.

Hallé F, Oldeman RA, Tomlinson PB (1978). *Tropical trees and forests: an architectural analysis*. Springer Verlag.

Heuret P, Meredieu C, Coudurier T, Courdier F, Barthélémy D (2006). Ontogenetic trends in the morphological features of main stem annual shoots of Pinus pinaster (Pinaceae). *American Journal of Botany* **93**, 1577-1587. https://doi.org/10.3732/ajb.93.11.1577.

Huang YS, Chen SS, Kuo-Huang LL, Lee CM (2005). Growth strain in the trunk and branches of *Chamaecyparis formosensis* and its influence on tree form. *Tree Physiol* **25**, 1119-1126. https://doi.org/10.1093/treephys/25.9.1119.

Huang YS, Hung LF, Kuo-Huang LL (2010). Biomechanical modeling of gravitropic response of branches: roles of asymmetric periphery growth strain versus self-weight bending effect. *Trees* **24**, 1151-1161. https://doi.org/10.1007/s00468-010-0491-0.

Hung LF, Tsai CC, Chen SJ, Huang YS, Kuo-Huang LL (2016). Study of tension wood in the artificially inclined seedlings of Koelreuteria henryi Dummer and its biomechanical function of negative gravitropism. *Trees* **30**, 609-625. https://doi.org/10.1007/s00468-015-1304-2.

Hung LF, Tsai CC, Chen SJ, Huang YS, Kuo-Huang LL (2017). Strain distribution, growth eccentricity, and tension wood distribution in the plagiotropic and orthotropic branches of Koelreuteria henryi Dummer. *Trees* **31**, 149-164. https://doi.org/10.1007/s00468-016-1464-8.

Kübler H (1959). Studien über Wachstumsspannungen des Holzes III. Längenänderungen bei der Wärmebehandlung frischen Holzes. Holz Rohst Werkst **17**(3), 77-86.

Nicholson JE (1971). A rapid method for estimating longitudinal growth stresses in logs. Wood Science and Technology **5**, 40-48. https://doi.org/10.1007/BF00363119.

Thibaut B (2019). *Three-dimensional printing, muscles, and skeleton: mechanical functions of living wood*. *Journal of Experimental Botany* **70**, 3453-3466. https://doi.org/10.1093/jxb/erz153.

Thibaut B, Gril J (2021). *Tree growth forces and wood properties*. Peer Community Journal **1**, e46. https://doi.org/10.24072/pcjournal.48.

Tsai CC, Hung LF, Chien CT, Chen SJ, Huang YS, Kuo-Huang LL (2012). *Biomechanical features of eccentric cambial growth and reaction wood formation in broadleaf tree branches*. *Trees* **26**, 1585-1595. https://doi.org/10.1007/s00468-012-0733-4.

Van Rooij A (2022). *Growth stress model*. Version V1. https://doi.org/10.57745/RECOBO.
Wang Y, Gril J, Sugiyama J (2009). Is the branch of Viburnum odoratissimum var. awabuki reaction wood? Unusual eccentric growth and various distributions of growth strain. In: 6th Plant Biomechanics Conference, pp. 328-334.

Yamamoto H, Yoshida M, Okuyama T (2002). Growth stress controls negative gravitropism in woody plant stems. Planta 216, 280-292. https://doi.org/10.1007/s00425-002-0846-x.

Yang JL, Baillères H, Okuyama T, Muneri A, Downes G (2005). Measurement methods for longitudinal surface strain in trees: a review. Australian Forestry 68, 34-43. https://doi.org/10.1080/00049158.2005.10676224.

Yoshida M, Okuyama: T (2002). Techniques for Measuring Growth Stress on the Xylem Surface Using Strain and Dial Gauges. Holzforschung 56, 461-467. https://doi.org/10.1515/HF.2002.071.

## Appendix A.

The calculation of integrals of system (3) requires preliminary elements. The situation of two consecutive rings is represented in Fig 12. Each position x in the geometrical reference frame is expressed with respect to the position x' in the pith reference frame according to the equation:

$$ (24) \qquad x = r \cos \theta = x' - \bar{e}R $$

with r the radius at time t and R the radius at the final time.

<figure><img src="image_18.png" /><figcaption>Figure 12 - Representation of two consecutive rings and the elements needed to calculate $\delta R(\theta)$</figcaption></figure>

The integrals of system (3) are computed as follows:

$$ \begin{aligned} \int_s \delta\sigma ds &= \int_s E [\delta a + (x + \bar{e}.R)\delta b] r \delta r d\theta \\ &= E\pi R^2 (\delta a + \bar{e}.R\delta b) \end{aligned} $$

$$ \begin{aligned} \int_s x' \delta\sigma ds &= \int_s [\delta a + (x + \bar{e}.R)\delta b] [x + \bar{e}.R] r \delta r d\theta \\ &= E\pi R^3 \left[ \bar{e}\delta a + R \left( \bar{e}^2 + \frac{1}{4} \right) \delta b \right] \end{aligned} $$

The tangential distribution of the radius increment $\delta R(\theta)$ is required to compute the maturation terms. Applying $\overrightarrow{O_R M} - \overrightarrow{O_{R+dR} M} = \overrightarrow{O_R O_{R+dR}}$ (Fig 12):
$$
\begin{array}{l}
(25a) \qquad \begin{cases} [R + \delta R(\theta)] \cos \theta - (R + \delta R) \cos (\theta + \delta\theta) = e_R \delta R \\ [R + \delta R(\theta)] \sin \theta - (R + \delta R) \sin (\theta + \delta\theta) = 0 \end{cases} \\
(25b)
\end{array}
$$

By setting $\delta\theta \to 0$, it comes:

$$
\begin{array}{l@{\qquad}l}
(26a) & \left\{ \begin{aligned} \cos(\theta + \delta\theta) &= \cos\theta - \sin\theta\delta\theta \\ \sin(\theta + \delta\theta) &= \sin\theta + \cos\theta\delta\theta \end{aligned} \right. \\
(26b) &
\end{array}
$$

Substituting (26) into (25), and combining (25a) and (25b), $\delta R(\theta)$ can finally be written as:

$$
(27) \qquad \boxed{\delta R(\theta) = \delta R [1 + e_R \cos \theta]}
$$

Then:

$$
\begin{align*}
\int_{\delta s} \sigma_0^i ds &= \int_{\delta s} \sigma_0^i(\theta) R \delta R(\theta) d\theta \\
&= \int_{\delta s} [\alpha + \beta \cos \theta] [1 + e \cos \theta] R \delta R(\theta) d\theta \\
&= \pi (2\alpha + e\beta) R \delta R
\end{align*}
$$

$$
\begin{align*}
\int_{\delta s} x' \sigma_0^i ds &= \int_{\delta s} \sigma_0^i(\theta) (x + e.R) R \delta R(\theta) d\theta \\
&= R^2 \delta R \pi (3\alpha e + \beta e^2 + \beta)
\end{align*}
$$

Appendix B.

The matrix system (7) becomes:

$$
\begin{array}{l@{\qquad}l}
(28a) & \left\{ \begin{aligned} \delta a &= \frac{\delta F_0 K_2 - \delta F_1 K_1}{K_0 K_2 - K_1^2} \\ \delta b &= \frac{\delta F_0 K_1 - \delta F_1 K_0}{K_1^2 - K_0 K_2} \end{aligned} \right. \\
(28b) &
\end{array}
$$

Then, numerators and denominators are calculated separately:

$$
K_0 K_2 - K_1^2 = E^2 \pi^2 R^6 \left(\bar{e}^2 + \frac{1}{4}\right) - E^2 \pi^2 R^6 \bar{e}^2 = \frac{(E \pi R^3)^2}{4}
$$

$$
\begin{align*}
\delta F_0 K_2 - \delta F_1 K_1 &= E \pi^2 R^5 \left[ -(2\alpha + e\beta) \left(\bar{e}^2 + \frac{1}{4}\right) + \bar{e} \left(3\alpha e + \beta e^2 + \beta\right) \right] \delta R \\
&\quad + E \pi R^3 \left[ R \delta N \left(\bar{e}^2 + \frac{1}{4}\right) + \bar{e} \delta M \right] \\
&= E \pi^2 R^5 \left[ \alpha \left(3e\bar{e} - 2\bar{e}^2 - \frac{1}{2}\right) + \beta \left(\bar{e}e^2 - e\bar{e}^2 + \bar{e} - \frac{e}{4}\right) \right] \delta R \\
&\quad + E \pi R^3 \left[ R \delta N \left(\bar{e}^2 + \frac{1}{4}\right) + \bar{e} \delta M \right]
\end{align*}
$$

$$
\begin{align*}
\delta F_0 K_1 - \delta F_1 K_0 &= E \pi^2 R^4 \left[ -\bar{e} (2\alpha + e\beta) + (3\alpha e + e^2\beta + \beta) \right] \delta R + E \pi R^2 [\bar{e} R \delta N + \delta M] \\
&= E \pi^2 R^4 \left[ \alpha (3e - 2\bar{e}) + \beta (1 + e^2 - e\bar{e}) \right] \delta R + E \pi R^2 [\bar{e} R \delta N + \delta M]
\end{align*}
$$

Putting the calculations together, system (28) becomes:
$$ 
\left\{ 
\begin{aligned}
\delta a &= \frac{4}{ER} \left[ \alpha \left( 3e\bar{e} - 2\bar{e}^2 - \frac{1}{2} \right) + \beta \left( \bar{e}e^2 - e\bar{e}^2 + \bar{e} - \frac{e}{4} \right) \right] \delta R \\
&\quad + \frac{4}{E\pi R^3} \left[ R\delta N \left( \bar{e}^2 + \frac{1}{4} \right) + \bar{e}\delta M \right] \\
\delta b &= \frac{-4}{ER^2} \left[ \alpha (3e - 2\bar{e}) + \beta (1 + e^2 - e\bar{e}) \right] \delta R + \frac{-4}{E\pi R^4} \left[ e\bar{R}\delta N + \delta M \right]
\end{aligned}
\right. 
$$

## Appendix C.

The following calculus is based on Fig 3.b. To get the vertical bending moment $M_y$ of unit n (eq 23), one need the calculation of each volume $V_n$ and center of gravity $G_n$. Let name $D(z)$ the diametral extension of the cone. It comes:

$$ (30) \qquad V_n = \int_0^{L_n} \frac{\pi D(z)^2}{4} dz $$

where $D(z) = D_n + \left(\frac{D_{n+1}-D_n}{L_n}\right) z$. One gives

$$ (31) \qquad O_n G_n = \frac{1}{V_n} \int_0^{L_n} \frac{\pi D(z)^2}{4} zdz $$

Setting $\gamma = \frac{D_{n+1}-D_n}{D_n}$ and $\xi = \frac{L_n}{z}$, equations (30) and (31) become:

$$ V_n = \frac{\pi D_n^2 L_n}{4} \int_0^1 (1 + \gamma \xi)^2 d\xi = \frac{\pi D_n^2 L_n}{4} \left( 1 + \gamma + \frac{\gamma^2}{3} \right) $$

$$ O_n G_n = \frac{1}{V_n} \frac{\pi D_n^2 L_n^2}{4} \left( \frac{1}{2} + \frac{2\gamma}{3} + \frac{\gamma^2}{4} \right) $$

So, finally, $O_n G_n$ can be written:

$$ (32) \qquad O_n G_n = \frac{L_n}{2} \left( \frac{1 + \frac{4}{3}\gamma + \frac{1}{2}\gamma^2}{1 + \gamma + \frac{1}{3}\gamma^2} \right) $$