ON THE HEAT IMPULSE METHOD FOR DEDUCING SAP FLOW 173 Now

Speed of sap flow in plants and trees is of interest to botanists and environmentalists because of its connection with the rate of utilisation of nutrients in the soil. An established method uses the transport of heat where an impulsive heat source is introduced along a radial line by a probe in the trunk sapwood. The temperature is monitored, upstream and downstream, and, by solving the heat flow equation in the moving fluid, the sap velocity may be deduced indirectly under some simplifying assumptions which chiefly render the method most useful when applied to trees of relatively large diameter. Transform methods are used to obtain the appropriate three- dimensional time-dependent solution in explicit form and values for the resulting sap velocity are compared with the existing two-dimensional theory.

Introduction Botanists and others concerned with forest management are interested in know- ing the speed at which sap rises in trees, since it is directly related to the rate at which nutrients in the soil can be utilised.Its variation can give valuable infor- mation about the health of tree species.Various methods of finding the sap speed are needed in practice including introduction of dyes, for example, but a common easily-applicable one in the field employs a probe which can introduce an impulsive heat source at time t = 0 into the sapwood of the trunk and then monitor the temperature 1 cm or so above and below its location.Information about the speed of flow, u, of sap in the tree trunk (see Edwards, et. al. (1996) for a proposed stan- dard nomenclature) can then be deduced.The method is therefore indirect and a satisfactory prediction of the speed relies on an accurate mathematical model for determining the temperature field and its dependence on u.Various assumptions have been made in the modelling process in order to simplify the mathematical so- lution of the problem and thereby render the determination of u as straightforward as possible.
If no account is taken of the curvature of the tree trunk, fully two-dimensional treatments are possible of which the simplest is that of Marshall (1958).In that work a solution for the temperature is used in which it is assumed that the heat source is an instantaneous line source along the whole z-axis directed normally to the trunk surface, which is itself assumed to be the xy plane in this approximation.Heat is thus convected by the sap and is conducted in the x direction, parallel to the trunk, and in the y direction.Clearly, in the Marshall model, there is no flow of heat in the z direction.
In order to take some account of variations, including heat flow, in the z direction, we here assume the thickness L of the xylem is finite with an instantaneous line source placed along the z-axis in 0 < z <_ L. A Newton's Law of Cooling condition is then taken to hold at the xylem-hardwood boundary.However,the assumption of a rectangular geometry is retained so the model would apply best to trees with trunks of relatively large radius.

