1. Introduction

Expanded knowledge of the carbohydrate pathway in trees is critical in agriculture, forestry and ecology. Improving on our understanding of translocation is particularly important in a climate change context so as to anticipate the effects of future environmental conditions on tree growth and carbon sequestration. Phloem transport is the process by which carbohydrates produced by photosynthesis in the leaves get distributed within a plant. Efficient transport ensures that the carbohydrate requirements of living tissues (respiration, growth) are met throughout the organism. In most trees, the phloem is a tissue layer located under the bark. Phloem is very thin (typically from 0.5 to 5 mm), i.e. at least 2 orders of magnitude less than its other dimensions. Because of that, it can be described as a three-dimensional manifold with a shape closely matching tree's external shape (minus the offset of bark's thickness). In a schematic view, sap flows in the phloem from leaves (sources) to roots (sinks). In reality, sinks are not only located at one end but also distributed all along the pathway. Even leaves can act as carbohydrate sinks when deciduous trees initiate new leaf growth at the beginning of the growing season. Phloem sap mostly consists in water carrying photosynthates by bulk. Inside the phloem, the sap moves through a network of elongated and interconnected cells. Those cells are referred to as sieve tube elements in angiosperms and as sieve cells in conifers. According to Münch [22], the osmotically generated hydrostatic phloem pressure is the force driving the long-distance transport of photoassimilates. Münch's hypothesis is generally accepted. For over 60 years, it has served as the basis for mathematical models of the phloem transport process [11, 2, 4, 7, 20, 29, 35, 15, 1].

We present the mathematical foundations and an implementation for a surface dynamic, anisotropic model of phloem transport. The purpose of this model is not to elaborate on the finer points of the physical process underlying phloem sap flow such as the contribution of diffusion [25], the radial leakage of solutes [1] or the presence of a relay mechanism [12]. Those aspects can be investigated individually. We adopt a standard mathematical model of coupled water-solute transport that is similar to Thomson and Holbrook's [35]. In this study, the model has been developed with large scale domains and arbitrary geometry in mind. In previous modelling attempts, phloem transport is treated as a one-dimensional process with the sap flowing though a file of sieve-tube elements connected end-to-end. Because plant stems are generally elongated with a high aspect ratio, they appear to be tubes or pipes. As a result, the 1D approach is also employed to simulate transport at the scale of the entire organism [13, 19]. However, a 1D model implicitly assumes that carbohydrates are in a common pool at any given position along stem's length. In trees, that assumption may not be valid.

Carbohydrates produced by a single branch are translocated along a downward, helicoidal pathway on the stem surface [9, 16]. That singular trajectory highlights two key points: i) carbohydrates predominantly follow the orientation of sieve elements with little lateral dispersion and ii) the direction of translocation does not correspond to the long axis of shoots and roots. In that context, most of the carbohydrates produced by a source are only available to sinks with a direct hydraulic connection to the source. Therefore, the difference in hydraulic properties along and across sieve elements as well as the lateral positioning of sources and sinks are essential to understand phloem transport in trees. Taking those features into account can be achieved by describing transport as a two-dimensional process [31]. The model we present is a surface model in that the flow of water and solute is neglected within the thickness of the phloem (i.e. intra-phloem radial transport). In mathematical terms, it does not present additional difficulty to extend the current model to the full 3D case. However, there are several reasons not to do so. Firstly, the macroscopic approach used here may not be applicable to a tissue less than a dozen cells wide [8]. Secondly, the fact that width is several orders of magnitude less than other conduit dimensions and that the computational domain must be defined for its smaller dimension would make the problem numerically untractable. Thirdly and more importantly, radial transport appears to follows a different cellular pathway, through ray parenchyma [33, 26], and Münch's hypothesis may not apply. On the other hand, using a surface model is not incompatible with simulating transport for three-dimensional surfaces and describing radial flow to adjacent tissues (as boundary conditions).

