We propose a hydrodynamic model of phloem transport, according to the Münch pressure flow hypothesis, with allowance for water and solute exchange along the pathway corresponding to the transport phloem (Van Bel, 2003) (Figure 1). The following assumptions are made:

The sieve tube is considered as a right circular cylinder of length L and radius R, such that R « L, limited by a porous membrane through which water and solutes fluxes occur. We use cylindrical coordinates where the axis along the pathway is denoted by z = [0, L], the radial coordinate is r = [0, R] and φ = [0, 2π ] is the azimuth. The tube is surrounded by unspecified tissue (for simplicity called apoplast in the following).

where α is the fraction of sieve plate area composed of sieve pore area; l is the length of sieve tube elements in between sieve plates; r_p and l_p define the sieve plate pore dimensions: radius and length respectively (Thompson and Holbrook, 2003a). As α = 1, l > l_p and R > r_p we have always β < 1. For the special case of no sieve plates the impedance factor is set to β = 1.

As in Phillips and Dungan (1993) the calculation starts from the steady state Navier-Stokes equation for an incompressible Newtonian fluid which, using the assumption of axisymmetrical flow, can be written in cylindrical coordinates as (Truskey et al., 2009) (2)ρ(ur∂ur∂r+uz∂ur∂z)=−∂p∂r+μ∗[∂∂r(1r∂(rur)∂r)+∂2ur∂z2 ]  (3)ρ(ur∂uz∂r+uz∂uz∂z)=−∂p∂z+μ*[1r∂∂r(r∂uz∂r)+∂2uz∂z2 ] where pressure p and velocities u_r and u_z are functions of axial coordinate z and radial coordinate r. In most cases we omit these dependences in our notation for simplicity. The effect of sieve plates is described introducing an effective viscosity μ^* = μ/β, thus having a more viscous fluid as β = 1 [see Equation (1)]. This approach treats the sieve tube as a conduit of essentially uniform resistivity, a close approximation to reality in view of the distribution of the sieve plates within sieve tubes (Weir, 1981; Henton et al., 2002). Pressure p includes the gravitational effect, i.e., p = p′ + ρgh; where p' is termed the hydrodynamic pressure inside the sieve tube, g the local acceleration of gravity and h is the vertical coordinate above a standard reference plane. Note that in contrast to Phillips and Dungan (1993) we use the Navier-Stokes equation including inertial forces [left hand side of Equations (2) and (3), respectively] to keep the formulation as general as possible. The actual impact of inertial forces is discussed in the Results section. The continuity equation for incompressible fluids in cylindrical coordinates is (Truskey et al., 2009) (4)1r∂(rur)∂r+∂uz∂z=0 and the convection-diffusion equation for dilute solutions in steady state is (Truskey et al., 2009) (5)ur∂c∂r+uz∂c∂z=D[1r∂∂r(r∂c∂r)+∂2c∂z2 ] where concentration c is a function of coordinates z and r in the same way as pressure and velocities. The axial and radial fluxes of solute inside the sieve tube are defined as sum of convective and diffusive components (6)jzs=uzc−D∂c∂z (7)jrs=urc−D∂c∂r

To allow for solute and water exchanges across the sieve tube membrane, the following boundary conditions are used: radial flux of solution through the membrane (at r = R) is given by Starling's equation (Nobel, 2009): (8)ur(R,z)=LP{p(R,z)−pout(z)−σRgT[c(R,z)−cout(z) ] } where the van't Hoff expression of osmotic pressure for dilute solution has been inserted. The reflection coefficient of the sieve tube membrane σ assumes a value between 0 (totally permeable) and 1 (semipermeable). Pressure, p_out, and solute concentration c_out, outside the tube do not to depend on r and are assumed to vary linearly with distance: (9)pout(z)=pout(0)+dpoutdzz and (10)cout(z)=cout(0)+dcoutdzz

Radial passive flux of solutes through the membrane on the one hand is given by Equation (7) for r = R, on the other hand it can be defined as the sum of convective flux (solute dragged by solvent) which is linear with the average solute concentration in the membrane, c(z), plus diffusive flux across the membrane that is described by the sieve tube membrane solute permeability, P_s, and the concentration difference across the membrane (Kedem and Katchalsky, 1958; Benedek and Villars, 2000).

