Modeling tree crown dynamics with 3D partial differential equations

We characterize a tree's spatial foliage distribution by the local leaf area density. Considering this spatially continuous variable allows to describe the spatiotemporal evolution of the tree crown by means of 3D partial differential equations. These offer a framework to rigorously take locally and adaptively acting effects into account, notably the growth toward light. Biomass production through photosynthesis and the allocation to foliage and wood are readily included in this model framework. The system of equations stands out due to its inherent dynamic property of self-organization and spontaneous adaptation, generating complex behavior from even only a few parameters. The density-based approach yields spatially structured tree crowns without relying on detailed geometry. We present the methodological fundamentals of such a modeling approach and discuss further prospects and applications.

1. Introduction

In terms of model scale, light sensitive functional-structural tree growth modeling has experienced the emergence of various trends. Organ-level approaches bring about a high precision of physiological processes, averting inaccuracies and effects of scale non-invariance, which may arise from simplifications in larger-scale approaches. Moreover, arbitrary small-scale biophysical or biochemical processes can in principle be readily induced. The LIGNUM model (Perttunen et al., 1996; Sievänen et al., 2008), the model by Sterck et al. (2005) or the L-Peach model (Allen et al., 2005) are examples of this model category. Local light interception (cf. also Chelle and Andrieu, 2007 for a methods review) determines the production and allocation of biomass. Their detail of physiological and morphological processes is at the same time the drawback of these models. On the one hand, the large number of organs implies high computational costs—all the more in competition scenarios of multiple trees. On the other hand, their detail can make these models susceptible to the propagation of errors, which could have been compensated for in an averaging rough scale approach.

Other organ-level models like Greenlab (Yan et al., 2004; Cournède et al., 2006) a priori focus on the topology in terms of the plant's structure. This implies the inability of easily taking physiological and structural responses to varying local light conditions into account. As one consequence, the approach cannot be straightforwardly applied to scenarios of competition for light: An additional yet non local competition index provides for this (Cournède et al., 2008). In another context of formal grammars used for tree growth simulation, Kurth and Sloboda (1999) present the 2D concept of the shadow-relevant cone of a shoot in order to take local light conditions into account, which in turn affect the rewriting system. Comparable methods have been used by Purves et al. (2007) and Takenaka (1994).

Sonntag (1996) presents a model in which the spatial motion and allocation of leaf area is based on heuristic rules on a 2D cellular space. Though quite different in terms of formalism, his approach bears conceptual resemblances to ours.

Models with a rougher scale use to impose certain characteristics on the crown shape. For instance, the Balance model Grote and Pretzsch (2002) and the model by Sorrensen-Cothern et al. (1993) describe the crown shape in terms of disk-like horizontal layers. This technique implies advantages with regard to the computational speed as well as a general robustness, compared to small-scale models. Yet it does so at the expense of a thorough plastic spatial crown structure. When applied to competition scenarios, these models often make use of empirically fitted competition indices (e.g., Pretzsch, 1992).

While being affiliated more with the latter forestry models in the attempt to describe crown structure and dynamics macroscopically for applications at the stand level, the present approach attempts a middle course between the fine organ-centered way and the rougher pre-imposition of the crown form. We characterize a tree's spatial foliage distribution via its local leaf area density. The focusing on this variable circumvents difficulties in terms of robustness in the geometrically detailed models accounting for individual leaf positions. At the same time, the locality allows for arbitrary spatial structures.

Applying Beer Lambert's law allows to express the local light conditions within a crown as a function of the local leaf area density. Aiming at an increased future light interception, local leaf area density is assumed to tend to move toward the light. This approach, which induces the spatial expansion of the crown, translates directly into a partial differential equation. Details are specified based on mass conservation and optimization considerations. The technique notably allows to account for a local and spontaneous adaptiveness with regard to changing environmental light conditions.

The mathematical approach of density-based partial differential equations has long established itself in spatial biology and ecology (cf. e.g., Okubo and Levin, 2002). In the context of macroscopic individual plant modeling, so far notably in the form of diffusion equations, it has proven applicable for root growth and proliferation (see Page and Gerwitz, 1974 for the original approach, Reddy and Pachepsky, 2001 for a review of later developments and Dupuy et al., 2010 for a current advance). A 2D diffusion approach for the foliage of crops with the objective to model competition in different field densities, without considering the vertical dimension, was presented by Beyer et al. (2014). Partial differential equations can generate, even from only a few terms, complex self-adapting dynamics, which is indeed present in biological systems. Attempts to reproduce this by simpler terms often requires a larger model framework and set of parameters. The latter, in turn, requires large data sets, which are often not available.

In this article we will exemplify the use of the leaf area density-based partial differential equation approach by means of a simplified model setup presented in section 2. The continuity of the approach suggests embedding into in the context of continuously growing trees sensu (Hallé et al., 1978), i.e., “with no marked endogenous cessation of extension [and] a more or less constant production of leaves and/or shoots throughout the year” (Barthélémy and Caraglio, 2007).