Other aspects relative to transport in trees have influenced the model's design. For instance, the model is based on the finite element method. Transport equations are integrated and solved numerically. The numerical approach provides the capability to carry out large-scale, long-time simulations. Any characteristic of the conduit (i.e. phloem's thickness or hydraulic conductivity) is defined on an element-by-element basis. In that manner, spatial heterogeneity can be easily represented. Finite elements also make it possible to solve transport equations for non-idealised geometries. Common geometrical irregularities such as nodal swelling, scars, burls, fluting and buttresses can be included in the model provided the biological shape can be characterized in the first place. Finally, carbohydrate unloading in the model can be represented as being time-and concentration-dependent, which allows combining the effects of sink dynamics on the patterns of carbon allocation in trees [21] with the effects associated to pathway's structure.

The paper is organised as follows. The first section is devoted to the analysis and qualitative properties of the model.The positivity conservation and the growth of the carbohydrate mass are proven. Although theoretical, this phase is necessary to ensure that the numerical approximated solutions preserve the properties of the biophysical process of transport. In the second section, we present numerical schemes to solve the highly nonlinear system of partial differential equations coupling carbohydrate transport and hydrostatic pressure in the phloem. Numerical simulations are presented in a the third section in order to evaluate the model and illustrate some of its capabilities. Those simulations include: a comparison and validation with an existing model [35] for the one-dimensional case; a parametric study; the application of the model to an existing tree; a preliminary investigation of the role played by sieve element orientation on carbohydrate distribution; simulating phloem transport on a branched, three-dimensional manifold.

---

2. Model description

In this section, we describe the model studied throughout the paper. The model is composed of a reaction-diffusion equation, coupled with a convection-reaction term. A schematic representation of a tree and of a phloem surface element are shown in figure 1.

---

2.1. Equations statement

The transport of a volume $\theta$ of water in a surface $\Omega$, a bounded domain of $\mathbb{R}^2$ of thickness $e$, is defined by the mass balance conservation

$\partial_t \theta + \nabla \cdot J_\theta + H_\theta = 0,$

where $J$, the water flux, is given by Darcy's law

$J_{\theta} = -\frac{e}{\mu} \left( k \nabla P \right).$

Here $P$ is phloem's pressure in Pa, $\mu$ is the sap viscosity in Pa\, s and depends on the concentration $C$ in solute, and $k$ in m$^2$ is the orthotropic permeability matrix

$k = \left( \begin{array}{cc} k_x & 0\\ 0 & k_y \end{array} \right).$

Typically, the $C$ value is 20% ($\approx$ 630 mol m$^{-3}$). At that C value, $\mu$ is twice that of pure water. $C$ values up to 1000 mol m$^{-3}$ and above are physiologically realistic. The dependency of $\mu$ on $C$ also theoretically affects flow efficiency [20]. It cannot be neglected beforehand.

$H_\theta$ denotes the radial water flux at the boundary between phloem and xylem. Here, the radial flux is function of the differential of water potential between phloem and xylem [35, 13]:

$H_\theta = -L_R(\psi - P + RT C) - V_s U,$

and $L_R$ the radial hydraulic conductivity in m Pa$^{-1}$ s$^{-1}$, $\psi$ the xylem hydrostatic pressure, $R$ the gas constant (in J mol$^{-1}$ K$^{-1}$), $T$ the temperature in K, $V_s$ the partial molal volume of sucrose in m$^3$ mol$^{-1}$ and $U$ the sucrose unloading (radial). It translates boundary conditions on entrance and exit as source and sink terms as follows