(11)jrs(R,z)=ur(R,z)c(R,z)−D∂c∂r(R,z)            =(1−σ)ur(R,z)c¯(z)+Ps[c(R,z)−cout(z) ]

The concentration of solute, c(z), within the sieve tube membrane is given by c¯(z)=c(R,z)+cout(z)2. The no-slip condition implies that (12)uz(R,z)=0 which means that the fluid velocity will be in the radial direction only, at the boundary.

Symmetry at the center of the tube implies that radial velocity as well as radial derivatives of axial velocity and concentration are zero at r = 0: u_r (0, z) = 0, ∂uz∂r(0,z)=0 and ∂c∂r(0,z)=0. According to our model assumptions, flow is already fully developed at the beginning of the transport phloem region considered in our model; loading of solutes and water has occurred in the source region, not being part of our model upstream. Therefore, the boundary conditions at the inlet of the tube z = 0 are 〈u_z(0)〉 = U, 〈c(0)〉 = c_i and 〈p(0)〉 = p_i where the bracket denotes the average over the cross section 〈c(z) 〉=2R2∫0Rc(r,z)rdr.

Note that the hydrodynamic equations (2) to (4) describing the solution movement inside the sieve tube do not depend on the solute concentration explicitly, and we take viscosity to be independent of concentration. The flow depends on concentration through the boundary condition Equation (8) only. Thus, we can treat independently the hydrodynamics problem that describes how pressure and velocity change, and the solute transport equation (5) which is solved once the velocity field is determined.

As with many problems in fluid dynamics and especially when applying the Navier-Stokes equation, the problem becomes more simple to solve using dimensional analysis (Regirer, 1960; Kundu and Cohen, 2008). With this method, it is possible to predict physical parameters that influence the sap flow and determine the relationships between several variables (pressure, velocity and concentration) when an exact functional relationship is unknown. This is not possible with a direct numerical solution of the governing equations (2) to (5). Following Regirer (1960), Phillips and Dungan (1993) and Thompson and Holbrook (2003b), we use system geometry and boundary conditions as scales to define dimensionless variables in the following way: (13)r^=rR,z^=zL,u^z=uzU  and  c^=cci

Inserting these definitions into the governing equations, the corresponding scaling of the remaining variables and parameters [unique except for numerical prefactors which were chosen according to Thompson and Holbrook (2003a)] becomes p^=pi^R28LUμ∗,u^r=urLRU,j^zs=jzsUci,j^rs=jrsLRUci,L^p=Lp16L2μ∗R3,P^s=Ps2LRU

Substituting the dimensionless variables into the governing equations, the dimensionless Navier-Stokes Equations (2) and (3) become: (14)Re ε3(u^r∂u^r∂r^+u^z∂u^r∂z^)=−8∂p^∂r^+ε2[∂∂r^(1r^∂(r^u^r)∂r^)+ε2∂2u^r∂z^2 ] (15)Reε(u^r∂u^z∂r^+u^z∂u^z∂z^)=−8∂p^∂z^+1r^∂∂r^(r^∂u^z∂r^)+ε2∂2u^z∂z^2

Here ε = R/L is an abbreviation for the aspect ratio and Re = ρRU/μ^* is the Reynolds number, the ratio of inertial to viscous forces. The continuity equation (4) now reads (16)1r^∂(r^u^r)∂r^+∂u^z∂z^=0

The convection diffusion equation (5) becomes (17)ε(u^r∂c^∂r^+u^z∂c^∂z^)=1Per[1r^∂∂r^(r^∂c^∂r^)+ε2∂2c^∂z^2 ] where Pe_r = RU/D is the Péclet number for radial flow. The components of the sieve tube solute flux are given by: (18)j^rs=−1Per∂c^∂r^+εu^rc^ (19)j^zs=−1Pez∂c^∂z^+u^zc^ where Pe_z = LU/D = Pe_r/ε is the Péclet number for axial flow. The Péclet number for any direction gives the ratio of the rate of convection of the solute in that direction by the flow of sap, to the corresponding rate of diffusion of that solute driven by a concentration gradient.

