Michel Haritopoulos
email: michel.haritopoulos@univ-orleans.fr

Cécile Capdessus
email: cecile.capdessus@univ-orleans.fr

Asoke K Nandi
email: a.nandi@liverpool.ac.uk

Foetal PQRST Extraction from ECG Recordings using Cyclostationarity-Based Source Separation Method

This work proposes a novel foetal electrocardiogram (FECG) extraction approach based on the cyclostationary properties of the signal of interest. The problem of FECG extraction can easily fit in a blind source separation (BSS) framework; taking into account specific statistical nature of the signal, that one wants to extract, leads to an algorithm able to estimate the FECG contribution to ECG recordings where the maternal ECG is predominant. We show that the proposed procedure provides estimates of the FECGs PQRST complexes without incorporating any prior knowledge concerning PQRST features. Discussions about foetal heart rate variability (HRV) estimation and future works conclude this paper.

I. INTRODUCTION

Given its numerous fields of application, the problem of source separation has been extensively investigated these last few decades [START_REF] Choi | Blind Source Separation and independent Component Analysis: A Review[END_REF], [START_REF] Comon | Handbook of Blind Source Separation -Independent Component Analysis and Applications[END_REF], [START_REF] Jutten | Advances in Blind Source Separation (BSS) and Independent Component Analysis (ICA) for non linear mixtures[END_REF], [START_REF] Comon | Independent Component Analysis, a New Concept[END_REF]. It copes with the possibility of retrieving one specific source from a set of sensors. The method is based on the assumption that the different sources reach the sensors through different propagation paths, so that the resulting mixtures are different. Different scenarios can be found, depending on the type of propagation and the hypotheses that can be made about the nature of the sources. Indeed, the sources can be either filtered (convolutive mixtures) by the propagation paths, or they can be delayed equally or not at all (instantaneous mixtures), which gives different mixture models. The classical hypothesis for the blind separation is independence of the sources, which leads to criteria based on second order and higher order statistics.

Here we address the problem of extracting one specific source rather than separating all the sources. This can be described by the classical equations given below. Let s(t) = [s 1 (t), s 2 (t), ..., s N (t)] † denote the source vector with † the transpose operator. These sources can be the mother's or the foetus' heartbeats, the breathing and other physical phenomena that generate some electromagnetic variations.

We suppose that at least as many sensors as sources, i.e. N , are available and the observations can be brought together in an observation vector

x(t) = [x 1 (t), x 2 (t), ..., x N (t)] † .
The attenuation due to the different propagation paths can be modelled as a N × N matrix A, called mixing matrix, which relates the observations to the sources:

x(t) = As(t). (1)
Let us suppose, without any loss of generality, that the source to be extracted, i.e. the foetus' heartbeats, is the source s 1 (t). The extraction then consists in estimating a vector

B = [b 1 , b 2 , ..., b N ] such that z(t) = Bx(t) is an estimate of s 1 (t).
The method we propose is simple and it is based on the measured value of the foetal heartbeat rate and on second order statistics. As it will be shown further, it leads to an accurate estimation.

II. DATA AND METHOD

A. Experimental Data

The dataset to which we applied the cyclostationary-based source separation method is the well-known DaISy database [START_REF]DaISy; Database for the Identification of Systems[END_REF]. It consists of eight electrocardiogram recordings (Fig. 1). Each of them represents a 5 sec length signal as captured by skin electrodes placed at different locations in the pregnant woman's body. The first five signals come from electrodes placed in the mother's abdominal region; the last three are obtained from mother's thoracic region. The sampling rate is 500 Hz and the eight observations are recorded simultaneously.

The mother's contribution to these ECGs is clearly visible in the last three signals contrary to the foetus' contribution which is dominated by the maternal ECG. This is mainly due to the much longer path taken by the signal generated from the foetus' heart until it arrives to the mother's chest located electrode. The electrical activity of the foetal heart is also disturbed by other noise signals, such as, the 50 Hz power line interference, random electronic noise due to the instrumentation and baseline wandering [START_REF] Assaleh | A Novel Technique for the Extraction of foetal ECG Using Polynomial Networks[END_REF].

The aim of the blind source separation technique we describe hereafter is to extract the signal of interest, i.e. the FECG, from a set of observations, also called mixtures.

B. Proposed Extraction Method

It has been shown in the literature that, under some assumptions concerning ECG's effective frequency and the distance between electrodes and the cardiac sources [START_REF] Comon | Handbook of Blind Source Separation -Independent Component Analysis and Applications[END_REF], [START_REF] Sameni | A Review of Foetal ECG Signal Processing; Issues and Promising Directions[END_REF], ECG recordings can be modelled as linear instantaneous additive mixtures. Moreover, artifacts introduced by breath and shift positions, and the propagation medium being timevarying, the stationary mixing is not always a good approximation [START_REF] Sameni | A Review of Foetal ECG Signal Processing; Issues and Promising Directions[END_REF].

Our algorithm is based on the assumption that the foetus' heartbeat, if not strictly periodic, is at least repetitive, i.e. cyclostationary, at a frequency α 0 that can be measured from the recordings. The extraction criterion is based on the hypothesis that no other source is cyclostationary at the same frequency, which is a realistic one since the mother's heartbeats and her breathing, for instance, are much slower phenomena. The criterion C(B) to minimise is the ratio between the power |P z | of the estimate and its cyclic power |P α0 z | at frequency α 0 :

C(B) = P z (B) P α0 z (B) . (2)
When a noiseless additive mixture model can be assumed, this criterion has been shown in [START_REF] Capdessus | Source Extraction Algorithm based on Cyclic Properties[END_REF] to exhibit a unique minimum over B which leads to the extraction of the cyclostationary source corresponding to the selected cyclic frequency. The source is retrieved with an indeterminacy over its amplitude, but its contribution to the observations could then be estimated by further processing.

The criterion (2) can be estimated from the observations and the matrix B by the expression:

C(B) = B R x (0)B † B R α0 x (0)B † (3)
where

R x (0) = x(t)x † (t) θ and R α0 x (0) = x(t)x † (t)e -2πjα0t
θ are the estimates of the covariance matrix of the observations and their cyclic covariance matrix at frequency α 0 , respectively, and . θ stands for the temporal averaging over θ seconds.

C. Cyclic Frequency Estimation

A survey of the literature on cyclostationarity can be found in [START_REF] Gardner | Cyclostationarity: Half a century of research[END_REF]. The method we used in order to estimate the cyclic frequency of the foetus' heartbeats is based on two assumptions. First, that the foetal heartbeat rate is greater than the mother's one and, second, that the dominant component of all the recordings is the mother's heartbeat one. Based on these two hypotheses, the foetal heartbeat rate can be measured from a plot of the envelope spectra of the observations.

The envelope spectrum of a signal s(t) is computed following the next three steps :

Apply the Hilbert transform to s(t), in order to generate the analytical signal s(t).

Compute the squared module of the analytical signal

|s(t)| 2 .
Apply Fourier transform to this squared analytical signal. Working on the analytical signal instead of s(t) itself, ensures that the high frequencies generated by squaring the signal do not lead to aliasing.

The envelope spectrum has been shown [START_REF] Randall | The Relationship between Spectral Correlation and Envelope Analysis in the Diagnostics of Bearing Faults and Other Cyclostationary Machine Signals[END_REF] to be a simple and efficient mean for detecting cyclostationary components, which produce spectral lines. If computed over the ECG signals measured on a pregnant woman's abdomen, one can expect to observe a strong family of spectral lines corresponding to the mother's heartbeats and a smaller one, at a higher rate, whose fundamental frequency is the foetal heartbeat rate.

The envelope spectrum of the first ECG channel is displayed on Fig. 2. The mother's heartbeat is the dominant phenomenon and produces a family of strong harmonic spectral lines whose fundamental frequency gives the average mother's heartbeat rate. A smaller spectral line appears between the first and second harmonics of the mother's component, which we suppose to correspond to the foetal heartbeat rate. We measured the frequency at its maximum and used it as the extraction cyclic frequency; we hence estimated α 0 around 4.49 Hz. We first applied the algorithm to the whole set of 8 sensors and then to a set of 4 sensors placed on the mother's abdomen, i.e. close to the foetus. What follows is a discussion about the obtained results.

III. RESULTS AND DISCUSSION

A. The Extracted Foetal ECG Signal

For our experiments, we used first the whole DaISy dataset and, then selected four recordings : x(t) = [x 1 (t), x 2 (t), x 3 (t), x 5 (t)] † . Each component of vector x(t) is a 1-dimensional signal of 2500 samples. We considered signals from channels 1, 2, 3 and 5 because they are placed on the abdomen and the contribution of the foetal heartbeats to these channels is significant. Information from channel 4 has not been used as our other experiments have shown that it introduces a high baseline wandering effect into the extracted signal. The set of observations is the input vector to our extraction algorithm which then minimises the criterion C(B) and outputs an estimated FECG signal denoted z(t). The only prior information the algorithm needs is the fundamental FECG cyclic frequency α 0 . Fig. 3 shows the normalised extracted FECG signal after application of our algorithm (y-axis scale is arbitrary) to the x(t) mixture vector (upper plot); for comparison purposes, we also plotted the normalised extracted FECG signal by using the whole recordings set (see section III-C for discussion). It exhibits 22 consecutive foetal heartbeat cycles. All these cycles are fairly cleared of interferences with mother's heartbeats or other noise sources. However, some artifacts still exist: a closer look near peaks around sample 2000 reveals the existence of lower height second peaks. One cannot say with certainty to which source of interference they do correspond. But the overall quality of the estimated FECG signal of interest is fairly good. Fig. 4 plots the spectra of the estimates; it can be seen that the mother's heartbeat (MHB) signal has been rejected from the 4 channel's estimate (top figure), except its first harmonic for the 8 channel's one (bottom figure).

Note, that, all BSS-based techniques estimate unknown source signals up to a sign term, permutation and amplitude scale. The advantage of the proposed technique is that we know, a priori, which signal is being estimated because we used its cyclostationary frequency to achieve our goal; so, we only extract the signal of interest. Also, expert eyes can easily decide which of the estimated signal z(t) and its opposite -z(t) is correct. Finally, the only indeterminacy that still persists is the one concerning the amplitude scale.

B. Foetal Heart Rate Variability