$U = \left\{ \begin{array}{rl} \widetilde{U} & {\rm{in}} \Omega_l {\rm{the loading area (source)}} \\ -\frac{\widetilde{U}}{C^*}C & {\rm{in}} \Omega_u {\rm{the unloading area (sink)}} \\ 0 & \mbox{ elsewhere}. \end{array} \right.$

Here $\widetilde{U}$ denotes a constant loading rate (mol m$^{-2}$ s$^{-1}$) and $C^*$ a reference sucrose concentration (mol m$^{-3}$). On the other hand, the variation of the volume $\theta$ depends on the pressure via the phloem thickness ($e$ in m) and phloem Young's modulus ($E$ in Pa) as

$\partial_t \theta = \frac{e}{E} \partial_t P.$

In other terms, the phloem is deformable and thickness depends on $P$. Concerning the amount of sucrose, the concentration $C$ is governed by

$e \partial_t C + \nabla \cdot J_C + H_C = 0,$

where $H_C = U$ and $J_C = C J_\theta.$ To sum up, the mass balance conservation is written, for $x\in \Omega$ and $t>0$

$\frac{e}{E} \partial_t P - \nabla \cdot \left( \frac{e}{\mu} \left( k \nabla P \right) \right) = L_R \left( \psi - P + RT C \right) + V_s U$

(1)

$e \partial_t C - \nabla \cdot \left( \frac{e}{\mu} C \left( k \nabla P \right) \right)= U,$

(2)

with initial data

$\label{eq5} P(0,x) = P_0(x) \mbox{ and } C(0,x) = C_0(x).$

(3)

For planar surfaces, we assume no flow of water or solute takes place on the vertical sides of the domain. We remind that in and out flows are imposed through the right hand side terms. We then impose Neumann boundary conditions:

$\label{eq4} \partial_n P(t,x) = \partial_n C(t,x) = 0.$

(4)

---

2.2. Theoretical qualitative study

In this section, the sap viscosity and the phloem thickness do not depend on the concentration of sucrose and the pressure, we assume in the following that $\mu$ and $e$ are also positive constants.

We first prove the conservation of the positivity.

Proposition 1. Let $(P_0,C_0)\in W^{2,\infty}(\Omega) \times W^{1,\infty}(\Omega)$ be positive initial data. Then, for all time $t\geq 0$ and a.e. $x\in \Omega$,

$P(t,x) \geq 0 \;and \; C(t,x) \geq 0.$

Proof. We have

$\partial_t P - \nabla \cdot \left( \frac{E}{\mu} ( k \nabla P) \right) = f(P,C),$

with

$f(P,C) := \frac{E L_R}{e}(\psi - P + RT C) + \frac{E V_s U}{e} \\ = (f(P,C) - f(0,C)) + f(0,C) = \frac{\partial f}{\partial P} (\theta P, C)P + f(0,C),$

where $0\leq \theta \leq 1.$ Consider $Q = \exp(\lambda t) P$, the equation becomes

$\partial_t Q - \nabla \cdot \left( \frac{E}{\mu} ( k \nabla Q) \right) - \left( \lambda + \frac{\partial f}{\partial P} (\theta P, C) \right) Q = \exp(-\lambda t) f(0,C).$

(5)

If there exists $(x^*,t^*)$ such that $Q(t^*,x^*) = \min_{t,x} Q(t,x) < 0$, then

$\partial_t Q(x^*,t^*) = 0, \; \nabla Q(x^*,t^*) = 0, \; \Delta Q(x^*,t^*) \geq 0,$

and $\lambda$ can be chosen such that the left-hand size of (5) and the right-hand size have opposite sign.

We deal similarly for $C$.

The phloem fills up and empties with respect to in and out flows. If $\nu ({\Omega _l}) = \nu (\Omega_u)$, then the total amount of sucrose $\int_\Omega C(t,x) dx$ is increasing as soon as $C < C^*$ in $\Omega_u$, whereas it becomes decreasing when $C > C^*$ in $\Omega_u$. More precisely, we have:

Proposition 2. For all time $t\geq 0$

$\frac{d}{dt} \int_\Omega C(t,x) dx = \frac{\widetilde{U}}{e}\left(\nu (\Omega_l) - \int_{\Omega_u} \frac{C(t,x)}{C^*}dx \right).$

Proof. Integrating (2) over $\Omega$ gives

$\frac{d}{dt} \int_\Omega C dx = \int_\Omega \nabla \cdot \left( \frac{1}{\mu} C \left( k \nabla P \right) \right) + \frac{U}{e} dx \\ = \int_\Omega \frac{U}{e} dx = \frac{\widetilde{U}}{e}\left(\nu (\Omega_l) - \int_{\Omega_u} \frac{C}{C^*}dx \right).$

Remark 1. Let $C^*>0$ be a given unloading strength. Then the steady state solution $C$ satisfies

$\int_{\Omega_u} C(x) dx = C^* \nu (\Omega_l).$

Remark 2. It seems difficult to prove the well-posedness of the system (1)-(2). A standard parabolic regularization $\varepsilon \Delta$ does not provide uniform estimates independent of $\varepsilon$ to evaluate the limit as $\varepsilon \to 0.$ We plan to study this issue in a future work.

---

3. Algorithm framework

We describe the numerical discretization used to approximate the pressure and the carbohydrate concentration given by the equations (1)-(2). To simulate large-time and large-space scales for realistic multi-dimensional phloem domain, the stability of the proposed scheme should not be too restrictive.

Let $n\in \mathbb{N}^*$. $P^n$, respectively $C^n$, denotes the approximation at time $t_n=n \Delta t$ of the pressure $P$, respectively the carbohydrate concentration $C$.

---

3.1. Splitting

To deal with the nonlinearity, the equation (2) is successively split [14] with a first order in time as, for $x\in \Omega$ and $t\in[t_n,t_{n+1}]$

$\partial_t \widetilde{C} - \left( \frac{k}{\mu} \nabla P \right) \cdot \nabla \widetilde{C} - \frac{U}{e} = 0 \mbox{ and } \partial_t C - C \, \nabla \cdot \left( \frac{k}{\mu} \nabla P \right) = 0,$

with $C(x,t_n) = \widetilde{C}(x,t_{n+1})$. Since $C(x,t)>0$ as soon as $C_0(x)>0$, the change of variables $\breve{C}=\log(C)$ is applied to the second equation to obtain

$\partial_t \breve{C} - \nabla \cdot \left( \frac{k}{\mu} \nabla P \right) = 0.$

Finally, the transformation $C = \exp(\breve{C})$ enforces the positivity.

---

3.2. Space and time discretization

Let $\varphi_1, \varphi_2, \varphi_3$ be test functions, with $\varphi_i, \nabla \varphi_i$ in $L^2$, for $1\leq i \leq 3$. Finite elements method in space is used, while the time discretization is obtained with a semi-implicit scheme.

This scheme is implemented using the software FreeFem₊₊[10]. The solution $P$ and $C$ are computed with $P_2$ and $P_1$ elements respectively whereas the pure convection terms are solved using the Characteristic-Galerkin method [27].

---

4. Numerical simulations

We present some numerical results obtained using Algorithm 1. The objective of the simulations is to illustrate the domain of applicability of the model. Geometries and parametrisation are case-specific; they are not properties of the model itself.

Algorithm 1 Semi-implicit scheme
Given $P_0, C_0$ and $N\in \mathbb{N}^*$.
For $n=1$ to $N$do
$\langle \frac{C^{n+1/2}-C^n}{\Delta t}, \varphi_1 \rangle - \langle \left( \frac{k}{\mu} \nabla P^n \right) \cdot \nabla C^n, \varphi_1 \rangle - \langle \frac{U^n}{e}, \varphi_1 \rangle = 0 \\ \breve{C}^{n+1/2}=\log(C^{n+1/2}) \\ \langle \frac{\breve{C}^{n+1}-\breve{C}^{n+1/2}}{\Delta t}, \varphi_2 \rangle + \langle \frac{k}{\mu} \nabla P^{n+1}, \nabla \varphi_2 \rangle + \int_{\partial \Omega} \frac{k}{\mu} \partial_n P^{n+1} \varphi_2 = 0 \\ \langle \frac{P^{n+1}-P^n}{\Delta t}, \varphi_3 \rangle + \langle \frac{E}{\mu} k \nabla P^{n+1}, \nabla \varphi_3 \rangle - \int_{\partial \Omega} \frac{E}{\mu} k \partial_n P^{n+1} \varphi_3 \\ - \langle \frac{E L_R}{e} (\psi - P^n + RT C^n, \varphi_3 \rangle - \langle \frac{V_s U^n}{e}, \varphi_3 \rangle = 0 \\ C^{n+1} = \exp(\breve{C}^{n+1}).$

---

4.1. Validation and comparison with an existing model

With no other 2D dynamic model of phloem transport available, we consider a case analogous to a 1D system so as to compare with the results of [35], a reference implementation of coupled water-solute transport. We simulate phloem transport in a 5 m-long, 10 cm-wide domain with a spatial step $\Delta x = 5$ mm. A 24-hour period is simulated with a time step $\Delta t = 1$ s. The simulation starts with initial pressure and carbohydrate concentration set to zero. Sap viscosity is a function of local solute concentration as a 15$^{th}$ order polynomial fitted on experimental viscosity for sucrose solutions at $T=293$K [34]. Parameter values are given in Table 1. They are chosen to be equivalent to those used in [35].

Figure 2 shows the predicted profiles of sucrose concentration (C), hydrostatic pressure (P) and axial water flux (J). All profiles are qualitatively very similar to those predicted in [35]. Peak positions, gradients magnitude and transitions over time are well reproduced. The profiles are also quantitatively comparable: the difference is within a few percent. It is only near steady-state, at $t=24$ h, that the variation becomes significant for $C$ (less than 10% underestimation) and $P$ (less than 20% underestimation). Although the behaviour at $t=24$ h is not particularly physical -simulated P is twice as high as the highest known measured value (2.4 MPa, [6]), loading and leaf dynamics are periodic, not constant [5] -both models should yield closer predictions. While the source of the discrepancy has not been identified, it could result from any of the following: differences in parametrisation due to domain geometry, numerical error, missing P derivative in eq. (2), implementation details, or the different $\mu = f(C)$ function.

The discrepancy originates from the present model predicting a flow that is slightly more efficient if compared to [35]. As a result less solute accumulates and less pressure builds up in the phloem. The difference is larger near steady-state because error accumulates at each time step. It is also more visible for C and P, which appeared very sensitive to small alterations of the flux.

In Figure 3, we show the numerical convergence of the scheme. Here different values of $\Delta t$ and $\Delta x$ are computed uniformly from 0.5 s to 10 s and from 1 mm to 100 mm respectively. We compute the maximum error by comparison with the approximated solution for $\Delta t = 0.1$ s, $\Delta x = 0.1$ mm after 12 hours. We also plot the lines $\Delta t$ and $\Delta x$ to demonstrate that the scheme is of order 1.

---

4.2. Parametric study

Several approaches are possible to simplify the governing equations (1) and (2). Sap viscosity can be treated as a constant instead of a function of sucrose concentration ($\mu = \mu_0$). This eliminates one term from eq. (2) and saves the computational cost of evaluating viscosity at all points of the grid. Alternatively, the effect of sucrose's partial molal volume can be neglected ($V_s = 0$). Also, phloem's thickness may also not be updated during the simulations ($e=e_0$). Each approach is considered individually as well as all at once and compared to results of the previous section.

Figure 4 shows how the axial water flux is affected by the proposed simplifications. Small changes in the J profile can have large effects on both C and P profiles (not shown here). With a constant sap viscosity, the flux is underestimated during the earlier stages and gets overestimated in the final stage. The position of peak flux also moves more slowly down the conduit with $\mu = \mu_0$. Ignoring the contribution of sucrose to sap volume ($V_s = 0$) causes the flux to be underestimated at any time. Initially, the effect is weak while the sucrose concentration is low but progressively increases in magnitude. Modelling the phloem thickness as being independent of pressure introduces only marginal differences. Although the longitudinal gradient of J is slightly higher between the loading zone and the front of the flow, the magnitude is comparable to that of the reference simulation. Because the simplification has only a minor influence on the J profile, it may seem advantageous to avoid computing phloem's deformation under flow. However, the simplification also has little consequences from a computational point of view; it does not modify the governing equations and all terms must still be evaluated. As expected [19], the flux predicted by combining all simplications is very close to the reference behaviour in near steady-state conditions. On the other hand, the flux in the early phases (0:10, 1:10) is strongly affected. Those simplifications are thus best avoided when attempting to describe phloem dynamics and rapid transitioning. Overall, none of the proposed simplifications appeared to be sufficiently beneficial to be included in the model.

---

4.3. Towards realistic designs

The objective of this application is to depict how the model can help in studying the behaviour of living trees. We simulate sap flow on the silhouette of an entire tree. In this application, a finite element mesh is created by Delaunay triangulation from a photograph of Te Matua Ngahere (see fig. 5). This kauri (Agathis australis) tree is the second largest in New Zealand. The height, diameter and volume of the tree trunk are equal to 10.21 m, 5.22 m and 208.1 m$^3$, respectively. There are obvious limitations to this approach. The surface is reconstructed from a photographic projection and flat. It means that the mesh is geometrically distorted compared to the actual 3D phloem layer and the results in this application get more inaccurate towards lateral edges. The main objective here is to show that the model can operate on a geometry that is not limited to rectangular channel but that has been retrieved from natural objects.

Figure 6 shows the tentative profiles for sucrose concentration and pressure in the phloem. Because of the aforementioned limitations, the profiles are not expected to be realistic. To improve on the quality of the results would require further steps such as reconstructing the surface of an existing tree stem in 3 dimensional space, determining the pattern of sieve cells' orientation on that surface, and identifying the strength of all carbon sources and sinks. Undertaking those steps is beyond the scope of this study.

---

4.4. Orthotropic transport

In this application, the potential for interaction between orthotropic (direction-dependent) transport and competing sinks is investigated. A simulation is carried out for a plate of dimensions $L=1$ m, $w=0.6$ m and $e=7.5$ $\mu$m. Three sources, denoted $r_i \; \forall i \in [1,3]$, are located near the top of the conduit. They are aligned diagonally to represent the spiral arrangement of branches (phyllotaxis). They have identical strength with a loading flux equal to $\widetilde{U}$. The botton region of the plate ($y < 0.1$ m) is divided into two sinks, $s_1$ on the left-hand side ($x < 0.3$ m) and $s_2$ on the right-hand side ($x>0.3$ m). The unloading at sinks is equal to $U_i = p_i \widetilde{U} \frac{C}{C^*}$ with $p_1 = 0.5$ and $p_2 = 1$. The lateral permeability is set to $k_T = 0.01 k_L$. The value is set arbitrarily to induce anisotropy in the system but not so high as to completely inhibit lateral transport. It represents the intermediate between isotropic transport ($k_T = k_L$) and the very high ratio ($k_T = 10^{-4} k_L$) suggested by [31] to recreate experimental conditions.

The distributions of sap pressure (P), sucrose concentration (C) and water flux (J) after 12 hours are shown in Figure 7. A vertical pressure ridge ($P> 0.95$ MPa) has formed on the left side of the conduit and a depression (0.7 MPa) has developed in the bottom right region corresponding to $s_2$. Like in the case of P, the C isocontours are slightly oblique, as would be expected from the difference of sink strength. Despite the fact that the three sources have the same strength, sucrose concentration at each source is not equal. Sucrose gets more concentrated as one moves leftward from one source to the next. The J pattern is opposite to those of P and C. The flow reaches higher values on the right side of the conduit because sources on that side are aligned with the stronger sink. In contrast, the sap flows less efficiently on the left side, which is aligned with the weaker sink. This causes a build-up in pressure and available sucrose. The role of sieve cell orientation on phloem transport is highlighted by strong longitudinal features in all profiles. Lateral transport is very limited. Nevertheless, the simulated distributions reveal a remarkable amount of interaction between conduit orientation and sink priorities.

---

4.5. Three-dimensional surfaces

There are limitations associated to representing the phloem of real trees as a planar surfaces. Tree's external surface must be figuratively cut and unrolled. It involves conformal mapping and setting periodic conditions at the lateral boundaries (for circumferential connectivity).

The alternative is to simulate transport directly on three-dimensional surfaces. Figure 8 shows a simulation carried out with the present model on a 3D domain. The mesh was created from an implicit surface using DELISO software [3]. Transport is simulated for a branched system, here a tree fork. Forks are common in deciduous, urban trees. Sap and carbohydrates flow from the primary branches and merge in the lower region, the main tree stem. Other geometrical features particular to trees (e.g. scars, buttresses, burls) can be included in simulations as long as a triangulated mesh can be generated to match the original geometrical configuration. Such a mesh was not available for this application.

---

5. Conclusion

The model presented here takes the mechanistic description of an osmotically-generated pressure flow and extends it to a phloem domain represented by a surface with anisotropic transport properties. The key features of the model were illustrated with applications. The model represents an important advance towards modelling the transport process for real, living trees. The transport equations are solved using a finite element method; each subdivision of the domain is assigned independent properties. The number of elements is only limited by the available computational ressource. With this approach, large-scale simulations become possible and the vast difference that exists between sieve cell dimensions and tree size can be bridged. The finite element approach is also particularly appropriate for describing a highly heterogenous biological material. As transport parameters are defined on an element-by-element basis, virtually any distribution of hydraulic properties and unloading rates can be represented in the model. Finally, the model can describe transport and source/sink dynamics with a fine time scale (under a second) for periods over several days. The effects on transport of heterogeneous distribution of hydraulic properties, periodic loading and interactions between sieve orientation and sink priority will all be investigated in future studies.

Major challenges have to be addressed before phloem transport can be computed for real tree configurations. The first challenge is to generate the external surface of an experimental tree, which could be done using 3D laser scanner or terrestrial LiDAR for trees of moderate size. Another challenge is to identify the distribution of diameter and orientation of sieve elements in a tree. The flow grain analogy [28] or models of auxin transport [17] could assist in approximating the pattern of fibre orientation. A third, very significant challenge is to evaluate the individual strength of all carbon sinks and sources within the tree over time. The emergence of experimental techniques relying on radiotracers and Positron Emission Tomography [23, 30] will be invaluable to monitor in vivo carbon dynamics and to inform models.

The presented model could be used to explore further additional aspects relative to the interaction between growth and transport. For example, it would be interesting to research optimality in the conflicting demands of carbohydrates for regulating transport (osmoregulation) and supplying living tissues (growth, respiration). From both a mathematical and biological point of view, the feedback loop between shape and transport is also of interest. On one hand, sap flow controls local growth by supplying the building material. On the other hand, tree stem geometry, derived from growth, impacts on where the sap flows. It has strong implications on the origin and the patterns of shape formation in the plant kingdom. Exploring the relation between tree shape and transport could be achieved by simulating phloem transport for virtual growing trees [31] generated using a Level Set method [32].

---

Acknowledgments

The numerical simulations were carried out on the High Parallel Computing platform MeCS with the support of the University of Picardie Jules Verne. This research was partially funded by the 'Growing Confidence in Future Forests' research programme (C04X1306). The authors would like to thank Jonathan Harrington, Dean Meason and the anonymous reviewers for their insightful comments on the first version of this manuscript. The authors express their gratitude to Tamal Dey for providing the DELISO software.

---