In dimensionless notation, the boundary conditions Equations (8), (11) and (12) at the sieve tube membrane are (20)u^r(1,z^)=L^P2{p^(1,z^)−p^out(z^)−σH^[c^(1,z^)−c^out(z^) ] } (21)1Per∂c^∂r^(1,z^)=εu^r(1,z^)[c^(1,z^)−(1−σ)c¯^(z^) ]−εP^s2 [c^(1,z^)                       −c^out(z^) ] (22)u^z(1,z^)=0

The dimensionless parameter H^=RgTciR28LUμ∗, defined by Π^=H^c^, gives the ratio of osmotic forces to viscous forces for a flow with a characteristic velocity U (Phillips and Dungan, 1993).

The boundary conditions at the center of the tube become û_r (0, ẑ) = 0, ∂u^z∂r^(0,z^)=0 and ∂c^∂r^(0,z^)=0. Boundary conditions at ẑ = 0 are 〈û_z〉 = 1, 〈ĉ〉 = 1 and 〈p^〉 = p^_i. Thus after non-dimensionalization the parameters governing the model behavior are L^p, Ĥ, P^_s, σ, ε, Re, Pe_r together with boundary condition values p^_i, p^_out (ẑ) and ĉ_out (ẑ).

Proceeding in the same way as Phillips and Dungan (1993) we expand the state variables, û_r, û_z, ĉ and p^ as power series of the small dimensionless aspect ratio ε:u^r=∑j=0∞εju^rj,u^z=∑j=0∞εju^zj,c^=∑j=0∞εjc^j and p^=∑j=0∞εjp^j. The accuracy of the expansions, i.e., the number of terms to include, will depend on the value of ε: the smaller ε, the more significant are the first terms compared to higher order terms. In the case of phloem transport, ε is smaller than 10⁻³ for typical sieve tube dimensions (Esau, 1969). For this reason we will only consider the first two terms (zeroth and first order) of the expansion to describe flow in a sieve tube. Higher order terms might be of interest when studying effects on a smaller length scale (e.g., single sieve tube elements).

Collecting terms of the appropriate order (ε⁰ and up to ε¹, respectively), using respective boundary conditions and inserting zeroth order results into the calculation of first order leads to the following expressions: axial velocity can be calculated from Equations (15) and (22) as (23)u^z0(r^,z^)​=​2∂p^0∂z^(r^2−1) (24)u^z1(r^,z^)​=​2∂p^1∂z^(r^2​−​1)+Re36∂2p^0∂z^2∂p^0∂z^(2r^6​−​9r^4​+​36r^2​−​29)  where pressures p^₀ and p^₁ are independent of radial coordinate r as a result of Equation (14). Using Equations (23) and (24), radial velocity follows from Equation (16) as (25)u^r0(r^,z^)=−12∂2p^0∂z^2(r^3−2r^) (26)u^r1(r^,z^)=−12∂2p^1∂z^2(r^3−2r^)−Re144∂∂z^(∂2p^0∂z^2∂p^0∂z^)                   (r^7−6r^5+36r^3−58r^)

Inserting Equation (20) into Equations (25) and (26) the relation of pressure and concentration becomes (27)∂2p^0∂z^2=L^P{p^0(z^)−p^out(z^)−σH^[c^0(1,z^)−c^out(z^) ] }        =:u^lat 0(z^) (28)∂2p^1∂z^2​=​L^P{p^1(z^)​−​σH^c^1(1,z^) }​−​3Re8∂∂z^(∂2p^0∂z^2∂p^0∂z^) 

The average of zeroth order axial velocity over the cross section is 〈u^z0 〉=−∂p^0∂z^. Thus in our model the pressure gradient at any distance from the origin is linear with the average axial velocity, but dependent both explicitly and implicitly on the radial flux of water at that location, given by the zeroth order function û_{lat 0} (ẑ), see Equation (27).

Evaluating the convection diffusion equation (17) for ε = 0 together with the boundary condition, we conclude ĉ₀ (ẑ) to be independent of radial coordinate: ∂c^0∂r^=0. Using this result and Equation (23) the first order expression of Equation (17) becomes (29)∂c^1∂r^=Per2∂p^0∂z^∂c^0∂z^(r^3−2r^)