Foetal heart rate variability (fHRV) is an important parameter that can provide early information to clinicians about the well-being of the foetus. It is defined as the normal irregular changes and fluctuations in the foetal heart rate. Short-as well as long-term fHRV monitoring may be useful for early detection of possible troubles concerning the foetus' health during pregnancy. Thus, we try to measure HRV from the extracted FECG using 4 recordings of Fig. 3 by detecting the peak of larger amplitude value, and hence, computing the fHRV from sets of 2 consecutive beats. Given that there are 22 foetus's heartbeat cycles to the estimated FECG, one can see in Fig. 5, 21 plots of the same star-shape marker that provide an estimation of the fHRV. The trend of foetal's instantaneous heart rate [START_REF] Barros | Heart Instantaneous Frequency (HIF): An Alternative Approach to Extract Heart Rate Variability[END_REF]

C. Related work

These last years, different source extraction and/or separation methods based on the periodic properties of the same ECG signals as the one we cope with, have been referenced in the literature. Jafari et al. [START_REF] Jafari | Sequential Blind Source Separation Based Exclusively on Second-Order Statistics Developed for a Class of Periodic Signals[END_REF] proved that, for periodic signals, the covariance matrix estimated for a time-lag equal to the fastest periodic component has properties that can be used for source separation. In [START_REF] Tsalaile | Sequential Blind Source Extraction For Quasi-Periodic Signals With Time-Varying Period[END_REF], [START_REF] Tsalaile | Blind Source Extraction of Heart Sound Signals From Lung Sound recordings Exploiting Periodicity Of The Heart Sound[END_REF], Tsalaile et al. proposed to extract a periodic signal by jointly diagonalising correlation matrices at time lags proportional to the period of the source to be extracted. Sameni et al. [START_REF] Sameni | Multichannel Electrocardiogram Decomposition Using Periodic Component Analysis[END_REF] proposed a method which consists in maximising the periodicity of the extracted signal at the frequency of the source of interest. The extraction is performed by jointly diagonalising two correlation matrices computed respectively at zero time lag and at a time-lag corresponding to the period of the source to be extracted. These two methods take into account the HRV and a prior measure of each period duration is required, which is easy on the mothers heartbeats but involves preprocessing in the case of the foetus' heartbeat. Sameni et al. [START_REF] Sameni | A Deflation Procedure for Subspace Decomposition[END_REF] presented an iterative algorithm that combines at each step both a diagonalisation/projection technique and a filtering technique in order to extract the components at the foetus' heartbeat rate. A periodicity measure is proposed in order to detect the mother's heartbeat interference on the estimate of the foetal one. This periodicity measure (PM) is given by :

P M = |E{x(t)x(t + τ t )}| E{x(t) 2 }E{x(t + τ t ) 2 } * 100%, (4)
where the period τ t depends on time, in order to take into account the heartbeat rate variability (PM ∈ [0, 100]).

The algorithm that we propose is simpler in that it requires only a measure of the mean cyclic frequency of the source to be extracted. In order to evaluate its performance we computed the PM proposed by Sameni et al. with the frequency of the mother's heartbeats. The results cannot be compared with the ones presented in [START_REF] Sameni | A Deflation Procedure for Subspace Decomposition[END_REF] since we used a mean value of the period, but it allows the evaluation of the performances depending on the set of sensors used.

When the signal contains no component at the tested periodicity, the PM should be 0%. The PM obtained with the 8 sensors is equal to 6.4%, while with the 4 sensors it reaches 0.3%. This can be explained by the fact that the sensors placed on the mother's chest bring few new information about the foetus heartbeats, while a better extraction can be performed using only the abdomen ones.

IV. CONCLUSIONS

The novel FECG extraction algorithm presented in this paper is based on the property of cyclostationarity that ECG signals exhibit at the frequency of the foetal heartbeat rate. The extraction process uses a BSS scheme. The algorithm is tested on real pregnant woman's non-invasive ECG recordings and proved to provide fairly denoised FECG signal estimation. Yet, the estimated heart rate variability of the foetus is not robust enough for clinicians' use.

The obtained results are promising. Present and future research work consists in estimating the contribution of the extracted FECG signal of interest on each acquisition channel, testing the algorithm's robustness against variations of the fundamental FECG cyclic frequency α 0 and, improving the extraction procedure in order to obtain significant fHRV characterisation, as well as a direct comparison with results provided by other methods (e.g., [START_REF] Zarzoso | Noninvasive Foetal Electrocardiogram Extraction: Blind Separation Versus Adaptive Noise Cancellation[END_REF]).

Fig. 1 .

 1 Fig. 1. Eight cutaneous electrode recordings from a pregnant woman

Fig. 2 .

 2 Fig.2. Zoom on the envelope spectrum of the first ECG channel: a mean of computing the foetus's heartbeats cyclic frequency α 0 .

Fig. 3 .

 3 Fig.3. The normalised extracted FECG signal, using four abdominal recordings (solid line) and all eight channels recordings (dash-dotted line).

Fig. 4 .

 4 Fig. 4. Spectra of the extracted FECG signals by using 4 channel recordings (top) and 8 channel recordings (bottom).

Fig. 5 .

 5 Fig. 5. Instantaneous heart rate from 22 consecutive foetus' heartbeat cycles of the extracted FECG signal

  in this figure is not very clear.

		4.7
		4.65
	Instantaneous heart rate of extracted FECG (Hz)	4.3 4.35 4.4 4.45 4.5 4.55 4.6
		4.25
		4.2	0	2	4	6	8	10	12	14	16	18	20	22
								PQRST complex number


Bernard Candelpergher
email: bernard.candelpergher@unice.fr

Marc-Antoine Coppo

A new class of identities involving Cauchy numbers, harmonic numbers and zeta values

Keywords: Mathematical Subject Classification (2000): 11B83, 11M41, 33B15, 40G99 Cauchy numbers, Bell polynomials, Harmonic numbers, Laplace-Borel transform, Mellin transform, Zeta values, Ramanujan summation, Hermite's formula

Improving an old idea of Hermite, we associate to each natural number k a modified zeta function of order k. The evaluation of the values of these functions F k at positive integers reveals a wide class of identities linking Cauchy numbers, harmonic numbers and zeta values.

Introduction

It has been well known since the second-half of the 19th century that the Riemann zeta function may be represented by the (normalized) Mellin transform (cf. [START_REF] Zeidler | Quantum Field Theory I : Basics in Mathematics and Physics[END_REF])

ζ(s) = 1 Γ(s) +∞ 0 t s-1 e -t
1 -e -t dt for (s) > 1 , and from late works of Hermite (cf. [START_REF] Hermite | Extrait de quelques lettres de M. Ch. Hermite à M. S. Pincherle[END_REF]) that one has also

ζ(s) - 1 s -1 = 1 Γ(s) +∞ 0 t s-1 e -t 1 -e -t ∞ n=1 λ n n! (1 -e -t ) n dt for (s) ≥ 1 ,
where λ 1 = 1 2 and λ n+1 = 1 0

x(1 -x) • • • (n -x) dx are the (non-alternating) Cauchy numbers 1 . Improving Hermite's idea, one may, more generally, consider Mellin transforms of type

F (s) = 1 Γ(s) +∞ 0 t s-1 e -t 1 -e -t f (1 -e -t ) dt with f (z) = ∞ n=1 ω n z n n k
for suitable sequences (ω n ) n≥1 of rational numbers. The simplest interesting case ω n = 1 corresponds to the Arakawa-Kaneko zeta function and has been studied extensively in [START_REF] Coppo | The Arakawa-Kaneko Zeta function[END_REF]. In this article, we investigate the case ω n = λ n n! , i.e., we study the function

F k (s) = 1 Γ(s) +∞ 0 t s-1 e -t 1 -e -t f k (1 -e -t ) dt with f k (z) = ∞ n=1 λ n n! z n n k (k = 0, 1, 2, . . . ),
which is a priori defined in the half-plane (s) ≥ 1 but analytically continues in the whole complex s-plane (Theorem 7). We call this function F k the modified zeta function of order k. An evaluation by two different methods of the values of F k at positive integers q leads to a new class of identities linking Cauchy numbers, harmonic numbers and zeta values. In the case k = 0, Hermite's formula for ζ (cf. [START_REF] Coppo | Nouvelles expressions des formules de Hasse et de Hermite pour la fonction zêta d'Hurwitz[END_REF]) is regained, i.e.,

F 0 (q) = ζ(q) - 1 q -1 = ∞ n=1
λ n n!n P q-1 (H (1) n , H (2) n , . . . , H (q-1)

n ) ,
where the polynomials P m are the modified Bell polynomials defined by the generating function exp

∞ k=1 x k z k k = ∞ m=0 P m (x 1 , • • • , x m ) z m ,
evaluated at harmonic numbers H (m) n = n j=1 1 j m . In the simplest higher case k = 1, this extension of Hermite's formula leads to the following new relation (Theorem 10):

F 1 (q) = ∞ n=1
λ n n!n 2 P q-1 (H n , H (2) n , . . . , H (q-1)

n ) = ∞ n=1 log(n + 1) n q + γζ(q) + ζ(q + 1) - ∞ n=1 H n n q - q-1 k=1 1 k ∞ n=1 1 (n + 1) k n q-k ,
where

H n = H (1)
n , and γ = lim n→∞ (H n -log n) is the Euler-Mascheroni constant.