Since merely selected dynamics are presented, without constructing a realistic model, we settle for illustrations of key qualitative properties of the approach instead of quantitative data comparison. Throughout the article, we will point out possibilities of extending and customizing process assumptions while preserving the advantages of the overall methodological framework. Extending future steps are discussed in section 3.

2. Model Framework

For clarity's sake all recurrently appearing variables and parameters and their definitions are listed in Table 1. We describe the spatial foliage distribution of a tree canopy by means of the (local) leaf area density α(x, t) ≥ 0 (in m² m⁻³) in a point x = (x₁, x₂, x₃) ∈ ℝ³ (each entry in m) at a time t ≥ 0 (in years), i.e., the spatial density of the total one-sided green leaf area (previously considered e.g., by Sinoquet et al., 2005). The map α(·, t): ℝ³ → ℝ, x ↦ α(x, t) is continuous for any t. For brevity, we will use the notion leaf density in place of leaf area density.

TABLE 1

Table 1. Key model variables and parameters.

We aim to describe the evolution t ↦ α(·, t). To this end we will first determine the biomass production B of a tree, which, along with the senescence S of old biomass, allows to describe the net biomass increment of a tree corresponding to a leaf density α(·, t) at a given time t:

$\frac{\partial }{\partial t}m(t)=B(t)-S(t)$

This, in turn, will be distributed among foliage and sapwood according to the pipe model theory (Shinozaki et al., 1964), specified in the subsequent paragraph. The coaction of foliage allocation and senescence in the case of a continuously growing tree induces what we abstractly interpret as a continuous motion of α(·, t), in particular directed toward the light, aiming at an increased future biomass production. This perspective leads to the description of the course of t ↦ α(·, t) essentially by means of a continuity equation.

Sapwood Associated to a Leaf Density: The pipe model theory by Shinozaki et al. (1964) allows to determine the sapwood mass corresponding to an arbitrary leaf density α(·, t). In the present context, the theory states that for any point x ∈ ℝ³ with α(x, t) > 0, a sapwood pipe, in charge of the transport of water and nutrients, leads from x down to the roots with a length denoted by ‖x‖_ϒ, its cross-sectional area being proportional to the leaf density α(x, t) at its tip via a constant P. The mass of the pipe leading to x then equals

$\underset{\text{cross-sectional area}}{\underbrace{\alpha (x,t)·P}}· \underset{\text{length}}{\underbrace{\Vert x{\Vert }_{{\rm Y}}}}·WD,$

WD denoting the wood density (in g m⁻³).

It is worth mentioning that, motivated by the limitations to the pipe model theory, pointed out e.g., by Tyree (1988), Pouderoux et al. (2001), and Deleuze and Houllier (2002), generalizations have been suggested: A noteworthy approach is the one by Bouchon et al. (1997), who, based on an allocation perspective, reason that the pipe does not necessarily have a constant cross-sectional area along its path, but more generally a one exponentially decreasing toward the stem base. This principle integrates in our context with only minor technical changes. Likewise, the approach by Letort et al. (2008), parametrically combining the pipe model approach and a uniform, common pool sapwood allocation, would be feasible.

As for the pipe's length ‖x‖_ϒ from x to the root tip, for the sake of simplicity, here we follow the multi-species approach used by Sonntag (1996), based on the branch architecture of coniferous species and the assumption that root length equals branch length, resulting in

$\begin{array}{cc}\Vert x{\Vert }_{{\rm Y}}:=\left|x_3\right|+2·\Vert (x_1,x_2)\Vert .& (1)\end{array}$

A more accurate choice for a particular species can be made by taking specific characteristics of its branching geometry (such as branching angles) and topology into account.

Finally, the sum of foliage and sapwood mass is given by

$\begin{array}{cc}\underset{\text{foliage mass}}{\underbrace{{\displaystyle {\int }_{{ℝ}^3}\frac{\alpha (x,t)}{SLA}} dx}}+\underset{\text{sapwood mass}}{\underbrace{{\displaystyle {\int }_{{ℝ}^3}\alpha }(x,t)·P·\Vert x{\Vert }_{{\rm Y}}·WDdx}}& (2)\end{array}$

where SLA denotes the specific leaf area (in m² g⁻¹).

2.1. Biomass Production

We determine the amount of biomass produced through photosynthetic activity by a given leaf density. To this, we take direct and diffuse radiation into account, cf. Fu and Rich (1999), using the horizontal celestial coordinate system, with ℝ² × {0} being the local horizon of a tree rooting in (0, 0, 0), and the vector (1, 0, 0) pointing north. The unit directional vector σ(t) ∈ S²₊: = {x ∈ ℝ³: ‖x‖ = 1, x₃ ≥ 0}, under which the sun is seen from the tree at daytime t reads σ(t) = (cos(−Az(t)) · cos(Alt(t)), sin(−Az(t)(t)) · cos(Alt(t)), sin(Alt(t))), where Alt(t) and Az(t) denote the time dependent altitude and azimuth, respectively.