The first order expression of the boundary condition Equation (21) is (30)1Per∂c^1∂r^(1,z^)=u^r0(1,z^)[c^0(1,z^)−(1−σ)c¯^0(z^) ]                        − P^s2[c^0(1,z^)−c^out(z^) ]

Combining Equations (29) and (30) and inserting Equations (25) and (27) leads to the following zeroth order equation for concentration (31)∂p^0∂z^∂c^0∂z^=−u^lat 0(c^0−(1−σ)c¯^0)+P^s(c^0−c^out)

The final first order expression for concentration follows from Equation (29) using the boundary conditions 〈ĉ₁ (z = 0)〉 = 0 and 〈û_z0 (z = 0)〉 = −∂p^₀(z = 0)/∂ẑ = 1 (32)c^1(r^,z^)=Per8∂p^0∂z^∂c^0∂z^(r^4−4r^2)−5Per24∂c^0(z=0)∂z^

As the velocity profile depends on pressure only [Equations (23) to (26)] and first order concentration follows from zeroth order expressions [Equation (32)] we need to know ĉ₀, p^₀, and p^₁ only in order to obtain the pressure and velocity profiles and thus describe phloem flow up to first order of the expansion in aspect ratio ε. These variables are obtained by solving the system of coupled differential equations (27), (28) and (31).

For a semipermeable membrane, i.e., σ = 1 and P^_s = 0, the zeroth order equations (23), (25), (27) and (31) are identical to findings of Phillips and Dungan (1993) except for different numerical prefactors stemming from slightly different definitions of the dimensionless quantities. Also in this case Equations (27) and (31) are equivalent to the governing equations of Thompson and Holbrook (2003b), the only difference being in the boundary conditions. In this light it is not surprising that the zeroth order results Equations (27) and (31) for a permeable membrane can also be obtained using an approach similar to Thompson and Holbrook (2003a), based on conservation equations and a local application of the Hagen-Poiseuille equation together with appropriate boundary conditions, see Appendix B.

As there is no analytical solution for the system of differential Equations (27), (28) and (29), they were transformed into a first order system of differential equations [Appendix C, Equations (C.3)–(C.7)] and solved with MATLAB (R2010b, MathWorks, Inc.) ode15s differential equations solver routine. The calculated zeroth order of concentration, ĉ₀, turgor pressure, p^₀, the first order of turgor pressure, p^₁, and their respective derivatives were then used to determine the average first order concentration 〈ĉ₁〉, Equation (C.8). Concentration profile and sieve tube turgor pressure were determined as sum of zeroth and first order. From these profiles, average over cross section of axial velocity Equation (C.9), radial flow at sieve tube membrane Equation (C.10), axial solute flux Equation (19) and radial solute flux Equation (18) were determined.

Unless otherwise specified, the parameter values of Table 1 were used (the basis for these values is given in Appendix A). The values of the non-dimensional model parameters related to the parameter values of Table 1 are shown in Table 2.

Figure 2 shows profiles of pressure, concentration and flow for different values of the reflection coefficient σ. In the case of a semipermeable membrane (σ = 1), concentration and pressure differences across the sieve tube membrane (c and p are higher than c_out and p_out) create a water potential difference, ΔΨ, across the sieve tube membrane that always draws water into the sieve tube (û_r < 0, Figure 2D), causing a pronounced non-linearity of pressure in the direction of flow (Figure 2A). In the case of a membrane that is permeable to solutes as well as water (σ < 1), the water influx depends on the value of the reflection coefficient σ and can even be zero (Figure 2D), leading to a pressure profile almost identical to Poiseuille flow (Figure 2A). Boundary condition Equation (20) shows that a water influx (û_r < 0) will occur as long as: (33)σ>p^(1,z)−p^out(z)H^[c^(1,z)−c^out(z) ]