For example, for q = 2, since P 1 (H n ) = H n and ∞ n=1 H n n 2 = 2ζ(3) (cf. [START_REF] Choi | Explicit evaluation of Euler and related sums[END_REF], [START_REF] Coppo | Nouvelles expressions des formules de Hasse et de Hermite pour la fonction zêta d'Hurwitz[END_REF]), then the previous relation may be written

F 1 (2) = ∞ n=1 λ n H n n!n 2 = ∞ n=1 log (n + 1) n 2 + γζ(2) -ζ(3) -1 ,
and this generalizes the known formula

F 0 (2) = ∞ n=1 λ n H n n!n = ζ(2) -1 .
The function F k also has an interesting interpretation in terms of Ramanujan summation (cf. [START_REF] Candelpergher | La sommation de Ramanujan[END_REF]) as underscored by Theorem 11. In particular, one shows the identity

F k (1) = ∞ n=1 λ n n! 1 n k+1 = R n≥1 P k (H n , H (2) n , . . . , H (k) n ) n ,
where, in the right member, R n≥1 denotes the sum (in the sense of Ramanujan) of the divergent series. This raises a kind of reciprocity between F k (1) and F 0 (k + 1).

Preliminaries

The non-alternating Cauchy numbers

Definition 1. The non-alternating Cauchy numbers (cf. [START_REF] Coppo | Nouvelles expressions des formules de Hasse et de Hermite pour la fonction zêta d'Hurwitz[END_REF], [START_REF] Merlini | The Cauchy numbers[END_REF]) are the sequence of (positive) rational numbers (λ n ) n≥1 defined by the exponential generating function

z log(1 -z) + 1 = n≥1 λ n n! z n . ( 1
)
Dividing by z and setting z = 1 -e -t and t > 0, this relation may be rewritten

1 1 -e -t - 1 t = ∞ n=1 λ n n! (1 -e -t ) n-1 . (2)
From (1), one may easily deduce the following recursive relation

n j=1 λ j j!(n -j + 1) - 1 n + 1 = 0 for n ≥ 1.
Example 1. The first non-alternating Cauchy numbers are

λ 1 = 1 2 , λ 2 = 1 6 , λ 3 = 1 4 , λ 4 = 19 30 , λ 5 = 9 4 .

The modified Bell polynomials evaluated at harmonic numbers

Definition 2. The modified Bell polynomials (cf. [START_REF] Chen | Dixon's F 2 (1)-series and identities involving harmonic numbers and the Riemann zeta function[END_REF], [START_REF] Coppo | Nouvelles expressions des formules de Hasse et de Hermite pour la fonction zêta d'Hurwitz[END_REF], [START_REF] Flajolet | Mellin Transforms and Asymptotics: finite differences and Rice's integrals[END_REF]) are the polynomials P m defined for all natural numbers m by P 0 = 1 and the generating function

exp   k≥1 x k z k k   = 1 + m≥1 P m (x 1 , ..., x m ) z m . ( 3
)
The general explicit expression for P m is

P m (x 1 , ..., x m ) = k 1 +2k 2 +3k 3 +•••=m 1 k 1 !k 2 !k 3 ! . . . x 1 1 k 1 x 2 2 k 2 x 3 3 k 3 . . .
One may also compute recursively the polynomials P m by means of the following relation

mP m (x 1 , . . . , x m ) = m k=1 x k P m-k (x 1 , . . . , x m-k ) (m ≥ 1) .
Proposition 1. For all natural numbers m, and each integer n ≥ 1,

+∞ 0 e -t (1 -e -t ) n-1 t m m! dt = P m (H n , . . . , H (m) n ) n , (4)
with

H (m) n = n j=1 1 j m and H n = H (1) n .
Proof. One starts from the classical Euler relation (cf. [START_REF] Zeidler | Quantum Field Theory I : Basics in Mathematics and Physics[END_REF])

B(a, b) = 1 0 u a-1 (1 -u) b-1 du = Γ(a)Γ(b) Γ(a + b) ,
and substitute u = e -t , a = 1 -z, and b = n + 1; then one obtains

+∞ 0 e -t (1 -e -t ) n e tz dt = n! (1 -z)(2 -z) . . . (n + 1 -z) .
Moreover, one has

n! (1 -z)(2 -z) . . . (n + 1 -z) = n! (n + 1)! × n j=0 (1 - z j + 1 ) -1 = 1 (n + 1) × exp(- n j=0 log(1 - z j + 1 )) = 1 (n + 1) × exp( n j=0 ∞ k=1 z k k(j + 1) k ) = 1 (n + 1) exp( ∞ k=1 H (k) n+1 z k k ) = ∞ m=0 P m (H (1) n+1 , . . . , H (m) n+1 ) n + 1 z m (by (3)).
Thus (4) results by identification of the term in z m .

Example 2. For small values of m, one has

P 1 (H n ) = H n , P 2 (H n , H (2) n ) = (H n ) 2 2 + H (2) n 2 , P 3 (H n , H (2) n , H (3) n ) = (H n ) 3 6 + H n H (2) n 2 + H (3) n 3 .

The Laplace-Borel transformation

We consider the vector space E of complex-valued functions f ∈ C 1 (]0, +∞[) such that for all ε > 0, there exists

C ε > 0 such that |f (t)| ≤ C ε e εt for all t ∈ ]0, +∞[ .
In particular, a function f ∈ E satisfies the following two properties: a) for all x with (x) > 0, t → e -xt f (t) is integrable on ]0, +∞[ , b) for all β with 0 < β < 1, t → |f (t)| 1 t β is integrable on ]0, 1[ . We recall now some basic properties (cf. [START_REF] Schiff | The Laplace transform : theory and applications[END_REF]) of the Laplace transformation in this frame which are appropriate for our purpose. Definition 5. Let f and g be two functions in E. The convolution product f * g of f and g is the function defined for all t > 0 by

(f * g)(t) = t 0 f (u)g(t -u) du .
Proposition 3 (cf. [START_REF] Schiff | The Laplace transform : theory and applications[END_REF]). If f ∈ E and g ∈ E, then f * g ∈ E and

L(f * g) = L(f ) L(g) . ( 5
)
Hence, if a ∈ E and b ∈ E then ab ∈ E since ab = L( a * b).

Theorem 1. Let a be a function in E. Then the series

n≥1 λ n n! +∞ 0 e -t (1 -e -t ) n-1 a(t) dt
converges and

∞ n=1 λ n n! +∞ 0 e -t (1 -e -t ) n-1 a(t)dt = +∞ 0 ( 1 1 -e -t - 1 t )e -t a(t) dt . ( 6
)
Proof. By (2)

+∞ 0 ( 1 1 -e -t - 1 t )e -t a(t)dt = +∞ 0 ∞ n=1 λ n n! (1 -e -t ) n-1 e -t a(t) dt .
In the right member, the order of +∞ 0 and ∞ n=1 may be interchanged since

+∞ 0 ∞ n=1 λ n n! (1 -e -t ) n-1 e -t a(t) dt = +∞ 0 ∞ n=1 λ n n! (1 -e -t ) n-1 e -t | a(t)| dt = +∞ 0 ( 1 1 -e -t - 1 t )e -t | a(t)| dt ,
and the convergence of this last integral follows from the assumption that a ∈ E.

Example 3. Let a(x) = 1 x s with (s) ≥ 1. Then a ∈ E and a(t) = t s-1 Γ(s) . Hence ∞ n=1 λ n n! +∞ 0 e -t (1 -e -t ) n-1 t s-1 Γ(s) dt = 1 Γ(s) +∞ 0 e -t ( 1 1 -e -t - 1 t )t s-1 dt =    γ if s = 1, ζ(s) - 1 s -1 if s = 1
where γ refers to the Euler constant. In particular, since

+∞ 0 e -t (1 -e -t ) n-1 dt = 1 n for each integer n ≥ 1 , then γ = ∞ n=1 λ n n! 1 n . 3 The operator D Proposition 4. If a ∈ E, then the integral +∞ 0 e -t (1 -e -t ) x-1 a(t)dt
converges for all x with (x) > 0.

Proof. If a ∈ E and (x) > 0, we may write for t ∈ ]0, +∞[,

e -t (1 -e -t ) x-1 a(t) ≤ e -t e (1-(x))(-log(1-e -t )) | a(t)| .
The convergence when t → +∞ results from the inequality

e -t e (1-(x))(-log(1-e -t )) | a(t)| ≤ e -t 1 -e -t | a(t)| ≤ 2e -t | a(t)| (for t ≥ log 2).
The convergence when t → 0 results from the inequality

e (1-(x))(-log(1-e -t )) ≤    1 if (x) ≥ 1, 1 (1-e -t ) (1-(x)) if 0 < (x) < 1 since the function t → e -t | a(t)| 1 (1 -e -t ) β is integrable at 0 for 0 < β < 1 by the definition of E (note that (1 -e -t ) -β ≤ (kt) -1 for small enough t).
Definition 6. Let a be a function in E. We call D(a) the function defined for all x with (x) > 0 by

D(a)(x) = +∞ 0 e -t (1 -e -t ) x-1 a(t)dt . ( 7
)
Remark 1. a) By Theorem 1, the series n≥1 λ n n! D(a)(n) converges and its sum is given by formula [START_REF] Choi | Explicit evaluation of Euler and related sums[END_REF].

b) The values of D(a) at positive integers may be computed directly without recourse to a. The development of (1 -e -t ) n by the binomial theorem gives

D(a)(n + 1) = n k=0 (-1) k n k a(k + 1) for all integer n ≥ 0. ( 8
) Definition 7. We call Λ the C 1 -diffeomorphism of R + defined by Λ(u) = -log(1-e -u ).
In particular, it is important to note that Λ is involutive:

Λ -1 = Λ .
Theorem 2. Let a be a function in E. Then the function D(a) ∈ E and, moreover, verifies the relation

D(a) = a(Λ) , ( 9
)
where a(Λ) denotes a • Λ.

Proof. The change of variables t = Λ(u) in [START_REF] Coppo | Nouvelles expressions des formules de Hasse et de Hermite pour la fonction zêta d'Hurwitz[END_REF] gives

D(a)(x) = +∞ 0 e -xu a(Λ(u)) du for (x) > 0.
Thus, D(a) = L( a(Λ)). It remains to prove that D(a) ∈ E. One has only to check that the function a(Λ) is in E. This function being in C 1 (]0, +∞[), it suffices to show that for all ε > 0, the function