For diffusive radiation, a uniform diffuse model (uniform overcast sky) is applied, in which incoming diffuse radiation is assumed to be the same from all sky directions. Let PAR_dir(t) and PAR_diff(t) (in MJ m^-2) denote the photosynthetically active direct and diffuse radiation at time t, respectively. Then the total radiation from direction v∈ S²₊ at time t is

with the indicator function _A(x) : = 1 if x ∈ A and 0 else.

2.1.1. Isolated tree

Incoming radiation is partly intercepted by the tree's foliage and partly passes through it. The fraction between 0 and 100% of radiation from direction v∈ S²₊ which actually reaches the point x ∈ ℝ³ can be determined using Beer-Lambert's law, where foliage characterized by leaf density acts as a light absorbing medium with locally varying α-concentration. This fraction reads

$\begin{array}{cc}\exp \text{}(-\Lambda ·{\displaystyle {\int }_{x+{ℝ}_{+}·v}\lambda }(\xi ,v)·\alpha (\xi ,t)d\xi )& (3)\end{array}$

where the extinction coefficient Λ ≤ 1 represents the mean light transmittance of foliage (Monteith, 1969; Nouvellon et al., 2000) and

$\lambda (x,v):=N(x)·v$

takes into account the angle between the sun ray and foliage in x, N(x) ∈ S²₊ denoting the unit normal to the plane in which foliage in x lies. N(x) can be chosen according to leaf angle distribution models without further ado (Wang et al., 2007) provide a review.

Assuming that local biomass production is proportional to the total locally intercepted ratiation via an energy efficiency μ (in g MJ⁻¹) (Monteith, 1977), the instantaneous local biomass production b(x, t) (in g m⁻³ s⁻¹) in x ∈ ℝ³ at time t for the leaf density α(·, t) reads

$\begin{array}{cc}\begin{array}{l}b(x,t)=\mu ·\hfill \\ {\displaystyle \underset{S_{+}^2}{\int }\underset{\text{radiation from direction }v\text{intercepted in }x}{\underbrace{\lambda (x,v)·\alpha (x,t)·\underset{\text{radiation reaching }x\text{from direction }v}{\underbrace{\exp (-\Lambda ·{\displaystyle {\int }_{x+{ℝ}_{+}·v}\lambda }(\xi ,v)·\alpha (\xi ,t)d\xi )·\text{PAR}(v,t)}}}}}dv. \end{array}& (4)\end{array}$

2.1.2. Population of trees

At a time t, let α₁(·, t), …, α_n(·, t) denote the leaf densities of n trees, shading each other and competing for light. Then the instantaneous local biomass production of tree i ∈ {1, …, n} in x generalizes to

$\begin{array}{l}{b}_i(x,t)=\mu ·{\displaystyle \underset{S_{+}^2}{\int } }\frac{{({\lambda }_i(x,v)·{\alpha }_i(x,t))}^2}{{\displaystyle \sum _{j=1}^n{\lambda }_j}(x,v)·{\alpha }_j(x,t)}·\hfill \\ \exp \left(-\Lambda ·{\displaystyle \underset{x+{ℝ}_{+}·v}{\int }{\displaystyle \sum _{j=1}^n{\lambda }_j}}(\xi ,v)·{\alpha }_j(\xi ,t) d\xi \right)\text{}·\text{}\text{PAR}\text{}(v,t)dv\hfill \end{array}$

The fraction $\frac{{\lambda }_i(x,v)·{\alpha }_i(x,t)}{{\displaystyle {\sum }_{j=1}^n{\lambda }_j}(x,v)·{\alpha }_j(x,t)}$ is the part of the incoming radiation in x that is attributed to tree i's foliage in x. In particular it reduces to 1 if two trees' crowns do not occupy common space.

2.2. Dynamics

2.2.1. Mass balance

We determine the instantaneous change in living mass (conductive sapwood and foliage mass) due to the production of new, and the senescence of old biomass. For convenience we assume that the senescence of leaves as well as the loss of conductivity of sapwood depend on time only. If sapwood and foliage that have existed for τ_W and τ_F years become nonconductive and senescent, respectively, then the living mass at time t reduces by

$\begin{array}{l}S(t):={\displaystyle {\int }_{{ℝ}^3}\frac{\frac{\partial }{\partial t}\alpha (x,t-{\tau }_F)}{SLA}}dx+{\displaystyle {\int }_{{ℝ}^3}\frac{\partial }{\partial t}}\alpha (x,t-{\tau }_W)·\hfill \\ P·\Vert x{\Vert }_{{\rm Y}}·WDdx.\hfill \end{array}$

At the same time it increases by the total instantaneous biomass production at t, i.e., B(t) : = ∫_ℝ³ b(x, t)dx. Thus we have

$\begin{array}{cc}\frac{\partial }{\partial t}m(t)=B(t)-S(t)& (5)\end{array}$

