A dynamic model of annual foliage growth and carbon uptake in trees

The growth of trees and other plants occurs through the interactive combination of photosynthesis and carbon (and other nutrient) assimilation. Photosynthesis enables the production of carbohydrate that can then be used in growing foliage, whereby photosynthesis is enabled. We construct a mathematical model of carbon uptake and storage, which allows the prediction of the growth dynamics of trees. We find that the simplest model allows uncontrolled foliage production through the positive feedback outlined above, but that leaf shading provides an automatic saturation to carbon assimilation, and hence to foliage production. The model explains the necessity for finite leaf area production at outbreak, and it explains why foliage density reaches a constant value during a growing season, while also non-leaf tissue also continues to grow. It also explains why trees will die when their carbon stores are depleted below a certain threshold, because the cost of foliage growth and maintenance exceeds the dynamic supply of carbon by photosynthesis.


INTRODUCTION
In the current ongoing debate about climate change, the role of the carbon cycle is of major importance (Bigg 2003). The net flux of carbon (as CO 2 ) to the atmosphere is some 2 Gt (gigatonnes or 10 12 kg), but the contribution from fossil fuel consumption is more than twice this, and deforestation is almost as important as that. In fact, the fluxes between the various reservoirs are much larger, and so an understanding of the way in which carbon is transferred is of some importance.
The growth of plants and trees is an interesting and important process in the discussion of climate, because of their huge importance as a sink of CO 2 in the atmosphere: the net removal by vegetation is some 5 Gt a year. Therefore the quantitative and predictive understanding of carbon uptake by plants is an issue of much current interest. On the other hand, the consequent effect of climate change on the dynamics of plant growth is of equal importance for the understanding of crop dynamics.
The interaction of plant growth with carbon uptake is a problem in plant physiology, which has attracted much attention in the botanical and ecological literature, but less so in the applied mathematical community. An understanding of the processes of development and growth is important both scientifically and practically, since the response of plants and trees to nutrient deficiency and attacks by viruses and pests has important implications for crop and forestry management. In this paper, we develop and analyse a mathematical model for the growth of trees, which considers the continuous evolution of two critical variables, the carbon store and the foliage density. Our model finds its origins in the pioneering work of Thornley (1976; see also Thornley & Johnson 1990), and indeed our model closely resembles the whole plant model developed in Thornley's book (see his eqns (7.13) and (7.14)). However, we go beyond Thornley's model in two respects. The first is that we give a derivation of the equations of the model from first principles, describing the processes of water and nutrient transport within the tree; and second, we analyse the model and show that, in its simplest form, it exhibits an unphysical indefinite growth of foliage. We then go on to include the effects of leaf shading and show that the resultant augmented model allows for growth to a stable steady state, consistent with annual foliage development in mature trees. In a subsequent paper, we will use this model as a constituent in studying the defoliative effects on mountain birch forests in Lapland of outbreaks of the autumnal moth Epirrita autumnata, and with this in mind we will parametrize the present model with values suitable for the growth of a mountain birch.
The outline of the paper is the following. In §2 we describe the necessary tree physiology to develop the model, and in §3 we develop the model from first principles. We show that the model allows unbounded growth, and thus in §4 we seek a modification that prevents this happening. This is at hand in the form of leaf shading (Thornley 1976, ch. 3), which causes the photosynthetic rate to become saturable. We show in §4 that the inclusion of this effect causes the model to behave in ways consistent with observation. Conclusions follow in §5.

TREE PHYSIOLOGY
We begin by discussing the physiology of tree growth, with specific reference where appropriate to mountain birch. Early growth of foliage in trees is thought to depend on carbohydrates accumulated during the previous growing season. Carbohydrates are the direct product of photosynthesis, and are involved in many different processes in plants, such as respiration and maintenance. The leaves use carbon dioxide, water and energy from sunlight to produce carbohydrates, which are then translocated to different parts of the tree, which need them. A significant fraction of carbohydrate produced from photosynthesis is used for the growth of plant tissue (new wood in the case of trees). Another large fraction is stored for use during times when there is no production from photosynthesis, and a large amount of this stored carbohydrate is then used for initial leaf growth in early spring. There is also a considerable amount of carbohydrate that is oxidized in respiration, which releases the energy for tree growth and metabolic processes.