For typical physiological conditions (Table 1), i.e., for the order of magnitude of pressure and concentration expected for sieve tubes, water influx into the sieve tube occurs only if σ > 0.7 (Figure 2D). If the permeability of the membrane, P_s, increases (with σ decreasing from σ = 1 to σ = 0.7; see Appendix A for the relationship of P_s and σ), the absolute value of the pressure gradient in the direction of flow becomes smaller with decreasing axial velocity (Figure 2C) as 〈u^z 〉≈−∂p^∂z^, see Equation (C.9). This trend arises from two factors that decrease the sieve tube solute concentration along the axis: first, the dilution created by water influx, dependent on the water potential difference across the sieve tube membrane; second, the passive efflux of solutes across the sieve tube membrane, dependent on the concentration difference between the sieve tube and the apoplast. The passive loss of solutes is favored by a higher sieve tube solute concentration compared with the surrounding apoplast. The decrease in concentration, as one moves further down the tube, means less water will enter because the water potential difference across the sieve tube membrane decreases, and more so for a more permeable membrane (Figure 2D). Consequently, due to the volume conservation, the increase in axial velocity is less for more permeable membranes (Figure 2C) and there is also a decrease in axial solute flux (Figure 2E). Due to the concentration difference across the sieve tube membrane there is solute efflux which increases for smaller values of the reflection coefficient, σ, corresponding to a leakier membrane (Figure 2F). However, the sieve tube membrane solute permeability, P_s, is in the order of 10⁻¹⁰ m s⁻¹ (see Appendix A), and too small to cause dramatic changes in the sieve tube concentration. If the permeability of the sieve tube membrane, P_s, increases (with a corresponding decrease in the reflection coefficient, σ) there is a smaller osmotic effect on water exchange, leading to less water entering the sieve tube (Figure 2D) and a smaller decrease of solute concentration (Figure 2B). Eventually, for a very permeable membrane there is a reversal and water moves out radially for σ < 0.7 (Figure 2D), and the axial velocity decreases along the axis (Figure 2C). Simultaneously, the pressure inside the sieve tube decreases less and less for smaller values of σ (Figure 2A). Axial velocity and axial solute flux at one point become null (Figures 2C,E). Numerical calculation has to stop at this point since, without specifying boundary conditions at the outlet of the tube, it cannot be decided whether there is inflow from the opposite side or the axial velocity remains zero along the tube.

The expected behavior of solute concentration 〈ĉ〉 for very permeable membranes (small values of σ) would be that 〈ĉ〉 decreases less with distance, approaching a value constant over distance for σ = 0, but not increasing over the initial value at z = 0. Obviously our numerical results of solute concentration (Figure 2B) deviate from this expected behavior. The steep increase in 〈ĉ〉 (Figure 2B) for very small values of σ stems from the fact that the pressure gradient on the left hand side of Equation (31) approaches zero much faster than radial solute flow on the right hand side. The extreme of this scenario is seen for a totally permeable membrane (σ = 0), in which case the radial water exchange is driven solely by the pressure difference across the sieve tube membrane, and the solute flux occurs predominantly through convection by the water moving across the sieve tube membrane. In this case Equation (31) becomes ∂p^0∂z^∂c^0∂z^=(−u^lat 0/2+P^s)(c^0−c^out), concentration inside the sieve tube c₀(z) should stay constant, and apoplast solute concentration c_out(z) would be expected to tend to c₀(z) as a consequence of the high permeability of the membrane. However, our model assumes c_out(z) to be independent of radial solute flux, with constant apoplastic concentration gradient over distance [model assumptions 8 and 9, see also Equation (10)], leading to a higher concentration difference between inside and outside than expected, and, as a mathematical consequence of Equation (31), an accumulation of solute inside the tube with distance. Thus the results for solute concentration (Figure 2B) and radial solute flow (Figure 2F), though consistent with the model assumptions, are not realistic for small values of σ.

In order to justify the truncation of the series expansion we compare zeroth and first order results. For velocity and pressure Equations (23) to (28) the first two orders structurally only differ by expressions proportional to the Reynolds number Re which is quite small in our case (Re = 8.94 × 10⁻⁵ for the values reported in Table 1). Thus we expect terms proportional to Re not to contribute significantly for the set of parameters given (Table 1). The remaining terms of first order in Equations (24), (26) and (28) should be of the same order of magnitude as zeroth order, depending on boundary conditions. Together with the fact that the aspect ratio ε is small (ε = 2 × 10⁻⁵ for the values of Table 1) we conclude that zeroth order strongly dominates the series expansions of velocity and pressure profiles.