u → e -εu | a(-log(1 -e -u ))| is bounded on ]0, +∞[. This results from the existence of C ε > 0 such that a(-log(1 -e -u )) ≤ C ε (1 -e -u ) ε for all u ∈ ]0, +∞[ . Example 4. Let a(x) = 1 x s with (s) ≥ 1. Then a(t) = t s-1 Γ(s)
. Thus, by [START_REF] Dilcher | Some q-series identities related to divisors functions[END_REF],

D( 1 x s ) = L Λ s-1 Γ(s) , ( 10
)
and if s = m + 1 with m a natural number and n ≥ 1, then by (4),

D( 1 x m+1 )(n) = P m (H n , . . . , H (m) n ) n . ( 11
)
Remark 2. Theorem 2 may be summarized in the following diagram

E D ----→ E   L -1  L E Λ ----→ E
where Λ ( a) = a(Λ). The algebraic properties of D are summed up in the following theorem. Proof. We can write D = LΛ L -1 and Λ is an automorphism of E which verifies Λ = (Λ ) -1 since Λ = Λ -1 . Furthermore,

D( 1 x ) = L(1) = 1 x .

The harmonic product

Our aim is to define the harmonic product of two functions a and b in E as being the unique function f of E such that

D(a) (x).D(b) (x) = D(f ) (x) .
Thus, we have to establish that such a function exists and is unique. In order to do this, we introduce first a Λ-convolution product of two functions in E. Proof. From the definition of the convolution product, one may write

The Λ-convolution product

a(Λ) * ( b(Λ) (t) = t 0 a(Λ(u)) b(Λ(t -u))du .
Now, for all ε > 0, there exists C ε > 0 and

D ε > 0 such that a(-log(1 -e -u )) ≤ C ε (1 -e -u ) ε and b(-log(1 -e -(t-u) )) ≤ D ε (1 -e -(t-u) ) ε for all u ∈ ]0, +∞[ . It follows that ( a(Λ) * b(Λ))(t) ≤ C ε D ε t 0 (1 -e -u ) ε (1 -e -(t-u) ) ε du .
One has also

t 0 (1 -e -u ) ε (1 -e -(t-u) ) ε du = 1 -e -t 1+2ε 1 0 u ε (1 -u) ε 1 (1 -(1 -e -t )u) ε+1 du ≤ 1 -e -t 1+2ε 1 0 1 (1 -(1 -e -t )u) ε+1 du ≤ 1 -e -t 1+2ε e tε -1 (1 -e -t )ε ≤ 1 -e -t 2ε e tε -1 ε ≤ e tε ε .
Hence, (

a(Λ) * b(Λ))(t) ≤ C ε D ε e tε
ε , which proves that this function belongs to E as required. Remark 3. The Λ-convolution product inherits the algebraic properties of the ordinary convolution product, i.e., bilinearity, commutativity, and associativity. This construction may be summarized in the following diagram:

The harmonic product

(a, b) ----→ ( a, b) ----→ ( a(Λ), b(Λ))       a b ← ----a b ← ----a(Λ) * b(Λ)
Remark 4. The harmonic product inherits the properties of the Λ-convolution product: it is bilinear, commutative and associative. Hence, if the numbers C k,l n are defined by

(X + Y -XY ) n = 0≤k≤n 0≤l≤n C k,l n X k Y l ,
then one has the following explicit formula:

(a b)(n + 1) = 0≤k≤n 0≤l≤n C k,l n a(k + 1)b(l + 1) ,
which can be rewritten in the following equivalent form:

(a b)(n + 1) = 0≤l≤k≤n (-1) k-l n k k l a(k + 1)b(n + 1 -l) (n ≥ 0) .
For small values of n, this enables one to compute 1) ,

(a b)(1) = a(1)b(
(a b)(2) = a(2)b(1) + a(1)b(2) -a(2)b(2) , (a b)(3) = a(3)b(1) + a(1)b(3) + 2a(2)b(2) -2a(3)b(2) -2a(2)b(3) + a(3)b(3) .
Theorem 5. Let

1 x k = 1 x 1 x • • • 1 x k (k = 1, 2, 3, • • • ) , where 1 x denotes (improperly) the function x → 1 x .
Then, for all natural numbers m ≥ 0,

1 x (m+1) = D( 1 x m+1 ) .
In particular, for all integers n ≥ 1,

1 x (m+1) (n) = P m (H n , . . . , H (m) n ) n . ( 14
)
Proof. By (13) we have

D( 1 x m+1 ) = D( 1 x . . . 1 x m+1 ) = D( 1 x ) (m+1) = 1 x (m+1) since D( 1 x ) = 1 x .
Thus, (14) results from [START_REF] Hermite | Extrait de quelques lettres de M. Ch. Hermite à M. S. Pincherle[END_REF].

The harmonic property

The following theorem explains the main reason why the harmonic product is called 'harmonic'.

Theorem 6. Let a ∈ E. Then

1 x a = A(x) x ,
where A denotes the function defined for (x) > 0 by

A(x) = +∞ 0 e -xt -1 e -t -1 e -t a(t)dt .
In particular, for each integer n ≥ 1,

1 x a (n) = A(n) n = 1 n n k=1 a(k) . ( 15
)
Proof. By the definition of the harmonic product, one has

1 x a = L(1 a) . Now (1 a)(Λ(u)) = (1 * a(Λ))(u) = u 0 a(Λ(v))dv = - Λ(u) +∞ a(t) e -t 1 -e -t dt
(by the change of variables t = Λ(v)). Hence,

(1 a)(u) = +∞ u a(t)
e -t 1 -e -t dt .

Thus, we have

1 x a = +∞ 0 e -xu +∞ u a(t) e -t 1 -e -t dt du = +∞ 0 t 0 e -xu du a(t) e -t 1 -e -t dt = 1 x +∞ 0 (1 -e -xt ) a(t) e -t 1 -e -t dt = A(x) x .
Furthermore, for each integer n ≥ 1, we have

A(n) = +∞ 0 e -nt -1 e -t -1 e -t a(t)dt = n k=1 a(k) .
Example 5.

1 x 1 x = D( 1 x 2 ) = L(Λ) = H(x)
x with H(x) = ψ(x + 1) + γ, ψ denoting the logarithmic derivative of Γ. In particular, for each integer n ≥ 1,

1 x 1 x (n) = H(n) n = H n n .
Example 6. For (s) ≥ 1, 1 x

1 x s = H (s) (x)
x , with

H (s) (x) = 1 Γ(s) +∞ 0 1 -e -xt 1 -e -t e -t t s-1 dt . For each integer n ≥ 1, 1 x 1 x s (n) = H (s) (n) n = H (s) n n = 1 n n m=1 1 m s .
From (15), by induction on k, we deduce the following important corollary.

Corollary 1. For each integer

k ≥ 2, 1 x k a (n) = 1 n   n≥n 1 ≥•••≥n k ≥1 a(n k ) n 1 . . . n k-1   . ( 16
)
Example 7. Applying (16) with a(x) = 1 x (and k = m), we get

1 x (m+1) (n) = 1 n   n≥n 1 ≥•••≥nm≥1 1 n 1 . . . n m   . ( 17
)
Hence, it follows from ( 14) and (17) that

P m (H n , H (2) n , . . . , H (m) n ) = n≥n 1 ≥•••≥nm≥1 1 n 1 . . . n m , ( 18
)
which is a nice reformulation of Dilcher's formula (cf. [START_REF] Boyadzhiev | Harmonic number identities via Euler's transform[END_REF], [START_REF] Dilcher | Some q-series identities related to divisors functions[END_REF]).

5 The modified zeta function F k

Integral representation

Definition 10. For all s ∈ C with (s) ≥ 1 and each natural number k, the modified zeta function of order k is defined by

F k (s) = 1 Γ(s) +∞ 0 t s-1 e -t 1 -e -t f k (1 -e -t ) dt with f k (z) = ∞ n=1 λ n n! z n n k . ( 19
)
Remark 6. By (2) and Example 3, one has

F 0 (s) = ζ(s) - 1 s -1 .
The fact that F k may be represented by a Mellin transform enables us to analytically continue this function outside its half-plane of definition by a standard analytic method (cf. [START_REF] Zeidler | Quantum Field Theory I : Basics in Mathematics and Physics[END_REF] section 6.7). λ n z n n! is equal to 1. Thus 1 is also the radius of convergence of the series ∞ n=1 λ n z n n!n k which defines an analytic function f k in the disc D(0, 1). Hence, the function

g k : t → f k (1 -e -t )
is analytic for all t ∈ C such that 1 -e -t ∈ D(0, 1). Since 1 -e 0 = 0, it follows that g k is analytic in a neighbourhood of 0. Since g k (0) = 0, the function t → g k (t) e -t 1 -e -t is itself analytic in a neighbourhood of 0. It follows that its Mellin transform analytically continues in the complex plane with simple poles at negative integers which are all cancelled by the poles of Γ. Theorem 8. For all s with (s) > 1 and each integer k ≥ 1,

F k (s) = ϑ(k)ζ(s)+ k j=1 (-1) j ϑ(k-j)Z j (s)+(-1) k 1 Γ(s) +∞ 0 t s-1 e -t 1 -e -t T k e -t -1 t dt (20) with ϑ(k) = ∞ n=1 λ n n! 1 n k , ( 21
)
Z j (s) = n>n 1 >n 2 >•••>n j >0 1 n s n 1 n 2 . . . n j , ( 22
) T f (t) = +∞ t e -u 1 -e -u f (u)du . ( 23
)
Proof. Formula (20) results from the integral representation (19) and the two following lemmas.

Lemma 1. For all t > 0,

f k (1 -e -t ) = k j=0 (-1) j ϑ(k -j) Λ j (t) j! + (-1) k T k ( e -t -1 t ) ,
where ϑ is defined by (21) and T is the operator defined by (23).

Proof. Let g k (t) = f k (1 -e -t
). The function g k verifies the recursive relation

g k (t) = e -t f k (1 -e -t ) = e -t 1 -e -t f k-1 (1 -e -t ) = e -t 1 -e -t g k-1 (t) . Thus g k (t) = t 0 e -u 1 -e -u g k-1 (u)du = g k (+∞) - +∞ t e -u 1 -e -u g k-1 (u)du with g k (+∞) = f k (1) = ϑ(k).
Thus, one has

g k (t) = ϑ(k) - +∞ t e -u 1 -e -u g k-1 (u)du = ϑ(k) -T (g k-1
) , and a repeated iteration k times of this relation gives

g k (t) = k-1 j=0 ϑ(k -j)(-1) j T j (1) + (-1) k T k (g 0 ) . Now, by (2), g 0 (t) = ∞ n=1 λ n (1 -e -t ) n n! = e -t -1 t + 1 ,
and thus

T k (g 0 ) = T k ( e -t -1 t ) + T k (1) .
Hence 1) and a tauberian theorem), one deduces that