for the change in living mass at time t. A priori, the possibility $\frac{\partial }{\partial t}m(t)<0$ is not excluded for arbitrary parameters. However, since m(t) > S(t), it follows that m(t) > 0 and thus B(t) > 0 for all t ≥ 0.

2.2.2. Leaf density dynamics

With this information on the global mass of the tree at hand, we consider the variable

$\widehat{\alpha }(x,t)=\frac{\alpha (x,t)}{{\displaystyle {\int }_{{ℝ}^3}\alpha }(x,t)dx},$

i.e., the leaf density modulo mass, or (mass-)relative leaf density, with the property that ∫_ℝ³$\widehat{\alpha }$(x, t)dx = 1 for all t, which is understood as an indicator of the spatial structure of the real leaf density α(·, t). Instead of describing the course of α(·, t) directly, we do so for $\widehat{\alpha }$(·, t) and deduce α(·, t) as α(·, t) = λ(t) · $\widehat{\alpha }$(·, t)}, where λ(t) is chosen such that the living mass corresponding to λ(t) · $\widehat{\alpha }$(·, t) (cf. (2)) equals indeed m(t) yielded by (5). Thus

$\begin{array}{cc}\lambda (t)=\frac{m(t)}{{\displaystyle {\int }_{{ℝ}^3}\frac{\widehat{\alpha }(x,t)}{SLA}} dx+{\displaystyle {\int }_{{ℝ}^3}\widehat{\alpha }}(x,t)·P·\Vert x{\Vert }_{{\rm Y}}·WDdx}.& (6)\end{array}$

The continuity of the growth process of the trees we consider suggests a description of the course of $\widehat{\alpha }$(·, t) in terms of a (mass-conserving) continuity equation,

$\begin{array}{cc}\frac{\partial }{\partial t}\widehat{\alpha }(x,t)={\nabla }_x·\varphi (x,t),& (7)\end{array}$

in which the relative leaf density is subject to a transport motion induced by a continuous flux ϕ: ℝ³ × ℝ₊ → ℝ³, itself determined by α(·, t). The idea of this transporting flux ϕ is that it incorporates both the effect of the allocation of new leaves and of the abscission of old ones on the spatial structure of foliage, which in combination, induces what we describe as an abstract motion of leaf density.

A predominant driver in the spatial dispersal of leaf density is the local expansion toward the light, aiming at an increase in future light interception. We formally embed this factor in the above framework, where it will take on the role of the flux ϕ. For some given leaf density α(·, t) let L: ℝ³ × ℝ₊ → ℝ be defined by

$\begin{array}{cc}\begin{array}{l}L(x,t)={\displaystyle {\int }_{S_{+}^2}\lambda (x,v)·\exp \left(-\Lambda ·{\displaystyle {\int }_{x+{ℝ}_{+}·v}\lambda (\xi ,v)·\alpha (\xi ,t)d\xi )}\right)}·\hfill \\ \text{PAR}(v,t)dv\hfill \end{array}& (8)\end{array}$

The function L measures the intercepted light in x per m² leaf area. The gradient ∇_xL points locally in the direction of the greatest rate of increase of intercepted light. In addition, similar to Beyer et al. (2014) we define the flux to correspond to the existing leaf density in x, so that finally we have

$\begin{array}{cc}\varphi (x,t)=k·\alpha (x,t)·{\nabla }_{x}L(x,t)& (9)\end{array}$

for a mobility constant k. This local gradient approach is motivated by “the observation that a tree is capable of acquiring also gradient information about its environment and that growth might be directed along these gradients (Schmidt and Wulff, 1993; Aphalo and Ballare, 1995)” (Sonntag, 1996).

The term ∇_x · ϕ describing a movement of leaf density toward the light contains spatial derivates of a function of integrals over α, which makes (7) a partial integro-differential equation.

2.3. Simulations

In this section we illustrate some structural properties of the model. For convenience we simplified radiation to be vertical only. Some details of the numerical implementation of the model are presented in the appendix.

Figures 1, 2 illustrate the evolution of leaf density in the course of time, as well as the vector field ϕ. The term is sensitive to any change in the local light conditions induced by shading; ϕ instantaneously adapts and points in the direction of the greatest light increase. Aside the general spatial expansion toward the light, we notably observe the predominant presence of foliage at the crown hull rather than its interior.

FIGURE 1

Figure 1. (x₁, x₃)-cross section of a leaf density-characterized crown in the course of time, subject to the present model dynamics, corresponding to t = 10 (A), 20 (B), 30 (C) and 40 (D) years with simulation parameters adopted from Letort et al. (2008). A darker color indicates a higher leaf density α. The arrows in (D) indicate the vector field ϕ.

FIGURE 2

Figure 2. 3D view of Figure 1D.

2.3.1. Population of trees