The Marshall Model
In order to describe the mathematical techniques employed in a simple special case we here derive the solution of the Marshall problem.In mathematical terms, the problem is described by the partial differential equation, OT cOT ( 02T c92T' 0-'-+ u-ffx _ + -y ) + Qh(x)5(y)5(t), in -oc < x, y, z < oo, t > 0, coupled with the conditions, T-+To as x--+-oc, +oc, y-+-cx, + oo, T= T0 at t=O.
Here c is the thermal diffusivity of the xylem and To is the constant ambient temperature.Q represents the heat source term whose dimensions are temperature x (length) 2.
It is convenient to use the temperature rise T-To W as a new dependent variable.Then W satisfies and W --+ 0 as x -c, +oo, y -oc, +cx, "[ We shall use transform methods to solve the initial-boundary value problems in this paper (see for example Duffy (1994)).
In this case we take the Laplace Transform of Equation (1) with respect to the time variable.Thus, OW t (1 where wt(,,) w(,,t)-dt.
On integrating this equation from y 0_ to y 0+, we obtain the jump condition, The solution of Equation (4) which satisfies Equation (5), is continuous at y 0 and which tends to zero at infinity, is where wLF and we take Re(/) > 0. Consider W LF for y > 0. The inverse Laplace Transform is in the usual notation for the Bromwich integral.If we replace s formally by a new variable of integration v, defined by, where i_ _i,,c_io f+ioo exp(_yv/-#)x/-#exp(cvt) dv In order to evaluate I we cut the complex v plane along the negative real axis as shown in Figure 1.The integrals round the curved paths can be shown to Im(v) Re(v) Figure 1.Contour in complex v plane for evaluating; I in (7) vanish in the limit as the radii of the large and small circular arcs tend to c and 0 respectively, so that On AB put v-.2eri, 0 < < cx', hence exp(-iy-a)t) The integral here may be evaluated as Re exp(iAy aA2t) dA exp(-y2/4at) { fo dX} 1_ v/r/at exp(_y/4at) 2 Hence, Y >0' >0.The same result holds for y < 0, since the whole solution is even in y, and hence the result quoted by Marshall (1958) for this problem is verified.Although the result is standard the method used will now be adapted to the more general problem discussed in Section 3.

Finite xylem thickness
Here, as mentioned earlier, we take account of the finite width of the sap region but, for simplicity, we retain a two-dimensional geometry, the boundaries of the xylem being represented by the planes z 0 and z = L.However, we consider the heat flow produced when an instantaneous point source of strength Q dzo/L, where dzo is an element of length, is placed at z z0, say, along the line x = y 0 (see Figure 2).
Thus, it is assumed that there is no flow of heat across the boundary at z 0 (the bark) but that at z = L, the boundary between sap wood and dry wood, the temperature satisfies Newton's Law of Cooling.The parameter K is not expected to be known.The method follows that of the previous section.The Laplace Transform in the t variable, WL, satisfies, Q dzo w + w (w + w + Wz) + ()()(-z0) (12) and the same boundary conditions (9), (10) and (11).On now taking the double Fourier Transform of (12) in the x and y variables we obtain, after some rearrange- ment, 5(z Zo) ( 2aL where W L (x, y, z, s) exp(-iwx iry) dxdy The boundary conditions on P-W LF':Fv are (10) and (11).
Taking first the region z < z0, the inverse Laplace Transform is wF= F, = 1 tic+ i 2rri !., c-ioo P e st ds 2rL /exp[(-aw2 ar2 iwu)t] /c,+, {x/ cosh[v/(L_ z0)] + If sinh[v/(L-Zo)]} cosh(v/-z) e tdv, c,-ioo x/ [x/ sinh(v/ L) + K cosh(x/ L)] (16 where the variable of integration v is here defined by s iwu 2 o.2 v--+--+w + Although the integrand contains v 1/2 it is actually a meromorphic function as an expansion in powers of v demonstrates.We find such an expansion to have the form l + K(L-zo) {l+O(v)} K so there is no pole at v 0, but poles of the integrand occur at values of v given by v/ sinh(v/-L)+ K cosh(/v L) O, (17) in the complex v plane in which no cut is required.
On writinga + ib, the real quantities a and b satisfy sinh(2aL) The assumption that a > 0, K > 0, L > 0 implies that the two sides of (18) have opposite sign since each denominator is positive, and a similar conclusion follows if a < 0. Hence a = 0 and the required solutions of Equation ( 19) correspond to pure imaginary values of x/.Equation (17) becomes tan(bL)-K/b (20) A sketch showing typical solutions is given in Figure 3.We see that solutions occur at bbn, (n 1,2, 3,...), the intersections of y-tan(bL) and y-K/b, where (n-1)Tr/L < b, < (2n-1)Tr/2L !-rd2L /l",J /' /' Figure 3. Sketches of y tan(x) and y const/x, whose intersections give the solutions x bl, b2, b3 Hence the singularities in the integrand of (16) occur at an infinite set of negative real values Vn -b2 and are, in fact, simple poles (justified a posteriori).The residue of the singularity at the nth such pole is {x/ cosh[x/n(L z0)] + K sinh[vr(L z0)]} cosh(vrn z) exp(avnt) lim ( The limit in (21), by L'H6pital's rule, is found to be (1 + KL) sinh(x/ L) + Lx/ cosh(v/n L)' the finiteness of which justifies the simple pole assumption.
The inversion of W F*F gives the same standard result as before so we obtain the temperature rise, Q dzo W-T-To= 2rLat exp{-[(x ut) 2 + y2]/4at} E f" (z, t; b, K, L) n--1 for z < z0.An analogous analysis can be carried out for the region z0 < z < L using Equation (15) and we obtain
The expression (24) is the temperature rise at any point (x, y, z) and time t, due to a point source of strength Qdzo/L situated at the point (0, 0, z0).On summing the temperature rises from all such sources placed continuously along the z axis between z = 0 and z = L, for the case where Q is constant, we integrate with respect to Zo and obtain the result for a uniform line source as, This solution is the basis for determining the sap velocity as described in Section 5. Its general behaviour is given in the next section and some properties and implications are discussed in Section 6.

Behaviour of the solution
The temperature rise is obtained by summing the series in Equation (25) for spec- ified values of the parameters a, K, L and u and at any point (x, y, z) at time t.In order to carry this out it is convenient to transform the problem to dimensionless variables" -x/L,y/L, -z/L, -t/t* where t* is a suitable time scale whose order of magnitude is that used in the probe.
Note that b,, is a function of e KL from (20), which becomes tan(b) e/b (27 with solutions bn, n 1,2, 3,..., satisfying (n-1)rr < bn < (2n-1)r/2.In summing the series in (26) the quantities bn must first be calculated by solving Equation (27).The convergence properties of the series are good because of the -2 exponential decay with bn, but at small values of , many terms may contribute to the sum.An important property in the present application is the near-independence of the sum on the parameter e, at least at values of cd typical of sap flows.
The behaviour of the solution (26) for the temperature rise is, in general terms, similar to that used by Marshall (1958).Thus, the heat pulse is convected with the sap velocity and simultaneously diffuses in the x, y and z directions.In the Marshall model there is no diffusion in the z direction because OT/Oz 0 identically.As an example of the present model we choose L 10 cm, c .0025 cm2/sec, e 1 (so that K 0.1 cm-1).
Figure 4 shows the temperature rise at the point x-1.5 cm, y=0 cm, z=5 cm, plotted as a function of time at various sap velocities u (in the range 0 to 60 cm/hr).
Similarly, for the same parameter values the temperature rise is plotted against x in Figure 5 for different times (in seconds) when u 30 cm/hr. Figure 4. Variation of temperature rise with time at x 1.5 cm, y 0, z 5 cm for different sap velocities.
These diagrams correspond to Figures 2 and 1 respectively presented by Marshall (1958) and are similar in shape.Thus, along the centre line y 0, z 5 cm of the xylem, the heat transfer behaves qualitatively the same as in the Marshall model.
In order to see the effect of the finite width of xylem, assumed in the present study, we present Figure 6, which shows the temperature rise distribution across the xylem for u = 30 cm/hr at different times along the line x = 1.5 cm, y 0. As expected on physical grounds, the solution yields a temperature rise which is uniform across most of the z-range, but which shows a fall-off to accommodate the new boundary condition at z = 10 cm.Thus, some of the heat generated by the line source will be transported into the hard wood and the question is how significant this is.x (cm) Figure 5. Variation of temperature rise W with x coordinate at y 0, z 5 cm, u 30 cm/hr at different times.

Determination of u from the solution
We assume that temperature measurements can be taken at y = 0, z = L/2 and z = -+-h at various times and that the thickness L of the zylem is known.Thus, if the instantaneous heat source Q is applied at t = 0, we measure the temperature rises W1 and W_I at x h and x -h, respectively, at t = /,From (24), we obtain, on reintroducing the dimensionless variables defined in Section 4, W1 exp{-[(l-tfi-l) 2 W-1 where/z-h/L, since all other terms in (25) are independent of x.Hence, ( relating the two unknown parameters and .Note that -0 in Equation (28) corresponds to W1 W-1.The remaining unknown parameter in the problem is We also measure temperature rises W2 and W3 at times 9 0 and, say, 5 0.5 as a suitable point in the sapwood.It is convenient to choose 2 21 and 3 3{1 so that [2 3 .T hen, again using Equation (26) and eliminating Q we obtain, where k(w) sin(w)cos(w/2)/[2w + sin(2w)], from which we deduce that }(,,) The equations (28), ( 29) and (30) are to be solved for fi, 5 and e.At small or moderate values of 6[, such as arise in the area of application of the sap velocity problem, it turns out that the sum of the series in (26) is very insensitive to changes in the parameter e. Very large variations in e can be substituted into the series without affecting the temperatures obtained on the centre line 0, 5 0.5.
This reflects the independence of these temperatures on the thermal properties of the hardwood since e is a dimensionless form of K in Equation (11).Consequently, we eliminate between (29) and (30) by subtraction to obtain 4W [ k(,) exp(-&2n)][P, k(, as the equation to solve for (i for a given e, say e 1 (see also Section 6).In fact this equation can be readily solved by root search and bisection and the resulting value of 0 can be used to find fi from Equation (28).The consistency can be checked by verifying that Equation (29), for example, is satisfied to a good approximation for the corresponding values of O and .0.8 0.9 Figure 6.Variation of temperature rise W with 5 at x 1.5 cm, y 0 when u 30 cm/hr at different times.

Discussion
It is of interest to determine whether results for the estimation of the sap velocity differ in the present method from those obtained using the earlier model of Marshall.Direct observational data is difficult to obtain so, as a test of the algorithm suggested in Section 5, the temperature rises W-l, Wz, W, W3 were computed for certain values of the parameters, viz.e 1, O .0025, fi-.08333333 at the points corresponding to -]-0.15, -0, 2-0.5, at {1 1, {2-2, [3-3. (These dimensionless values correspond to a length scale L 10 cm thickness of the xylem and t* 100 sees, so c 0.0025 cm2/sec and u 0.008333 cm/sec 30 cm/hr.)The following temperature rises were calculated: W_ 0.00432024 Wz 0.64118037 _/Tr W2 0.49310355 Q/rr W3 0.23884349 (32) (Since only ratios of W's are used the constant factor Q/r is not significant.) Eight decimal places, obtained with double precision calculations on the com- puter, are retained in (32) so as to provide adequately for loss of significance in the inverse calculation, in which these values of the W's are substituted into Equation (31) and the equation solved as described in Section 5.With a 0 as initial guess, the values of a and u are readily recovered; these solutions are unchanged to 6 decimal places over a wide range of e.However, it turns out that the Marshall model gives the same values of a and u for the temperature rises quoted in (32).
This prompted a closer look at the relationship between the two solutions.We can write our solution (25) for the temperature rise as W exp{-[(x-ut) 2 + y]/4at} S, Lrrat where E sin(b,L)cos(bnz)exp (-abnt) 2b, L + sin(2b, L) n=l and depends on z and as well as the parameters a and e.The method used for finding a from this expression involves only ratios of W values, for example as quoted in (31).Now, if the series S in (33), over the relevant range of parameters, is essentially constant (i.e.independent of z or t or e), then the ratio of W values at different times reduces to the values which would be obtained from the Marshall model.Thus we are led to consider in more detail the series, S(5, :, c) E sin(n) _cos(nS) e_xp(-Cn) 2bn + sin (2bn) where cand , tan(,) c, n 1,2, 3, It is shown in the Appendix that S(5, c, O) 1/4, for all 5, 0 < 5 < 1, and e > 0, although the rate of convergence of the series is slow.Hence, if e is small and positive we would expect that S(5, e, c) 1/4, with an increased rate of convergence (since for sufficiently large c, S is essentially determined by its first term).It follows that, for small c, the ratio of two such series is unity to a good approximation so that the Marshall solution applies.In the observations, if we assume typical maximum values of { are about 5 and of & about 0.0025, we see that c is indeed small; Table 1 gives values of S(5, e, c) for different c and e at three values of 5.
For the biological application, with c having a maximum value of about 0.01, the relevant part of the table is towards the top row.There the temperature rise is almost independent of e and 5 and the general conclusion seems to be that the Marshall theory gives an accurate prediction of W except possibly near the xylem-hardwood boundary.Situations where the fall-off calculated for W near this boundary might play a more important role could arise in investigations into effects of possible non-uniformities in sap velocity.There is evidence of a reduction in sap velocity near the internal xylem-boundary but to study this would require a more detailed investigation of the heat flow in the vicinity of this boundary.

Table 1 .
Values of S(5, e, c).The three entries in each position of the table refer to 2 0.2, 0.5, 0.8 from top to bottom respectively.