g k (t) = k-1 j=0 ϑ(k -j)(-1) j T j (1) + (-1) k T k (1) + (-1) k T k ( e -t -1 t ) . Since ϑ(0) = ∞ n=1 λ n n! = 1 (by (
g k (t) = k j=0 ϑ(k -j)(-1) j T j (1) + (-1) k T k ( e -t -1 t ) ,
and, now, it remains to prove that

Λ j (t) j! = T j (1) ,
which follows from the recursive relation

Λ j (t) j! = - t +∞ e -u 1 -e -u Λ j-1 (u) (j -1)! du = T Λ j-1 (j -1)! .
Lemma 2. Let Z j (s) be defined by ( 22). Then, for all s ∈ C with (s) > 1,

Z j (s) = 1 Γ(s) +∞ 0 t s-1 e -t
1 -e -t Λ j (t) j! dt .

Proof. From the recursive relation

∂ Λ j (t) j! = Λ j-1 (t) (j -1)! ∂Λ(t) = - e -t 1 -e -t Λ j-1 (t) (j -1)! = - m>0 e -mt Λ j-1 (t) (j -1)! ,
and Λ(t) = n>0 e -nt n , one may check by induction on j that

Λ j (t) j! = n 1 >n 2 >...>n j >0 e -n 1 t n 1 1 n 2 • • • 1 n j .
Furthermore, one has

1 Γ(s) +∞ 0 t s-1 e -N t e -t 1 -e -t dt = n>N 1 n s (for (s) > 1) . Hence 1 Γ(s) +∞ 0 t s-1 e -t 1 -e -t Λ j (t) j! dt = n>n 1 >n 2 >•••>n j >0 1 n s 1 n 1 1 n 2 • • • 1 n j = Z j (s) .

Values of F k at integers

Theorem 9. For all s in C with (s) ≥ 1 and each natural number k, then

F k (s) = ∞ n=1 λ n n!n k D 1 x s (n) . ( 24
)
In particular, for all natural numbers m,

F k (m + 1) = ∞ n=1 λ n n! P m (H n , H (2)
n , . . . , H

(m) n ) n k+1 . ( 25
)
Proof. The change of variables t = Λ(u) in (19) enables to write

F k (s) = 1 Γ(s) +∞ 0 f k (e -u )(Λ(u)) s-1 du . Since D( 1 x s ) = L Λ s-1 Γ(s)
, we deduce (24) from this last expression of F k (s). Moreover, by [START_REF] Hermite | Extrait de quelques lettres de M. Ch. Hermite à M. S. Pincherle[END_REF], one also has D(

1 x m+1 )(n) = P m (H n , . . . , H (m) n ) n
, which proves (25).

Corollary 2. Let ϑ(s) be the Dirichlet series defined for (s) > 0 by

ϑ(s) = ∞ n=1 λ n n! 1 n s .
Then for each natural number k,

ϑ(k + 1) = F k (1) . ( 26
)
Example 8.

F 0 (1) = ∞ n=1 λ n n!n = γ = ϑ(1) , F 0 (2) = ∞ n=1 λ n H n n!n = ζ(2) -1 , F 0 (3) = 1 2 ∞ n=1 λ n H 2 n n!n + 1 2 ∞ n=1 λ n H (2) n n!n = ζ(3) - 1 2 , F 1 (1) = ∞ n=1 λ n n!n 2 = ϑ(2) , F 1 (2) = ∞ n=1 λ n H n n!n 2 , F 1 (3) = 1 2 ∞ n=1 λ n H 2 n n!n 2 + 1 2 ∞ n=1 λ n H (2) n n!n 2 .

Identities linking Cauchy numbers, harmonic numbers and zeta values

Theorem 10. For all integers q ≥ 2,

F 1 (q) = ∞ n=1 λ n n!n 2 P q-1 (H n , H (2) n , . . . , H (q-1) n ) = ∞ n=1 log(n + 1) n q + γζ(q) + ζ(q + 1) - ∞ n=1 H n n q - q-1 k=1 1 k ∞ n=1 1 (n + 1) k n q-k . (27)
Proof. By (20) and (25), one may write (2) n , . . . , H (q-1)

F k (q) = ∞ n=1 λ n n!n k+1 P q-1 (H n , H
n ) = ϑ(k)ζ(q) + k j=1 (-1) j ϑ(k -j)Z j (q) + (-1) k 1 Γ(q) +∞ 0 t q-1 e -t
1 -e -t T k e -t -1 t dt .

(28)

We apply now (28) with k = 1. This gives

F 1 (q) = γζ(q) - n≥1 H n-1 n q + 1 Γ(q) +∞ 0 t q-1 e -t 1 -e -t E 1 (t) dt , with E 1 (t) = -Ei(-t) = +∞ t e -u u du. Thus F 1 (q) = γζ(q) - n≥1 H n n q + ζ(q + 1) + I(q) ,
where

I(q) = 1 Γ(q) +∞ 0 t q-1 e -t 1 -e -t E 1 (t) dt = 1 Γ(q) ∞ n=1 +∞ 0 e -nt t q-1 E 1 (t) dt . Since E 1 (t) = -γ -log t + ∞ n=1 (-1) n-1 n t n n! ,
and -γ -log t = log x x (cf. [START_REF] Schiff | The Laplace transform : theory and applications[END_REF]), then E 1 = log(x+1)

x .

Thus

+∞ 0 e -nt t q-1 E 1 (t) dt = (-1) q-1 log(x + 1)

x (q-1) (n) .
Hence, by a calculation of the (q -1)th derivative, we get

I(q) = (-1) q-1 (q -1)! ∞ n=1 log(x + 1) x (q-1) (n) = ∞ n=1 log(n + 1) n q - q-1 k=1 1 k ∞ n=1 1 (n + 1) k n q-k .
Remark 7. 1) We recall Euler's formula (cf. [START_REF] Choi | Explicit evaluation of Euler and related sums[END_REF])

∞ n=1 H n n q = 1 2 (q + 2)ζ(q + 1) -1 2 q-2 k=1 ζ(k + 1)ζ(q -k) for q > 2, 2ζ (3)
for q = 2.

2) From

∞ n=1 1 (n + 1)n = 1, and the decomposition

1 (n + 1) k n q-k = 1 (n + 1) k-1 n q-k - 1 (n + 1) k n q-k-1 (0 < k < q) ,
the sum of the series ∞ n=1 1 (n + 1) k n q-k may be expressed as a linear combination of zeta values and integers.

Example 9.

∞ n=1 log (n + 1)

n 2 + γζ(2) -ζ(3) -1 = ∞ n=1 λ n H n n!n 2 , ∞ n=1 log (n + 1) n 3 + γζ(3) - 1 10 ζ(2) 2 - 1 2 ζ(2) = 1 2 ∞ n=1 λ n H 2 n n!n 2 + 1 2 ∞ n=1 λ n H (2) n n!n 2 , ∞ n=1 log (n + 1) n 4 + γζ(4) -2ζ(5) + ζ(2)ζ(3) - 2 3 ζ(3) + 1 3 ζ(2) - 1 2 = 1 6 ∞ n=1 λ n H 3 n n!n 2 + 1 2 ∞ n=1 λ n H n H (2) n n!n 2 + 1 3 ∞ n=1 λ n H (3) n n!n 2 .

Link with the Ramanujan summation

The function F k has strong connections with the Ramanujan summation (cf. [START_REF] Candelpergher | La sommation de Ramanujan[END_REF], [START_REF] Candelpergher | Ramanujan summation and the exponential generating function ∞ k=0 z k k! ζ (-k)[END_REF]).

If a ∈ E, then the series n≥1 a(n) may be written However, this last integral may be divergent at 0. Nevertheless we can renormalize it by removing the singularity at zero. This may be done merely by subtracting the polar part 1 t of 1 1-e -t . From Theorem 1, we know that

+∞ 0 ( 1 1 -e -t - 1 t )e -t a(t)dt = ∞ n=1 λ n n! +∞ 0 e -t (1 -e -t ) n-1 a(t)dt = ∞ n=1 λ n n! D(a) (n) .
This justifies the following definition.

Definition 11. Let a be a function in E = L(E). The Ramanujan sum of the series

n≥1 a(n) is defined by R n≥1 a(n) = +∞ 0 ( 1 1 -e -t - 1 t )e -t a(t)dt = ∞ n=1 λ n n! D(a) (n) . ( 29
) Lemma 3. Let a and b in E. Then R n≥1 (a b)(n) = ∞ n=1 λ n n! D(a)(n) D(b)(n) . ( 30
)
Proof. This results directly from ( 12) and (29).

Theorem 11. for all s ∈ C with (s) ≥ 1, one has

F 0 (s) = R n≥1 1 n s and F k (s) = R n≥1 1 x k 1 x s (n) for k ≥ 1 . ( 31
)
Proof. By ( 24) and (30), taking into account the invariance of 1

x by D, one may write

R n≥1 1 x k 1 x s (n) = ∞ n=1 λ n n! D 1 x k (n)D 1 x s (n) = ∞ n=1 λ n n! 1 x k (n)D 1 x s (n) = ∞ n=1 λ n n!n k D 1 x s (n) = F k (s) .
In particular, by ( 14), one deduces from (31) the following identity.

Corollary 3.

For each natural number k,

F k (1) = ϑ(k + 1) = ∞ n=1 λ n n! 1 n k+1 = R n≥1 P k (H n , H (2) n , . . . , H (k) n ) n . ( 32
)
Example 10.

ϑ(1) = ∞ n=1 λ n n!n = R n≥1 1 n = γ , ϑ(2) = ∞ n=1 λ n n!n 2 = R n≥1 H n n , ϑ(3) = ∞ n=1 λ n n!n 3 = 1 2 R n≥1 H 2 n n + 1 2 R n≥1 H (2) n n .
Remark 8. Comparing (32) with

F 0 (k + 1) = ∞ n=1 λ n n!n P k (H n , H (2) n , . . . , H (k) n ) ,
one may observe a kind of reciprocity between F k (1) and F 0 (k + 1). This results from the fact that D = D -1 .

Remark 9. In the case q = 1, ( 27) is meaningless since both the series n≥1 log(n + 1) n and n≥1

H n n diverge. However, since

log(x + 1) -(ψ(x + 1) + γ) = +∞ 0 (e -xu -1)( 1 1 -e -u - 1 u )e -u du , it follows that   log(x + 1) x - ψ(x + 1) + γ x   (t) = +∞ t 1 1 -e -u - 1 u e -u du ,
and then one may easily deduce from (29) the relation

R n≥1 log(n + 1) n - R n≥1 H n n = - γ 2 2 ,
which may be rewritten in the following form (cf. Example 10):

R n≥1 log(n + 1) n = ϑ(2) - 1 2 ϑ(1) 2 .