Together with the adjustments in terms of biomass production addressed at the end of section 2.1, the approach generalizes to competition scenarios when α(·, t) is replaced by ∑_{i ≥ 1} α_i(·, t) in (8), alongside the different initial states α₁(·, 0), α₂(·, 0), …. We illustrate the dynamic effects to which competition gives rise by means of a simplified scenario: Consider a sufficiently large stand, in which the trees' stem bases, as a point set in ℝ² × {0}, generate a Voronoi-tesselation which is regular. If radiation is assumed to be radially symmetric (as done e.g., by Perttunen et al., 1998), the analysis of all competing trees reduces to that of a single one for which periodic boundary conditions for (7) are added on the boundary of the tree's 3D cell, i.e., the extension of the appropriate 2D Voronoi cell in the x₃-dimension.

Periodic boundary conditions induce that the light conditions on the other side of the boundary are considered identical to those within, accounting for another tree growing in equal measure and shading its environment.

This implies that, when, for simplicity, further assuming the essential light incidence to be vertical, periodic boundary conditions reduce to no flux conditions on the boundary: There, due to the identical light conditions on the other side of the boundary, the light gradient ∇_xL, governing the flux ϕ, changes from pointing further outwards to a zero flux.

Figure 3 shows the different stages of this scenario for an underlying square tesselation.

FIGURE 3

Figure 3. (x₁, x₃)-cross section and 3D view of the stages in a competition scenario at t = 10 (A), 20 (B), 30 (C) and 40 (D) years with periodic bounday conditions at 3.5 m. Due to properties of self-symmetry, visible in Figure 1, these stages are in fact qualitatively similar no matter the cell size.

The overly sharp boundary between two crowns is a consequence of the unrestrained mobility of foliage in this simplified approach, and would fade when additional features corresponding to a branch structure are included, cf. section 3.

The conspicuous concentration of leaf density at the upper edges of the canopy in the competition case (cf. Figure 3C) results from the two factors of (i) the regular expansion of leaf density uninfluenced by the boundary condition and (ii) the immigration of leaf density from lower regions whose horizontal expansion had abruptly turned into a vertical one after reaching the Voronoi cell boundary. We are unable to provide evidence or counterevidence to determine whether this phenomenon corresponds to reality. In any case, as observable in Figure 3D, this is only an intermediate state, before in the long-term a homogeneously distributed high concentration of leaf density in the top canopy establishes itself.

Wheras an isolated tree grows radially symmetric, this symmetry is eventually broken for a non-circular Voronoi cell, such as the square used just now. Figure 4 illustrates this.

FIGURE 4

Figure 4. Cross-section of (A) a diagonal plane (B) the (x₁, x₃)-plane at t = 15 years. (C) 3D view.

3. Discussion and Future Prospects

The aim of this article was to show how local leaf area density, a concept opposing geometrically detailed individual leaf configurations, can be used to approach macroscopic tree crown dynamics. Its integration into the formal framework of partial differential equations allowed to rigorously formulate the growth toward light. In the simulations we observed the generation of self-organization and adaptiveness that come along with this modeling approach. The simplistic model framework was meant to draw attention to the key mechanisms and their dynamic effects.

Foliage dynamics are by nature coupled to the tree's branch structure, which has not been taken into account on a topological or geometric level in this article. Future work, with the aim of introducing more spatial heterogeneity to the approach, begins here. Taking merely the stem and the most vigorous primary branches into account while leaving the finer structures to the leaf area density concept may already suffice to tackle crown plasticity satisfactorily. While in the present simplified model the motion of local leaf area density is governed by light only and otherwise unrestrained, a simple branch architecture can add sort of directional inertia to that motion, channeling local leaf area density in several major directions, resulting in heterogeneous foliage clustering, representing individual branches. In particular, this includes the incorporation of genetically predetermined branching angle spectra. Introducing a branch structure, even if only a rough one, is moreover accompanied by a refinement of the pipe length term (1), governing the distribution of mass between foliage and wood according to the pipe model theory, thus determining secondary growth.

Taking into account the organization of growth units, i.e., as weakening our assumptions on neoformation and polycyclism in section 2.1.1, or considering immediate vs. delayed bud outbreak, would bring the model closer to actual tree architecture dynamics.

Alongside the phototropism considered in this article, more biomechanical constraints that have feedback influences on the growth, such as hydraulic aspects sensu (Ryan and Yoder, 1997; Tyree and Zimmermann, 2002), in particular in the context of growth limitation, the avoidance of interlocked growth due to mechanical stress of touching branches (Oliver and Larson, 1990) or gravitropism represent perspectives for model extensions.

The application and validation of a refined model based on the theoretical framework presented in this paper will benefit from empirical data on local leaf area density. Conceivable ways to obtain this include the following three, which are currently being practically explored: Firstly, from the direct recordings of local light intensities at various positions {x⁽¹⁾, …, x⁽ⁿ⁾} ⊂ ℝ³ within a canopy, the map

$\begin{array}{l}\{x^{(1)},\dots ,x^{(n)}\}\to {ℝ}_{\ge 0}\hfill \\ x\mapsto \alpha (x)\hfill \end{array}$