Vascular system
Nearly all living plants have a vascular system that is crucial for the transport of water, minerals, food and other substances to the different parts of the plant. This vascular system consists of two main pathways: the xylem and the phloem. The main difference between these is that the phloem is used to transport substances down the plant, whereas the xylem is used to transport substances (water and other minerals) up the plant. The phloem is considered to be the main pathway for carbohydrate transport, and generally very little carbohydrate is thought to travel in the xylem. It is thought that most of the reserve carbohydrates are stored in parenchyma cells in the xylem of the trunk for birch species. This was found by Olofinboba (1969) while comparing storage carbohydrates in African and European trees. These storage cells are thought to constitute 10 per cent of the xylem; the rest of the pathway is made up of dead cells, i.e. their only function is structural, and they do not store carbohydrate. Xylem is the main pathway for the transport of water and nutrients up the tree; the phloem is the main pathway for the products of photosynthesis down the tree. Transport down the phloem is driven largely by osmotic pressure differences; the flow of water up the xylem is driven by a pressure gradient that is induced by transpiration at the leaves, such that the water pressure in the xylem is below the saturation vapour pressure. This is known as the cohesion-tension theory ( Dixon & Joly 1895) and, although well accepted, it is still controversial ( Tyree 1997).
Xylem is made up of microscopic tubular vessels, which are distributed across the trunk in small clusters. Takahashi (1987) measured the diameters and distributions of xylem vessels in six birch species that grow in Siberia. For a birch species with characteristics similar to those of mountain birch, he measured an average vessel diameter of 47 mm, and a distribution of approximately 100 per mm 2 , with an average spacing of 23 mm. The vessels are never much longer than 4 mm, but the transported water can easily pass from one vessel to another by passing through small pits before entering another xylem vessel. For simplicity, we will suppose the vessels run continuously up the tree trunk.
Mineral nutrients and carbohydrates are transported in the vascular system as solutes in water. We suppose that stored carbohydrate in parenchyma cells is released to the water in the system in a diffusive fashion.

Photosynthesis
The production of carbohydrates occurs through photosynthesis in the leaves of the plant. Photosynthesis is the mechanism that reduces atmospheric CO 2 to carbohydrates using the energy from light and the release of oxygen from water. The basic chemical reaction can be described by and occurs at chloroplast cells within the leaves. These lie below the epidermis of the leaves, in which there are small openings called stomata. The stomata open and close, depending on climatic conditions such as temperature and wind strength. Atmospheric CO 2 is absorbed through the stomata, where it reacts as in (2.1) with water derived through the xylem from the roots. The excess water produced is then released by evaporation to the atmosphere through the stomata, this process being called transpiration, while the products of photosynthesis, the photosynthates, are transported away from the leaf through the phloem.
The carbohydrates thus produced are removed to the parts of the tree where they are needed, and to the storage cells. The amount of carbohydrate produced depends on how much water, CO 2 and light are available, and on how much foliage is present. Many experiments have been performed to measure how much CO 2 is absorbed in plants and then transformed into carbohydrate, and the results can vary considerably between species ( James 1973;Kozlowski & Pallardy 1997;Newell et al. 2002).