Link with the Arakawa-Kaneko zeta function

For (s) ≥ 1 and k ≥ 1, one can define in an algebraic fashion the function ξ k by

ξ k (s) = ∞ n=1 D 1 x k 1 x s (n) = ∞ n=1 1 n k D 1 x s (n). ( 33
)
In particular, for all natural numbers m, one has (cf. [START_REF] Coppo | The Arakawa-Kaneko Zeta function[END_REF], Corollary 1) 1 -e -t Li k (1 -e -t ) dt , which is the analogue of (19) (with Li k in place of f k ) and also the original definition of the Arakawa-Kaneko zeta function (cf. [START_REF] Arakawa | Multiple zeta values, Poly-Bernoulli numbers and related zeta functions[END_REF], [START_REF] Coppo | The Arakawa-Kaneko Zeta function[END_REF]). Thus, taking into account the facts that ξ k (1) = ζ(k + 1) and Li 1 (1 -e -t ) = t, and following the same process as in the proof of Theorem 8, one obtains the following analogue of (20): ), we have for positive integers m the nice formula

ξ k (m + 1) = ∞ n=1 1 n k D 1 x m+1 (n) =
F k,ω (m + 1) = ∞ n=1 ω n n k D 1 x m+1 (n) = ∞ n=1 ω n P m (H n , H (2)
n , . . . , H

(m) n ) n k+1
, which extends (25). However, this formula is more theoretical than practical because of the fast increase in the size of polynomials P m : the number of monomials in P m is equal to the number p(m) of partitions of m, as shown by the explicit expression of the mth modified Bell polynomial.

Definition 3 .Definition 4 .

 34 Let f be a function in E. The Laplace transform L(f ) of f is defined byL(f )(x) = +∞ 0 e -xt f (t) dt for (x) > 0 .Proposition 2 (cf.[START_REF] Schiff | The Laplace transform : theory and applications[END_REF]). Let E = L(E) be the image ofE under L. If a is a function in E, then a) a is an analytic function of x in the half-plane (x) > 0 , b) a(x) → 0 when (x) → +∞ , c) L : E → E is an isomorphism. Let a ∈ E. TheBorel transform of a is the unique function a ∈ E such that a = L( a). One has the two reciprocal formulas a(t) = 1 2iπ c+i∞ c-i∞ e zt a(z) dz for all c > 0 and t > 0 , and a(x) = +∞ 0 e -xt a(t) dt for (x) > 0 .

Theorem 3 .

 3 The operator D is an automorphism of E which verifies D = D -1 and lets the function x → 1 x invariant.

Proposition 5 .

 5 If a and b are in E, then a(Λ) * b(Λ) ∈ E.

Definition 8 .

 8 Let a and b be two functions in E. The Λ-convolution product a b of a and b is defined by a b = Λ (Λ ( a) * Λ ( b)) , or equivalently (since Λ = (Λ ) -1 ) ( a b)(Λ) = a(Λ) * b(Λ) .

Definition 9 .

 9 Let a and b two functions in E. The harmonic product a b of a and b is defined by a b = L( a b) ∈ E .

Theorem 4 .Remark 5 .

 45 Let a and b be in E. Then, D(a b) = D(a) D(b) , (12) and D(ab) = D(a) D(b) . (13) Proof. One knows from Theorem 2 that D = LΛ L -1 . Hence D(a b) = LΛ L -1 (a b) = LΛ ( a b) = L(Λ ( a) * Λ ( b)) , and it follows from (5) and (9) that L(Λ ( a) * Λ ( b)) = L(Λ ( a))L(Λ ( b)) = D(a) D(b) which proves (12). Moreover, (12) enables us to write D(D(a) D(b)) = D 2 (a) D 2 (b) = ab (since D = D -1 ), and so D(a b) = D 2 (D(a) D(b)) = D(a) D(b)which proves[START_REF] Schiff | The Laplace transform : theory and applications[END_REF]. The values of (a b)(n) may be computed without recourse to a and b. By elementary transformations, it can be shown that (a b)(n + 1) = t-s )(e -t + e -s -e -t e -s ) n a(t) b(s)dtds .

Theorem 7 .

 7 The function F k analytically continues in the whole complex plane as an entire function.Proof. The function z → 1 log(1 -z) + 1 zbeing analytic in the disc D(0, 1) with a singularity at 1, we deduce from (1) that the radius of convergence of the series ∞ n=1