for a discrete, but arbitrary fine domain can be obtained by applying Beer-Lambert's law in the reverse way. Secondly, high-definition multi-directional 3D terrestrial laser scan data and appropriate skeletonization algorithms allow to relate to a leafless tree a set of cylinders representing branch segments (Raumonen et al., 2013). Automatically removing all scan points corresponding to cylinder diameters (i.e., branch thicknesses) above a certain threshold allows to spatially isolate recent shoots, for which then a relation to foliage can be assumed. Thirdly, scanning a tree at the beginning of spring during the process of bud opening, when—in contrast to the case of a fully developed canopy—laser rays still reasonably penetrate the crown, and deducting from this image the point cloud yielded by a scan of the completely leafless tree, may allow to obtain a local bud density, from which the the leaf density can be deduced. Thus, assessing the spatial foliage distribution in terms of local leaf area density in a functional-structural crown dynamics model can make good use of this data type for model calibration and validation.

The property of locally and spontaneously adapting to changing light conditions suggests that the present partial differential equation approach can be applied to competition scenarios both in pure as in mixed tree groups. Empirical findings, based on laser scans, about the plasticity Bayer et al. (2013) may thus be approached from a functional-structural modeling point of view. These perspectives are currently being explored.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors would like to thank Pascal Laurent-Gengoux and two anonymous reviewers for their helpful remarks. This work was supported by a doctoral scholarship of the the Heinrich Böll Foundation.

References

Allen, M., Prusinkiewicz, P., and Dejong, T. (2005). Using L-Systems for Modeling the Architecture and Physiology of Growing Trees: The L-PEACH Model. New Phytol. 166, 869–880. doi: 10.1111/j.1469-8137.2005.01348.x

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Aphalo, P., and Ballare, C. (1995). On the importance of information-acquiring systems in plant-plant interactions. Funct. Ecol. 9, 5–14. doi: 10.2307/2390084

CrossRef Full Text

Barthélémy, D., and Caraglio, Y. (2007). Plant architecture: a dynamic, multilevel and comprehensive approach to plant form, structure and ontogeny. Ann. Bot. 99, 375–407. doi: 10.1093/aob/mcl260

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bayer, D., Seifert, S., and Pretzsch, H. (2013). Structural crown properties of Norway spruce (Picea abies [L.] Karst.) and European beech (Fagus sylvatica [L.]) in mixed versus pure stands revealed by terrestrial laser scanning. Trees 27, 1035–1047. doi: 10.1007/s00468-013-0854-4

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Beyer, R., Etard, O., Cournède, P.-H., and Laurent-Gengoux, P. (2014). Modeling spatial competition for light in plant populations with the porous medium equation. J. Math. Biol. doi: 10.1007/s00285-014-0763-1. [Epub ahead of print].

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bouchon, J., de Reffye, P., and Barthélémy, P. (1997). Modélisation et Simulation de L'architecture des Végétaux. Montpellier: INRA.

Chelle, M., and Andrieu, B. (2007). Modelling the light environment of virtual crop canopies. Funct. Struct. Plant Model. Crop Product. 22, 75–89. doi: 10.1007/1-4020-6034-3_7

CrossRef Full Text

Cournède, P.-H., Kang, M., Mathieu, A., Barczi, J.-F., Yan, H., Hu, B., et al. (2006). Structural factorization of plants to compute their functional and architectural growth. Simulation 82, 427–438. doi: 10.1177/0037549706069341

CrossRef Full Text

Cournède, P.-H., Mathieu, A., Houllier, F., Barthélémy, D., and de Reffye, P. (2008). Computing competition for light in the GREENLAB model of plant growth: a contribution to the study of the effects of density on resource acquisition and architectural development. Ann. Bot. 101, 1207–1219. doi: 10.1093/aob/mcm272

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Deleuze, C., and Houllier, F. (2002). A flexible radial increment taper equation derived from a process-based carbon partitioning model. Ann. Forest Sci. 59, 141–154. doi: 10.1051/forest:2002001

CrossRef Full Text

Dupuy, L., Gregory, P., and Bengough, A. (2010). Root growth models: towards a new generation of continuous approaches. J. Exp. Bot. 61, 2131–2143. doi: 10.1093/jxb/erp389

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Fu, P., and Rich, P. (1999). “Design and implementation of the Solar Analyst: an ArcView extension for modeling solar radiation at landscape scales,” in Proceedings of the 19th Annual ESRI User Conference, (San Diego).

Grote, R., and Pretzsch, H. (2002). A model for individual tree development based on physiological processes. Plant Biol. 4, 167–180. doi: 10.1055/s-2002-25743

CrossRef Full Text

Hallé, F., Oldeman, R., and Tomlinson, P. (1978). Tropical Trees and Forests. Berlin: Springer-Verlag. doi: 10.1007/978-3-642-81190-6

CrossRef Full Text