The comparison is less obvious for the concentration profile where first order Equation (32) is proportional to the Péclet number Pe_r. Since the Péclet numbers in our case are Pe_r = 3.62 and Pe_z = 1.81 × 10⁵, it is obvious from Equation (19) that the diffusion component of the solute flux in the axial direction is negligible compared to the convection component. However this situation does not occur for the radial direction, Equation (18). Clearly, diffusive processes (in this case passive) can be important for solute exchange between the sieve tube and its surroundings, and especially if the radial convection flow is small. Figure 3 shows a comparison of zeroth and first order concentration, i.e., the numerical solutions of Equations (31) and (C.8) for the parameter values of Table 1 and different values of the reflection coefficient σ. First order concentration never exceeds the order of magnitude of zeroth order. Together with the small aspect ratio ε it is thus clear that first order in general does not contribute significantly for this set of parameters (Table 1). This implies that out of the set of dimensionless parameters governing the model behavior together with the three boundary condition values p^_i, p^_out (ẑ) and ĉ_out (ẑ), there are only three independent significant ones: L^p, Ĥ and the reflection coefficient σ. Permeability P^_s depends on σ, see Equation (A.2). Aspect ratio ε, Reynolds number Re (i.e., the effect of inertial forces) and Peclét number Pe_r play a role only when the first order of the series expansion becomes relevant, such as in wider sieve tubes for which ε is bigger than we have considered (Table A1).

For a semipermeable membrane (σ = 1), the exploration of the model behavior for a wide range of possible values of the dimensionless parameters has been provided by Thompson and Holbrook (2003b). Extending this approach to semipermeable membranes would be laborious due to the additional parameters involved. Instead we performed a sensitivity analysis of our model in order to compare the significance of the model parameters, with particular focus on the role the reflection coefficient plays compared to the other parameters, varying each model parameter and boundary condition stepwise from −50 to 50% in relation to the parameter values of Table 1 [which are somewhat similar to the standard parameter set of Thompson and Holbrook (2003b)]. Figure 4 shows the relative change of the state variables pressure, average solute concentration, average axial velocity and radial velocity at the sieve tube membrane. Variation of the parameters causes the dimensionless values of pressure at the inlet of the tube p^(z = 0) to differ even if initial pressure p_i is kept constant. To avoid this effect we converted dimensionless results back into original dimensions for the sensitivity analysis. Since the values of the state variables (except for radial velocity) are fixed by boundary conditions at the inlet of the tube, see Figure 2, changes at the end of the tube (i.e., z = L) are taken as indicator of model sensitivity. Boundary conditions p_i, c_i, and U correspond to state variables p, 〈c〉, and 〈u_z〉, respectively. In cases where these boundary conditions vary, curves of related state variables are not shown in Figure 4 because the variation of initial values dominates the changes of these variables.

Variation of parameters and boundary conditions is straightforward except for the reflection coefficient σ ∈ [0, 1] which is the only parameter limited by a lower as well as an upper boundary. We chose σ = 1 as reference and show negative variation only (such that −50% corresponds to σ = 0.5). Doing so, σ overall appears as one of the most sensitive model parameter, see Figure 4. Given the difficulty of comparing results for varying σ with other (only partially limited) parameters, this result should be treated cautiously regarding a direct numerical comparison of the fractional sensitivities. Still it shows that even small variations of sieve tube solute permeability have a strong effect on the resulting system behavior. Thus there remains no reason to ignore radial solute exchange as it apparently plays in the same league as hydraulic conductivity, geometry and other system properties.

Another obvious result from Figure 4 is that apoplast pressure, concentration and respective gradients (Table 1) do not contribute to the system behavior as strongly as other parameters and boundary conditions. The apoplast surrounding the phloem has been considered to have a constant water potential in most phloem transport models; very few authors (Tyree et al., 1974; Weir, 1981; Thompson and Holbrook, 2003b) considered the more realistic situation of a water potential varying with height. Our result suggests that a water potential gradient can be disregarded, since it constitutes the least influential parameter for all state variables. Nevertheless, in a branched architecture the apoplastic water potential of alternative sinks may well be important (Kaufman and Kramer, 1967; Lang and Düring, 1991).

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.