e 1 -

 1 -nt a(t) dt , and a formal permutation of n≥1 and +∞ 0 e -t e -t a(t) dt.

  Li k (e -u )(Λ(u)) s-1 du , and the change of variables t = Λ(u) leads to the integral representation ξ k (s

ξ 1 j=1(- 1 )m 2 , 6 Conclusion

 1126 k+1 (s) = ζ(k+1)ζ(s)+ kj ζ(k+1-j)Z j (s)+(-1) k 1 Γ(s) +∞ 0 t s-1 e -t 1 -e -t T k (t)dt .(34) In particular, in the simplest case k = 1, sinceT (t) = +∞ t e -u1 -e -u udu = (34) again gives the formulaξ 2 (s) = ζ(2)ζ(s) -s Arakawa and Kaneko (cf. [1] Theorem 6 (ii)).Most of the general results given for the modified zeta function F k , especially Theorem 7, Theorem 8, and Theorem 9, also apply (with minor adaptations) to a wide class of functions including the Arakawa-Kaneko zeta function ξ k , specifically to the class of functions represented by normalized Mellin transforms of typeF k,ω (s) = 1 Γ(s) +∞ 0 t s-1 e -t 1 -e -t f k,ω (1 -e -t ) dt with ω = (ω n ) n≥1 and f k,ω (z) = ∞ n=1 ω n n k z n .In particular, under the assumption that|ω n | n k = O( 1 n


Laurie Boithias
email: l.boithias@gmail.com

Supat Frederic C Do

S Siltecho

S Isarangkool

Na Ayutthaya

J Junjittakarn

C Hammecker

Transpiration, growth

Keywords:

INTRODUCTION

As a consequence of increasing demand for natural rubber and competition with oil palm in the traditional areas, the extension of Hevea brasiliensis cultivation is accelerating in drought-prone areas, such as in the southern part of northeast Thailand. Although growth and latex production are two distinct physiological phenomena in rubber trees [START_REF] Rao | Studies on the relationship between yield and meteorological parameters of para rubber tree (Hevea brasiliensis)[END_REF], they are both strongly related to evaporative demand and water availability. The importance of the water relations for latex production is known [START_REF] Pakianathan | Water relations and latex flow[END_REF], and some authors have studied the impact of water constraints on tree water status, girth increment and latex production [START_REF] Rao | Influence of soil, plant and meteorological factors on water relations and yield in Hevea brasiliensis[END_REF][START_REF] Chandrashekar | An analysis of growth and drought tolerance in rubber during the immature phase in dry subhumid climate[END_REF]. Water stress generally results from an imbalance between transpiration, driven by evaporative demand, and root-water-uptake allowed by soil water availability. Drought stress occurs whenever soil water availability drops below a certain threshold, thus inducing restrictions to transpiration and growth. Frequently, but not invariably, soil dryness is coupled with strong evaporation caused by air dryness. Conversely, an atmospheric drought may occur despite water being available in the soil. A model that predicts tree water use, tree growth and latex production model can be a useful tool to analyze key points of improvement in rubber tree cultivation in water-limited areas. Growth models have become a tool for making rational land use decisions and for monitoring sustainable agricultural systems if their capacity to predict crop response and identify management options is well established (e.g. [START_REF] Verdoodt | Daily simulation of potential dry matter production of annual field crops in tropical environments[END_REF]. They can help to identify environmental factors limiting growth and resource use and to assess productivity and profitability at various time scales (e.g. [START_REF] Shamudzarira | Simulating response of maize to nitrogen fertilizer in semi-arid Zimbabwe[END_REF]. For perennial crops much less work on production models has been carried out, probably due to data limitation, relatively high research costs and the difficulties of accumulated errors in long term simulations [START_REF] Zuidema | A physiological production model for cocoa (Theobroma cacao): model presentation, validation and application[END_REF].

The WaNulCAS model (Water Nutrients and Light Capture in Agroforestry Systems) was developed by the International Centre for Research on Agroforestry (ICRAF, Indonesia) as a research tool with the major objective to synthesize existing knowledge and hypotheses on above-and below-ground resource-use by trees and crops at the patch-scale [START_REF] Van Noordwijk | WaNuLCAS, a model of water, nutrient and light capture in agroforestry systems[END_REF]. Although it can undertake economic evaluations, WaNuLCAS is primarily a biophysical model that also provides a great range of outputs, including soil water, transpiration, growth of trees and crops for a given environment. This model was chosen for this study since it has already been parameterised for pure rubber tree stands and it included a specific module of latex production. However, it has mainly been used for long term (multi-annual) predictions of rubber tree growth in pure stand [START_REF] Yahya | Estimating growth (girth), yield and aboveground biomass of Hevea brasiliensis: Validation of the WaNuLCAS model[END_REF] or in agro-forestry systems [START_REF] Pinto | Simulation of agroforestry systems with sugarcane in Piracicaba, Brazil[END_REF]. Moreover, to our knowledge, this model has not been tested for rubber tree growing in drought prone area. Hence, this study had three aims: firstly, to test the possibility of using the WaNuLCAS model as an exploratory tool to analyze relationships between water use, growth and latex production on an intra-annual, daily basis and in a pure stand. The second was to evaluate the relative accuracy of the simulations of soil water content, tree transpiration, cambial growth and latex production. Finally, the third objective was to identify key features in the various parameterizations which may be adapted to improve overall prediction quality in drought prone areas. To address these issues, simulations were compared to field measurements in a mature rubber tree stand growing in the water-limited area of northeast Thailand. The period of analysis concerned seven months of full foliation and latex tapping from May to November. Four subperiods were distinguished, combining contrasted magnitudes of evaporative demand and soil water availability.

MATERIALS AND METHODS

The WaNuLCAS model

Version 3.2 of the WaNuLCAS model is used in this study. WaNuLCAS was developed as a collection of modules in which processes are computed at a daily time step. A short description of the modules used for the simulation is given below. For more theoretical background on the modules we refer to the original model documentation of [START_REF] Van Noordwijk | WaNuLCAS, a model of water, nutrient and light capture in agroforestry systems[END_REF].

The soil is divided into four layers in the vertical and four spatial zones in the horizontal, and each layer is characterised by different soil properties (Figure 1). The soil hydraulic module includes pedotransfer functions (Wösten et al.1998) that are used to derive the constants of the Van Genuchten (1980) equation depending on soil texture, bulk density and soil organic matter content. These constants determine the relationship between soil water content, soil water potential, hydraulic conductivity and potential tree water uptake in each cell. To calculate water infiltration, the model estimates a layerspecific field capacity. Two definitions of field capacity are used: (i) field capacity is the soil water content at which downward drainage would become less than a critical value of conductivity (Kcrit) and would effectively stop, or (ii) field capacity is the soil water content that is in hydrostatic equilibrium with a water table at a distance defined from the bottom of layer 4. The highest of these two values in any cell is used. A saturated hydraulic conductivity (Ksat) is generated from the pedotransfer functions. The water balance includes vertical and horizontal transport. It incorporates the water balance inflows and outflows described in Table 1 at daily time steps.

Surface run-off occurs if (i) daily rainfall exceeds daily maximum infiltration or if (ii) daily rainfall exceeds the soil potential water storage. Allocation to surface runoff or infiltration in the first soil layer depends on a ratio, called SSI in this study, that is the ratio between a reference Ksat (for degraded soil) and the Ksat generated by the pedotransfer function. Infiltration then follows the tipping bucket principle for wetting subsequent layers of soil, filling a cascade of soil layers up to their field capacity.

Tree water uptake, i.e. transpiration, is driven by evaporative demand, with the possibilities determined by tree root length density and soil water availability in the various cells to which a plant has access. A "potential transpirational demand" is estimated from LAI and ET 0 . Hence, the transpiration depends on a water demand reduction factor or stress factor:

0 ) ( τ θ τ f =
where τ is the tree water demand, f(θ) is the water demand reduction factor (where θ is the soil water content) and τ 0 is the potential transpirational demand. The reduction factor depends on the "plant water potential" which is estimated from the "soil water potential".

Canopy light capture is determined by the LAI, the canopy height and the Beer's law extinction coefficient specific to rubber trees. The growth reserve pool varies with the tree potential growth and the minimum stress factor regarding light, water and nitrogen. Part of the growth reserves are allocated to growth. Allometric equations are used to relate tree girth to tree biomass. They allocate predicted biomass to the different tree components: below-ground and total, above-ground biomass, which in turn is partitioned into roots, leaves and twigs, as well as branches. As for tree growth, part of the "tree growth reserves pool" is allocated to "daily latex production". Tapped latex depends on the "fraction of latex stock" that can be tapped every day ("tapping fraction") and on existing "Brown Bast", a physiological disease that induces a reduction of latex production [START_REF] Paardekooper | Exploitation of the rubber tree[END_REF]. The version of WaNuLCAS used only predicts latex production using a d/2 tapping pattern (one day of tapping -one day of rest).

Data measurements

Experiment site and forestry system: The experiment site was located in Baan Sila (N15° 16' 23.6" E103° 04' 51.3", altitude = 150 m) between Satuk and Khu Muang in the Buriram province of northeast Thailand. The field was a monoculture stand of RRIM600 rubber clones planted in 1995 (i.e. 12 years old in 2007) at a 7mx2.5m tree spacing over 5.5ha, in a deep sandy soil. Average annual rainfall in this area is 1176 mm (Khu Muang station). The wet season lasts from May to October but intermittent dry spells often occur in June and July. Rubber trees have been tapped since June 2003. The stand was modelled as a half alley system as shown in Figure 1. The 2-D unit represents a flat area 3.5m wide and 1.8m deep. The width was arbitrarily split into four zones of 1, 1, 1 and 0.5 m. Layer bottoms were set according to the main variations of soil properties in the profile, at depths of 0.25, 0.5, 1 and 1.8m. Trees were grown in zone 1 at a density of 571 trees.ha -1 (Table 2 (a)).

Water balance

Soil physical properties: the main soil properties [START_REF] Hartmann | Rubber tree Trunk Phloem Necrosis in North-East Thailand: 1-Investigations on soil heterogeneities and disease distribution, IRRDB Meeting[END_REF] of each modelled layer are provided in Table 2 ( rainfall was approximately 965mm (879mm from May to November), 18% less than the long-term annual average. The wet season is often interspersed with a short dry and hot period (scarce rainfall and high evaporative demand) in June and July, a feature which occurred in 2007. Dry and cool conditions prevailed again from November. Climate data inputs were set equal in the four spatial zones. In the rubber stand, volumetric soil water content was measured once or twice each month with a neutron probe (Troxler 3300, USA) calibrated for the experimental soil. Calibration was carried out separately for the upper (0-0.2m depth) and lower (0.2-1.8m depth) layers. Measurements were made every 0.2m from the ground down to a depth of 1.8m. Twelve access tubes of 2.0m in length were set up in pairs: one tube in the planting line between two trees, and the other one in the middle of the inter-row. Average measured moisture from the 12 tubes were used as observational data. Initial soil water contents (i.e. in May) are given in Table 2 (b). The Kcrit value was calculated in order to get the field capacity of layer 1 equal to the soil water content at a soil water potential of -100 cm. No water table was observed below the stand in Baan Sila. A value of 10m depth was set as input for the distance from the bottom of layer 4 to the supposed water table. Thus, the value of field capacity computed from definition (ii) was insignificant. The SSI parameter (influencing the ratio between run-off and infiltration) was calibrated in order to obtain soil water stock values during periods of high rainfall (August-September-October) similar to observed values. Daily tree water use, i.e. transpiration, was estimated by the daily total sap flow in trunk xylem, ignoring changes in tree water storage [START_REF] Na Ayutthaya | Transient thermal dissipation method of xylem sap flow measurement: multi-species calibration and field evaluation[END_REF]. The measurements of xylem Sap Flow Density (SFD) on 11 trees in Baan Sila were made using the Transient Thermal Dissipation (TTD) system [START_REF] Do | Influence of natural temperature gradients on measurements of xylem sap flow with thermal dissipation probes. 2. Advantages and calibration of a non-continuous heating system[END_REF][START_REF] Na Ayutthaya | Transient thermal dissipation method of xylem sap flow measurement: multi-species calibration and field evaluation[END_REF]. A cycle of 10min heating and 20min cooling was used to measure sap flow density every 30min. The zero flux signal was determined every night assuming that sap flow was negligible at the end of the night (see details in Isarangkool Na [START_REF] Na Ayutthaya | Transient thermal dissipation method of xylem sap flow measurement: multi-species calibration and field evaluation[END_REF]. Probes were inserted into the trunks at a height of 1.8m above the soil. At this height, average sapwood area was estimated to be 2.04 dm 2 (SD = 0.47). Three probes were inserted into each trunk to take circumferential variability into account. An average radial profile of sap flux density was also taken into account [START_REF] Na Ayutthaya | Transient thermal dissipation method of xylem sap flow measurement: multi-species calibration and field evaluation[END_REF]. All probes were connected to a data logger (CR10X, Campbell Scientific, Leicester, U.K.). Hourly SFD was cumulated over a 24h period to calculate daily SFD and multiplied by sapwood area. Average daily total flow of the representative trees was divided by the soil surface theoretically available for a tree (17.5m²) in order to estimate the observed tree water uptake in mm. Isarangkool Na Ayutthaya et al. ( 2011) provide a detailed analysis of a data subset of these field measurements of water use.

Tree growth and light capture: the potential growth (0.01kg.m -2 ) was taken from WaNuLCAS default values. As inputs, the model required initial stem and canopy above-ground biomasses (Table 2 (a)). Total above-ground biomass was calculated using allometric relation between girth and weight. Constants for RRIM600 rubber trees were set according to [START_REF] Chantuma | Carbon sequestration: direct rubber biomass measurement and visual wood weight[END_REF]:

5623 . 2 0082 . 0 G W * =
where W is the above ground tree weight in kg and G the tree girth in cm, measured at 170cm above the ground. The simulated girth increments were compared to observations. The girths of 232 trees were measured using a metre rule once or twice per month in 2007 in Baan Sila. The maximum LAI (Table 2 (a Latex production: Observed latex production was estimated from the amount of rubber sheets produced each tapping day on the whole plot, knowing the number of tapped trees and an average weight per rubber sheet of 1200g. The "tapping fraction" was calibrated to fit the range of expected latex production at the start of tapping period in May, taking into account the fact that 30% of the trees of the plot were suffering from Brown Bast and that the tapping system in the model could only be d/2 instead of the observed 2d/3. The used parameters in the latex module are given in Table 2 (

Evaluation of model predictions

The simulation covers the period from May the 1 st to November the 30 th in 2007. The simulation assumed no limitations due to nutrients (nitrogen or phosphorus), weeds and pests. Predicted data regarding soil water content, transpiration, girth increment and latex production were compared to field measurements. The analysis used determination coefficients (R 2 ) and relative root mean square errors RMSE rel :

n O P O RMSE n i i i rel ∑ = - = 1 2 ) ( 100
where O and P represent observed and predicted values, respectively, and O represents the observed value average. The number of observations is n. R² and RMSE rel for transpiration estimates were compared over the seven months of study and for the four climatic sub-periods (Figure 2): one combining high evaporative demand and low soil water constraint (May); two periods combining high evaporative demand and high soil water constraint (June-July and November); and one combining low evaporative demand and low soil water constraint (August-September-October).

RESULTS

Water balance

Soil water content: Observations ranged from 5 to 25% (Figure 3). Soil moisture decreased in each layer during the short dry period of June and July. Model results followed the same pattern but they were often overestimated, particularly in layer 3. R² of layers 1, 2, 3 and 4 were 0.218, 0.307, 0.277 and 0.496, respectively, and RMSE rel values were 34.7, 39.8, 68.7 and 22.3%, respectively. For the total profile, R² and RMSE rel values were 0.461 and 34.9%, respectively.

Tree transpiration: Measured values ranged from 0.3 to 2.3mm.d -1 (Figure 4).

Transpiration dramatically decreased during the intermittent droughts in June and July. It also progressively decreased in November at the beginning of the long dry period.

Predicted and observed values were in the same range during periods without soil water constraint and with low or moderate evaporative demand (ET 0 < 2mm.d -1 , September-October). In contrast, model-predicted values were significantly overestimated in periods of high evaporative demand (ET 0 > 2mm.d -1 : start of May, June and July, November).

Moreover, the dramatic reduction of transpiration in June and July due to soil drought was not simulated. Finally, R² and RMSE rel values equalled 0.104 and 94 %, respectively.

Tree growth

Observed tree girth increased at the beginning of the rainy season (May-June) and then stopped and decreased in July in accordance with the decline in transpiration (Figure 5).

Girth increased again from August to October. The total increment reached 1.3cm over the seven-month period. The model significantly overestimated girth increment. Final R² and RMSE rel values were 0.916 and 208 %, respectively.

Latex prediction

Observed values ranged from 1.3 and 19kg.ha -1 .d -1 and the pattern was not correlated to variation in transpiration (Figure 6). As for growth, the model did not predict any reduction due to water stress and the latex production is largely overestimated. Total predicted production (2888 kg.ha -1 ) was twice the measured production (1416kg.ha -1 ) over the simulated period. Finally, R² and RMSE rel values were 0.423 and 169 %, respectively.

DISCUSSION

WaNuLCAS model as an exploratory tool

The test of the possibility of using WaNuLCAS model as an exploratory tool for a pure rubber tree stands appeared to be successful. It simulated daily changes in transpiration, girth increment and latex production within plausible ranges [START_REF] Silpi | Effects of tapping activity on the dynamics of radial growth of Hevea brasiliensis trees[END_REF][START_REF] Pakianathan | Water relations and latex flow[END_REF][START_REF] Watson | Climate and soil[END_REF]. In addition, the model runs were easily implemented and several shortcuts were available to test the schematic inputs of rainfall, ET 0 , etc.

Relative accuracy of simulations

With the current model version and parameters, results showed large inaccuracies for certain model predictions: soil water content (R² = 0.461, RMSE rel = 34.9%), tree transpiration (R² = 0.104, RMSE rel = 94 %), tree girth increment (R² = 0.916, RMSE rel = 208.4%) and latex production (R² = 0.423, RMSE rel = 169.5%). However, as the water shortage in June and July was not predicted by the model, the significant overestimations in girth increment and latex production appear logical and can not be discussed in further depth. Hence, the focus of this investigation is foremost on the soil water balance and tree water use. Model estimations of tree transpiration were correct in conditions without soil and atmospheric droughts. Major inaccuracies resulted from (i) the non-prediction of soil water shortage and (ii) the overestimation of tree transpiration under high evaporative demand and non-limiting conditions of soil water.

The soil water shortage was not predicted because the soil water content was overestimated in all layers except in the lowest one. In WaNuLCAS, as in all reservoir models, infiltration water fills each layer of soil until their field capacity is reached and the excess moisture wets the subsequent soil layer. The speed of the process depends on the infiltration rate. When the deepest layer reaches its field capacity, excess water is then lost by deep drainage. A small error in predicting soil water content in upper layers may thus be amplified in deeper layers. Bias in predicting soil water content can come from errors in predicting (i) the layers field capacity and (ii) the infiltration rate of each layer, particularly in the upper-most layer which determines the infiltration/runoff ratio. Aware of these points, measured Ksat and a calibrated ratio between infiltration and runoff (SSI) were used. The default values of Ksat and SSI provided by the pedotransfer functions of WaNuLCAS were even higher than our estimated values and induced even greater overestimates of soil water content. Other authors reported inadequacies of the pedotransfer functions and predicted infiltration by WaNuLCAS for tropical soils [START_REF] Hodnett | Marked differences between van Genuchten soil water-retentions parameters for temperate and tropical soils: a new waterretention pedotransfer function developed for tropical soils[END_REF]Walker et al., 2007;[START_REF] Pansak | Assessing soil conservation strategies for upland cropping in Northeast Thailand with the WaNuLCAS model[END_REF]. To explain the overestimate of soil water content, despite the use of measured Ksat and calibrated SSI values, several hypotheses are proposed. Firstly, measured Ksat is still overestimated, particularly in the upper layer. This is possible because the variability between field measurements of Ksat was very large (Table 2 (b)). Moreover, significant surface flooding and run-off were observed for each important rainfall event. Secondly, the infiltration of water is not homogeneous between the tree row and in the inter-row space.

Surface flooding and run-off was of greater importance in the inter-row spaces which are at a lower level than the tree rows. In addition, preferential vertical flow may occur in the inter-row spaces, which may induce deep, localised drainage.

In addition, the model overestimated maximum transpiration because it simulated tree transpiration that follows the evaporative demand without any stomatal regulation due to atmospheric conditions. The only water-regulation included in the model originates from soil water shortage, similar to the majority of models. However, observed data showed that effective regulation of transpiration occurred above ET 0 equal to 2.0mm (Fig. 4). This is, however, a very recent result for rubber trees (Isarangkool Na Ayyuthaya et al., 2011), which may differ according to the clone (here RRIM 600). The regulation of transpiration at high evaporative demand, whatever the soil water availability, has been quoted for several species and environments [START_REF] David | Constrains on transpiration from an evergreen oak tree in southern Portugal[END_REF][START_REF] Bovard | Environmental controls on sap flow in a northern hardwood forest[END_REF][START_REF] Oguntunde | Measurement and modeling of transpiration of a rain-fed citrus orchard under subhumid tropical conditions[END_REF][START_REF] Bush | Wood anatomy constrains stomatal response to atmospheric vapor pressure deficit in irrigated, urban trees[END_REF]. Moreover, the overestimate of transpiration in the simulation should further decrease the soil water reserve, which emphasizes the current overestimation of soil water content by the model.

The current version of the WaNuLCAS model has, however, already predicted correct values of rubber tree growth and production over long periods in other regions [START_REF] Pinto | Simulation of agroforestry systems with sugarcane in Piracicaba, Brazil[END_REF][START_REF] Yahya | Estimating growth (girth), yield and aboveground biomass of Hevea brasiliensis: Validation of the WaNuLCAS model[END_REF]. We assume that these results likely correspond to environmental conditions without high evaporative demand.

Key points of parameterization and improvement

Under the conditions of degraded soil and high evaporative demand investigated, the key point of parameterization and improvement are the soil water balance and the regulation of tree transpiration under high evaporative demand. It is likely that the ability to separate soil properties between horizontal zones in the model will facilitate the improvement of the simulated soil water balance. This is particularly the case in rubber tree plantations, where several years of different soil-, weed-and litter-management practices between tree row and inter-row spaces may influence the topography and soil properties related to soil water infiltration and balance. Secondly, the regulation of transpiration as a function of evaporative demand requires the introduction of a new equation or reduction factor.

The transpiration model in [START_REF] Granier | A generic model of forest canopy conductance dependent on climate, soil water availability and leaf area index[END_REF] provides an example of a multiple regulation of canopy conductance where atmospheric drought is taken into account whatever the soil water conditions.

CONCLUSIONS

In conclusion, the WaNuLCAS model appears to be a useful exploratory tool for the simulation of water use, growth and latex production in a pure stand of rubber tree on a daily basis. Results do, however, show that under the conditions of degraded soil and high evaporative demand investigated, the modules of soil water balance and tree transpiration require particular parameterizations and improvement in order to more accurately represent these processes.

  b). Textures were assumed homogeneous in the four zones. The soil was degraded and poorly structured. Bulk density was high in layer 2: this compaction resulted from previous cassava cultivation. Saturated hydraulic conductivities (Ksat) were measured with four replications. The magnitude of the measured standard deviation suggested significant spatial variability of this parameter. Soil water storage, infiltration and evaporation: WaNuLCAS requires the daily inputs of two climatic variables: rainfall and reference evapotranspiration. In Baan Sila, climatic data (air temperature, relative humidity, global radiation, wind speed and rainfall) were automatically recorded every 30 minutes with a Minimet station (Skye Instruments Ltd, U.K.) in an open field located to the east of the rubber stand, 50m away from any tree. Recorded climate data were then aggregated as maximal/minimal (temperature, humidity), average (wind speed) or cumulated (radiation, rainfall) daily values. Daily values of reference potential evapotranspiration were calculated using the Penman-Monteith equation. Figure 2 shows the input values of climatic data. 2007

  )) for the year 2007 was measured by collecting fallen leaves in one square meter open boxes from November 2007 till the complete defoliation at the end of January 2008. Root biomass and root length density were assumed constant for the simulation period. The average profile of root length density was taken from measurements in the site (Table 2 (b), Pierret et al., unpublished report). Default model values were used for other physiological parameters.

  c). Other parameters have been set following WaNuLCAS default values.

Figure 1 .

 1 Figure 1. General layout of soil layers and spatial zones in the WaNuLCAS model. Modelled trees were planted in Zone 1. Adapted from Van Noordwijk and Lusiana (1999).

474 25Figure 2 .

 4742 Figure 2. Daily climate data from May to December 2007: rainfall and reference evapotranspiration.

Figure 3 .

 3 Figure 3. Observed and predicted soil water content: (a) and (b) Layer 1; (c) and (d) Layer 2; (e) and (f)

  27

Figure 4 .

 4 Figure 4. (a) Observed and predicted tree transpiration; (b) Predicted vs. observed. Diagonal line indicates the 1-1 line.

Figure 5 .

 5 Figure 5. (a) Observed and predicted cumulated stem girth increment; (b) Predicted vs. observed. Diagonal line indicates the 1-1 line.

Figure 6 .

 6 Figure 6. (a) Observed and predicted latex production; (b) Predicted vs. observed. Diagonal line indicates the 1-1 line.

Acknowledgments. The authors thank all French and Thai counterparts from the Institute of Research for Development (IRD), the Khon Kaen University and the Land Development Department. We particularly thank the Land Development Department of Bangkok and Dr. Darunee Chairod for the assistance in neutron probe measurements. We are grateful to Dr Henri Robain for his support as coordinator of the agreement program between IRD and the Rubber French Institute. Were also are very grateful to Mr. Chaipat Sirichaiboonwat who welcomed us so kindly in his rubber tree plantation. Finally, we sincerely thank Ms Ni'matul Khasanah (ICRAF, Indonesia), Dr Fiona Tummon (CNRS, France) and Dr Erwan Motte (STARLAB, Spain) for their kind assistance in achieving this work.

Funding. This Thai-French research was funded by the French Institute of Research for Development (IRD), the Rubber French Institute (IFC) and by Michelin/Socfinco/SIPH Plantations companies.