Kurth, W., and Sloboda, B. (1999). Tree and stand architecture and growth described by formal grammars. II. Sensitive trees and competition. J. Forest Sci. 45, 53–63.

Letort, V., Cournède, P.-H., Mathieu, A., de Reffye, P., and Constant, T. (2008). Parametric identification of a functional-structural tree growth model and application to beech trees (Fagus sylvatica). Funct. Plant Biol. 35, 951–963. doi: 10.1071/FP08065

CrossRef Full Text

Monteith, J. (1969). “Light interception and radiative exchange in crop stands,” in Physiological Aspects of Crop Yield, eds J. Eastin, F. Haskins, C. Sullivan, and C. van Bavel (Madison, WI: American Society of Agronomy and Crop Science Society of America), 89–111.

Monteith, J. (1977). Climate and the efficiency of crop production in britain. Proc. R. Soc. Lond B 281, 277–294. doi: 10.1098/rstb.1977.0140

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Nouvellon, Y., Begue, A., Moran, M., Seen, D., Rambal, S., Luquet, D., et al. (2000). PAR extinction in shortgrass ecosystems: effects of clumping, sky conditions and soil albedo. Agric. Forest Meteorol. 105, 21–41. doi: 10.1016/S0168-1923(00)00194-5

CrossRef Full Text

Okubo, A., and Levin, S. (2002). Diffusion and Ecological Problems, Modern Perspectives, 2nd Edn. Berlin; Heidelberg; New York, NY: Springer.

Oliver, C., and Larson, B. (1990). Forest stand dynamies. New York, NY: F. McGraw-Hill.

Page, E., and Gerwitz, A. (1974). Mathematical models, based on diffusion equations, to describe root systems of isolated plants, row crops, and swards. Plant Soil 41, 243–254. doi: 10.1007/BF00017252

CrossRef Full Text

Perttunen, J., Sievänen, R., and Nikinmaa, E. (1998). LIGNUM: a model combining the structure and the functioning of trees. Ecol. Modell. 108, 189–198. doi: 10.1016/S0304-3800(98)00028-3

CrossRef Full Text

Perttunen, J., Sievänen, R., Nikinmaa, E., Salminen, H., Saarenmaa, H., and Vakeva, J. (1996). LIGNUM: a tree model based on simple structural units. Ann. Bot. 77, 87–98. doi: 10.1006/anbo.1996.0011

CrossRef Full Text

Pouderoux, S., Deleuze, C., and Dhôte, J. (2001). Analyse du rendement des houppiers dans un essai d'éclaircie de hètre grace à un modèle à base écophysiologique. Ann. Forest. Sci. 58, 261–275. doi: 10.1051/forest:2001125

CrossRef Full Text

Pretzsch, H. (1992). Modellierung der kronenkonkurrenz von fichte und buche in rein- und mischbeständen. Allgemeine Forst- und Jagdzeitung 163, 203–213.

Purves, D., Lichstein, J., and Pacala, S. (2007). Crown plasticity and competition for canopy space: a new spatially implicit model parameterized for 250 north american tree species. PLoS ONE 2:e870. doi: 10.1371/journal.pone.0000870

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Raumonen, P., Kaasalainen, M., Akerblom, M., Kaasalainen, S., Kaartinen, H., Vastaranta, M., et al. (2013). Fast automatic precision tree models from terrestrial laser scanner data. Remote Sens. 5, 491–520. doi: 10.3390/rs5020491

CrossRef Full Text

Reddy, V., and Pachepsky, Y. (2001). Testing a convective-dispersive model of two-dimensional root growth and proliferation in a greenhouse experiment with maize plants. Ann. Bot. 87, 759–768. doi: 10.1006/anbo.2001.1409

CrossRef Full Text

Ryan, M., and Yoder, B. (1997). Hydraulic limits to tree height and tree growth. Bioscience 47, 235–242. doi: 10.2307/1313077

CrossRef Full Text

Schmidt, J., and Wulff, R. (1993). Light Spectral Quality, Phytochrome and Plant Competition. Trends Ecol. Evol. 8, 46–51.

Pubmed Abstract | Pubmed Full Text

Shinozaki, K., Yoda, K., Hozumi, K., and Kira, T. (1964). A quantitative analysis of plant form - the pipe model theory I. Basic Analysis. Japan. J. Ecol. 14, 97–105.

Sievänen, R., Perttunen, J., Nikinmaa, E., and Kaitaniemi, P. (2008). Toward extension of a single tree functional-structural model of scots pine to stand level: effect of the canopy of randomly distributed, identical trees on development of tree structure. Funct. Plant Biol. 35, 964–975. doi: 10.1071/FP08077

CrossRef Full Text

Sinoquet, H., Sonohat, G., Phattaralerphong, J., and Godin, C. (2005). Foliage randomness and light interception in 3-D digitized trees: an analysis from multiscale discretization of the canopy. Plant Cell Environ. 28, 1158–1170. doi: 10.1111/j.1365-3040.2005.01353.x