A MATHEMATICAL MODEL
There are three different types of carbohydrates (monosaccharides, oligosaccharides and polysaccharides) present in most plants (Barbaroux et al. 2003), but we will lump these together in our model. One of the variables of the model will thus be the carbon store, denoted C, which is measured in units of weight of carbohydrate, kilograms of carbohydrate (kgC). A second variable is the amount of foliage, denoted F, which is measured in units of weight of foliage, kilograms of foliage (kgF). The final variable is the woody biomass, denoted W, which is measured in units of weight of woody biomass, kilograms of woody biomass (kgW). More realistic models would include further compartments, but we avoid such complications here for purposes of simplicity.
We model the continuous variation of foliage, woody biomass and carbon by three ordinary differential equations, which take the form dC dt In these equations, P denotes the rate of production of carbohydrate by photosynthesis; S 1 denotes the rate of consumption of carbohydrate in foliage and woody tissue production; MR C denotes the rate of consumption of carbohydrate in supporting metabolic processes; S 2 is the rate of production of foliage; MR F represents foliage respiration; S 3 is the rate of production of woody biomass; and MR W is the rate of woody biomass respiration. The production of carbohydrate by photosynthesis is assumed to be proportional to leaf surface area, and we therefore take it to be proportional to F, thus P Z aF: ð3:2Þ Thornley (1976) suggested a similar relationship. The coefficient a will depend on received photosynthetically active radiation and water supply. A discussion of this (assuming wet conditions) is given by Thornley & Johnson (1990, p. 224 ff.). In the present notation, we write where s l is the leaf areal density, in units of kgF m K2 , and P l (Thornley & Johnson's P n ) is the specific photosynthetic rate, in units of kgCO 2 m K2 d K1 . The molecular weights M CH 2 O of organic carbon and M CO 2 of CO 2 provide the relevant conversion from carbon dioxide weight to carbohydrate weight. P l is described by a diffusive flux from the atmosphere, thus where r s and r b are stomatal and boundary-layer resistances, respectively, measured in units of d m K1 , and C a and C l are free air and leaf air concentrations of CO 2 , in units of kg m K3 . Thornley & Johnson suppose a specific leaf respiration rate of R d , which is the local production rate of CO 2 . The leaf CO 2 concentration is then determined by equating the net production of CO 2 to its removal by photosynthesis, described by a Michaelis-Menten-like pair of equations; this leads to where I is the leaf irradiance, and a p and K p are constants associated with the photosynthetic reactions. From these, we can eliminate C l (as the unique positive solution of a quadratic equation), and thus obtain an expression for a in terms of I, C a and r s Cr b . In the following, we assume it to be constant, as is commonly done (e.g. Magnani et al. 2000). The metabolic rates MR C , MR F and MR W are simply taken to be proportional to C, F and W, respectively, thus and then the specific leaf respiration rate R d is given by whereâ (kgC kgF K1 ) is the carbohydrate fraction of leaf tissue. It is clear that the rates in (3.6) must go to zero as C, F and W go to zero, respectively. Since MR F represents a decay rate, a linear dependence seems appropriate. Carbon use must saturate at high levels of carbon, since the tree has a finite capacity for carbon uptake. We ignore this here, partly because our later concern will be with trees under attack, where it is more likely that a linear dependence of metabolism on C is appropriate. The sink and source terms S 1 , S 2 and S 3 are more complicated. If S 1 is known, then we assume that where r and r 0 are conversion factors that determine foliage and woody biomass production in terms of carbohydrate consumed. These are simplifying assumptions, since they suppose a direct connection between foliage production and carbon release from storage. Thornley (1976) proposes a similar connection between the sink and source terms. It then remains to prescribe S 1 , and we discuss this in some detail. First, we find in appendix A that the water flux q is defined implicitly by the equation where m Z r a c p Lgðr s C r b Þr w s l ; ð3:10Þ and the various constants are defined in appendix A, as well as in table 1. This defines q implicitly in terms of F. Using values from table 1, we find p w z133 MPa, which compares with the typical measured leaf values of hydraulic potential of j l zK0.7 MPa. Because we expect jj l j/ p w , it follows from (A 5) that p l zp sat , which thus defines q through (A 4) and (A 3), and then (A 1) determines the value of j l , given j r . Specifically, we have note that this is independent of tree height h.
To calculate the carbon transport to the foliage, we denote by c w the concentration of carbon in the water in the xylem, measured as the weight of carbohydrate equivalent per unit volume of water (kgC m K3 ), so that S 1 Zqc w is the flux to the leaves. In appendix B, we show that The flux q(F ) is implicitly defined by (3.9) and, we anticipate, by (3.11).

Non-dimensionalization
To simplify the model, we scale the variables as We anticipate thatq/ q 0 C 0 , and an appropriate balance of terms then suggests that we choose   Tyree & Ewers (1991) where the dimensionless parameters are defined by The dimensionless version of (3.9) is then

Phase plane analysis
A first comment concerns the woody biomass equation for w. Because d 0 / 1, we see that w changes little over the course of a growing season and can be taken as approximately constant. This of course simply represents the fact that trees grow over a much longer time scale than a year. In any case, the equation for w uncouples from those for c and f, unless account is taken of the variation of stomatal conductance 1/r s , and thus m, and thus a and b (via C 0 ), on tree height (and hence w) (Ryan et al. 2000). This dependence would cause a and b to decrease somewhat with w, but, since w itself changes slowly, this coupling has little effect on a seasonal time scale. We thus formally put d 0 Z0 and consider only the first two equations in (3.19). We begin by ignoring the very small terms in 3, d r and d a . We then have Q zf ð1KH Þ ð 3:26Þ (which corresponds to (3.11)), and the equations ( where, for convenience, we now drop the asterisk on the dimensionless time t Ã . These two ordinary differential equations are easily studied in the (c,f ) phase plane. The solutions should represent the evolution of carbon store and foliage density during a growing season. We might expect that f would increase towards a steady state, and carbon store might do so as well or gradually increase. Carbon stores are depleted during winter when photosynthesis is inoperative.
There are two steady states, cZfZ0, and which is in the positive quadrant if, as we assume, aOb (and the relative humidity H!1). On cZ0, _ cZ af O 0, while the c-axis is invariant. Thus, if f and c are initially positive, they remain so. The nullclines are given in the positive quadrant by cZ 1=ð1KH Þ and fZ0 for _ f Z 0, and f Z dc=ðaKbð1KH ÞcÞ for _ cZ 0. These intersect at (c Ã , f Ã ), which is a saddle point. Trajectories lying below the two stable separatrices that approach (c Ã , f Ã ) tend towards the origin, which is a stable node, while those above asymptote towards the upper unstable separatrix from (c Ã , f Ã ), and both c and f grow unboundedly. The phase portrait is shown in figure 1.
The phase plane in figure 1 has a good point as well as a bad one. It suggests that if the tree is initially depleted of carbon ðc( 1Þ, then, if there is insufficient initial foliage, the carbon store will be depleted and the foliage will die: the tree will not survive. If cT 1, then the same fate can occur, but only if the initial foliage is very small. Note that the initial foliage must be positive for the tree to survive; this is consistent with the fact that, at budbreak, the tree does indeed produce an initial positive foliage density, depleting its carbon store to do so, as indicated by the trajectories in cT 1. However, the model then predicts unbounded growth in both carbon and foliage through an uncontrolled positive feedback between the two. It is this unwanted behaviour that we seek to redress in §4.

LEAF SHADING
One simple way in which we might suppose that foliage growth is curtailed is through the same mechanism that limits the growth of trees themselves. Although other suggestions have been made, Ryan & Yoder (1997) plausibly suggested that tree growth becomes limited in association with the height of the tree h. The mechanism of this hydraulic limitation hypothesis is that, as h increases, the tension in the water column in the xylem must increase (j l decreases in (A 1)) to drive the same water flow, making the xylem more prone to cavitation and thus damage. To prevent this, leaf stomatal conductance (1/r s ) is reduced, and this decreases photosynthetic rate and thus also the carbon available for woody (and perhaps foliage) growth. The hypothesis has been supported by Schäfer et al. (2000), Delzon et al. (2004) and Addington et al. (2006), for example, although not conclusively, and others do not agree (e.g. Becker et al. 2000). In the present model, hydraulic limitation does not seem a likely cause for reducing foliage growth, since our description is essentially appropriate for mature trees, and because there is actually no serious dynamic effect of reducing water flow in the model (3.19): reduced stomatal conductance 1/r s causes reduced m in (3.10), and thus increased C 0 in (3.18). From (3.20), this reduces a and b proportionally, so that the critical inequality aOb is unaltered. The phase plane of figure 1 is unaltered, and the problem of unbounded foliage growth remains. We may note that the choice of respiratory terms in (3.14) is actually consistent with the respiration hypothesis (Ryan & Yoder 1997), although we should also emphasize that the issue of tree height is not central to our thesis.
There are two covert assumptions in the model equation (3.1), which we might alter in seeking a resolution to the prediction of uncontrolled foliage growth. One was mentioned after (3.6); that is, the carbohydrate metabolization term M C must saturate at large C. Since in the dimensionless model this term is small, we do not expect that its modification at a large c will have any serious effect on the solution behaviour.
The other assumption is that the rate of photosynthesis is proportional to leaf area and thus foliage density. This assumption ignores the fact that sunlight reaching the surface of a leaf will depend on the degree to which it lies in shade, and leaves deeper within the canopy will receive less sunlight. This leads to the assumption that photosynthesis is a saturable function of foliage F, such as described by McMurtrie & Wolf (1983) and West (1987), who supposed that P in (3.1) was given by where s l is the leaf areal density introduced before; A c is the tree canopy area (assumed constant); and k l is a number, taken illustratively as 0.6 by McMurtrie & Wolf (1983). Introduction of (4.1) into the dimensionless model yields the corresponding dimensionless equations (in the approximated form analogous to (3.27)) where a retains its definition in (3.20) provided we define a Z k l P 0 s l ; ð4:3Þ and the parameter k is given by With the values in table 1, we find that aZ0.1 (as estimated independently) for the value s l Z0.08 kgF m K2 . Since, in fact, we hope that leaf shading will provide a limiting mechanism for carbon and foliage growth, it is natural to assume that kwO(1), and we can finally use this to choose F 0 so that kZ1, i.e.
With the values in table 1, we find F 0 Z0.4 kgF, as previously assumed, and therefore the values of a, b and d retain their previous estimated values.
Taking kZ1, the model (4.2) takes the form This has a fixed point as before where fwd, and the local structure of the phase plane is the same as in figure 1. In addition, there is a further fixed point where cZc 0 Z1/(1KH ), and fZf 0 ZO(1) is the larger root of the equation which has such an O(1) solution for aT b. The fact that a and b are large (and d/ 1) enables us to describe the dynamics of the solution when fZO(1). First, c rapidly approaches a quasi-steady state given by The right-hand side of this expression is a monotonically decreasing function of f, which equals (a/b)O1 when fZ0 and tends to zero as f/N. The migration of f along this curve is given by The right-hand side of (4.9) is a unimodal (onehumped) function with a unique positive steady state at fZf 0 (the positive solution of (4.7)) which is globally stable. From this we understand that the steady state (c 0 , f 0 ) is a stable node. The phase portrait for the system is shown in figure 2. The dynamics for small f are still described by figure 1, which is thus a close-up of figure 2 near the unsteady saddle (not readily discernible in figure 2).

CONCLUSIONS
In this paper, we have presented a relatively simple model, similar to others in the literature, in which we describe the dynamics of carbon assimilation, storage and use in mature trees, with a view to describing how the carbon storage evolves over the growing season. We have focused our attention on the mountain birch, because in a succeeding paper we will describe how this dynamic evolution is affected by the presence of the larvae of the moth E. autumnata. Our model presents two advances over previous work. First, the dynamics of carbon uptake and transport within the xylem is described physically, so that the nonlinearity of the model is based on coherent physical principles. Second, we study the global dynamic behaviour of the model in a way that has not been done before. In particular, we find that the model provides a descriptive reason for why at budbreak leaves form with a small but finite capacity for photosynthesis. In addition, we show that leaf shading provides a practical mechanism by which foliage growth is limited and the carbon storage reaches a steady state.
Both of these features provide insight into the mechanisms of carbon regulation during tree growth. The stable steady state shown in figure 2 at cZ5, fz0.56 corresponds to a carbon store of 0.65 kgC and a foliage quantity of 0.22 kgF. During the winter the carbon store is depleted, so that in spring the initial values (in dimensionless units) will be c!5 and f small. For example, the separatrix in figure 2 from cZ0, fZ0.257 goes approximately through the values cZ3, fZ0.1. So if cZ3 at the end of winter (corresponding to 0.39 kgC), then for an initial foliage less than fZ0.1 (corresponding to 0.04 kgF), carbon will first increase, but then irreversibly decrease, and tree death will ensue.
In reality, we would expect the carbon store to be less depleted. Figures 3 and 4 show the time series of foliage and carbon quantities as a function of time (in dimensional units) starting from dimensionless values cZ4.8, fZ0.2, corresponding to CZ0.62 kgC and FZ0.08 kgF. There is an initial slight dip in foliage quantity before it rises monotonically towards its final steady state. The carbon store initially rises as photosynthesis is initiated and then is depleted by the growth of foliage. Depending on the parameters that are chosen, the final approach to the steady state can be somewhat oscillatory, which is probably unphysical. Figures 3 and 4 do highlight a shortcoming of this model. Although we have attempted to use the best estimates available for the parameters, we see that the time scale for approach to equilibrium is much too long, of the order of 1000 days. This is too long by a factor of perhaps 30. The time scale is given by t 0 Z1/m F , where m F is the specific respiration rate. In this simple model, somewhat analogous to simple first-order rate laws used in epidemiology, we are actually assuming an exponential distribution of foliage lifetime with a mean of 1/m F Z67 days; evidently this is much too long, and the model ought to be refined, perhaps by means of specific inclusion of foliage lifetime as well as foliage respiration.
A simpler alternative is just to shorten the time scale by increasing m F . There is then a consequent effect on a, b and d. If one does increase m F , then indeed the time scale is shortened, but the dynamics becomes increasingly oscillatory, which also seems implausible.  A dynamic model of carbon uptake in trees A. C. Fowler et al. 1093 We think that this shortcoming represents an issue concerning the veracity of the model which needs to be addressed. For example, specific respiration rates decrease rapidly with foliage age ( Kozlowski & Pallardy 1997), and this would have a compensating effect. However, we also think that the robust qualitative conclusions of the present model provide a useful basis on which more elaborate and realistic models can be built.

APPENDIX A
In this appendix, we derive an expression for the water flux q through the xylem. This is given by where h is the tree height; j r and j l are the hydraulic potentials in roots and leaves, respectively; and the hydraulic conductivity k is supposed to be given by a formula based on the assumption of Poiseuille flow through n xylem vessels of diameter d, where h is the viscosity of water ( Tyree & Ewers 1991). The flux of fluid to the leaves is essentially balanced by evaporation from the foliage, and we suppose that this is proportional to the quantity of foliage, thus where E is the evaporation rate measured in kg m K2 d K1 (1 dZ1 dayZ8.64!10 4 s), and r w is the density of water. Tan et al. (1978) give an expression for the evaporation rate which can be written in the form E Z r a c p ðp l K p a Þ Lgðr s C r b Þ ; ðA 4Þ where p l and p a are the water vapour pressures at the leaf and in the atmosphere, respectively; r a is the density of moist air; c p is the specific heat of moist air; L is the latent heat of vaporization of water; and g is the so-called psychrometric constant. The stomatal conductance 1/r s is a decreasing function of the vapour pressure deficit p l Kp a (Oren et al. 1999), so that the evaporation rate E in (A 4) grows more slowly than linearly and may saturate at high vapour pressure deficit. The derivation of equation (A 4) is given in appendix C. In order to proceed, we must relate the hydraulic potential at the leaf, j l , to the water vapour pressure p l . This is done in appendix D, where we find the relationship j l Z p w ln p l p sat where p sat is the saturation vapour pressure; r w is the density of water; R is the universal gas constant; T is the absolute temperature; and M v is the molecular weight of water. We eliminate j l and p l from (A 1) and (A 3) using (A 4) and (A 5) to obtain where m Z r a c p Lgðr s C r b Þr w s l : ðA 7Þ

APPENDIX B
We assume that the carbon transport to the foliage is equal to the flux from the carbon store. If V s is the volume of the carbon store, then C/V s is the effective concentration of stored carbon. We suppose that stored carbon enters the xylem diffusively, so that the relationship equating the flux to the foliage with that to the xylem (kgC s K1 ) is where D is a diffusion coefficient and w is the width of xylem vessel walls. It follows from this that the xylem concentration of carbon is where q Ã Z nphdD w : ðB 3Þ Finally, we allow for the fact that the number of xylem vessels is related to the carbon store, most simply because the carbon is stored in parenchyma cells in the xylem. Put another way, if r c is the density of structural carbohydrate in the wood, then the total mass of carbon in the trunk is approximately r c hA, where A is the trunk cross-sectional area. If the carbon reserves available for growth are a fraction s of this, then C Z sr c hA; ðB 4Þ while an estimate for the number of xylem vessels is where l is the xylem spacing. This leads to n Z C sr c hl 2 : ðB 6Þ From (A 2) and (B 3), we have