CrossRef Full Text

Sonntag, M. (1996). Effect of morphological plasticity on leaf area distribution, single tree, and forest stand dynamics. Bayreuther Forum Ökologie 52, 205–222.

Sorrensen-Cothern, K., Ford, E., and Sprugel, D. (1993). A model of competition incorporating plasticity through modular foliage and crown development. Ecol. Monogr. 63, 277–304. doi: 10.2307/2937102

CrossRef Full Text

Sterck, F., Schieving, F., Lemmens, A., and Pons, T. (2005). Performance of trees in forest canopies: explorations with a bottom-up functional-structural plant growth model. New Phytol. 166, 827–843. doi: 10.1111/j.1469-8137.2005.01342.x

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Takenaka, A. (1994). A simulation model of tree architecture development based on growth response to local light environment. J. Plant Res. 107, 321–330. doi: 10.1007/BF02344260

CrossRef Full Text

Toro, E. (2009). Riemann Solvers and Numerical Methods for Fluid Dynamics. Berlin; Heidelberg; New York, NY: Springer. doi: 10.1007/b79761

CrossRef Full Text

Tyree, M. (1988). A dynamic model for water flow in a single tree: evidence that models must account for hydraulic architecture. Tree Physiol. 4, 195–217. doi: 10.1093/treephys/4.3.195

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Tyree, M., and Zimmermann, M. (2002). Xylem Structure and the Ascent of Sap, 2nd Edn. Berlin; Heidelberg; New York, NY: Springer. doi: 10.1007/978-3-662-04931-0

CrossRef Full Text

Wang, W., Li, Z., and Su, H. (2007). Comparison of leaf angle distribution functions: effects on extinction coefficient and fraction of sunlit foliage. Agric. Forest Meteorol. 143, 106–122. doi: 10.1016/j.agrformet.2006.12.003

CrossRef Full Text

Yan, H., Kang, M., de Reffye, P., and Dingkuhn, M. (2004). A dynamic, architectural plant model simulating resource-dependent growth. Ann. Bot. 93, 591–602. doi: 10.1093/aob/mch078

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

4. Appendix

Numerical Implementation

A finite volume scheme Toro (2009) is used, in which we consider ${\overline{\alpha }}_{ijk}:=\frac{1}{{\Delta }^2}{\displaystyle {\int }_{{\Sigma }_{ijk}}\alpha }(x,t)dx$ (with a similar notation for the other variables) on a regular mesh with cells Σ_ijk = [iΔx, (i + 1) Δx[× [jΔx,(j + 1) Δx[×[kΔx,(k + 1)Δx] for i, j, k∈ ℤ and Δx « 1.

The local light incidence (8) (and similarly the local biomass production (4)) for vertical radiation, v = e₃, is computed as

$\begin{array}{l}{\overline{L}}_{ijk}(t)={\overline{\lambda }}_{ijk}(e_3)·\exp (-\Lambda ·\Delta x{\displaystyle \sum _{\kappa \ge k}{\overline{\lambda }}_{ijk}(e_3)·{\overline{\alpha }}_{ijk}(t)·}\hfill \\ \text{PAR}(e_3,t)dv)\hfill \end{array}$

The general case for v∈ S²₊ is conceptually similar, yet technically more extensive: For a given v, the above sum reaches over all cells Σ_i'j'k' whose intersection with the line through the center of the cell Σ_i'j'k' and pointing in the direction v is non-empty, taking into accout the individual length ℓ_i′j′k′(v) ≤ Δx of this intersection by means of the additional coefficient $\frac{{\ell }_{i^{\prime }{j}^{\prime }{k}^{\prime }}(v)}{\Delta x}$ in the sum term.

As for the PDE, let

${\overline{\varphi }}_{1,ijk}(t)=k·\frac{{\overline{\alpha }}_{i+1jk}(t)+{\overline{\alpha }}_{ijk}(t)}{2}·\frac{{\overline{L}}_{i+1jk}(t)-{\overline{L}}_{ijk}(t)}{\Delta x}$

and ϕ_2,ijk, ϕ_3,ijk be defined accordingly for the other spatial entries. The standard discretization of the divergence is given by

$\begin{array}{l}{\overline{\text{div}}}_{ijk}(t)=\frac{{\overline{\varphi }}_{1,ijk}(t)-{\overline{\varphi }}_{1,i-1jk}(t)}{\Delta x}+\frac{{\overline{\varphi }}_{2,ijk}(t))-{\overline{\varphi }}_{2,ij-1k}(t)}{\Delta x}\hfill \\ +\frac{{\overline{\varphi }}_{3,ijk}(t))-{\overline{\varphi }}_{3,ijk-1}(t))}{\Delta x}\hfill \end{array}$

and finally, using the Euler method for Δt « 1, we obtain

${\overline{\alpha }}_{ijk}(t+\Delta t)={\overline{\alpha }}_{ijk}(t)+\Delta t·{\overline{\text{div}}}_{ijk}(t).$